function h = ConfidenceEllipse(X,  confidence, conc, lw,lc)
%
% X - n x 2 matrix. It is a sample of size n from bivariate normal
% distribution with mean mu=(mu1, mu2) and covariance Sigma.
%
% confidence - confidence level. Default is 95% CE.
% conc -  0 for confidence ellipse (for the mean mu=(mu1, mu2))
% and 1 for concentration ellipse. Default 0.
% lw - ellipse line width. Default 4.
% lc - ellipse line color. Default black =[0 0 0].
% Output h -- handle
% Example of Use:
%     >  mu = [2 3];
%     > SIGMA = [1 -0.5; -0.5 3];
%     > X = mvnrnd(mu,SIGMA,100);
%     > plot(X(:,1),X(:,2),'r.')
%     > hold on
%     > ConfidenceEllipse(X, 0.95, 0, 4, [1 0 0])
%             %    95% Confidence ellipse for (mu1, m2).
%     > ConfidenceEllipse(X, 0.95, 1, 2, [0 1 0])
%             %    95% Concentration ellipse for data X.
% ==========================================
% BVidakovic 2013
 

if nargin == 1
    confidence=0.95; conc=0; lw=4; lc=[0 0 0]'
end
if nargin ==2
    conc=0; lw=4; lc=[0 0 0];
end
if nargin ==3
   lw=4; lc=[0 0 0];
end
if nargin ==4
   lc=[0 0 0];
end

[n p] = size(X);
Xbar = mean(X);

numPtsOnEllipse = 361; %point at each degree
theta = linspace(0,2*pi,numPtsOnEllipse)';
C = finv(confidence,p,n-p)*p*(n-1)/(n*(n-p));
if conc==1
    C=C*n;  %the concentration ellipse
end
[pcomp  scores  latent] = princomp(X);

K = pcomp * diag(sqrt(C*latent));
parametriceq = [cos(theta),sin(theta)]*K' + ...
           repmat(Xbar,length(theta),1);
xtheta =parametriceq (:,1);  %x(theta)
ytheta =parametriceq (:,2);  %y(theta)

% Add ellipse to current plot
hand = line(xtheta, ytheta,'Clipping','off', 'LineWidth',lw , ...
'Color', lc);
if nargout > 0
    h  = hand;
end
