|
|
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
% 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
|
|
