nirscript.m

PCA and PLS aims:to get some insight into the bilinear factor models Principal Component Analysis

%Plot the design points:
clear all
close all
load NIRdata
plotNIRDesign(conc)
text(conc(:,1),conc(:,2),conc(:,3),ObjLab)
xlabel('ethanol')
ylabel('water')
zlabel('iso-propanol')

%Fig1

%Plot the pure components for the three different temperatures:
figure
subplot(3,1,1)
plot(wavelen,[ethanol30; ethanol50; ethanol70])
legend('30C','50C','70C',2)
title('ethanol')
ylabel('A.U.') 
grid
subplot(3,1,2)
plot(wavelen,[water30; water50; water70])
title('water')
ylabel('A.U.')
grid 
subplot(3,1,3)
plot(wavelen,[isopropanol30; isopropanol50; isopropanol70])
xlabel('wavelength (nm)')
ylabel('A.U.')
grid
title('iso-propanol')  

%Fig2

%Plot the different mixture-spectra
figure
subplot(3,1,1)
plot(wavelen,spec30,'b')
ylabel('A.U.')
grid
subplot(3,1,2)
plot(wavelen,spec50,'g')
ylabel('A.U.')
grid 
subplot(3,1,3)
plot(wavelen,spec70,'r')
xlabel('wavelength (nm)')
ylabel('A.U.')
grid 

%Fig3

% A PCA-script to analyse the 30C mixture-set could go as follows
[Xm30,mx30] = meanc(spec30);
[T30,P30,ssq30] = mypca(Xm30,5);
figure
subplot(2,1,1) 
plot(1:5,ssq30*100)
grid
xlabel('factor')
ylabel('X-var. expl. (%)')
title('PCA mix. 30C')
subplot(2,1,2)
plot(wavelen,P30(:,1:3))
grid
xlabel('wavelength (nm)')
legend('p_1','p_2','p_3',2)
figure
plotNIRDesign(T30(:,1:3));
text(T30(:,1),T30(:,2),T30(:,3),ObjLab);

%Fig4
%Fig5

% Try the same for the mixture-sets recorded 50C and 70C
[T70] = mypca(meanc(spec70),5);
figure
plot3(T30(:,1),T30(:,2),T30(:,3),'.b',T70(:,1),T70(:,2),T70(:,3),'ob')
hold on
for a=1:length(ObjLab), plot3([T30(a,1) T70(a,1)],[T30(a,2) T70(a,2)],[T30(a,3) T70(a,3)],'b'); end
hold off
title('Mix. (.) PCA scores 30C; (o) PCA scores 70C')
grid
[T30c,T70c,diagnos] = procrusND(T30,T70);
figure
plot3(T30c(:,1),T30c(:,2),T30c(:,3),'.b',T70c(:,1),T70c(:,2),T70c(:,3),'ob')
hold on
for a=1:length(ObjLab), plot3([T30c(a,1) T70c(a,1)],[T30c(a,2) T70c(a,2)],[T30c(a,3) T70c(a,3)],'b'); end
hold off
text(T30(:,1),T30(:,2),T30(:,3),ObjLab);
title('Mix. (.) PCA scores 30C; (o) PCA scores 70C after proc. rotation')
grid

[T30p,P30p,W30p,U30p,Q30p,B30p,ssq3p] = mypls(Xm30,conc,5);
figure
plot3(T30(:,1),T30(:,2),T30(:,3),'.b',T30p(:,1),T30p(:,2),T30p(:,3),'ob')
hold on
for a=1:length(ObjLab), plot3([T30(a,1) T30p(a,1)],[T30(a,2) T30p(a,2)],[T30(a,3) T30p(a,3)],'b'); end
hold off
title('Mix. 30 (.) PCA scores; (o) PLS scores')
grid
[T30c,T30pc,diagnos] = procrusND(T30,T30p);
figure
plot3(T30c(:,1),T30c(:,2),T30c(:,3),'.b',T30pc(:,1),T30pc(:,2),T30pc(:,3),'ob')
title('Mix. 30 (.) PCA scores; (o) PLS scores after proc. rotation')
grid