stepwise逐步回归

clear clc data=xlsread('week10_eg9_data'); [n,m]=size(data); y=data(:,m); z=(y.^0.07-1)/0.07; X=cell(m-1,1); for i=1:m-1 X{i}=data(:,i); end aE=0.05; aD=0.05; %-------------------------------------------------------------------------- fprintf('Statistics for Entry:Step 1\n'); %-------------------------------------------------------------------------- p=2; b_2=cell(4,1); s_2=ones(4,1); F=ones(4,1); Prob=ones(4,1); for j=1:m-1 x=[ones(n,1),X{j}]; [b_2{j},bint,r,rint,stats]=regress(z,x,0.05); zHat=x*b_2{j}; aver_z=mean(z); sigma_2=stats(4); C=inv(x'*x); s_2(j)=sqrt(sigma_2*C(2,2)); F(j)=(b_2{j}(2)/s_2(j))^2; Prob(j)=1-fcdf(F(j),1,n-2); end fprintf('Variable|F |Prob>F \n'); V={' X1 ',' X2 ',' X3 ',' X4 '}; for i=1:4 fprintf('%s|%-7.4f|%-10.4f',V{i},F(i),Prob(i)); fprintf('\n'); end f=1; for k=1:3 if F(k+1)>F(f)&&Prob(k+1)<aE f=k+1; end end if Prob(f)<aE u=1; else u=0; end if u==1 fprintf('添加X%d\n',f); else fprintf('不添加\n'); end %-------------------------------------------------------------------------- x=[ones(n,1),X{f}]; [b_f,bint,r,rint,stats]=regress(z,x,0.05); PE=b_f; C=inv(x'*x); SE=ones(2,1); F_f=ones(2,1); Prob_f=ones(2,1); sigma_2=stats(4); for l=1:2 SE(l)=sqrt(sigma_2*C(l,l)); F_f(l)=(b_f(l)/SE(l))^2; Prob_f(l)=1-fcdf(F_f(l),1,n-2); end fprintf('回归模型参数：\n') fprintf(' Variable |Parameter Estimate |Standard Error |F |Prob>F \n'); fprintf(' INTERCEP |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(1),SE(1),F_f(1),Prob_f(1)); fprintf(' X4 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(2),SE(2),F_f(2),Prob_f(2)); if Prob_f(2)>aD fprintf('出错，aE应小于或等于aD'); end %-------------------------------------------------------------------------- fprintf('Statistics for Entry:Step 2\n'); %-------------------------------------------------------------------------- p=3; b_3=cell(3,1); s_3=ones(3,1); F=ones(3,1); Prob=ones(3,1); for j=1:m-2 x=[ones(n,1),X{j},X{4}]; [b_3{j},bint,r,rint,stats]=regress(z,x,0.05); zHat=x*b_3{j}; aver_z=mean(z); sigma_2=stats(4); C=inv(x'*x); s_3(j)=sqrt(sigma_2*C(2,2)); F(j)=(b_3{j}(2)/s_3(j))^2; Prob(j)=1-fcdf(F(j),1,n-3); end fprintf('Variable|F |Prob>F \n'); V={' X1 ',' X2 ',' X3 '}; for i=1:3 fprintf('%s|%-7.4f|%-10.4f',V{i},F(i),Prob(i)); fprintf('\n'); end f=1; for k=1:2 if F(k+1)>F(f)&&Prob(k+1)<aE f=k+1; end end if Prob(f)<aE u=1; else u=0; end if u==1 fprintf('添加X%d\n',f); else fprintf('不添加\n'); end %-------------------------------------------------------------------------- x=[ones(n,1),X{f},X{4}]; [b_f,bint,r,rint,stats]=regress(z,x,0.05); PE=b_f; C=inv(x'*x); SE=ones(3,1); F_f=ones(3,1); Prob_f=ones(3,1); sigma_2=stats(4); for l=1:3 SE(l)=sqrt(sigma_2*C(l,l)); F_f(l)=(b_f(l)/SE(l))^2; Prob_f(l)=1-fcdf(F_f(l),1,n-3); end fprintf('回归模型参数：\n'); fprintf(' Variable |Parameter Estimate |Standard Error |F |Prob>F \n'); fprintf(' INTERCEP |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(1),SE(1),F_f(1),Prob_f(1)); fprintf(' X3 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(2),SE(2),F_f(2),Prob_f(2)); fprintf(' X4 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(3),SE(3),F_f(3),Prob_f(3)); for i=2:3 if Prob_f(i)>aD fprintf('剔除X%d\n',i-1); end end %-------------------------------------------------------------------------- fprintf('Statistics for Entry:Step 3\n'); %-------------------------------------------------------------------------- p=4; b_4=cell(2,1); s_4=ones(2,1); F=ones(2,1); Prob=ones(2,1); for j=1:m-3 x=[ones(n,1),X{j},X{3},X{4}]; [b_4{j},bint,r,rint,stats]=regress(z,x,0.05); zHat=x*b_4{j}; aver_z=mean(z); sigma_2=stats(4); C=inv(x'*x); s_4(j)=sqrt(sigma_2*C(2,2)); F(j)=(b_4{j}(2)/s_4(j))^2; Prob(j)=1-fcdf(F(j),1,n-4); end fprintf('Variable|F |Prob>F \n'); V={' X1 ',' X2 '}; for i=1:2 fprintf('%s|%-7.4f|%-10.4f',V{i},F(i),Prob(i)); fprintf('\n'); end f=1; k=1; if F(k+1)>F(f)&&Prob(k+1)<aE f=k+1; end if Prob(f)<aE u=1; else u=0; end if u==1 fprintf('添加X%d\n',f); else fprintf('不添加\n'); end %-------------------------------------------------------------------------- x=[ones(n,1),X{f},X{3},X{4}]; [b_f,bint,r,rint,stats]=regress(z,x,0.05); PE=b_f; C=inv(x'*x); SE=ones(4,1); F_f=ones(4,1); Prob_f=ones(4,1); sigma_2=stats(4); for l=1:4 SE(l)=sqrt(sigma_2*C(l,l)); F_f(l)=(b_f(l)/SE(l))^2; Prob_f(l)=1-fcdf(F_f(l),1,n-4); end fprintf('回归模型参数：\n'); fprintf(' Variable |Parameter Estimate |Standard Error |F |Prob>F \n'); fprintf(' INTERCEP |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(1),SE(1),F_f(1),Prob_f(1)); fprintf(' X2 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(2),SE(2),F_f(2),Prob_f(2)); fprintf(' X3 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(3),SE(3),F_f(3),Prob_f(3)); fprintf(' X4 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(4),SE(4),F_f(4),Prob_f(4)); for i=2:4 if Prob_f(i)>aD fprintf('剔除X%d\n',i-1); end end %-------------------------------------------------------------------------- fprintf('Statistics for Entry:Step 4\n'); %-------------------------------------------------------------------------- p=5; x=[ones(n,1),X{1},X{2},X{3},X{4}]; [b_5,bint,r,rint,stats]=regress(z,x,0.05); zHat=x*b_5; aver_z=mean(z); sigma_2=stats(4); C=inv(x'*x); s_5=sqrt(sigma_2*C(2,2)); F=(b_5(2)/s_5)^2; Prob=1-fcdf(F,1,n-5); fprintf('Variable|F |Prob>F \n'); fprintf(' X1 |%-8.4f|%-10.4f',F,Prob); fprintf('\n'); if Prob<aE f=1; else f=0; end if f==1 fprintf('添加X1\n'); else fprintf('不添加\n'); end %-------------------------------------------------------------------------- x=[ones(n,1),X{1},X{2},X{3},X{4}]; [b_f,bint,r,rint,stats]=regress(z,x,0.05); PE=b_f; C=inv(x'*x); SE=ones(5,1); F_f=ones(5,1); Prob_f=ones(5,1); sigma_2=stats(4); for l=1:5 SE(l)=sqrt(sigma_2*C(l,l)); F_f(l)=(b_f(l)/SE(l))^2; Prob_f(l)=1-fcdf(F_f(l),1,n-5); end fprintf('回归模型参数：\n'); fprintf(' Variable |Parameter Estimate |Standard Error |F |Prob>F \n'); fprintf(' INTERCEP |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(1),SE(1),F_f(1),Prob_f(1)); fprintf(' X1 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(2),SE(2),F_f(2),Prob_f(2)); fprintf(' X2 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(3),SE(3),F_f(3),Prob_f(3)); fprintf(' X3 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(4),SE(4),F_f(4),Prob_f(4)); fprintf(' X4 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(5),SE(5),F_f(5),Prob_f(5)); for i=2:5 if Prob_f(i)>aD fprintf('剔除X%d\n',i-1); end end %-------------------------------------------------------------------------- x=[ones(n,1),X{1},X{2},X{3}]; [b_f,bint,r,rint,stats]=regress(z,x,0.05); PE=b_f; C=inv(x'*x); SE=ones(4,1); F_f=ones(4,1); Prob_f=ones(4,1); for l=1:4 SE(l)=sqrt(sigma_2*C(l,l)); F_f(l)=(b_f(l)/SE(l))^2; Prob_f(l)=1-fcdf(F_f(l),1,n-4); end fprintf('回归模型参数：\n'); fprintf(' Variable |Parameter Estimate |Standard Error |F |Prob>F \n'); fprintf(' INTERCEP |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(1),SE(1),F_f(1),Prob_f(1)); fprintf(' X1 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(2),SE(2),F_f(2),Prob_f(2)); fprintf(' X2 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(3),SE(3),F_f(3),Prob_f(3)); fprintf(' X3 |%-20.4f|%-16.4f|%-9.4f|%-7.4f\n',b_f(4),SE(4),F_f(4),Prob_f(4)); %-------------------------------------------------------------------------- fprintf('Statistics for Entry:Step 5\n'); %-------------------------------------------------------------------------- p=5; x=[ones(n,1),X{1},X{2},X{3},X{4}]; [b_5,bint,r,rint,stats]=regress(