[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: how can I do Principal Components Analysis with octave?
From: |
Vic Norton |
Subject: |
Re: how can I do Principal Components Analysis with octave? |
Date: |
Wed, 12 May 2004 09:31:01 -0400 |
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
-------------------------------------------------------------