偏最小二乘法PLS回归NIPALS算法的Matlab程序及例子.doc

偏最小二乘法PLS回归NIPALS算法的Matlab程序及例子function T,P,W,Wstar,U,b,C,B_pls,. Bpls_star,Xori_rec,Yori_rec,. R2_X,R2_Y=PLS_nipals(X,Y,nfactor)% USAGE: T,P,W,Wstar,U,b,C,Bpls,Bpls_star,Xhat,Yhat,R2X,R2Y=PLS_nipals(X,Y,nfact)% PLS regression NIPALS algorithm PLS回归NIPALS算法% Compute the PLS regression coefficients PLS回归系数的计算% X=T*P Y=T*B*C=X*Bpls X and Y being Z-scores% B=diag(b)% Y=X*Bpls_star with X being augmented with a col of ones% and Y and X having their original units% T*T=I (NB normalization SAS)% W*W=I% Test for PLS regression% Herve Abdi November 2002/rev November 2004% Version with T, W, and C being unit normalized% U, P are not% nfact=number of latent variables to keep 保持潜在变量的数量% default = rank(X)X_ori=X;Y_ori=Y;if exist(nfactor)=1;nfactor=rank(X);endM_X=mean(X);M_Y=mean(Y);S_X=std(X);S_Y=std(Y);X=zscore(X);Y=zscore(Y);nn,np=size(X) ;n,nq=size(Y) ;if nn= n; error(Incompatible # of rows for X and Y);end% Precision for convergenceepsilon=eps;% # of components kepts% Initialistion% The Y setU=zeros(n,nfactor);C=zeros(nq,nfactor);% The X setT=zeros(n,nfactor);P=zeros(np,nfactor);W=zeros(np,nfactor);b=zeros(1,nfactor);R2_X=zeros(1,nfactor);R2_Y=zeros(1,nfactor);Xres=X;Yres=Y;SS_X=sum(sum(X.2);SS_Y=sum(sum(Y.2);for l=1:nfactort=normaliz(Yres(:,1);t0=normaliz(rand(n,1)*10);u=t;nstep=0;maxstep=100;while ( ( (t0-t)*(t0-t) epsilon/2) & (nstep This is your title.% - Change it to fit your dataze_title=PLS. Wines. ;% *% This is the data matrix of the Predictors (IV)% - Change it for your exampleX=. 7 7 13 7 4 3 14 7 10 5 12 5 16 7 11 3 13 3 10 3 ;% get the # of rows and columns %nI,nJ=size(X);%*% - These are your predictors names.% - Change them to fit your data% You need as many names as mcX has columnsn=0;%n=n+1;nom_x(n)= Price;n=n+1;nom_x(n)= Sugar;n=n+1;nom_x(n)= Alcohol;n=n+1;nom_x(n)= Acidity;% Test the # of colums namesif nJ=n; erreur(Error - (Wrong # of column names)!);end%*% - These are your observations names.% - Change them to fit your datal=0;l=l+1;nom_rl=W_1;l=l+1;nom_rl=W_2;l=l+1;nom_rl=W_3;l=l+1;nom_rl=W_4;l=l+1;nom_rl=W_5;if nI=l; erreur(Error - (Wrong # of row names)!);end% *% This is the data matrix of the Dependant Variables (DV)% - Change it for your exampleY=. 14 7 8 10 7 6 8 5 5 2 4 7 6 2 4 ;% get the # of rows and columns %nI2,nK=size(Y);%*% - These are your predictors names.% - Change them to fit your data% You need as many names as mcX has columnsm=0;%m=m+1;nom_y(m)= Hedonic;m=m+1;nom_y(m)= Meat;m=m+1;nom_y(m)= Dessert;% Test the # of colums namesif nK=m; erreur(Error - (Wrong # of column names)!);end%*% Call PLS_nipals%nfact2keep=2 ;% nfact gives the number of latent variables% the default is equal to the rank of XT,P,W,Wstar,U,b,C,Bpls,Bpls_star,Xhat,Yhat,R2X,R2Y=. PLS_nipals(X,Y,nfact2keep)% Plot the Observations (T vectors)%percent_of_inertiaX=100*R2X;percent_of_inertiaY=100*R2Y;% The axes to keep for the plotsaxe_horizontal=1;axe_vertical=2;%le_taux=. tau_X_,int2str(axe_horizontal),=,. int2str(percent_of_inertiaX(axe_horizontal),%, . tau_X_,int2str(axe_vertical),=,. int2str(percent_of_inertiaX(axe_vertical),%, . tau_Y_,int2str(axe_horizontal),=,. int2str(percent_of_inertiaY(axe_horizontal),%, . tau_Y_,int2str(axe_vertical),=,. int2str(percent_of_inertiaY(axe_vertical),%;% Plot here %figure(1);clf% Now plot the Observations Scores T%ze_tRC=ze_title, Observations (T).;titre=ze_tRC, le_taux;plotxyha(T,1,2,titre,nom_r);axis(equal) ;% Now plot the X Scores W%figure(2);clfze_tRC=ze_title, Predictors (W).;titre=ze_tRC, le_taux;plotxyha(W,1,2,titre,nom_x);axis(equal) ;% Now plot the Y Scores U%figure(3);clfze_tRC=ze_title, DV (C - Non Ortho).;titre=ze_tRC, le_taux;plotxyha(C,1,2,titre,nom_y);% axis(equal) ;%作图函数plotxy定义function plotxy(F,axe1,axe2,titre,noms);% % * This is a Test Version *% July 1998 Herve Abdi% Usage plotxy(F,axe1,axe2,title,names);% plotxy plots a MDS or PCA or CA graph of component #Axe1 vs #Axe2% F is a matrix of coordinates% axe1 is the Horizontal Axis% axe2 is the Vertical Axis% title will be the title of the graph% Axes are labelled Principal Component number% names give the names of the points to plot (def=numbers)nomdim=Dimension;nrow,ncol=size(F);if exist(noms)=0; nomsnrow,1=; for k=1:nrow;nomsk,1=int2str(k);endendminx=min(F(:,axe1);maxx=max(F(:,axe1);miny=min(F(:,axe2);maxy=max(F(:,axe2);hold off; clf;hold on;rangex=maxx-minx;epx=rangex/10;rangey=maxy-miny;epy=rangey/10; axis(equal);axis(minx-epx,maxx+epx,miny-epy,maxy+epy) ;%axis(equal);%axis;plot ( F(:,axe1),F(:,axe2 ),.);la