phuigui.m

最小二乘法可以用来处理一组数据, 可以从一组测定的数据中寻求变量之间的依赖关系, 这种函数关系称为经验公式

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的变量个数
aver_x=mean(x);          %x矩阵的均值，是个行向量
std_x=std(x);            %x矩阵的标准差，是个行向量
for j=1:m3                %对x进行标准化
    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的变量个数
aver_y=mean(y);                %y矩阵的均值，是个行向量
std_y=std(y);                   %y矩阵的标准差，是个行向量
for j=1:m4                      %对y进行标准化
    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);          %标准化后x的秩(好象没用上)
h=0;                 %用来统计循环次数的
flag=0;              %是一个只能取0和1的符号，用来控制循环，跳出循环
Xsh=struct('t',[],'u',[],'w',[],'v',[],'a',[],'b',[]);        %是一个结构体，用来记录每次循环并且把循环的值放到同一个变量中
  for i=1:1000
   while flag==0       %当flag==0时候才能执行
       h=h+1;          %用来统计循环次数的
    M=x0'*y0*y0'*x0;   %根据书上的公式
    [V,D]=eig(M);       %求M的特征向量V和特征值D
     p=sqrt(max(max(D)));     %求M特征值里面最大的值
    [n,m]=size(D);      %以下步骤到k，求最大特征值所在的列，然后求出最大特征值所对应的特征向量
    max_d=D(1,1);        %是求出一个矩阵里面那个值最大，并且记录最大值的行和列
    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的权重，把最大的特征值对应的特征向量附给结构体中的w
    
  Xsh(h).v=y0'*x0*Xsh(h).w/p;   %v为y的权重，并写到构体中
  Xsh(h).t=x0*Xsh(h).w;           %t为x的得分向量，并写到构体中
  Xsh(h).u=y0*Xsh(h).v;         %u为y的得分向量，并写到构体中
  Xsh(h).a=(inv(Xsh(h).t'*Xsh(h).t)*Xsh(h).t'*x0)';   %a为x的载荷量，就是回归方程中参数向量
  Xsh(h).b=(inv(Xsh(h).t'*Xsh(h).t)*Xsh(h).t'*y0)';   %b为y的载荷量，就是回归方程中参数向量
  Xsh(k);       %总结构体
  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进行叠代
           y0=F1;                %用残差阵F1代替y0进行叠代
      else
          flag=1;              %达到终止条件，用flag==1来表示
          %disp('将计就计')
      end
   
  end
end
h               %显示出叠代次数
[n1,n2]=size(x0);         %以下步骤是通过循环就出结构体中个参数的每次叠代的值
for i=1:n2           %行，表示x0的变量个数
    for j=1:h        %列,表示循环的次数
    W(j,i)=Xsh(j).w(i);
    end
end
disp('---------X模型权重为(行向量)----------')
W
for i=1:m4            %行，表示y0的变量个数
    for j=1:h          %列表示循环的次数
    V(j,i)=Xsh(j).v(i);
    end
end
disp('---------Y模型权重为（行向量)----------')
V 
for i=1:m4           %行，表示y0的变量个数
    for j=1:h          %列表示循环的次数
    B(j,i)=Xsh(j).b(i);
    end
end
disp('---------Y负荷量为（行向量)----------')
B

for i=1:n2            %行，表示x0的变量个数
    for j=1:h         %列表示循环的次数
    A(j,i)=Xsh(j).a(i);
    end
end
disp('---------X负荷量为（行向量)----------')
A
disp('-----------标准化方程系数为(行向量)----------')
z=B'*W               %求出标准化方程系数
disp('-----------原始系数(行向量)-------')
D=inv(diag(std_x))*z'*diag(std_y);          %原始变量还原，diag是以行向量为准生成一个对角矩阵，inv是求逆
D=D'                       %把列向量转为行向量
D0=aver_y-aver_x*inv(diag(std_x))*z'*diag(std_y);     %求出y=d0+d1x+……中的d0,是个行向量
D0=D0'
%综合指标的计算
h =[   41.6959
   24.3422
   18.0158
   15.9461];%各个指标贡献率，已求
ZC=D'*pc