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