Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Correcting confidence elipse scripts
Replies: 0

 Daniel Robbins Posts: 27 Registered: 6/12/12
Correcting confidence elipse scripts
Posted: Jul 18, 2012 9:39 AM

Hi,

I am trying to work out how to generate confidence ellipses for datasets.

I have tried the following code:

%# generate data
num = 50;
X = [ mvnrnd([0.5 1.5], [0.025 0.03 ; 0.03 0.16], num) ; ...
mvnrnd([1 1], [0.09 -0.01 ; -0.01 0.08], num) ];
G = [1*ones(num,1) ; 2*ones(num,1)];

gscatter(X(:,1), X(:,2), G)
hold on
axis equal

for k=1:2
%# indices of points in this group
idx = ( G == k );

%# substract mean
Mu = mean( X(idx,:) );
X0 = bsxfun(@minus, X(idx,:), Mu);

%# eigen decomposition [sorted by eigen values]
[V D] = eig( X0'*X0 ./ (sum(idx)-1) ); %#' cov(X0)
[D order] = sort(diag(D), 'descend');
D = diag(D);
V = V(:, order);

t = linspace(0,2*pi,100);
e = [cos(t) ; sin(t)]; %# unit circle
VV = V*sqrt(D); %# scale eigenvectors
e = bsxfun(@plus, VV*e, Mu'); %#' project circle back to orig space

%# plot cov and major/minor axes
plot(e(1,:), e(2,:), 'Color','k');
%quiver(Mu(1),Mu(2), VV(1,1),VV(2,1), 'Color','k')
%quiver(Mu(1),Mu(2), VV(1,2),VV(2,2), 'Color','k')

end

However this seems to drastically underestimate the confidence ellipse. I also tried this function I found in the file exchange:

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.

% written by Douglas M. Schwarz
% schwarz@kodak.com

n = size(xy,1);
mxy = mean(xy);

numPts = 181; % The number of points in the ellipse.
th = linspace(0,2*pi,numPts)';

p = 2; % Dimensionality of the data, 2-D in this case.

k = finv(conf,p,n-p)*p*(n-1)/(n-p);
% Comment out line above and uncomment line below to use ftest toolbox.
% k = fdistinv(p,n-p,1-conf)*p*(n-1)/(n-p);

[pc,score,lat] = princomp(xy);
% Comment out line above and uncomment 3 lines below to use ftest toolbox.
% xyp = (xy - repmat(mxy,n,1))/sqrt(n - 1);
% [u,lat,pc] = svd(xyp,0);
% lat = diag(lat).^2;

ab = diag(sqrt(k*lat));
exy = [cos(th),sin(th)]*ab*pc' + repmat(mxy,numPts,1);

% Add ellipse to current plot
h = line(exy(:,1),exy(:,2),'Clipping','off');
if nargout > 0
hh = h;
end

However this seems to drastically over estimate the ellipse.

Can anyone help me work out which code is nearer to being correct and how to amend the code?