Home > fcn > ppr.m

ppr

PURPOSE ^

PPR Projection-Pursuit non-linear Regression

SYNOPSIS ^

function [Y_fit,params]=ppr(X,Y,varargin);

DESCRIPTION ^

 PPR Projection-Pursuit non-linear Regression

 Y_FIT=PPR(X,Y) fits a PPR model to the data Y_FIT = Sum_i{ f_i(w_i'*X) }
 minimizing the sum of error squares sum(sum(abs(Y_FIT-Y).^2));
 X and Y are a predictor and output matrices, respectively (columns are observations)
 The functions f_i are estimated using piece-wise polynomial fitting (see PPF for additional options)

 optional arguments:
   'components'    : number of components extracted (default size(X,1))
   'maxiter'       : number of iterations (default 50*size(X,1))
   'disp'          : display on/off (defaults 1=on)
   'control'       : PPF number of control points (default 5)
   'order'         : PPF order of polynomial fit (default 3)
   'continuity'    : PPF order of continuity (default 2)

 [Y_FIT,PARAMS]=PPR(...) also returns the parameters of the fit
 Y_NEW=PPR(X_NEW,PARAMS) estimates new data using the previous fit

 example:
 X=randn(6,1e3); 
 n=3; m=5; 
 Y=.01*randn(1,size(X,2)); for n1=1:m, idx=ceil(size(X,1)*rand(ceil(rand*min(n,size(X,1))),1)); Y=Y+(prod(X(idx,:))); end
 [Y_fit,params]=ppr(X,Y);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 
0002 function [Y_fit,params]=ppr(X,Y,varargin);
0003 % PPR Projection-Pursuit non-linear Regression
0004 %
0005 % Y_FIT=PPR(X,Y) fits a PPR model to the data Y_FIT = Sum_i{ f_i(w_i'*X) }
0006 % minimizing the sum of error squares sum(sum(abs(Y_FIT-Y).^2));
0007 % X and Y are a predictor and output matrices, respectively (columns are observations)
0008 % The functions f_i are estimated using piece-wise polynomial fitting (see PPF for additional options)
0009 %
0010 % optional arguments:
0011 %   'components'    : number of components extracted (default size(X,1))
0012 %   'maxiter'       : number of iterations (default 50*size(X,1))
0013 %   'disp'          : display on/off (defaults 1=on)
0014 %   'control'       : PPF number of control points (default 5)
0015 %   'order'         : PPF order of polynomial fit (default 3)
0016 %   'continuity'    : PPF order of continuity (default 2)
0017 %
0018 % [Y_FIT,PARAMS]=PPR(...) also returns the parameters of the fit
0019 % Y_NEW=PPR(X_NEW,PARAMS) estimates new data using the previous fit
0020 %
0021 % example:
0022 % X=randn(6,1e3);
0023 % n=3; m=5;
0024 % Y=.01*randn(1,size(X,2)); for n1=1:m, idx=ceil(size(X,1)*rand(ceil(rand*min(n,size(X,1))),1)); Y=Y+(prod(X(idx,:))); end
0025 % [Y_fit,params]=ppr(X,Y);
0026 %
0027 
0028 
0029 
0030 % alfnie@bu.edu
0031 % 05/06
0032 
0033 fields={'components','size(X,1)',...    % Number of components
0034         'maxiter','50*size(X,1)',...    % Number of iterations for the estimation of each component
0035         'decay','((params.maxiter:-1:1)/params.maxiter)/3',... % S.D. of random steps during estimation of each component
0036         'tol',.01,...                   % When reaching R2=1-tol the algorithm stops
0037         'resort',.90,...                % If resort*R2 of a new component is larger than R2 of the previous component it will resort them
0038         'disp',1,...                    % Graphic display on/off
0039         'control','5',...               % PPF parameter
0040         'order',[],...                  % PPF parameter
0041         'continuity',[],...             % PPF parameter
0042         'weight',[],...                 % PPF parameter
0043         'type',[]};                     % PPF parameter
0044 params=[]; for n1=1:2:nargin-2, params=setfield(params,lower(varargin{n1}),varargin{n1+1}); end
0045 for n1=1:2:length(fields), if ~isfield(params,fields{n1}) | isempty(getfield(params,fields{n1})), if isstr(fields{n1+1}), params=setfield(params,fields{n1},eval(fields{n1+1})); else, params=setfield(params,fields{n1},fields{n1+1}); end; end; end
0046 
0047 if ~isstruct(Y), % fit
0048     sX=size(X);
0049     params.comp0=mean(Y,2);
0050     Y=Y-repmat(params.comp0,[1,size(Y,2)]);
0051     e0=sum(sum(abs(Y).^2));
0052 
0053     params.comp=[];
0054     Y0=Y;
0055     rnow=1;
0056     ncomp=1;
0057     while ncomp<=params.components & rnow>params.tol,
0058         if params.disp, cumdisp; end
0059         if length(params.comp)>=ncomp, % if already computed this component
0060             w=params.comp(ncomp).w;
0061             [Y_fit,Pfit]=ppf(w'*X,params.comp(ncomp).ppf);
0062             err=sum(sum(abs(Y_fit-Y).^2));
0063             if params.disp, subplot(2,params.components,params.components+ncomp); plot(w'*X,Y,'b.',w'*X,Y_fit,'r.'); set(gca,'xtick',[]); end
0064             if params.disp, cumdisp(['C',num2str(ncomp),'=',num2str(100*(rnow-err/e0),'%04.1f'),'%% ']); end
0065         else,
0066             % projection pursuit
0067             w=zeros(sX(1),1); 
0068             err=inf;
0069             ERR=nan+zeros(1,params.maxiter);
0070             for niter=1:params.maxiter, 
0071                 w_new=w+params.decay(niter)*randn(sX(1),1); w_new=w_new/(eps+sqrt(sum(abs(w_new).^2)));
0072                 [Y_fit_new,Pfit_new]=ppf(w_new'*X,Y,params.control,params.order,params.continuity,params.weight,params.type);
0073                 err_new=sum(sum(abs(Y_fit_new-Y).^2));
0074                 if err_new<err, 
0075                     err=err_new; w=w_new; Y_fit=Y_fit_new; Pfit=Pfit_new;
0076                     if params.disp>1, idxplot=ceil(rand*size(Y,1)); subplot(2,params.components,params.components+ncomp); plot(w'*X,Y(idxplot,:),'b.',w'*X,Y_fit(idxplot,:),'r.'); set(gca,'color','y','xtick',[]); end
0077                     if params.disp, cumdisp(['C',num2str(ncomp),'=',num2str(100*(rnow-err/e0),'%04.1f'),'%% ']); end
0078                 end
0079                 ERR(niter)=err_new/e0;
0080                 if params.disp>1, subplot(2,params.components,ncomp); plot(ERR(1:niter),'.-'); set(gca,'color','y'); drawnow; end
0081             end
0082             if params.disp, subplot(2,params.components,ncomp); plot(ERR(1:params.maxiter),'.-'); if ncomp==1, ylabel('error'); end; set(gca,'ylim',[0,1]); box off; end
0083         end
0084         params.comp(ncomp).ppf=Pfit;
0085         params.comp(ncomp).w=w;
0086         params.comp(ncomp).r=rnow-err/e0;
0087         if params.disp, idxplot=ceil(rand*size(Y,1)); subplot(2,params.components,params.components+ncomp); plot(w'*X,Y(idxplot,:),'b.',w'*X,Y_fit(idxplot,:),'r.'); if ncomp==1, h=ylabel('Y'); set(h,'rotation',0); end; h=xlabel(['z_{',num2str(ncomp),'}']); set(gca,'xtick',[]); h=title(['C',num2str(ncomp),' = ',num2str(100*params.comp(ncomp).r,'%2.1f'),'%']); set(h,'units','norm','position',[.4,.9,0],'fontweight','bold'); set(gcf,'color','w','doublebuffer','on'); drawnow; end;
0088         rnow=err/e0;
0089         Y=Y-Y_fit;
0090         if params.resort && ncomp>1,
0091             rall=cat(1,params.comp(1:ncomp).r);
0092             idx=find((.01+rall(1:end-1))./rall(end)<params.resort);
0093             if ~isempty(idx),
0094                 params.comp=params.comp([1:idx(1)-1,ncomp]);
0095                 Y=Y0;
0096                 rnow=1;
0097                 ncomp=0;
0098                 if params.disp, fprintf('\n'); end
0099             end
0100         end
0101         ncomp=ncomp+1;
0102     end
0103     params.components=ncomp-1;
0104     Y_fit=repmat(params.comp0,[1,size(Y,2)])+Y0-Y;
0105 else,   % estimate
0106     params=Y;
0107     Y_fit=repmat(params.comp0,[1,size(X,2)]);
0108     for ncomp=1:params.components,
0109         Y_fit=Y_fit+ppf(params.comp(ncomp).w'*X,params.comp(ncomp).ppf);
0110     end
0111 end
0112

Generated on Tue 27-Mar-2007 12:06:24 by m2html © 2003