Re: how can I do Principal Components Analysis with octave?

[Vic Norton]

Hi Mario,


I'll give you my slant on Principal Components Analysis. I prefer
singular value decomposition to the covariance-eigen-value route. I
understand the former better, and it is more efficient.


Start with an m x n data matrix X and an m x 1 vector of weights w. We
require w(i) > 0 for i = 1,..,m and sum(w) = 1. In the normal
application w(i) = 1/m for i = 1,..,m.


Then
   E = sum(diag(w) * X)
is the 1 x n matrix of expected values of the columns of X,
   V = (X - ones(m, 1) * E)' * diag(w) * (X - ones(m, 1) * E)
is the n x n covariance matrix corresponding of X
(V = ((m - 1)/m) * cov(X) in the normal case),
   S = sqrt(diag(V))'
is the 1 x n matrix of standard deviations of the columns of X, and
   C = diag(1 ./ S) * V * diag(1 ./ S)
is the correlation matrix. (C = corrcoef(X) in the normal case.)


Now let's geometrize these computations. First reassign
   w = sqrt(w)
and put
   X = diag(w) * X.
Then the expected values of the columns of the original X are
   E = w' * X.
One can think of E as the perpendicular projection of the new X
onto w-axis (more or less) with
   F = X - w * E
the complementary projection. Now
   S = sqrt(sum(F .^ 2))
is the 1 x n matrix of standard deviations of the columns of the
original X, and
   V = F' * F
is the covariance matrix above. To normalize the situation set
   F1 = F * diag(1 ./ S).
Then
   C = F1' * F1
is the correlation matrix above.


OK, let's get to principal component analysis. We work with F or
F1. It is not necessary to compute the covariance or correlation
matrices.


If the units of measurement are important, do
   [U, Sig, W] = svd(F, 1)
otherwise do
   [U1, Sig1, W1] = svd(F1, 1).
The singular values, the coefficients of the diagonal matrices Sig and
Sig1, measure the significance of the various factors:
   Sig(1) >= Sig(2) >= ... >= 0
      with  sum(diag(Sig) .^ 2) = trace(V) = total variance
and
   Sig1(1) >= Sig1(2) >= ... >= 0
      with  sum(diag(Sig1) .^ 2) = trace(C) = n.
These singular values are, in effect, principal standard deviations.
Their squares are the eigenvalues of the covariance and correlation
matrices.


(Note in passing. Put  Y = Sig * W' and Y1 = Sig1 * W1'.
Then Y and Y1 are n x n and  V = Y' * Y, C = Y1' * Y1.)


For two dimensional plots, plot the n columns of
   Y([1, 2], :) = diag( [Sig(1), Sig(2)] ) * W(:, [1, 2])'
or
   Y1([1, 2], :) = diag( [Sig1(1), Sig1(2)] ) * W1(:, [1, 2])',
point by point. These points represent your data. The horizontal, 1-axis
is the most significant here. It is the axis of maximum variance
(correlation?) in the data.



Regards,

Vic


At 9:58 AM +0100 5/11/04, Ted Harding wrote:
>  On 11-May-04 rino mailing wrote:
>  > I'd like to do Principal Components Analysis with octave
>  >
>  > What are the command I ave to write?
>  >
>  > How to plot the result?
>
>  You find the principal components by first finding the eigenvalues
>  and eigenvectors of the covariance or correlation matrix (depending
>  on whether you want the results to be independent of the units in
>  which the variables are measured).
>
>  So let C be (say) the covariance matrix. Then
>
>    [V,D] = eig(C)
>
>  gives you the eigenvectors as columns of V, and the eigenvalues
>  as the diagonal elements of D, in increasing order of magnitude.
>
>  So, therefore, ifr X is the N by k matrix of data, C=cov(X),
>  [V,D]=eig(C), then X*V(:,k) gives you the first principal component,
>  X*V(:,k-1) the second, and so on.
>
>  I'm not sure what you mean by "plot the result". Some people want
>  to see a graph which shows the proportions of variance accounted
>  for by the successive principal components; for this you can plot
>  the diagonal elements diag(D) of D in reverse order, perhaps in the
>  form 100*diag(D)/sum(diag(D)). Or the amount of variance accounted
>  for by the first, first+second, first+second+third, etc; for this
>  you can plot the (reverse-order) cumsum of diag(D).
>
>  Or you may want to say plot the first principal component against
>  the second, so: PC1=X*V[:,k];PC2=X*V[:,k-1];plot(PC1,PC2)
>
>  This may be useful in exploring possible clustering; or you may
>  already have an index say Grp which identifies each row of X
>  as belonging to a group by Grp==1, Grp==2, Grp==3, ... , in which
>  case you can plot these subsets of (PC1,PC2) in different colours,
>  to see where the groups fall in the plot.
>
>  And so on ...
>
>  Hoping this helps,
>  Ted.
>
>
>  --------------------------------------------------------------------
>  E-Mail: (Ted Harding) <address@hidden>
>  Fax-to-email: +44 (0)870 167 1972
>  Date: 11-May-04                                       Time: 09:58:49



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------