www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xreglinear/stepwise.m

    function [m,OK,Stats,B]=stepwise(m,StepTerms)
%STEPWISE Perform stepwise term reduction on the model
%
%  [M, OK, PRESS, B] = STEPWISE(M, STEPTERMS) performs a stepwise operation
%  on the model.  If STEPTERMS is a list of term indices, those terms will
%  be toggled in/out of the current model state in as optimal a fashion as
%  possible.  If STEPTERMS is a logical vector the same size as the number
%  of terms in the model, the model will have all of its terms "stepped"
%  out as indicated.  Terms that have been set as being forced in or out of
%  the model will retain these settings.
%
%  The model store values (eg a new Q and R) and the model coefficients are
%  all updated during a stepwise operation.

%  Copyright 2000-2008 The MathWorks, Inc. and Ford Global Technologies, Inc.


OK= 1;
Stats=[];
B=[];
Nobs = size(m.Store.y,1);
FX = m.Store.X;
x = m.Store.D;
if nargin >= 2
    StepTerms= StepTerms(:);
    
    % Flag that indicates that a full QR decomposition needs to be
    % performed
    FULL_REFRESH = false;
    
    if islogical(StepTerms) && length(StepTerms) == length(m.TermsOut)
        % StepTerms specifies terms to be excluded
        StepTerms = m.TermStatus==3 & StepTerms;
        m.TermsOut = StepTerms;
        StepTerms = find(StepTerms);
        FULL_REFRESH = true;
    else
        % Toggle state of specified terms
        m.TermsOut(StepTerms)= ~m.TermsOut(StepTerms);
        
        if islogical(StepTerms)
            % An edge usage case is a logical vector not the length of the
            % term vector.
            StepTerms = find(StepTerms);
        end
        
        if numel(StepTerms)~=1 || any(~m.TermsOut(StepTerms))
            FULL_REFRESH = true;
        end
    end


    m.TermsOut(m.TermStatus==1)= false;
    m.TermsOut(m.TermStatus==2)= true;
    FX(:,m.TermsOut)=[];

    if ~FULL_REFRESH && Nobs>100
        % the problem is large enough and we are just removing a single term

        % find term being deleted
        pos= sum(~m.TermsOut(1:StepTerms-1))+1;
        [Q,R,df,ri]= qrdelete(m,pos);
    else
        % need to compute the new Q and R
        [Q,R,OK,df,ri]= qrdecomp(m,x,[],FX);
    end

    if ~OK
        return
    end


    
    % Store qr for use by stats and further stepwise
    m.Store.Q=Q;
    m.Store.R=R;
    % pad the y (only has an effect for ridge)
    y = [m.Store.y; zeros(size(Q,1) - Nobs,1)];
    if ~all(m.TermsOut)
        H = sum(Q.*Q,2);
        %truncate H for ridge
        % residuals
        r= y- Q*(Q'*y);
    else
        H= zeros(size(Q,1),1);
        r= y;
    end
    m.Store.H= H(1:Nobs,1);
    % truncate for ridge
    r = r(1:Nobs, 1);
    % calculate MSE
    if df>0
        mse = sum(r.*r)/df;
    else
        mse= 0;
    end

    if OK
        % recalc coeffs
        m.Beta = calcBeta(m);
        % store variance info
        m= var(m,ri*sqrt(mse),mse,df);
    else
        m= var(m,[],0,Inf);
    end
else
    Q= m.Store.Q;
end

if nargout>2
    % PRESS
    Stats =press(m);
end
if nargout>3
    % find next terms
    FX = m.Store.X;

    beta  = m.Beta;
    NextPress= NaN(size(beta));

    % pad the y (only has an effect for ridge)
    y = [m.Store.y; zeros(size(Q,1) - Nobs,1)];

    % R2x = collinear(m);

    [ri,mse,df] = var(m);    % df(SSE)= n-p

    df = repmat(df,size(beta));
    Bstd = zeros(size(beta));
    if ~all(m.TermsOut)
        %                     sqrt(diag(ri*ri'));
        Bstd(~m.TermsOut,:) = sqrt(sum(ri.*ri,2));
    end
    
    TempTermsOut= m.TermsOut;
    % Add/Delete Terms to Model one at a time
    TermsIn= Terms(m);
    % determine cases when m.TermsOut needs to be updated before qrdecomp
    % is called
	UpdateTermsRequired =  true;
    
    for i = find(m.TermStatus==3)'
        % toggle term 'i' in or out of model
        TermsIn(i) = ~TermsIn(i);
        if UpdateTermsRequired
            % need to change termsout so qrdecomp does not have dimension
            % problems
            % need to change termsout so diag(lam)-FX'*FX works
            m.TermsOut(i)= ~TermsIn(i);
        end
        if ~TermsIn(i) && Nobs>100
            % use qrdelete for 'large' problems
            j= sum(TermsIn(1:i-1))+1;
            [q,r,dfi]= qrdelete(m,j);
            qrok= true;
        else
            % recalculate qr is quicker for small problems
            Xt= FX(:,TermsIn);
            [q,r,qrok,dfi]= qrdecomp(m,x,[],Xt);
        end
        % Check for rank deficient cases
        if qrok
            if ~isempty(q);
                H= sum(q.*q,2);
                % truncate H for ridge
                H = H(1:Nobs,1);
            else
                H= 0;
            end

            res= y - q*(q'*y);
            % truncate res for ridge
            res = res(1:Nobs,1);
            if any(H==1)
                NextPress(i)= NaN;
            else
                NextPress(i) = sum( (res./(1-H)).^2 );
            end
            if TermsIn(i)
                % i.e. term was out and we are adding it back

                % Calculate STD of coeff.

                Xpos= find(TermsIn)==i;
                if Xpos>100
                    b = r(Xpos:end,Xpos:end)\(q(:,Xpos:end)'*y);
                    beta(i) = b(1);
                else
                    b = r\(q'*y);
                    beta(i) = b(Xpos);
                end

                df(i) = dfi;

                ei= zeros(1,size(r,1));
                ei(Xpos)= 1;
                if strcmp(m.qr,'ols')
                    r1= ei/r;
                else
                    r1= ((ei/r)/r')*Xt';
                end

                if df(i)
                    s2 = sum(res.*res)/df(i);
                    Bstd(i) = sqrt( s2 * r1*r1' );
                else
                    Bstd(i)= NaN;
                end

            end
        else
            % Rank deficient X matrix
            beta(i)     = NaN;
            Bstd(i)     = NaN;
            NextPress(i) = NaN;
        end
        % return term 'i' to starting state
        TermsIn(i) = ~TermsIn(i);
        if UpdateTermsRequired
            m.TermsOut(i)= ~TermsIn(i);
        end
    end

    B= [beta,Bstd,df,NextPress];
    m.TermsOut= TempTermsOut;
end