www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregrbf/iteraterols.m
function [om,OK]=iteraterols(m) % XREGRBF/ITERATEROLS iterates rols to select reguarization parameter % % Lambda selection algorithm for RBFs. Implemented as an xregoptmgr. This % is normally a sub xregoptmgr (of e.g trialwidths) % % Copyright 2000-2014 The MathWorks, Inc. and Ford Global Technologies, Inc. omNestCent = rols(m); % make the xregoptmgr for rols omNestCent = setAltMgrs(omNestCent,{@rols});% give a choice of center selection algorithms om= contextimplementation(xregoptmgr,m,@i_iteraterols,[],'IterateRols',@iteraterols); om = setAltMgrs(om,{@iterateridge,@iteraterols,@stepitrols}); om = AddOption(om,'CenterSelectionAlg',omNestCent,'xregoptmgr', 'Center selection algorithm',2);% % fit parameters om = AddOption(om,'MaxNumIter',10,{'int',[1,100]}, 'Maximum number of iterations');%number of iterations om= AddOption(om,'Tolerance',0.005,{'numeric',[eps 1]}, 'Minimum change in log10(GCV)');% stopping criterion om= AddOption(om,'nTestValues',5,{'int',[0,100]},'Number of initial test values for lambda');%number of test values of lambda om= AddOption(om,'CheapMode',1,'boolean', 'Do not reselect centers for new width'); om= AddOption(om,'PlotFlag',0,'boolean', 'Display');%plot flag om= AddOption(om,'cost',Inf,'numeric',[],false);% field to store cost (GCV), not gui-settable OK = 1; function [m,cost,OK,varargout]=i_iteraterols(m,om,~,x,y,varargin) % iteraterols will iterate the lambda parameter to minimise % a model critierion such as GCV or MML % a different set of centers will be chosen for each value of lambda % Inputs: % m - rbf object with centers % x - matrix of data points % y - target values % Outputs % m - new rbf object varargout = cell(1,max(0,nargout-3)); %parameters %%%%%%%%%%%%%%%%%%%%% plotflag = get(om,'PlotFlag'); % whether to plot or not - set to 0 if using uOptRols or OptRols with this routine ntest = get(om,'nTestValues'); % number of test values of lambda niter = get(om,'MaxNumIter'); % number of iterations tol = get(om,'Tolerance'); cheapflag = get(om,'CheapMode'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% centalg = get(om,'CenterSelectionAlg'); lambda= get(m,'lambda'); m2= m; m2.centers= x; PHI= x2fx(m2,x); if cheapflag try CenterTol= get(centalg,'CenterTol'); catch CenterTol= 0; end if CenterTol>0 candcenters = curvefitlib.internal.uniqueWithinTol(x,CenterTol*ones(1,size(x,2))); m2.centers= x; PHI= x2fx(m2,candcenters); end J = CalcJacob(m,x); [Q,R] = xregqr(J); n = size(R,1); W= Q*sparse(1:n,1:n,diag(R),n,n); GCV= log10GCV(m,x,y); else % re-choose the centers for the new width [m,GCV] = run(centalg,m,[],x,y, PHI); end lambdaRange = [1e-10,100]; if ntest > 0 % test several lambda values to determine a rough start-point (based on GCV) lamx = logspace(log10(lambdaRange(1)),log10(lambdaRange(2)),ntest); log10GCVy = zeros(1,length(lamx)); for i = 1:length(lamx); m = set(m,'lambda',lamx(i)); if cheapflag % calculate GCV for current centers log10GCVy(i)= log10GCV(m,x,y); else % rechoose the centers [~,log10GCVy(i)] = run(centalg,m,[],x,y,PHI); end end [~,index] = min(log10GCVy); newlam = lamx(index); % best lambda from the vector else newlam = max(get(m,'lambda'),1e-12); end oldGCV = Inf; i=1; diff = Inf; while i <= niter && diff> tol && newlam >= lambdaRange(1) && newlam <= lambdaRange(2) % search for new lambda lambda(i) = newlam; m = set(m,'lambda',newlam); if cheapflag % calculate GCV for current centers [cost,newlam] = iCalcGCV(W,y,newlam); else % rechoose the centers [m,cost,~,newlam] = run(centalg,m,[],x,y,PHI); end GCV(i)= cost; % [GCV(i), newlam] = i_calcGCV(W,y,lambda(i)); diff = abs(oldGCV-GCV(i)); oldGCV = GCV(i); i = i+1; end [bestGCV,ind] = min(GCV); bestlambda =lambda(ind); % update the centers for the new lambda [m,cost] = run(centalg,m,[],x,y,PHI); OK =1; if plotflag == 1 plothand =figure('MenuBar','none',... 'ToolBar','none',... 'NumberTitle','off',... 'Name','Results of IterateRols',... 'Tag','mbcfitplot',... 'Color',get(0,'DefaultUicontrolBackgroundColor')); a = axes('Parent',plothand,... 'XScale','log',... 'ButtonDownFcn', @(h,evt) mv_zoom(h)); line('Parent',a,'XData',lambda,'YData',(GCV),'Color','r','Marker','o','LineStyle','-'); line('Parent',a,'XData',lambda(1),'YData',(GCV(1)),'Color','k','Marker','o','MarkerFaceColor','k','LineStyle','none'); line('Parent',a,'XData',bestlambda,'YData',(bestGCV),'Color','g','Marker','*','LineStyle','none'); mbctitle(a,['Results of IterateRols for Width = ' num2str(m.width)]); mbcxlabel(a,'lambda'); mbcylabel(a,'log10(GCV)'); end function [GCV,newlam] = iCalcGCV(W,y,lam) dw= sum(W.*W,1); N= length(y); ee = sum(y.*y) - sum( (((y'*W).^2)./(lam +dw)).*((2*lam + dw)./(lam + dw))); % calculate df here rather than in qrdecomp df = N - sum(dw./(dw+lam'));% N - gamma nu = sum(dw./(dw+lam).^2); wAw = sum( ((y'*W).^2) ./ ((dw+lam).^3) ); if abs(wAw) > eps && abs(df) > sqrt(eps) GCV = N*ee/(df)^2; GCV = log10( GCV ); newlam = (nu/df)*ee/wAw; else GCV = Inf; newlam = 0; end return