stddemo.m

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

echo off
%STDDEMO Demos STDSSLCT and STDGEN for instrument standardization
echo on
% This script demonstrates the functions STDSSLCT and
% STDGEN for doing instrument standardization.

echo off
% Copyright
% Barry M. Wise
% 1994
echo on

% We're going to use some NIR spectrometer data from two instruments.
% Let's load the data and plot it up.
%pause

echo off
load nir_data
plot(lamda,spec1,'-r',lamda,spec2,'-b')
title('NIR Spectra')
xlabel('Frequency')
ylabel('Absorbance')
text(825,1.5,'Instrument number 1 spectras shown in red')
text(850,1.25,'Instrument number 2 spectras shown in blue')
echo on

% This plot shows both of the sets of spectra.  Here they don't look
% all that different, so lets plot the difference.
pause

echo off
plot(lamda,spec1-spec2,'-r')
title('Difference between NIR Spectra from Instruments 1 and 2')
xlabel('Frequency')
ylabel('Absorbance Difference')
echo on

% Now the difference is much more apparent.
pause

% The next thing we're going to do is get PLS calibration models that
% predict the concentration of each of the 5 analytes in the mixtures.
% We'll make a model based on the known concentrations and the spectra
% from the first instrument.  Later, we will use the models along with the
% transformed spectra from instrument 2.  We won't attempt to optimize
% the models for this excercise.  Instead, we'll just use 5 LVs in each
% one, since we know the spectra should have a true rank of 5 beause
% there are 5 analytes.

% This will take a couple minutes. Hang on.
pause

echo off 
[mspec1,mns1] = mncn(spec1);
[mspec2,mns2] = mncn(spec2);
[mconc,mnsc] = mncn(conc);
[P1,Q1,W1,T1,U1,b1,ssqdif1] = pls(mspec1,mconc(:,1),5);
c1p = plspred(mspec1,b1,P1,Q1,W1,5);
[P2,Q2,W2,T2,U2,b2,ssqdif2] = pls(mspec1,mconc(:,2),5);
c2p = plspred(mspec1,b2,P2,Q2,W2,5);
[P3,Q3,W3,T3,U3,b3,ssqdif3] = pls(mspec1,mconc(:,3),5);
c3p = plspred(mspec1,b3,P3,Q3,W3,5);
[P4,Q4,W4,T4,U4,b4,ssqdif4] = pls(mspec1,mconc(:,4),5);
c4p = plspred(mspec1,b4,P4,Q4,W4,5);
[P5,Q5,W5,T5,U5,b5,ssqdif5] = pls(mspec1,mconc(:,5),5);
c5p = plspred(mspec1,b5,P5,Q5,W5,5);
echo on

% Now that we have the models, we can look at the actual concentrations
% and the fit based on the PLS models.

echo off
plot(conc,rescale([c1p c2p c3p c4p c5p],mnsc),'o'), hold on
zax = axis;
plot([0 zax(2)],[0 zax(2)],'-r'), hold off
title('Actual versus Fit Concentrations Based on Instrument 1')
xlabel('Actual Analyte Concentration')
ylabel('Analyte Concentration Fit to Model')
text(5,50,'Each analyte shown as different color')
pause
echo on

% You can see that our PLS models, all based on 5 LVs, fit the data
% quite well.  Now we'll use the STDSSLCT function for selecting
% samples with high leverage to use for calibration transfer. We
% will select 5 samples out of the 30 available for calculating
% the transform between instruments.

[specsub,specnos] = stdsslct(spec1,5);
pause

% We can now use this subset of the samples to calculate the
% transform. (In this case, all the samples have already
% been measured on both instruments. If only the instrument 1
% samples had been measured, STDSSLCT would tell you which samples
% should be measured on the second instrument.)

% Now we can use the STDGEN function to obtain a transform
% that converts the response of intrument 2 to look like that of
% instrument 1.  STDGEN can be used to generate direct or
% piecewise direct standardizations with or without additive
% background correction. In this case we'll use both direct
% and picewise direct standardization with additive background
% correction and compare results. For the piecewise direct
% standardization we will use a window of 5 channels.

pause

[stdmatp,stdvectp] = stdgen(spec1(specnos,:),spec2(specnos,:),5);
[stdmatd,stdvectd] = stdgen(spec1(specnos,:),spec2(specnos,:),0);

% Now we can convert the second spectra by multiplying by the 
% transform matrices, and adding the background correction
% as follows.

cspec2p = spec2*stdmatp + ones(30,1)*stdvectp;
cspec2d = spec2*stdmatd + ones(30,1)*stdvectd;

% Lets look at the difference between instrument 1 spectra
% and the corrected instrument 2 spectra.

echo off
plot(lamda,spec1-spec2,'-r',lamda,spec1-cspec2p,'-b')
hold on
plot(lamda,spec1-cspec2d,'-g'), hold off
title('Difference between NIR Spectra from Instruments 1 and 2')
xlabel('Frequency')
ylabel('Absorbance Difference')
text(825,-.05,'Difference before correction shown in red')
text(850,-.075,'Difference after piecewise correction shown in blue')
text(875,-.1,'Difference after direct correction shown in green')
pause
echo on

% You can see that the differences are much smaller.  Now lets
% see how the predictions look based on the instrument 1 models 
% and the standardized instrument 2 spectra. This will take a 
% minute to make all the new predictions.

echo off
cmspec2p = scale(cspec2p,mns1);
cmspec2d = scale(cspec2d,mns1);
mspec2 = scale(spec2,mns1);
c1p2cp = plspred(cmspec2p,b1,P1,Q1,W1,5);
c2p2cp = plspred(cmspec2p,b2,P2,Q2,W2,5);
c3p2cp = plspred(cmspec2p,b3,P3,Q3,W3,5);
c4p2cp = plspred(cmspec2p,b4,P4,Q4,W4,5);
c5p2cp = plspred(cmspec2p,b5,P5,Q5,W5,5);
c1p2cd = plspred(cmspec2d,b1,P1,Q1,W1,5);
c2p2cd = plspred(cmspec2d,b2,P2,Q2,W2,5);
c3p2cd = plspred(cmspec2d,b3,P3,Q3,W3,5);
c4p2cd = plspred(cmspec2d,b4,P4,Q4,W4,5);
c5p2cd = plspred(cmspec2d,b5,P5,Q5,W5,5);
c1p2 = plspred(mspec2,b1,P1,Q1,W1,5);
c2p2 = plspred(mspec2,b2,P2,Q2,W2,5);
c3p2 = plspred(mspec2,b3,P3,Q3,W3,5);
c4p2 = plspred(mspec2,b4,P4,Q4,W4,5);
c5p2 = plspred(mspec2,b5,P5,Q5,W5,5);
plot(conc,rescale([c1p c2p c3p c4p c5p],mnsc),'o'), hold on
zax = axis;
plot([0 zax(2)],[0 zax(2)],'-r'), hold off
title('Actual versus Fit Concentrations Based on Instrument 1')
xlabel('Actual Analyte Concentration')
ylabel('Analyte Concentration Fit to Model')
text(5,45,'Each analyte shown as different color')
pause
echo on

% I put that last plot back up just so you'd remember how good
% the fit based on instrument 1 is.  Lets look at some 
% predictions based on the unstandardized instrument 2 spectra 
% and the instrument 1 models.

echo off  
plot(conc,rescale([c1p2 c2p2 c3p2 c4p2 c5p2],mnsc),'o'), hold on
title('Actual Concentrations vs. Predictions Based on Instrument 2');
zax = axis;
plot([0 zax(2)],[0 zax(2)],'-r'), hold off
xlabel('Actual Analyte Concentration')
ylabel('Predicted Analyte Concentration')
text(5,45,'Each analyte shown as different color')
pause
echo on

% As you can see, this doesn't look real great.  Now we can now
% look at the predictions based on the standardized spectra.

echo off 
plot(conc,rescale([c1p2cp c2p2cp c3p2cp c4p2cp c5p2cp],mnsc),'o')
hold on
plot(conc,rescale([c1p2cd c2p2cd c3p2cd c4p2cd c5p2cd],mnsc),'*')
zax = axis;
plot([0 zax(2)],[0 zax(2)],'-r'), hold off
title('Actual Concentrations vs. Predictions Based on Standardized Instrument 2');
xlabel('Actual Analyte Concentration')
ylabel('Predicted Analyte Concentration')
text(5,50,'Each analyte shown as different color')
text(7,47,'o Piecewise direct standardized samples')
text(9,44,'* Direct standardized samples')
hold off
pause
echo on

% I think you will agree that the predictions from the piecewise
% direct method are pretty good. The predictions using direct
% standardization are not as good.

% Let's put some numbers on the difference by calculating the 
% sum of squares errors.
echo off
ssq1 = sum(sum((conc-rescale([c1p c2p c3p c4p c5p],mnsc)).^2));
ssq2 = sum(sum((conc-rescale([c1p2 c2p2 c3p2 c4p2 c5p2],mnsc)).^2));
ssq2p = sum(sum((conc-rescale([c1p2cp c2p2cp c3p2cp c4p2cp c5p2cp],mnsc)).^2));
ssq2d = sum(sum((conc-rescale([c1p2cd c2p2cd c3p2cd c4p2cd c5p2cd],mnsc)).^2));
disp('  ')
disp('Sum of Squares Error for')
s = sprintf('Instrument 1 Fit Error              = %g',ssq1);
disp(s)
s = sprintf('Unstandardized Instrument 2         = %g',ssq2);
disp(s)
s = sprintf('Piecewise Standardized Instrument 2 = %g',ssq2p);
disp(s)
s = sprintf('Direct Standardized Instrument 2    = %g',ssq2d);
disp(s)
disp('  ')
pause
echo on

% So things did get quite a bit better with standardization!
% How much better?
echo off
disp('  ')
s = sprintf('By a factor of %g for piecewise and',ssq2/ssq2p);
disp(s)
s = sprintf('by %g for direct standardization.',ssq2/ssq2d);
disp(s)
echo off