fc948aeb1ae31907536485958bc03af9d352b9d6
[scilab.git] / 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