crcvrnd.m

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

function [press,fiterr,minlvp,b] = crcvrnd(x,y,split,iter,lv,powers,ss,mc)
%CRCVRND Cross validation for CR models using SDEP
%  Inputs are the matrix of predictor variables (x), matrix
%  of predicted variables (y), number of divisions of the data
%  (split), number of iterations of the cross validation to
%  perform (iter), maximum number of latent variables to calculate (lv),
%  the vector of continuum parameters to be considered (powers),
%  an optional variable (ss) which can be used to make the
%  routine use contiguous blocks of data when set to 1. Outputs are 
%  the prediction residual error sum of squares matrix (press),
%  the calibration fit error sum of squares matrix (fiterr), 
%  number of latent variables and power at minimum PRESS (minlvp),
%  and the final regression vector (b) at minimum PRESS. 
%  
%  Note: This cross validation routine is based on the SDEP 
%  or Standard Deviation of Error of Prediction method.
%  The cross validation procedure is repeated after re-ordering
%  the data and the final cumulative PRESS is reported.
%
%  I/O format is: 
%  [press,fiterr,minlvp,b] = crcvrnd(x,y,split,iter,lv,powers,ss,mc);


%  Copyright
%  Barry M. Wise
%  1991
%  Modified February 1994

[mx,nx] = size(x);
[my,ny] = size(y);
if mx ~= my
  error('Number of samples must be the same in both blocks')
end
if nargin < 7
  ss = 0;
end
if nargin < 8
  mc = 1;
end
[mp,np] = size(powers);
press = zeros(np,lv);
clc
ind = ones(split,2);
for i = 1:split
  ind(i,2) = round(i*mx/split);
end 
for i = 1:split-1
  ind(i+1,1) = ind(i,2) +1;
end
for kk = 1:iter
  if ss == 0
    [garb,inds] = sort(randn(1,mx));
	x = x(inds,:);
	y = y(inds,:);
  else
	if kk ~= 1
	  [garb,inds] = max(randn(1,mx));
	  x = [x(inds+1:mx,:); x(1:inds,:)];
	  y = [y(inds+1:mx,:); y(1:inds,:)];
	end
  end
  for i = 1:split
    home
    s = sprintf('Now working on test set %g of %g for iteration %g',i,split,kk);
    disp(s)

    if mc ~= 0
      [calx,mnsx] = mncn([x(1:ind(i,1)-1,:); x(ind(i,2)+1:mx,:)]);
      testx = scale(x(ind(i,1):ind(i,2),:),mnsx);
      [caly,mnsy] = mncn([y(1:ind(i,1)-1,:); y(ind(i,2)+1:mx,:)]);
      testy = scale(y(ind(i,1):ind(i,2),:),mnsy);
    else
      calx = [x(1:ind(i,1)-1,:); x(ind(i,2)+1:mx,:)];
      testx = x(ind(i,1):ind(i,2),:);
      caly = [y(1:ind(i,1)-1,:); y(ind(i,2)+1:mx,:)];
      testy = y(ind(i,1):ind(i,2),:);
    end
    bbr = cr(calx,caly,lv,powers);
	c = 0;
    for j = 1:np
      for k = 1:lv
	    c = c + 1;
        ypred = testx*bbr((c-1)*ny+1:c*ny,:)';
        press(j,k) = press(j,k) + sum(sum((ypred-testy).^2));
	  end
    end
  end
  surf(powers,1:lv,press')
  set(gca,'Xdir','reverse')
  set(gca,'Xscale','log')
  set(gca,'Ydir','reverse')
  title('Cumulative PRESS vs. Number of LVs and Power')
  ylabel('Number of Latent Variables')
  zlabel('PRESS')
  drawnow
end
[m1,i1] = find(press==min(min(press)));
minlvp = [i1 powers(m1)];
t = sprintf('Minimum PRESS is at %g LVs and Power %g',minlvp(1),minlvp(2));
disp(t)
disp('  ')
disp('Now working on final CR model')
b = cr(x,y,lv,powers);
c = 0;
for j = 1:np
  for k = 1:lv
    c = c + 1;
    ypred = x*b((c-1)*ny+1:c*ny,:)';
    fiterr(j,k) = sum(sum((ypred-y).^2));
  end
end
c = (m1-1)*lv + i1;
b = b((c-1)*ny+1:c*ny,:)';
plot(b), hold on, plot(b,'o'), plot(zeros(nx,1),'-g'), hold off
s = sprintf('Regression Coefficients Best PRESS Model %g LVs Power %g',...
minlvp(1), minlvp(2));
title(s)
xlabel('Variable Number')
ylabel('Regression Coefficient')