83.100 	88.200 	4067.000 	2858.900 
87.000 	90.500 	6109.000 	3221.200 
88.200 	92.800 	4865.000 	2942.500 
83.500 	85.100 	4135.000 	2603.100 
84.300 	87.600 	4163.000 	2840.600 
86.400 	93.700 	3507.000 	3040.900 
82.000 	89.300 	4067.000 	3062.000 
82.900 	90.900 	3948.000 	3594.000 
88.000 	91.300 	4182.000 	2622.200 
84.000 	94.700 	7058.000 	2894.300 
85.700 	93.100 	7159.000 	3005.000 
83.100 	86.800 	3452.000 	2834.500 
84.600 	86.500 	5115.000 	2963.000 
89.900 	90.500 	4135.000 	3056.500 
85.500 	84.700 	4650.000 	2881.700 
83.900 	93.600 	6767.000 	3158.500 
83.800 	92.300 	3948.000 	3202.500 
84.900 	87.200 	3950.000 	3980.600 
92.300 	89.700 	4400.000 	2679.500 
81.600 	86.100 	4163.000 	3069.700 
85.600 	89.800 	3569.000 	3192.400 
82.800 	87.800 	4515.000 	2818.100 
81.000 	89.800 	1661.000 	2906.400 
82.200 	92.700 	7583.000 	2990.100 
90.100 	95.500 	4351.000 	3917.400 
84.500 	91.100 	5727.000 	2534.500 
83.700 	89.700 	4981.000 	2670.400 
85.300 	91.300 	5159.000 	2709.300 
81.500 	89.700 	3548.000 	2774.700 
85.600 	88.200 	4135.000 	2938.900 
85.600 	88.200 	4135.000 	2938.900 
86.900 	92.200 	6576.000 	2998.400 
86.500 	87.200 	7248.000 	3123.100 
80.600 	98.900 	4050.000 	3265.400 
83.600 	93.800 	5274.000 	3191.200 
77.100 	86.000 	1736.000 	2282.600 
81.000 	87.700 	1682.000 	3026.000 
86.500 	91.400 	5059.000 	3205.300 
82.200 	90.100 	5833.000 	3748.600 
86.300 	89.200 	3948.000 	2628.300 
86.900 	93.900 	4650.000 	2973.600 
90.000 	89.500 	4328.000 	3062.000 
85.300 	85.900 	1134.000 	3223.000 
84.300 	94.000 	5624.000 	3180.800 
84.000 	83.800 	1809.000 	2730.000 
85.700 	91.900 	4213.000 	3173.400 
84.300 	89.500 	4650.000 	2818.900 
81.400 	85.300 	1524.000 	2291.400 
81.600 	89.400 	4086.000 	2864.700 
81.600 	89.400 	4086.000 	2864.700 
82.800 	91.300 	6183.000 	3000.300 
85.500 	92.700 	4865.000 	2843.200 
];
[n3,m3]=size(x);        %求x的行和列，其中m3為X的變量個數(shù)
aver_x=mean(x);          %x矩陣的均值，是個行向量
std_x=std(x);            %x矩陣的標(biāo)準(zhǔn)差，是個行向量
for j=1:m3                %對x進(jìn)行標(biāo)準(zhǔn)化
    for i=1:n3
        sdd_x(i,j)=(x(i,j)-aver_x(j))/std_x(j);
    end
end
x0=sdd_x

    [n4,m4]=size(y);           %求y的行和列，其中m4為y的變量個數(shù)
aver_y=mean(y);                %y矩陣的均值，是個行向量
std_y=std(y);                   %y矩陣的標(biāo)準(zhǔn)差，是個行向量
for j=1:m4                      %對y進(jìn)行標(biāo)準(zhǔn)化
    for i=1:n4
        sdd_y(i,j)=(y(i,j)-aver_y(j))/std_y(j);
    end
end
y0=sdd_y

r=rank(x0);          %標(biāo)準(zhǔn)化后x的秩(好象沒用上)
h=0;                 %用來統(tǒng)計循環(huán)次數(shù)的
flag=0;              %是一個只能取0和1的符號，用來控制循環(huán)，跳出循環(huán)
Xsh=struct('t',[],'u',[],'w',[],'v',[],'a',[],'b',[]);        %是一個結(jié)構(gòu)體，用來記錄每次循環(huán)并且把循環(huán)的值放到同一個變量中
  for i=1:1000
   while flag==0       %當(dāng)flag==0時候才能執(zhí)行
       h=h+1;          %用來統(tǒng)計循環(huán)次數(shù)的
    M=x0'*y0*y0'*x0;   %根據(jù)書上的公式
    [V,D]=eig(M);       %求M的特征向量V和特征值D
     p=sqrt(max(max(D)));     %求M特征值里面最大的值
    [n,m]=size(D);      %以下步驟到k，求最大特征值所在的列，然后求出最大特征值所對應(yīng)的特征向量
    max_d=D(1,1);        %是求出一個矩陣?yán)锩婺莻€值最大，并且記錄最大值的行和列
    k=1;
     for i=1:n 
        for j=1:m
           if D(i,j)>max_d
               max_d=D(i,j);
               k=j;
           end
        end
      end
max_d;
k;
  Xsh(h).w=V(:,k);      %w為x的權(quán)重，把最大的特征值對應(yīng)的特征向量附給結(jié)構(gòu)體中的w
    
  Xsh(h).v=y0'*x0*Xsh(h).w/p;   %v為y的權(quán)重，并寫到構(gòu)體中
  Xsh(h).t=x0*Xsh(h).w;           %t為x的得分向量，并寫到構(gòu)體中
  Xsh(h).u=y0*Xsh(h).v;         %u為y的得分向量，并寫到構(gòu)體中
  Xsh(h).a=(inv(Xsh(h).t'*Xsh(h).t)*Xsh(h).t'*x0)';   %a為x的載荷量，就是回歸方程中參數(shù)向量
  Xsh(h).b=(inv(Xsh(h).t'*Xsh(h).t)*Xsh(h).t'*y0)';   %b為y的載荷量，就是回歸方程中參數(shù)向量
  Xsh(k);       %總結(jié)構(gòu)體
  x0_=Xsh(h).t*Xsh(h).a';       %x0的估計值
  y0_=Xsh(h).t*Xsh(h).b';       %y0的估計值
  E1=x0-x0_;                    %x殘差值
  F1=y0-y0_;                    %y殘差值

      if norm(E1,inf)>0.000001        %如果殘差陣F1（本程序用E1來代替拉）中元素的絕對值近似為0，則停止疊代
           x0=E1;                %用殘差陣E1代替x0進(jìn)行疊代
           y0=F1;                %用殘差陣F1代替y0進(jìn)行疊代
      else
          flag=1;              %達(dá)到終止條件，用flag==1來表示
          %disp('將計就計')
      end
   
  end
end
h               %顯示出疊代次數(shù)
[n1,n2]=size(x0);         %以下步驟是通過循環(huán)就出結(jié)構(gòu)體中個參數(shù)的每次疊代的值
for i=1:n2           %行，表示x0的變量個數(shù)
    for j=1:h        %列,表示循環(huán)的次數(shù)
    W(j,i)=Xsh(j).w(i);
    end
end
disp('---------X模型權(quán)重為(行向量)----------')
W
for i=1:m4            %行，表示y0的變量個數(shù)
    for j=1:h          %列表示循環(huán)的次數(shù)
    V(j,i)=Xsh(j).v(i);
    end
end
disp('---------Y模型權(quán)重為（行向量)----------')
V 
for i=1:m4           %行，表示y0的變量個數(shù)
    for j=1:h          %列表示循環(huán)的次數(shù)
    B(j,i)=Xsh(j).b(i);
    end
end
disp('---------Y負(fù)荷量為（行向量)----------')
B

for i=1:n2            %行，表示x0的變量個數(shù)
    for j=1:h         %列表示循環(huán)的次數(shù)
    A(j,i)=Xsh(j).a(i);
    end
end
disp('---------X負(fù)荷量為（行向量)----------')
A
disp('-----------標(biāo)準(zhǔn)化方程系數(shù)為(行向量)----------')
z=B'*W               %求出標(biāo)準(zhǔn)化方程系數(shù)
disp('-----------原始系數(shù)(行向量)-------')
D=inv(diag(std_x))*z'*diag(std_y);          %原始變量還原，diag是以行向量為準(zhǔn)生成一個對角矩陣，inv是求逆
D=D'                       %把列向量轉(zhuǎn)為行向量
D0=aver_y-aver_x*inv(diag(std_x))*z'*diag(std_y);     %求出y=d0+d1x+……中的d0,是個行向量
D0=D0'
%綜合指標(biāo)的計算
h =[   41.6959
   24.3422
   18.0158
   15.9461];%各個指標(biāo)貢獻(xiàn)率，已求
ZC=D'*pc