gusucode.com > 一个很好用得基于MATLAB的偏最小二乘回归得程序 > pls/pls.m

    % 非常好用得PLS程序

if exist('opt') 
disp(' ')
else
lv=input(' How many latent variables should be calculated (default 3!) ');if isempty(lv),lv=3;end;
Rx=input(' What is the order of X (default 3) ');if isempty(Rx),Rx=3;end;
if Rx==2
disp(' ')
disp(' Well, a tri-linear model will be made, but with on variable in the third order')
disp(' The only difference between ordinary and this bi-PLS is that no P loadings are')
disp(' introduced ')
disp(' ')
disp(' Hit any key to continue'),disp(' '),pause,end
Ry=input(' What is the order of Y (default 2 or 1) ');if isempty(Ry),Ry=2;end;

Xidx=['I';'J';'K';'L';'M';'N'];
Yidx=['Iy';'Jy';'Ky'];

[I,Jx]=size(X);[I,Jyy]=size(y);
if Rx==2, J=Jx;K=1;end
if Rx>2
for rx=2:Rx
if exist(Xidx(rx));rrx=eval(Xidx(rx));
if isempty(rrx),rrx=0;end,else,rrx=0;end;
str=([Xidx(rx),'=input('' What is the dimension of the ' , num2str(rx) , ' order of X (default ', num2str(rrx), ') '');']);
eval(str)
if isempty(eval(Xidx(rx))),str=([Xidx(rx),'=rrx;']);eval(str);end;
end,else,J=Jx;
end

if Ry>2
for ry=2:Ry
if exist(Yidx(ry,:));
rrx=eval(Yidx(ry,:));
if isempty(rrx),rrx=0;end
else,rrx=0;end
str=([Yidx(ry,:),'=input('' What is the dimension of the ' , num2str(ry) , ' order of Y  (default ', num2str(rrx), ') '');']);
eval(str)
if isempty(eval(Yidx(ry,:))),str=([Yidx(ry,:),'=rrx;']);eval(str),end;
end
else
Jy=Jyy;
end,clear rrx


end % if isempty(opt)


yres=y;
Xres=X;
ypred=zeros(size(y));
xmodel=zeros(I,Jx);
T=zeros(I,lv);
Wj=zeros(J,lv);
Wk=zeros(K,lv);
if Rx>3,Wl=zeros(L,lv);if Rx>4,Wm=zeros(M,lv);end,end
B=zeros(lv,lv);
Q=zeros(lv,Jyy);
Qj=zeros(Jy,lv);
if Ry>2,Qk=zeros(Ky,lv);end
U=zeros(I,lv);


	sakX=ssq(Xres); saky=ssq(y);
	
	for f=1:lv						% #2
	[u,bbb] = n_pca(yres,1,2);clear bbb
	maxit=250; it=0; ugl=u*2;;

		while (norm(u-ugl)/norm(u))>1e-8		% 		% #3		
		ugl=u;it=it+1;

%______________CALCULATE T

	if Rx<4   % (meaning Rx=3)
	kovxy=reshape((Xres'*u),J,K);
	[wj,wk] = n_pca(kovxy,1,2);
	wj=wj/norm(wj);
	wk=wk/norm(wk);w=reshape(wj*wk',J*K,1)';
	T(:,f)=Xres*w';
	Wj(:,f)=wj;
	Wk(:,f)=wk;
	end

	if Rx==4
	kovxy=reshape((Xres'*u),J,K*L);
	[wj,wk,wl]=n_pca(kovxy,1,3,K,L);
	wj=wj/norm(wj);	wk=wk/norm(wk);wl=wl/norm(wl);
	w=kron(wl,kron(wk,wj))';
	T(:,f)=Xres*w';
	Wj(:,f)=wj;
	Wk(:,f)=wk;
	Wl(:,f)=wl;end

	if Rx==5
	kovxy=reshape((Xres'*u),J,K*L*M);
	[wj,wk,wl,wm]=n_pca(kovxy,1,4,K,L,M);
	wj=wj/norm(wj);	wk=wk/norm(wk); wl=wl/norm(wl); wm=wm/norm(wm);
	w=kron(wm,kron(wl,kron(wk,wj)))';
	T(:,f)=Xres*w';
	Wj(:,f)=wj;
	Wk(:,f)=wk;
	Wl(:,f)=wl;
	Wm(:,f)=wm;
	end


%_______________CALCULATE Q & U

	if Ry<3
	q=T(:,f)'*yres;q=q'/norm(q);Qj(:,f)=q;
	u=yres*q;
	U(:,f)=u;end


	if Ry==3
	q=reshape((yres'*T(:,f)),Jy,Ky);[qj,qk] = n_pca(q,1,2);qj=qj/norm(qj);qk=qk/norm(qk);
	for i=1:I,yi=reshape(yres(i,:),Jy,Ky);
	u(i,1)=(yi*qk)'*qj;end							% #3
	Qj(:,f)=qj;Qk(:,f)=qk;U(:,f)=u;end

end


%_______________CALCULATE B

	B(1:f,f)=inv(T(:,1:f)'*T(:,1:f))*T(:,1:f)'*U(:,f); %y

%_______________CALCULATE PREDICTIONS

	if Ry<3
	ypred=T(:,1:f)*B(1:f,1:f)*Qj(:,1:f)';end
	
	if Ry==3, load=[]; for a=1:Ky load=[load qj'*qk(a)];end;Q(f,:)=load;
	ypred=T(:,1:f)*B(1:f,1:f)*Q(1:f,:);end
	

%________________CALCULATE RESIDUALS

	if Jy > 1,fprintf('number of iterations: %g',it);disp(' '),disp(' '),end
	xmodel=xmodel+T(:,f)*w;
	Xres=X-xmodel;
	yres=y-ypred;
	fprintf('Explained part of X: %g % ',(1-ssq(Xres)/sakX)*100);
	disp(' ')
	fprintf('Explained part of y: %g ',(1-ssq(y-ypred)/saky)*100);
	disp(' ')
	disp(' ')

end					% #2


if it==maxit
('Algoritmen failed')
end
clear w str sakX sakY rx ry maxit kovxy k i k a b c Xidx Yidx f l q u it Jyy ugl wj wk wl wm yres Xres