function ellipsplot(center,covar,data,dist,nid,labx,laby)

%ELLIPSPLOT plots the 97.5% tolerance ellipse of the bivariate data set
% (data). The ellipse is defined by those data points whose distance (dist) 
% is equal to the squareroot of the 97.5% chisquare quantile with 2 degrees of
% freedom. 
%  
% Required input arguments:
%   center  : estimate of the center of the data set
%    covar  : estimate of the covariance matrix of the data set 
%     data  : the two-dimensional data matrix
%     dist  : distance needed to flag data points outside the ellipse
%
% Optional input arguments:
%      nid  : number of data points with largest distance
%             to be identified (default value: 3)
%     labx  : a label for the x axis (default value: 'X1')
%     laby  : a label for the y axis (default value: 'X2')
%
% I/O: ellipsplot(center,covar,data,dist,nid,labx,laby)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at: 
%              http://wis.kuleuven.be/stat/robust.html
%
% Last update: 23/10/2003

set(gcf,'Name', '97.5% Tolerance ellipse', 'NumberTitle', 'off');
n=length(dist);
if size(data,2)~=2
     disp('The tolerance ellipse is only drawn for two-dimensional data sets')
else
    if nargin==4
        nid=3;
        labx='X1';
        laby='X2';
    elseif nargin==5
        labx='X1';
        laby='X2';
    elseif nargin==6
        laby='X2';
    end
     deter=covar(1,1)*covar(2,2)-covar(1,2)^2;
     ylimit=sqrt(7.37776*covar(2,2));
     y=-ylimit:0.005*ylimit:ylimit;
     sqtdi=sqrt(deter*(ylimit^2-y.^2))/covar(2,2);
     sqtdi([1,end])=0;
     b=center(1)+covar(1,2)/covar(2,2)*y;
     x1=b-sqtdi;
     x2=b+sqtdi;
     y=center(2)+y;
     ellip=[x1,x2([end:-1:1]);y,y([end:-1:1])]';
     xmin=min([data(:,1);ellip(:,1)]);
     xmax=max([data(:,1);ellip(:,1)]);
     ymin=min([data(:,2);ellip(:,2)]);
     ymax=max([data(:,2);ellip(:,2)]);
     xmarg=0.05*abs(xmax-xmin);
     ymarg=0.05*abs(ymax-ymin);
     xmin=xmin-xmarg;
     xmax=xmax+xmarg;
     ymin=ymin-ymarg;
     ymax=ymax+ymarg;
     x1=data(:,1)';
     x2=data(:,2)';
     plot(x1,x2,'o')
     xlabel(labx)
     ylabel(laby)
     xlim([xmin,xmax]);
     ylim([ymin,ymax]);
     box on
     [dist,ind]=sort(dist);
     ind=ind(n-nid+1:n)';
     text(x1(ind),x2(ind),int2str(ind));   
     title('Tolerance ellipse (97.5%)');         
     line(ellip(:,1),ellip(:,2),'Color','r');
 end