[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Plotting a principal component analysis with contours
From: |
Gerrit J. Kiers |
Subject: |
Plotting a principal component analysis with contours |
Date: |
Fri, 4 May 2001 00:24:48 +0200 |
Hello all,
I want to use Octave to plot a principal component analysis, including
contour ellipses around the data for say 75%, 90% and 95% probabilities.
I use sources that have been published on this list in early 1999 (see
130.txt, 131.txt and 143.txt in the 1999-archive). A contour is plotted
using the 2x2 covariance matrix with function file 'ellipse.m' (contribution
of John W. Eaton).
I lack the mathematical background required for determining the right values
for 'level' for the different probability contours. Level is a multiplier
(variable) in 'ellipse.m'. It is multiplied with the square root of the
eigenvalues.
I assume that the probability along the two principal axis of the ellipse
can
be calculated with:
> normal_pdf(*somethingtodowitheigenvalue*,Mean_t(1,1),Var_t(1,1) )
and
> normal_pdf(*somethingtodowitheigenvalue*,Mean_t(1,2),Var_t(1,2) )
But I cannot get my finger behind it. Frankly: I have no clue about
Eigenvalue. And I've not found an statistical book or handbook on multi
variate analysis that I can access easily to solve my questions. Please
help me on this.
Then I still would need to use an inverse of the normal_pdf(x,m,v), so a
given probability x results in a value for level. Does such an inverse
exist? Does anyone know of another valid approach?
And finally I wonder if it is correct to use only 1 value for level. I
expect that since mean and variance are 1x2 matrices, level should also be a
1x2 matrix.
It would generate these lines
a = level(1,2) * sqrt (l(1,1));
b = level(2,2) * sqrt (l(2,2));
(as opposite to the current
a = level * sqrt (l(1,1));
b = level * sqrt (l(2,2)); )
Am I correct here? Or does the Eigenvalue do the trick here?
Any help is highly appreciated.
Kind regards,
Gerrit Kiers
2 files:
pca_plot.m
#########################################################
# This function file is in development.
t = [rand()+randn(100,1), rand()+randn(100,1)];
Covariance_matrix = cov(t)
[V,v] = eig(cov(t));
PCA = V
Eigenvalue = v
Mean_t = mean(t)
Var_t = var(t)
level = 1; # BEWARE THIS IS A NONSENSE VALUE
[x, y] = ellipse (cov(t), level);
x75 = x; y75 = y;
level = 2; # BEWARE THIS IS A NONSENSE VALUE
[x, y] = ellipse (cov(t), level);
x90 = x; y90 = y;
level = 3; # BEWARE THIS IS A NONSENSE VALUE
[x, y, major, minor] = ellipse (cov(t), level);
x95 = x; y95 = y;
gset size ratio 1;
gset nokey
plot (t(:,1), t(:,2), "*", minor(:,1)+mean(t(:,1)), minor(:,2)+mean(t(:,2)),
"-0", major(:,1)+mean(t(:,1)), major(:,2)+mean(t(:,2)),"-0",
x75+mean(t(:,1)), y75+mean(t(:,2)), "-b", x90+mean(t(:,1)),
y90+mean(t(:,2)), "-b", x95+mean(t(:,1)), y95+mean(t(:,2)), "-b");
#########################################################
ellipse.m
#########################################################
## Plotting an ellipse with an aspect ratio not equal to one can be
## mighty confusing.
gset size ratio 1;
## [x, y, major, minor, bbox] = ellipse (amat, level, n)
##
## Given a 2x2 matrix, generate ellipse data for plotting.
function [x, y, major, minor, bbox] = ellipse (amat, level, n)
if (nargin < 3)
n = 181;
endif
if (nargin > 1)
[v, l] = eig (amat);
dl = diag(l);
if (any (imag (dl)) || any (dl <= 0))
error ("ellipse: amat must be positive definite");
endif
## Generate contour data.
a = level * sqrt (l(1,1));
b = level * sqrt (l(2,2));
t = linspace (0, 2*pi, n)';
xt = a * cos (t);
yt = b * sin (t);
## Rotate the contours.
ra = atan2 (v(2,1), v(1,1));
cos_ra = cos (ra);
sin_ra = sin (ra);
x = xt * cos_ra - yt * sin_ra;
y = xt * sin_ra + yt * cos_ra;
## Endpoints of the major and minor axes.
major = minor = (v * diag ([a, b]))';
major (2,:) = -major (1,:);
minor (1,:) = -minor (2,:);
## Bounding box for the ellipse using magic formula.
ainv = inv (amat);
xbox = sqrt (level * ainv(1,1));
ybox = sqrt (level * ainv(2,2));
bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];
else
usage ("ellipse (amat, level, n)");
endif
endfunction
#########################################################
-------------------------------------------------------------
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
-------------------------------------------------------------
- Plotting a principal component analysis with contours,
Gerrit J. Kiers <=