scilab/modules/statistics/macros/pca.sci

   1 function [lambda,facpr,comprinc]=pca(x,N)
   2 //
   3 //This  function performs  several  computations known  as
   4 //"principal component  analysis". It includes  drawing of
   5 //"correlations circle",  i.e. in the  horizontal axis the
   6 //correlation   values  r(c1;xj)   and  in   the  vertical
   7 //r(c2;xj). It is an extension of the pca function.
   8 //
   9 //The  idea  behind this  method  is  to  represent in  an
  10 //approximative  manner a  cluster of  n individuals  in a
  11 //smaller  dimensional subspace.  In order  to do  that it
  12 //projects the cluster onto a subspace.  The choice of the
  13 //k-dimensional projection subspace is  made in such a way
  14 //that  the distances  in  the projection  have a  minimal
  15 //deformation: we are looking for a k-dimensional subspace
  16 //such that the squares of the distances in the projection
  17 //is  as  big  as  possible  (in  fact  in  a  projection,
  18 //distances can only stretch).  In other words, inertia of
  19 //the projection  onto the k dimensional  subspace must be
  20 //maximal.
  21 //
  22 //x is a nxp (n  individuals, p variables) real matrix.
  23 //
  24 //lambda is a px2 numerical matrix. In the first column we
  25 //find the  eigenvalues of V,  where V is  the correlation
  26 //pxp matrix  and in the  second column are the  ratios of
  27 //the   corresponding   eigenvalue   over   the   sum   of
  28 //eigenvalues.
  29 //
  30 //facpr  are  the principal  factors:  eigenvectors of  V.
  31 //Each  column is an  eigenvector element  of the  dual of
  32 //R^p. Is an orthogonal matrix.
  33 //
  34 //comprinc  are  the  principal components.   Each  column
  35 //(c_i=Xu_i)  of  this  nxn  matrix  is  the  M-orthogonal
  36 //projection of individuals onto principal axis.  Each one
  37 //of this columns is a linear combination of the variables
  38 //x1,  ...,xp   with  maximum  variance   under  condition
  39 //u'_iM^(-1)u_i=1.
  40 //
  41 //Verification: comprinc*facpr=x
  42 //
  43 //References: Saporta, Gilbert, Probabilites,  Analyse des
  44 //Donnees et Statistique, Editions Technip, Paris, 1990.
  45 //
  46 //author: carlos klimann
  47 //
  48 //date: 2002-02-05
  49 //commentary fixed 2003-19-24
  50 //
  51   [lhs,rhs]=argn(0)
  52   if  rhs==0 then
  53     error('pca must have at least one parameter'),
  54   end
  55   if x==[] then
  56     lambda=%nan
  57     facpr=%nan
  58     comprinc=%nan,
  59     return,
  60   end
  61   [rowx colx]=size(x)
  62   if rhs==2 then
  63     if size(N,'*')<>2 then error('Second parameter has bad dimensions'), end,
  64     if max(N)>colx then disp('Graph demand out of bounds'), return, end
  65   else
  66     N=[1 2],
  67   end
  68   xbar=sum(x,'r')/rowx
  69   y=x-ones(rowx,1)*xbar
  70   std=(sum(y .^2,'r')) .^ .5
  71   V=(y'*y) ./ (std'*std)
  72   [lambda facpr]=bdiag(V,(1/%eps))
  73   [lambda order]=sort(diag(lambda))
  74   lambda(:,2)=lambda(:,1)/sum(lambda(:,1))
  75   facpr=facpr(:,order)
  76   comprinc=x*facpr
  77
  78
  79   w=winsid();if w==[] then w=0, else w=max(w)+1,end
  80   fig1=scf(w)
  81
  82   rc= (ones(colx,1)* sqrt((lambda(N,1))')) .* facpr(:,order(N))
  83
  84   rango=rank(V)
  85   ra=[1:rango]'
  86   if rango<=1 then return, end
  87   plot2d(-rc(ra,1),rc(ra,2),style=-10)
  88   legend('(r(c1,xj),r(c2,xj)')
  89   ax=gca();ax.x_location="middle";ax.y_location = "middle";
  90   blue=color('blue')
  91   for k=1:rango,
  92     xstring(rc(k,1),rc(k,2),'X'+string(k)),
  93     e=gce();e.foreground=blue;
  94   end
  95   title(' -Correlations Circle- ')
  96   fig2=scf(w+1);
  97   plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0])
  98   plot2d(ra,lambda(ra,2),style=9)
  99   ax=gca(); ax.grid=[31 31]
 100   plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0])
 101   plot2d(ra,lambda(ra,2),style=9)
 102   for k=1:rango,
 103     xstring(k,0,'l'+string(k)),
 104     e=gce();e.font_style=1
 105   end
 106   title(' -Eigenvalues (in percent)- ')
 107   ylabel('%')
 108
 109 endfunction
 110
 111