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