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

    function [L,B,Wc,sigma2,J]= gls(L,X,Y,B,Wc)
%GLS Generalized least squares (with weights)

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


%L.covmodel = covmodel('power','');
% hard-coded for the moment

global STOP


DATA= [X,Y];

if any(cellfun('isempty',Wc))
    DisplayFit(L,'OLS for initial parameter estimates');

    OK= all(isfinite(B),1);

    B= B(:,OK);

    [B,yhat,res,J]= gls_fitB(L,B,DATA);


    for i=1:size(Y,3)
        Xs= X(:,:,i);
        if ~isempty(L.covmodel)
            Wc{i}= choltinv(L.covmodel,yhat{i},X{i});
        end
    end

end


df= (size(Y,1)- size(Y,3)*size(L,1));

[res,J]= gls_costB(B,L,DATA,Wc);
Cost= sqrt(sum(double(res).^2)/df);

% calculate initial model fit, residuals etc
[res,Jr,yhat]= gls_costB(B,L,DATA);

if ~isempty(L.covmodel)
    for i=1:size(Y,3)
        if size(Wc{i},1)~= length(yhat{i})
            Wc{i}= choltinv(L.covmodel,yhat{i},X{i});
        end
    end
end



iter=0;
RUNNING=~isempty(L.covmodel);
STOP=0;

if RUNNING
    str= sprintf('   Iter %9s %9s %9s %9s','sigma','norm(B)','norm(C)','cparam');
    DisplayFit(L,str);
    DisplayFit(L, @() i_CreateDisplayString(0, Cost, B, L.covmodel));
end
while RUNNING && ~STOP
    OldCost= Cost;
    oldB= B;
    oldW= double(L.covmodel);

    iter= iter+1;

    [L.covmodel,Wc]= gls_FitW(L.covmodel,yhat,res,Jr,X);

    % update parameter estimates
    % run lsqnonlin
    [B,yhat,res,J]= gls_fitB(L,B,DATA,Wc);
    Cost= sqrt(sum(double(res).^2)/df);

    % need raw residuals for next step
    [res,Jr,yhat]= gls_costB(B,L,DATA);

    DisplayFit(L, @() i_CreateDisplayString(iter, Cost, B, L.covmodel));

    deltaC= norm(Cost-OldCost);
    deltaB= norm(B-oldB);
    deltaW= norm(oldW-double(L.covmodel));
    if (deltaC/Cost < L.FitOptions.TolSigma) || ...
            (deltaB < L.FitOptions.TolParams) || ...
            (deltaW < L.FitOptions.TolCov) || ...
            (iter >= L.FitOptions.MaxIter)
        RUNNING=0;
    end
end

sigma2= Cost.^2;


function str = i_CreateDisplayString(Iter, Cost, B, cm)
str = [sprintf('GLS:%3d ',Iter), ...
    sprintf('%9.5g ', Cost, i_SafeNorm(B), i_SafeNorm(double(cm)), double(cm))];


%--------------------------------------------------------------------------
function n = i_SafeNorm( x )
% NORM protected against INF and NAN
%
% If x constrain Inf or NaN than norm x errors with "NaN or Inf prevents
% convergence". As we know how to respond to these cases, we look from them
% and return an appropriate value.

if any( isnan( x(:) ) ),
    % If there are NANs then the norm must also be NAN
    n = NaN;
elseif any( isinf( x(:) ) ),
    % If there are no NANs but there are INFs then any sensible norm should
    % return INF.
    n = Inf;
else
    % if no NAN and no INF, we should safely be able to call NORM
    n = norm( x );
end