Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Correcting confidence elipse scripts

Subject: Correcting confidence elipse scripts

From: Daniel Robbins

Date: 18 Jul, 2012 13:39:35

Message: 1 of 1

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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us