www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xreglinear/prune.m
function [om,ok]= prune(m) %PRUNE Prune least squares model from last term % % [OM, OK] = PRUNE(M) creates an optimisation manager that is set up to % prune steppable terms from the model starting at the final term. The % model configuration with the best criteria value is returned. % Copyright 2000-2010 The MathWorks, Inc. and Ford Global Technologies, Inc. om= contextimplementation(xregoptmgr,m,@i_prune,[],'Prune',@prune); [s,list]= summary(m); str= sprintf('%s|',list{:}); om= AddOption(om,'Criteria',list{1},str(1:end-1)); % Summary Statistic Criteria om= AddOption(om,'MinTerms',0,{'int',[0 Inf]},'Minimum number of terms');% minimum number of terms in model om= AddOption(om,'Tolerance',0,{'numeric',[0 1000]});% tolerance for prune om= AddOption(om,'IncludeAll',0,'boolean','Include all terms before prune');% include all terms om= AddOption(om,'isinitialised',0,'boolean',[],false);% flag to skip initial fit om= AddOption(om,'Display',0,'boolean',[]);% plot results at end of fit ok=1; %-------------------------------------------------------------------------- function [m,cost,OK,s]= i_prune(m,om,x0,x,y,varargin) if ~get(om,'isInitialised') m= leastsq(m,x,y); end Criteria= get(om,'Criteria'); if get(om,'IncludeAll') % All Terms included in model InModel = false(size(m)); m = stepwise(m,InModel); end MinTerms = get(om,'minterms'); % Convert MinTerms to refer to the minimum number of steppable terms required MinTerms = max(0, MinTerms - sum(m.TermStatus==1)); % get the stats with all the terms in TermsIn = Terms(m); % Take terms out one at a time. tord = termorder(m); indices = fliplr(tord(:)'); ok = ismember(indices,find(m.TermStatus==3)); indices = indices(ok); % Make sure that we keep the minimum number of terms required if MinTerms<=length(indices) indices = indices(1:(end-MinTerms)); else % There are no spare terms to step out. indices = []; end s = zeros(length(indices)+1,1); % first row is current model s2full= mse(m); s(1)=summary(m,Criteria,m.Store.D,m.Store.y,s2full); np= s; np(1)= numParams(m); for i = 1:length(indices); [m,OK] = stepwise(m,indices(i));%toggle the state and get the fit np(i+1)= numParams(m); if OK s(i+1)= summary(m,Criteria,m.Store.D,m.Store.y,s2full); else s(i+1)= NaN; end end Tol= get(om,'tolerance'); [minCrit,list]= summary(m); critint= strcmp(list,Criteria); if minCrit(critint) [ms,ind]= min(s(:,1)); if Tol>0 tgt= ms + abs(ms)*Tol/100; pos= find(s(:,1)<tgt); ind= pos(end); end else [ms,ind]= max(s(:,1)); if Tol>0 tgt= ms - abs(ms)*Tol/100; pos= find(s(:,1)>tgt); ind= pos(end); end end % fit selected model TermsIn(indices(1:ind-1))= false; [m,OK] = stepwise(m,~TermsIn); cost= s(ind); if get(om,'Display'); figure plot(np,s(:,1),np(ind),s(ind,1),'*r') ylabel(Criteria) xlabel('Number of Parameters') end