function [c,m,v]=norm1(C) 
%對數(shù)據(jù)進行標準化處理 
[n,s]=size(C); 
for i=1:n 
for j=1:s 
c(i,j)=(C(i,j)-mean(C(:,j)))/sqrt(cov(C(:,j))); 
end 
end 
m=mean(C); 
for j=1:s 
v(1,j)=sqrt(cov(C(:,j))); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function w=maxdet(A) 
%求矩陣的最大特征值 
[v,d]=eig(A); 
[n,p]=size(d); 
d1=d*ones(p,1); 
d2=max(d1); 
i=find(d1==d2); 
w=v(:,i); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,q,w,wh,f0,FF]=fun717(px,py,C) 
% px自變量的輸入個數(shù) 
% py因變量的輸入個數(shù)。 
% C輸入的自變量和因變量組成的矩陣 
% t提取的主成分 
% q為回歸系數(shù)。 
% w最大特征值所對應的特征向量。 
% wh處理后的特征向量 
% f0回歸的標準化的方程系數(shù) 
% FF原始變量的回歸方程的系數(shù) 
c=norm1(C); %norm1為標準化函數(shù) 
%y=c(:,px+1:px+py); %截取標準化的因變量 
E0=c(:,1:px); %截取標準化的自變量
F0=c(:,px+1:px+py); %截取標準化的因變量
A=E0'*F0*F0'*E0; %關(guān)聯(lián)矩陣
w(:,1)=maxdet(A); % maxdet為求取最大特征值的函數(shù)
t(:,1)=E0*w(:,1); %第一主成分 
p(:,1:px)=(E0'*t(:,1)/(t(:,1)'*t(:,1)))'; % p參
E(:,1:px)=E0-t(:,1)*(E0'*t(:,1)/(t(:,1)'*t(:,1)))'; % 殘差矩陣
%進入循環(huán)
for i=0:px-2 
B(:,px*i+1:px*i+px)=E(:,px*i+1:px*i+px)'*F0*F0'*E(:,px*i+1:px*i+px); 
w(:,i+2)=maxdet(B(:,px*i+1:px*i+px));  
t(:,i+2)=E(:,px*i+1:px*i+px)*w(:,i+2); 
p(:,px*i+px+1:px*i+2*px)=(E(:,px*i+1:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t(:,i+2)))'; 
E(:,px*i+px+1:px*i+2*px)=E(:,px*i+1:px*i+px)-t(:,i+2)*(E(:,px*i+1:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t(:,i+2)))'; 
end 

for s=1:px 
q(:,s)=p(1,px*(s-1)+1:px*s)'; 
end 

[n,d]=size(q); 

for h=1:px 
iw=eye(d); 
for j=1:h-1 
iw=iw*(eye(d)-w(:,j)*q(:,j)'); 
end 
wh(:,h)=iw*w(:,h); 
end

for j=1:py 
zr(j,:)=(regress(F0(:,j),t))';
end 

for j=1:px 
for i=1:py 
w1=wh(:,1:j); 
zr1=(zr(i,1:j))'; 
f0(i,:,j)=(w1*zr1)'; %生成標準化變量的方程的系數(shù)矩陣 
end 
[c,m,v]=norm1(C); 
ccxx=ones(py,1)*m(1,1:px); 
ccy=(v(1,px+1:px+py))'*ones(1,px); 
ccx=ones(py,1)*(v(1,1:px)); 
ff=ccy.*f0(:,:,j)./ccx; 
fff=-(sum((ccy.*ccxx.*f0(:,:,j)./ccx)')-m(1,px+1:px+py))'; 
FF(:,:,j)=[fff,ff]; %生成原始變量方程的常數(shù)項和系數(shù)矩陣 
end 


實驗結(jié)果：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t,q,w,wh,f0,FF]=fun717(3,3,C)

t =

   -0.6429   -0.5914   -0.1305
   -0.7697   -0.1667    0.1339
   -0.9074    0.5212    0.0477
    0.6884    0.6800    0.3458
   -0.4867   -1.1328   -0.1818
   -0.2291    0.0717    0.0247
   -1.4037    0.0767   -0.5718
    0.7436    0.2106   -0.0322
    1.7151    0.6549   -1.5572
    1.1626   -0.1668    0.3329
    0.3645   -0.7007    0.2011
    0.7433   -0.6983    0.0023
    1.1867    0.7570    0.3362
   -4.3898    0.7600    0.2554
   -0.8232   -0.9738   -0.0826
   -0.7490    0.5211   -0.6672
   -0.3929    0.2034    0.5642
    1.1993   -0.7827    0.0924
    1.0485   -0.3729    0.3190
    1.9424    1.1294    0.5677


q =

   -0.6659    0.0198   -0.6575
   -0.6760    0.3547    0.2871
    0.3589    1.1942   -0.6967


w =

   -0.5899   -0.4688   -0.6575
   -0.7713    0.5680    0.2871
    0.2389    0.6765   -0.6967


wh =

   -0.5899   -0.3679   -0.9346
   -0.7713    0.6999    0.8023
    0.2389    0.6356   -0.2228


f0(:,:,1) =

   -0.2015   -0.2635    0.0816
   -0.2454   -0.3209    0.0994
   -0.0843   -0.1103    0.0342


f0(:,:,2) =

   -0.0778   -0.4989   -0.1322
   -0.1385   -0.5244   -0.0854                       %取兩個成分的標準數(shù)據(jù)回歸系數(shù)
   -0.0604   -0.1559   -0.0073


f0(:,:,3) =

    0.3683   -0.8818   -0.0258
    0.2872   -0.8898    0.0161
   -0.2590    0.0146   -0.0546


FF(:,:,1) =

   29.2002   -0.0431   -0.4350    0.0598
  430.2511   -0.6220   -6.2712    0.8625
  150.4807   -0.1752   -1.7662    0.2429


FF(:,:,2) =

   47.0197   -0.0167   -0.8237   -0.0969
  612.5671   -0.3509  -10.2477   -0.7412            %取兩個成分的原始數(shù)據(jù)回歸系數(shù)
  183.9849   -0.1253   -2.4969   -0.0518


FF(:,:,3) =

   47.9684    0.0788   -1.4558   -0.0190
  623.2817    0.7277  -17.3872    0.1393
  179.8868   -0.5379    0.2338   -0.3886


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%主成分分析程序
px=3;
py=1;
c=norm1(A);
E0=c(:,1:px);
F0=c(:,px+py);
f0=pinv(E0'*E0)*E0'*F0;
disp(['y1=' num2str(f0(1)) 'x1+(' num2str(f0(2)) ')x2+(' num2str(f0(3)) ')x3'])
[c,m,v]=norm1(A); 
%待改進
for j=1:py 
ccxx=ones(py,1)*m(1,1:px); 
ccy=(v(1,px+1:px+py))'*ones(1,px); 
ccx=ones(py,1)*(v(1,1:px)); 
ff=ccy.*f0(:,:,j)./ccx; 
fff=-(sum((ccy.*ccxx.*f0(:,:,j)./ccx)')-m(1,px+1:px+py))'; 
FF(:,:,j)=[fff,ff];
end 
FF