Plotting a principal component analysis with contours

[Gerrit J. Kiers]

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
-------------------------------------------------------------