PCA for confidence ellipses

25 views (last 30 days)
david  nickerson
david nickerson on 19 Dec 2011
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 Comments
bym
bym on 20 Dec 2011
I would say, if you don't understand it - don't use it
david  nickerson
david nickerson on 20 Dec 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.

Sign in to comment.

Answers (1)

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

Community Treasure Hunt

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

Start Hunting!