function daplot(x,group,center,covar,classic,method)
%DAPLOT plots 97.5% tolerances ellipses of the bivariate data set x, which
% consists of several groups defined by the input argument 'group'.
% Center and covar are estimates of the center and covariance matrix of each group,
% obtained with a classical ('CDA') or a robust ('RDA') discriminant analysis. See cda.m and rda.m
% Method indicates whether linear ('linear') or quadratic ('quadratic')
% classification is used.
%
% For technical reasons, 6 different groups can be plotted (with different symbols).
% In case there are more groups, please adapt lines 22 and 23.
%
% I/O: daplot(x,group,center,cov,'CDA','linear')
%
% Written by S. Verboven on 28/07/03
% Last update: 17/02/2004
set(gcf,'Name', 'Discriminant analysis', 'NumberTitle', 'off');
n = size(x,1);
[m,p] = size(center);
lev=unique(group);
marking={'bx';'ro';'kdiamond';'y+';'g.';'m*'};
linecolor={'b-','r-','k-','y-','g-','m-'};
if p~=2
disp('Warning: Tolerance ellipses are only drawn for two-dimensional data sets.')
return
elseif length(lev)>6
error(['Only 6 groups can be drawn with different symbols. Please adapt the code',...
' of daplot.m if you want to draw more groups.'])
else
legstr=[];
for i=1:m
data{i}=x(group==lev(i),:);
plot(data{i}(:,1),data{i}(:,2),marking{i});
hold on
legstr=[legstr; sprintf(['Group',num2str(i)])];
end
for i=1:m
if strcmp(method,'linear')
Z{i} = ellipse(center(i,:),covar);
elseif strcmp(method,'quadratic')
Z{i} = ellipse(center(i,:),covar{i});
end
plot(Z{i}(:, 1), Z{i}(:, 2), linecolor{i})
hold on
end
end
if strmatch('CDA',classic,'exact') & strmatch('linear',method,'exact')
title('Classical linear discriminant analysis')
elseif strmatch('CDA',classic,'exact') & strmatch('quadratic',method,'exact')
title('Classical quadratic discriminant analysis')
elseif strmatch('RDA',classic,'exact') & strmatch('linear',method,'exact')
title('Robust linear discriminant analysis')
elseif strmatch('RDA',classic,'exact') & strmatch('quadratic',method,'exact')
title('Robust quadratic discriminant analysis')
end
hold off
legend(legstr,0)
%--------------
function out = ellipse(loc,covar)
A=covar;
detA = A(1, 1)*A(2, 2) - A(1, 2)^2;
dist = sqrt(chi2inv(0.975, 2));
ylimit = sqrt(A(2, 2)) * dist;
y = - ylimit:0.01*ylimit:ylimit;
sqroot.discr = sqrt(detA/A(2, 2).^2 * (A(2, 2) * dist.^2 - y.^2));
sqroot.discr([1, length(sqroot.discr)]) = 0;
b = loc(1) + A(1, 2)/A(2, 2) * y;
x1 = b - sqroot.discr;
x2 = b + sqroot.discr;
y = loc(2) + y;
out= [[x1' y']; flipud([x2' y'])];