0001
0002 function [Y_fit,params]=ppr(X,Y,varargin);
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 fields={'components','size(X,1)',...
0034 'maxiter','50*size(X,1)',...
0035 'decay','((params.maxiter:-1:1)/params.maxiter)/3',...
0036 'tol',.01,...
0037 'resort',.90,...
0038 'disp',1,...
0039 'control','5',...
0040 'order',[],...
0041 'continuity',[],...
0042 'weight',[],...
0043 'type',[]};
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),
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,
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
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,
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