pcr1.m

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

function [t,p,b] = pcr1(x,y,pc)
%PCR1 Principal components regression for univariate y.
%  The inputs are the matrix of predictor variables (x),
%  vector of predicted variable (y), and maximum number
%  of principal components to consider (pc). The outputs
%  are the x-block scores (t), the x-block loadings (p) and
%  the matrix of regression coefficients (b) for each number 
%  of principal components, where each row corresponds to the
%  PCR model for that number of principal components.
%  The I/O format is: [t,p,b] = pcr1(x,y,pc);

%  Copyright
%  Barry M. Wise
%  1992
%  Modified by B.M. Wise, November 1993

[mx,nx] = size(x);
[my,ny] = size(y);
if mx ~= my
  error('Number of samples in x- and y-blocks not equal.')
elseif pc > nx
  error('pc must be <= number of x-block variables.')
end
if nx < mx
  cov = (x'*x)/(mx-1);
  [u,s,v] = svd(cov);
  p = v(:,1:pc);
else
  cov = (x*x')/(mx-1);
  [u,s,v] = svd(cov);
  v = x'*v;
  for i = 1:pc
    v(:,i) = v(:,i)/norm(v(:,i));
  end
  p = v(:,1:pc);
end
t = x*p;
b = zeros(pc,nx);
ssqty = zeros(pc,1);
ssqy = y'*y;
for i = 1:pc
  r = t(:,1:i)\y;
  b(i,:) = (p(:,1:i)*r)';
  dif = y - t(:,1:i)*r;
  ssqty(i,1) =  ((ssqy - dif'*dif)/ssqy)*100;
end
temp = diag(s)*100/(sum(diag(s)));
temp = temp(1:pc);
ssqy = zeros(pc,1);
ssqy(1,1) = ssqty(1,1);
for i = 1:pc-1
  ssqy(i+1,1) = ssqty(i+1,1)-ssqty(i,1);
end
ssq = [(1:pc)' temp cumsum(temp) ssqy ssqty];
disp('  ')
disp('       Percent Variance Captured by PCR Model   ')
disp('  ')
disp('           -----X-Block-----    -----Y-Block-----')
disp('   LV #    This PC    Total     This PC    Total ')
disp('   ----    -------   -------    -------   -------')
format = '   %3.0f     %6.2f    %6.2f     %6.2f    %6.2f';
for i = 1:pc
  tab = sprintf(format,ssq(i,:)); disp(tab)
end
disp('  ')