声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1989|回复: 1

[应用数学] 逐步回归算法程序应用问题

[复制链接]
发表于 2015-9-17 09:59 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我采用下面的逐步回归算法,为什么只运行一步?只选择出一个变量之后就跳出循环?请大神赐教



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))
                 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               
                 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                    
                     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

%第二个函数 xiaoqu.m
function R=xiaoqu(R,k)
%多元逐步回归分析
%对R作消去变换
%作者:唐世星
%2006.11.20
G=1/R(k,k);
m=length(R);
for i=1:m
     for j=1:m
         if i~=k&j~=k
             R(i,j)=R(i,j)-R(i,k)*R(k,j)*G;
         end
     end
end
for i=1:m
     R(k,i)=R(k,i)*G;
     R(i,k)=-R(i,k)*G;
end
R(k,k)=G;
回复
分享到:

使用道具 举报

发表于 2015-9-17 10:55 | 显示全部楼层
好长的代码,有什么提示信息吗?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 11:40 , Processed in 0.061297 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表