Soft Independent Modeling of Class Analogy (SIMCA)

by

 

15 Mar 2011 (Updated )

M-files for classes modeling and prediction by SIMCA

simca.m
function simca_model=simca(x,classes)
% 
% Model for Soft Independent Modeling of Class Analogy (SIMCA)
% 
% simca_model=simca(x,classes)
%
% input:
% x (samples x descriptors)	for calibration
% classes (samples x 1)     classes numbers (must be >0)
% 
% output:
% simca_model struct with:
%     models (1 x No. classes) struct with: (one for each class)
%         x (samples x descriptors)   input for this class
%         clas (samples x 1)          numbers for this class
%         S (samples x pc)            scores for this model
%         L (descriptors x pc)        loadings for this model
%         ev (pc x 1)                 eigenvalues for this model
%         T2 (samples x 1)            T for samples of this model
%         Q2 (samples x 1)            Q for samples of this model
%         T2lim (1 x 1)               T limit for this model
%         Q2lim (1 x 1)               Q limit for this model
%     stats (1 x 1) struct with:
%         classes (No. classes x 1)	  class types
%         T2 (samples x 1)            T for calibration samples
%         Q2 (samples x 1)            Q for calibration samples
%         class_c (samples x 1)       Predicted classes for calibration samples
%         RMSEP (1 x 1)               Root Mean Square Error for Calibration
%         R2 (1 x 1)                  Correlation Coefficient for calibration
%         suc (1 x 1)                 Success (%) of classification for calibration samples
% 
% See also simcapred
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil
%  

%---classes---
y=classes;
yu=unique(y);
for i=1:length(yu)
ncl=yu(i);
ncli=find(y==ncl);
simca_model.models(i).x=x(ncli,:);
simca_model.models(i).clas=y(ncli,:);
end

%---CV---

for inc=1:length(yu)
    
    wb=waitbar(0,'Please wait...','Name','Performing cross-validation');
    
    x0=simca_model.models(inc).x;
    y0=simca_model.models(inc).clas;
    if min(size(x0))<=20;
        ncp=(min(size(x0)))-1;
    else
        ncp=20;
    end
    
    for incp=1:ncp
                
        waitbar(incp/ncp);

        for iam=1:size(x0,1)
            xv=x0(iam,:);
            xc=x0;
            xc(iam,:)=[];
            
            yv=y(iam,:);
            yc=y;
            yc(iam,:)=[];
            
            [E L varpc eive]=pcar(xc);
            ns=xv*L(:,1:incp);
            r=sum((ns*(L(:,1:incp))'-xv).^2);
                        
            rCP(iam,incp)=r;
            vCP(incp)=sum(varpc(1:incp));
        end
    end
    
    close(wb);
    
       press=sum(rCP);
       
       subplot(1,2,1);
       plot(press,'-o');
       t = sprintf('PRESS for class %g',yu(inc));
       title(t)
       xlabel('No. of PC on model');
       ylabel('PRESS');
       
       subplot(1,2,2);
       plot(vCP,'-o');
       t = sprintf('Explained variance for class %g',yu(inc));
       title(t)
       xlabel('No. of PC on model');
       ylabel('Explained variance (%)');

       cpret = input('How many PC do you want to retain?   ');
       
       close;
       clear rCP;
       clear vCP;
%---PCA---       
       [E L varpc eive]=pcar(x0);
       L=L(:,1:cpret);
       E=x0*L;
       
       
%---T---
if cpret>1
  t  = E*inv(diag(eive(1:cpret)));
  T  = sum((t.^2)')';
else
  T  = E.^2/eive(1,1);
end

%---Q---

res = (x0-E*L')';
Q   = (sum(res.^2))';

%---Tlim---
[m,n]=size(x0);

if m > 300
  Tlim = (cpret*(m-1)/(m-cpret))*finv(0.95,cpret,300);
elseif m == cpret;
  Tlim = 0;
else
  Tlim = (cpret*(m-1)/(m-cpret))*finv(0.95,cpret,m-cpret);
end

%---Qlim---
me=size(eive,1);
if cpret < min([m n])
  cl=2*95-100;
  t1 = sum(eive(cpret+1:me,1));
  t2 = sum(eive(cpret+1:me,1).^2);
  t3 = sum(eive(cpret+1:me,1).^3);
  h0 = 1-2*t1*t3/3/(t2^2);
  if h0<0.001
     h0=0.001;
  end
  ca = sqrt(2)*erfinv(cl/100);
  h1 = ca*sqrt(2*t2*h0^2)/t1;
  h2 = t2*h0*(h0-1)/(t1^2);
  Qlim = t1*(1+h1+h2)^(1/h0);
else
    Qlim = 0;

end

       simca_model.models(inc).S=E;
       simca_model.models(inc).L=L;
       simca_model.models(inc).ev=eive(1:cpret);
       simca_model.models(inc).T2=T;
       simca_model.models(inc).Q2=Q;
       simca_model.models(inc).T2lim=Tlim;
       simca_model.models(inc).Q2lim=Qlim;
        
end
simca_model.stats.classes=yu;

%---pred calib---
pred=simcapred(x,simca_model,classes);

simca_model.stats.T2=pred.T2;
simca_model.stats.Q2=pred.Q2;
simca_model.stats.class_c=pred.class_p;
simca_model.stats.RMSEC=pred.RMSEP;
simca_model.stats.R2=pred.R2;
simca_model.stats.suc=pred.suc;




%--- PCA function ---
function [E L varpc eive]=pcar(x)

m=size(x,1);
[u s L]=svd(x,'econ');
E=u*s;
eive=diag(s);
var=diag(s).^2;
varpc=var/sum(var)*100;

Contact us