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) <***@nessie.mcc.ac.uk>
> 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
-------------------------------------------------------------