cr.m

偏最小二乘算法在MATLAB中的实现

function b = cr(x,y,lv,powers);
%CR Continuum Regression for multivariate y
%  The inputs are the matrix of predictor variables (x), matrix
%  or vector of predicted variables (y), number of latent variables 
%  to consider (lv) and a vector of powers to consider (powers).
%  The output is a matrix of regression vectors or matrices (b) where
%  each row, or block of ny rows (where ny is the number of predicted 
%  variables), of b correspond to a regression vector or matrix for a 
%  particular power and number of latent variables. The regression vectors
%  or are grouped by power and ordered by number of latent variables.
%  For example, if there are 2 y variables and 10 lvs considered, the first
%  20 rows of b are the models for the first power considered, the first
%  2 rows are the regression matrix for the first power and 1 latent variable.
%  
%  The I/O syntax is b = cr(x,y,lv,powers);
%
%  This routine is based on one developed by Sijmen de Jong and uses
%  the continuum power method in canonical space. The routine was
%  tested and modifed by Barry M. Wise.
[mx,nx] = size(x);
[my,ny] = size(y);
if mx >= nx
  [u,s,v] = svd(x,0);  
else 
  [v,s,u] = svd(x',0);
end
s1 = s(1);
s = diag(s)/s1;
rnk = sum(s>1e-14);
s = s(1:rnk);
u = u(:,1:rnk);
v = v(:,1:rnk);
if lv > rnk
 disp('Number of Latent Variables > Rank of x-block')
 disp(sprintf('Changing number of LVs to Rank x-block = %g',rnk))
 lv = rnk;
end
np = length(powers);
uty = u'*y;
b = zeros(rnk,np*lv*ny);
k = 0;
for j = 1:np
  sp2 = s.^(2*powers(j));
  f = uty;
  zz = zeros(rnk,lv);
  for a = 1:lv
    if ny == 1
      z = sp2.*f;
	else
	  ss = diag(sp2)*f;
	  [qa,covsq] = eig(ss'*f);
	  z = ss*qa(:,find(diag(covsq)==max(diag(covsq))));
	end
    z = z-zz(:,1:max(1,a-1))*(z'*zz(:,1:max(1,a-1)))';
    z = z/sqrt(z'*z);
    zz(:,a) = z;
    c = z'*f;
    f = f-z*c;
    k = k+1;
	lastb = b(:,ny*(k-(a>1)-1)+1:ny*(k-(a>1)));
    b(:,ny*(k-1)+1:ny*k) = lastb + (z./s)*(c/s1);
  end
end
b = (v*b)';