dyzbhg.m

逐步回归算法程序

function dyzbhg(xy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                多元逐步回归分析                    %   
%        作者：唐世星           2006.11.22           % 
%   xy为待输入的原始数据，按照先x后y按列排列的数组   %
%   如：x1 x2 x3 x4 y等等                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%计算离差阵R(m,m)
[n,m]=size(xy);

xy_aver=mean(xy);%求均值
disp(['均值为：' num2str(xy_aver)])
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=');
    if isempty(F1)==1
        F1=1;
    end
    if isempty(F2)==1
        F2=1;
    end
    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        
        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;
        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)])

    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(['各自变量的方差贡献:' num2str(V)])
    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