www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregunispline/costknot.m
function [s,g]= costknot(knots,m,X,Y,nk,c0,alpha); %COSTKNOT % Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc. persistent FAC % First reconstruct the knot sequence. h=diff([-1;knots(:);1]); h(h<sqrt(eps))=sqrt(eps); penalty= 1-alpha*sum(log((nk+1)/2*h)); m3= set(m.mv3xspline,'knots',knots); [m3,OK]= leastsq(m3,X,Y); yhat= eval(m3,X); if OK r= Y-yhat; else r= zeros(size(Y)); r(:)= 10/sqrt(length(Y)); end sseJ=penalty*r; s= sum(sseJ.^2); if ~isempty(c0) s= s/c0; else FAC=[]; c0=1; end if nargout>1 m.mv3xspline=m3; % numerical jacobian thresh= max(abs(knots)*sqrt(eps),sqrt(eps)); [Jy,FAC]= numjac('fbnleval',m,knots,yhat,thresh,FAC,0,[],[],X); % penalty jacobian dhdk= [diag(ones(nk,1));zeros(1,nk)] - [zeros(1,nk);diag(ones(nk,1))]; dk= -alpha* ((1./h)'*dhdk); % penalised jacobian J= del ( pen*yhat) /delk J= -penalty*Jy + r*dk; g= 2*sseJ'*J/c0; end