function [x,y,xr1,yr1,xr2,yr2] = ellipse_ci(ybar,Sigma,n,gamma,npts); 
% 
% Construct 100*gamma % confidence ellipsoid for mu
% assuming Y is bivariate normal with mean E(Y) = mu and var(ybar) = Sigma. 
% Output: x and y are vectors (npts x 1) giving points on the 
% ellipsoid. 
% 
% If Sigma is known, then input -n.  Otherwise input n = sample size.
% 
if n< 0
  crit=-chi2inv(gamma,2)/n;
else
  crit=((n-1)*2/(n*(n-2)))*finv(gamma,2,n-2);
end
n=npts+1; 
if npts==4, 
     n=4; 
end; 
[U,L,V]=svd(inv(Sigma)); 
L1=L(1,1); 
L2=L(2,2); 
Sql=[sqrt(L1) 0; 0 sqrt(L2)];
f=zeros(2,1); 
g=zeros(2,1); 
del=(2*pi)/npts; 
x=zeros(n,1); 
y=zeros(n,1); 
t=-del; 
for j=1:n; 
t=t+del; 
r=crit/(L1*cos(t)^2+L2*sin(t)^2); 
r=sqrt(r); 
f(1)=r*cos(t); 
f(2)=r*sin(t); 
g=U*f; 
x(j)=ybar(1)-g(1); 
y(j)=ybar(2)-g(2);
end;
 h=sqrt(crit)*[-1 0;1 0;0 1;0 -1]; 
 h=1.25*h*inv(Sql)*U';
 xr1=h(1:2,1)+ybar(1);
 xr2=h(3:4,1)+ybar(1);
 yr1=h(1:2,2)+ybar(2);
 yr2=h(3:4,2)+ybar(2);