Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



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 nby2 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, 2D in this case.
k = finv(conf,p,np)*p*(n1)/(np); % Comment out line above and uncomment line below to use ftest toolbox. % k = fdistinv(p,np,1conf)*p*(n1)/(np);
[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



