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