1 function [lambda,facpr,comprinc]=pca(varargin)
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 // update disable graphics output Allan CORNET 2008
51 // request user
52   [lhs,rhs]=argn(0)
54   x = [];
55   N = [];
56   no_graphics = %f;
58   select rhs,
59         case 1 then
60           x = varargin(1);
61           N=[1 2];
62         case 2 then
63           x = varargin(1);
64           N = varargin(2);
65         case 3 then
66           x = varargin(1);
67           N = varargin(2);
68           no_graphics = varargin(3);
69         else
70         error(sprintf("%s: Wrong number of input argument(s).\n","pca"));
71         end
73   if ( type(no_graphics) <> 4 ) then
74     error(sprintf("%s: Wrong type for third input argument: Boolean expected.\n" ,"pca"));
75   end
77   if x==[] then
78     lambda=%nan;
79     facpr=%nan;
80     comprinc=%nan;
81     return;
82   end
84   [rowx colx]=size(x)
85   if size(N,'*')<>2 then error('Second parameter has bad dimensions'), end,
86   if max(N)>colx then disp('Graph demand out of bounds'), return, end
88   xbar=sum(x,'r')/rowx
89   y=x-ones(rowx,1)*xbar
90   std=(sum(y .^2,'r')) .^ .5
91   V=(y'*y) ./ (std'*std)
92   [lambda facpr]=bdiag(V,(1/%eps))
93   [lambda order]=sort(diag(lambda))
94   lambda(:,2)=lambda(:,1)/sum(lambda(:,1))
95   facpr=facpr(:,order)
96   comprinc=x*facpr
98   if ~no_graphics then
99     w = winsid();
100     if (w == []) then
101       w = 0;
102     else
103       w = max(w)+1;
104     end
105     fig1 = scf(w);
107     rc = (ones(colx,1)* sqrt((lambda(N,1))')) .* facpr(:,order(N)) ;
108     rango = rank(V);
109     ra = [1:rango]';
110     if ( rango <= 1 ) then return, end
111     plot2d(-rc(ra,1),rc(ra,2),style=-10);
112     legend('(r(c1,xj),r(c2,xj)');
113     ax=gca();
114     ax.x_location="middle";
115     ax.y_location = "middle";
116     blue=color('blue')
117     for k=1:rango,
118       xstring(rc(k,1),rc(k,2),'X'+string(k));
119       e=gce();
120       e.foreground=blue;
121     end
122     title(' -Correlations Circle- ');
123     fig2 = scf(w+1);
124     plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0]);
125     plot2d(ra,lambda(ra,2),style=9);
126     ax=gca();
127     ax.grid=[31 31];
128     plot2d3([0;ra;rango+1]',[0; lambda(ra,2);0])
129     plot2d(ra,lambda(ra,2),style=9)
130     for k=1:rango,
131       xstring(k,0,'l'+string(k)),
132       e=gce();e.font_style=1
133     end
134     title(' -Eigenvalues (in percent)- ')
135     ylabel('%')
136   end
137 endfunction