bhg.rar_matlab逐步回归_stepwise regression_逐步回归_逐步回归 matlab_逐步回归算法

function dyzbhg(xy) %多元逐步回归分析 %作者：唐世星 %2006.11.20 %xy为待输入的原始数据，按照先x后y按列排列的数组 %如：x1 x2 x3 x4 y等等 %clc;%clear all; %计算离差阵R(m,m) [n,m]=size(xy); %F1=0;F2=0; %disp('均值为:') xy_aver=mean(xy)%求均值 for i=1:m for j=1:i R(i,j)=0; for k=1:n R(i,j)=R(i,j)+(xy(k,i)-xy_aver(i))*(xy(k,j)-xy_aver(j)); end R(j,i)=R(i,j); end SR(i)=sqrt(R(i,i));%计算对角线元素的平方根 end %disp('************ Deviation Matrix & Value of SR (离差阵R&SR) ***********') %输出离差阵R,及SR %[R SR'] %计算相关系数R(m,m) for i=1:m for j=1:i R(i,j)=R(i,j)/(SR(i)*SR(j)); R(j,i)=R(i,j); end end %disp('********** Correlation Coefficient Matrix (相关系数阵R) **********')%输出相关系数阵R %R flag=1;%是否重复进行逐步回归的标志 while(flag) disp('******** Stepwise Regression Analysis Start *************') F1=input('剔除门坎值:F1='); F2=input('引入门坎值:F2='); S=0;%计算步数 L=0;%引入方程的自变量个数 FQ=n-1;%残差平方和的自由度 disp('************** Discriminant Value of Contribution V **************') Imin(1)=0;Imax=1:m-1;%定义已引入(最小)和未引入(最大)变量的序号 inn=0;outt=0;%引入和剔除的变量的顺序号 while(1) % pause VN=1E+08;%已引入方程的自变量贡献的最小值 VX=0;%未引入方程的自变量贡献的最大值 IN=0;%贡献最小的已引入的自变量序号 IX=0;%贡献最大的未引入的自变量序号 S=S+1; disp(['--------- step = ' int2str(S) '------------'])%输出步骤数 for i=1:m-1 if R(i,i)<1E-08 continue end % disp(['VMAX=' int2str(VX) '; IMAX=' int2str(IX)]) %输出Vmax=VX;Imax=IX; V(i)=R(m,i)^2/R(i,i);%计算已引入的变量的方差贡献 if V(i)>=0 if V(i)>VX %寻找未引入变量方差贡献的最大值 for in=1:length(Imax) if i==Imax(in) VX=V(i);IX=i; end end end end if abs(V(i))<VN %寻找已引入变量方差贡献的最小值 for out=1:length(Imin) if i==Imin(out) VN=abs(V(i));IN=i; end end end %disp(['方差贡献:V=' num2str(V(i)) 'VX=' num2str(VX) 'IX=' int2str(IX) 'VN=' num2str(VN) 'IN=' int2str(IN)]) end % Imax(inn+1)=IX;inn=inn+1; t=find(Imax==IX); Imax(t)=[]; disp(['******** 方差贡献V **********' num2str(V)]) disp(['VMAX=' num2str(VX) '; IMAX=' int2str(IX)]) %输出Vmax=VX;Imax=IX; % disp(['VMIN=' num2str(VN) '; IMIN=' int2str(IN)]) %输出Vmin=VN;Imin=IN; if S==1 disp(['S=' int2str(S)]) %输出S=1 else disp(['VMIN=' num2str(VN) '; IMIN=' int2str(IN)]) %输出Vmin=VN;Imin=IN; end if S==1%||S==2||S==3 FE=VX*(n-L-2)/(R(m,m)-VX); disp(['FE=' num2str(FE)]) %输出 FE if FE<F1 if L~=0 disp('Neither Delete Out Nor Select In!') else disp('May Be Smaller F1 And F2') disp('The Stepwise Regression Analysis End!') break;%程序结束 end else L=L+1;FQ=FQ-1;K=IX; disp(['X' int2str(K) ' Be Selected In']) Imin(outt+1)=IX;outt=outt+1; disp(['L = ' int2str(L) ]) R=xiaoqu(R,K) %调用子函数,执行消去变换 if L~=m-1 continue; end disp('Already Selecting End') break; end else %计算剔除变量的F检验值 FT=VN*(n-L-1)/R(m,m); disp(['剔除变量的F检验值' num2str(FT)]) if FT>=F2 FE=VX*(n-L-2)/(R(m,m)-VX); disp(['***FE=' num2str(FE)]) %输出 FE if FE<F1 if L~=0 disp('Neither Delete Out Nor Select In!') disp('The Stepwise Regression Analysis End!') break;%程序结束 else disp('May Be Smaller F1 And F2') disp('The Stepwise Regression Analysis End!') break;%程序结束 end else L=L+1;FQ=FQ-1;K=IX; disp(['X' int2str(K) ' Be Selected In']) disp(['L = ' int2str(L) ]) Imin(outt+1)=IX;outt=outt+1; R=xiaoqu(R,K) %调用子函数,执行消去变换 if L~=m-1 continue; end disp('Already Selecting End') break; end else L=L-1;FQ=FQ+1;K=IN; disp(['X' int2str(K) ' Be Deleted Out']) disp(['L = ' int2str(L) ' (No. of Variable Selected)']) R=xiaoqu(R,K) %调用子函数 continue end end end %输出相应的计算结果 for i=1:m-1 kk=R(i,m)*R(m,i); if kk<0 B(i)=R(i,m)*SR(m)/SR(i); else B(i)=0; end end B0=xy_aver(m); for i=1:m-1 B0=B0-B(i)*xy_aver(i); end disp(['回归系数为:' num2str(B0) ' ' num2str(B)]) disp('回归方程为：') disp(['Y=' num2str(B0)]) for i=1:m-1 if B(i)~=0 if B(i)>0 disp(['+' num2str(B(i)) 'X' int2str(i)]); else disp([num2str(B(i)) 'X' int2str(i)]); end end end Q=SR(m)^2*R(m,m);%残差平方和 disp(['Sum of SQuares of Residual Error(残差平方和) Q = ' num2str(Q)]) S=SR(m)*sqrt(R(m,m)/FQ);%剩余标准差 disp(['Standard Deviation(剩余标准差，即模型误差的均方根) S = ' num2str(S)]) RR=sqrt(1-R(m,m));%复相关系数 disp(['Multiple Correlation Coefficient(复相关系数) R = ' num2str(RR)]) FF=FQ*(1-R(m,m))/(L*R(m,m));%回归方程显著性检验的F值 disp(['F Value for Test of Regression(回归方程显著性检验，即回归模型的统计量) F = ' num2str(FF)]) %F=SH*(m-n-1)/(SX*n);%F-统计量 %PROB = 1 - fcdf(FF,m,n-length(Imin)-1)%与统计量F对应的概率P for i=1:m-1 CC=R(i,i)*R(m,m); T(i)=R(i,m)/sqrt(CC/FQ);%各回归系数的t检验值 R1(i)=R(i,m)/sqrt(CC+R(i,m)^2);%各自变量的偏相关系数 end disp(['t Test Value of Argument(各回归系数的t检验值):' num2str(T)]) disp(['Partial Corre.Coeffi.Ofargu.(各自变量的偏相关系数):' num2str(R1)]) %for i=1:n % y(i)=B0; % for j=1:m-1 % y(i)=y(i)+B(j)*xy(i,j); % end % E(i)=xy(i,m)-y(i); % PC(i)=E(i)/xy(1,m)*100; %end %x=1:length(xy); %disp(' No. 回归值 误差 误差百分比%') %[x' y' E' PC'] flag=input('是否重新进行逐步回归分析(1:是;0:否)：'); end