function PLS=pls(x0,y0,max); 
%[co,tt,pp,rr,uu,wstar,ww,cc]=plswhw(x0,y0,max) 
%programmed according to Wang Huiwen's book by Hongdong Li,Oct. 2006. 
%max is the interation epoch.co is the regression coefficients linking x0 and  
y0 
%model: x=t*p'   y=t*r'=u*q' 

xstd=standard(x0); 
ystd=standard(y0); 
x=xstd;y=ystd; 
%main 
matrix1=x'*y*y'*x; 
matrix2=y'*x*x'*y; 
for i=1:max 
w=fff(matrix1); ww(:,i)=w; 
c=fff(matrix2); cc(:,i)=c; 
t=x*w;   tt(:,i)=t; 
p=x'*t/norm(t)^2;pp(:,i)=p; 
u=y*c;uu(:,i)=u; 
r=y'*t/norm(t)^2;rr(:,i)=r; 
x=x-t*p'; 
y=y-t*r'; 
matrix1=x'*y*y'*x; 
matrix2=y'*x*x'*y; 
end 
yre=y; 
%get the regression weights********************* 
[m,n]=size(x); 
unit=eye(n); 
wstar(:,1)=ww(:,1); 
for i=2:max 
   wnew=eye(n,n);     
   for j=1:i-1 
       wnew=wnew*(unit-ww(:,j)*pp(:,j)'); 
   end 
   wstar(:,i)=wnew*ww(:,i);    
end 
%regression coefficients************************* 
coefficient=0; 
for j=1:max 
coefficient=coefficient+wstar(:,j)*rr(:,j)'; 
end 
%get regression coefficient linking x0 and y0**************** 
[xs,zx]=standard(x0); 
covx=cov(x0); 
[m2,n2]=size(x0);[m3,n3]=size(y0); 
co=zeros(n2+1,n3); 
for k=1:n3 
  constant=0; 
  for i=1:n2 
    co(i,k)=coefficient(i,k)*sqrt(cov(y0(:,k)))/sqrt(covx(i,i)); 
    constant=constant-zx(i)*coefficient(i,k)*sqrt(cov(y0(:,k)))/sqrt(covx(i,i) 
); 
  end 
  co(n2+1,k)=constant+mean(y0(:,k)); 
end 
%******************************************** 
x_expand=[x0 ones(size(x0,1),1)]; 
y_estimated=x_expand*co; 
error=y_estimated-y0; 
%******************************************** 
SST=sum((y0-mean(y0)).^2);SSR=sum((y_estimated-mean(y0)).^2); 
SSE=sum((y0-y_estimated).^2);R2=1-SSE/SST; 
%******************************************** 
 vip=vipp(xstd,ystd,tt,ww); 
%Output**************************************  
PLS.coef=co;PLS.W=wstar;PLS.xPC=tt;PLS.yPC=uu;PLS.eigvec_X=ww; 
PLS.var_imp=vip;PLS.y_est=y_estimated;PLS.res=error; 
PLS.SST=SST;PLS.SSR=SSR;PLS.SSE=SSE;PLS.RMSEP=sqrt(SSE/length(y0));PLS.R2=R2; 

%subfunctions******************************** 
function y=fff(x);   %find the eigenvector with maximum eigenvalue 
[m,n]=size(x);[a,b]=eig(x); 
diagb=diag(b); 
maxb=find(diagb==max(diagb)); 
y=a(:,maxb); 
%******************************************** 
function [xs,zx1]=standard(x);   %standarization 
[m,n]=size(x);zx1=mean(x); 
zx=repmat(zx1,m,1);x1=x-zx; 
covx=cov(x); 
for i=1:n 
    x1(:,i)=x1(:,i)/sqrt(covx(i,i)); 
end 
xs=x1; 
%******************************************* 
function vip=vipp(x,y,t,w); 
%to calculate the vip for each variable to the response 
%vip=sqrt(p*q/s) 
%initializing 
[m,p]=size(x); 
[m,h]=size(t); 
[p,h]=size(w); 
%calculate s 
for i=1:h 
    corr=corrcoef(y,t(:,h)); 
    co(i,1)=corr(1,2)^2; 
end 
s=sum(co); 
for i=1:p 
    for j=1:h 
        d(j,1)=co(j,1)*w(i,j)^2; 
    end 
    q=sum(d); 
    vip(i,1)=sqrt(p*q/s); 
end 
%*************************************************