% 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;                %方差陣