mainfactoranalysis.m

matlab《主成分分析》的源代码,是化学计量学的模式识别的主要程序。

function [t,p1,eigv]=pcaana(x)
    [n,p]=size(x);
    minum=n;
    ml=n;
    if p<minum
        minum=p;
    else
        ml=p;
    end
    %DIM e1(1,1),EIGV(1,1),RATIO(1,1),RSD(1,1),ind(1,1)
    t=zeros(n,minum);
    p1=zeros(p,minum);
    e1=zeros(n,p);
    eigv=zeros(1,minum);
    ratio=zeros(1,minum);
    rsd=zeros(1,minum-1);
    ind=zeros(1,minum-1);
    %   mean-centred
    %   DIM x1(1),x2(1)
    x1=zeros(p);
    x2=zeros(p);
    for i=1:p
         x1(i)=0;
         x2(i)=0;
         for j=1:n
             x1(i)=x1(i)+x(j,i);
             x2(i)=x2(i)+x(j,i)*x(j,i);
         end
         lxx=x2(i)-x1(i)*x1(i)/n;
         x1(i)=x1(i)/n;
         x2(i)=sqrt(lxx/(n-1));
    end
    for i=1:p
        for j=1:n
            x(j,i)=(x(j,i)-x1(i))/x2(i)/(n-1);
        end
    end
    %   NIPALS algorithm
    %   calculating score and load matrix
    %    DIM e(1,1),t0(1,1),tc(1,1),pc(1,1)
    %   DIM tct(1,1),pct(1,1)
    %   DIM tcttc(1,1),pctpc(1,1)
    e=x;
    tc=zeros(n,1);
    pc=zeros(p,1);
    t0=zeros(n,1);
    for z=1:minum
        for j=1:n
            tc(j,1)=e(j,1);
        end
        eeand=10000;
        while eeand>=1e-04
           t0=tc;
           tct=tc';
           tcttc=tct*tc;
           ttt=1/tcttc(1,1);
           pct=tct*e;
           pct=ttt*pct;
           pc=pct';
           pctpc=pct*pc;
           ppp=1/sqrt(pctpc(1,1));
           pct=ppp*pct;
           pc=pct';
           pctpc=pct*pc;
           ppp=1/pctpc(1,1);
           tc=e*pc;
           tc=ppp*tc;
           eeand=0;
           for j=1:n
                eeand=eeand+abs(t0(j,1)-tc(j,1));
           end
       end
         eigv(1,z)=0;
        for j=1:n
             eigv(1,z)=eigv(1,z)+tc(j,1)*tc(j,1);
        end
        for j=1:n
             t(j,z)=tc(j,1);
        end
        for j=1:p
             p1(j,z)=pc(j,1);
        end
        for j=1:n
            for k=1:p
                 e1(j,k)=t(j,z)*p1(k,z);
            end
        end
        e=e-e1;
    end
    %   calculation four principal component parameters
    total =0;
    for i=1:minum
         total=total+eigv(1,i);
    end
    for i=1:minum
         typ=0;
        for j=1:i
             typ=typ+eigv(1,j);
        end
         ratio(1,i)=typ/total;
    end
    for i=1:minum-1
         tot=0;
        for j=i+1:minum
             tot=tot+eigv(1,j);
        end
         rsd(1,i)=sqrt(tot/ml/(minum-i));
         ind(1,i)=rsd(1,i)/(minum-i)/(minum-i);
    end
    subplot(2,2,1)
    plot(rsd,'.-')
    subplot(2,2,2)
    plot(ind,'.-')
    subplot(2,2,3)
    plot(ratio,'.-')
    subplot(2,2,4)
    plot(eigv,'.-')