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

    function [B,yhat,res,J,OKall]= gls_fitB(poly,B,DATA,Wc)
% POLYNOM/GLS_FITB least-squares estimation of localpspline
%
% [B,res,J,yhat]= gls_fitB(ps,B,DATA,Wc)
%   ps    localpspline object
%   B     initial parameter matrix (cols= sweeps)
%   DATA  sweepset of data to fit [X,Y]
%   Wc    optional weights Wc'*Wc= inv(covmatrix)

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



Ns= size(DATA,3);

if isa(DATA,'sweepset')
   CDATA= cell(DATA);
elseif ~isa(DATA,'cell') 
   CDATA= {DATA};
else
   CDATA= DATA;
end
if nargin < 4 | isempty(Wc)
   Wc=cell(Ns,1);
end

% setup Jacobian pattern
Jcell= cell(1,Ns);
ts= cellfun('size',CDATA,1);
for i=1:Ns;
   Jcell{i}= ones(ts(i),1);
end
JPattern= spblkdiag(Jcell{:});

p0 = datum(poly,0);
OKall= zeros(1,Ns);

for i=1:Ns
   d= CDATA{i};
   Xs= d(:,1:end-1);
   Ys= d(:,end);
   
   tp=0;
   [poly,OK]= leastsq(p0,Xs,Ys,Wc{i});
   
   
   p= double(poly);
   p= [zeros(size(B,1)-length(p),1); p];
   B(:,i)= p;
   if nargout>1
      % this makes sure that the polynomial is the right size
      % and sets the datum
      poly= update(poly,p);
      [res{i},J{i},yhat{i}]= lsqcost(poly,Xs-datum(poly),Ys,Wc{i});
   end
   OKall(i)= OK;
end
% not sure if I want to do this in here or whether this is best post GLS 

if isa(DATA,'sweepset')
   res= cell2sweeps(DATA(:,1),res);
else
   res= cat(1,res{:});
end