pwrdemo.m

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

echo on
%PWRDEMO Demonstrates POWERPLS for continuum regression (CR)
% This script illustrates the Wise Ricker method for
% continuum regression using the POWERPLS function in
% the PLS_Toolbox.

echo off
%  Copyright
%  Barry M. Wise
%  1992
%  Modified May 1994
echo on

% The data we are going to work with is from a Liquid-Fed
% Ceramic Melter process. The variables are temperatures
% in the molten glass tank and the tank level.
%
% Lets start by loading and plotting the data.  Hit a key when
% you are ready.
pause

load plsdata

echo off
plot(xblock1)
title('X-block data for Continuum Regression Demo')
xlabel('Sample Number')
ylabel('Temperature')
pause
echo on

% You can probably already see that there is a very regualar
% variation in the data.  When you see the level plotted, you
% will know why.
pause

echo off
plot(yblock1)
title('Y-Block data for Continuum Regression Demo')
xlabel('Sample Number')
ylabel('Level')
pause
echo on

% Yes that's right, the level correlates with temperature.  Now
% lets use that fact to build a PLS model that uses temperature
% to predict level.  Lets start by mean centering the data.

[mx,mnx] = mncn(xblock1);
[my,mny] = mncn(yblock1);

% Now that the data is scaled we can use the POWERPLS routine 
% to develop a calibration. We'll also make some models using 
% PCR, PLS and MLR and compare them to the CR models.
% We will cross validate over the region of powers of 1/2 to
% 1/8. (From previous tests, I've observed that the best models
% for this problem lie in this range.) Because this is time
% series data, we will use block cross validation, i.e. 
% contiguous blocks of data will be used as test sets.
% We will start by calculating the CR models, using a maximum
% of 5 latent variables. While we're at it, we'll also determine
% the cross validation error on an MLR model.
pause

echo off
powers = [1/16 1/32 1/64 1/128 1/256 1/512 ...
 1/1024 1/2048 1/4096].^0.25; 
errmat = zeros(5,9);
ssq = zeros(1,45);
mlrssq = 0;
for i = 1:4
  e = i*75;
  b = e-74;
  testx = mx(b:e,:);
  testy = my(b:e,:);
  calx = [mx(1:b-1,:); mx(e+1:300,:)];
  caly = [my(1:b-1,:); my(e+1:300,:)];
  m = powerpls(calx,caly,5,powers);
  bmlr = calx\caly;
  ypreds = testx*m';
  ypredmlr = testx*bmlr;
  yact = testy(:,ones(45,1));
  ydif = ypreds - yact;
  ydifmlr = ypredmlr - testy;
  ssq = ssq + sum(ydif.^2);
  mlrssq = mlrssq + sum(ydifmlr.^2);
end
disp('Now working on final CR model')
m = powerpls(mx,my,5,powers);
bmlr = mx\my;
[mincr,indcr] = min(ssq);
mcr = m(indcr,:);
for i = 1:9, e = i*5; b = e-4; errmat(:,i) = ssq(b:e)'; end
echo on

% It is interesting to look at the table of cross-validation
% errors displayed below
echo off
disp('              Number of Latent Variables')
disp('    Power')
disp([0 1:5; powers' errmat'])
pause
echo on

% Note that the best model is found at a power of 0.1768
% and one latent variable. 


% We can now use the PCRCVBLK and PLSCVBLK routines to 
% develop PCR and PLS models.
pause

[pcrp,pcrcp,minlvpls,bpls] = pcrcvblk(mx,my,4,10,1);
[plsp,plscp,minlvpcr,bpcr] = plscvblk(mx,my,4,10,1);

% Note the errors of cross validation for the different
% techniques as displayed below:
echo off
disp('           Cross Validation Error     ')
disp('      CR       PCR       PLS       MLR')
disp([min(ssq) min(pcrcp) min(plscp) mlrssq])
pause
echo on

% Based on this we would expect that the model from CR
% would have slightly better predictive abilites compared
% to the PCR and PLS models and significantly better 
% than the MLR model.