PCA for confidence ellipses

I have a cloud of two dimensional data (catesian or polar coordinates, don't mind which) and want to plot a confidence ellipse based on a principle components analysis. i have some code to do this (see below), but i also want to get out all the information i can about the orientation of the elipses/ relative sizes of the principle axes. any explanation would be very helpful, as part of this code was given to me and i don't really understand it, though i have done my best to comment it appropriately.
function hh = confellipse2(xy,conf)
%CONFELLIPSE2 Draws a confidence ellipse.
% CONFELLIPSE2(XY,CONF) draws a confidence ellipse on the current axes
% which is calculated from the n-by-2 matrix XY and encloses the
% fraction CONF (e.g., 0.95 for a 95% confidence ellipse).
% H = CONFELLIPSE2(...) returns a handle to the line.
n = size(xy,1);
mxy = mean(xy);
numPts = 200; % The number of points in the ellipse.
th = linspace(0,2*pi,numPts)';
%dimensionality of the data
p = 2;
%convert confidence rating (eg 0.95) into z score - relative to size of
%sample(n) and the dimensionality of the data, n-p is therefore the degrees
%of freedom.
k = finv(conf,p,n-p)*p*(n-1)/(n-p);
% principle components analysis, lat gives eigenvalues
[pc,score,lat] = princomp(xy);
ab = diag(sqrt(k*lat));
exy = [cos(th),sin(th)]*ab*pc' + repmat(mxy,numPts,1);
% Add ellipse to current plot
%plot(x,y)
hold on
h = line(exy(:,1),exy(:,2),'Clipping','off');
if nargout > 0
hh = h;
end

2 comentarios

bym
bym el 20 de Dic. de 2011
I would say, if you don't understand it - don't use it
david  nickerson
david nickerson el 20 de Dic. de 2011
ok, in that case proecsm could you please suggest an alternative method? how would you go about solving the problem? any help would be much appreciated.

Iniciar sesión para comentar.

Respuestas (1)

Mohammad Abdolrahmani
Mohammad Abdolrahmani el 8 de Jun. de 2017

1 voto

This code was written by Douglas M. Schwarz, and you have just removed his name and claim it's yours. Not fair!

Categorías

Preguntada:

el 19 de Dic. de 2011

Respondida:

el 8 de Jun. de 2017

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by