Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: Correcting confidence elipse scripts
Replies: 0  

Advanced Search

Back to Topic List Back to Topic List  
Daniel Robbins

Posts: 27
Registered: 6/12/12
Correcting confidence elipse scripts
Posted: Jul 18, 2012 9:39 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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
% last modified: 12 June 1998

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?

Thanks in advance
Dan



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.