www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@localtruncps/costknot.m
function [res,B,J,YHAT,PS]= costknot(knots,ts,DATA,Wc,c0) % localpspline/COSTB % % [res,B,J,YHAT,PS]= costknot(knots,ps,DATA,Wc) % Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc. if nargin<5 c0=1; end Ns= size(DATA,3); if nargout>1 B= zeros(size(ts,1),Ns); J= cell(1,Ns); PS= J; YHAT= cell(1,Ns); end nk= length(ts.knots); knots= reshape(knots,nk,Ns); if nargin<4 | isempty(Wc) Wc=cell(1,Ns); end res= DATA(:,1); for i= 1:Ns % extract data for current sweep d= DATA{i}; x= d(:,1:end-1); y= d(:,end); ny = length(y); % Residual Calculation ts.knots= knots(:,i); % possibly expand data [ts,OK]= leastsq(ts,x,y,Wc{i}); if ~OK % this occurs because there is no data above the knots X= x2fx(ts,x); X= X(:,Terms(ts)); stat= getstatus(ts); all0= all(X(:,ts.order+1:end)==0); ts= setstatus(ts,find(all0)+ts.order+1,0); [ts,OK]= leastsq(ts,x,y,Wc{i}); ts.xreglinear= set(ts.xreglinear,'status',stat); end yhat= eval(ts,x); % do we have to weight the conditional problem ? wci=[]; if ~isempty(Wc{i}) wci= Wc{i}; end r = y-yhat; % weighting if ~isempty(wci) r= wci*r; end % should be in a sweepset res{i}= r; if nargout > 1 % Model object setup % set up polynomial objects in x-k Ji= jacobian(ts,x); if ~isempty(wci) Ji= wci*Ji; end % what length should these be ? ni YHAT{i}= yhat; PS{i}= ts; % parameter matrix B(:,i)= double(ts); J{i}= Ji; end end if nargout==1 res= double(res)/c0; end