pca.m

利用奇异值法（SVD）对载入数据进行主成分分析

% s 环境阵，p 特征向量阵 ,A 特征值阵;并假定方差阵 R 为单位阵,X 遥外测差向量
% S 是6n*m的矩阵 
% S'*S为实对称阵
 %输入S以及向量X;
load('S.mat');
load('X.mat');
 s=S';
 c=s*S;
[p,A]=schur(c);         % A为15阶
p1=p';
for j=1:26             %主成分特征值排序
    for i=1:(27-j)
        if A(i,i)<A(i+1,i+1)
            n=A(i,i); 
            A(i,i)=A(i+1,i+1);
            A(i+1,i+1)=n;
            m=p1(i,:);
            p1(i,:)=p1(i+1,:);
            p1(i+1,:)=m;
        end 
    end
end                  %  解出新的排序后的特征向量阵
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                          %选择主成分
 sum=A(1,1);
 for i=2:27
     sum=sum+A(i,i);
 end
SUM=A(1,1);
for i=2:27
    if (1-SUM/sum)<0.0000001
        break
    else SUM=SUM+A(i,i);
        N=i;
    end 
end
P2=p1(1:N,:);  %主成分矩阵
Pa=P2';
Aa=A(1:N,1:N);   %主成分特征值矩阵
Ab=pinv(Aa);    %求逆 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D=Pa*Ab*P2*s*X;          %主成分估计值  
Ed=Pa*Ab*P2;                %方差阵