function result=rda(x,group,varargin)
%RDA performs linear and quadratic robust discriminant analysis
% on the data matrix x with known group structure. It is based on the
% MCD estimator (see mcdcov.m), hence it has to be applied to
% low-dimensional data.
%
% The Robust Discriminant method is described in:
% Hubert, M., Van Driessen, K. (2004),
% "Fast and Robust Discriminant Analysis,"
% Computational Statistics and Data Analysis, 45, 301-320.
%
% Required input arguments:
% x : training data set (matrix of size n by p).
% group : column vector containing the group numbers of the training
% set x. For the group numbers, any strict positive integer is
% allowed assuming that the first group is the one with the smallest group number.
%
% Optional input arguments:
% alpha : (1-alpha) measures the fraction of outliers the MCD-algorithm should
% resist. Any value between 0.5 and 1 may be specified. (default = 0.75)
% method : String which indicates whether a 'linear' (default) or 'quadratic'
% discriminant rule should be applied
% misclassif : String which indicates how to estimate the probability of
% misclassification. It can be based on the
% training data ('training'), a validation set ('valid'),
% or cross-validation ('cv'). Default is 'training'.
% membershipprob : Vector which contains the membership probability of each
% group (sorted by increasing group number). If no priors are given, they are estimated as the
% proportions of regular observations in the training set.
% valid : If misclassif was set to 'valid', this field should contain
% the validation set (a matrix of size m by p).
% groupvalid : If misclassif was set to 'valid', this field should contain the group numbers
% of the validation set (a column vector).
% predictset : Contains a new data set (a matrix of size mp by p) from which the
% class memberships are unknown and should be predicted.
% plots : If equal to 1, one figure is created with the training data and the
% MCD tolerance ellipses for each group. This plot is
% only available for bivariate data sets. For technical reasons, a maximum
% of 6 groups is allowed. Default is one.
% classic : If equal to one, classical linear or quadratic discriminant analysis will be performed
% (see also cda.m). (default = 0)
% compare : If equal to one, the classical CDA analysis will be performed
% with the same weights and the same priors as the robust analysis
% has been performed. This is especially useful to compare the robust
% and classical result on the same data with the same priors. (default = 0)
%
% I/O: result=rda('alpha',0.5,'plots',0,'misclassif','training','method','linear',...
% 'membershipprob',proportions,'valid',y,'groupvalid',groupy,'classic',0);
% The user should only give the input arguments that have to change their default value.
% The name of the input arguments needs to be followed by their value.
% The order of the input arguments is of no importance.
%
% Examples: out=rda(x,group,'method','linear')
% out=rda(x,group,'plots',0)
% out=rda(x,group,'valid',y,'groupvalid',groupy)
%
% The output is a structure containing the following fields:
%
% result.assignedgroup : If there is a validation set, this vector contains the assigned group numbers
% for the observations of the validation set. Otherwise it contains the
% assigned group numbers of the original observations based on the discriminant rules.
% result.scores : If there is a validation set, this columnvector of size m contains the maximal discriminant
% scores for each observation from the validation set. Otherwise it is a columnvector of size n
% containing the maximal discriminant scores of the training set.
% result.method : String containing the method used to obtain the discriminant rules (either 'linear' or 'quadratic').
% This is the same as the input argument method.
% result.cov : If method equals 'linear', this is a matrix containing the estimated common covariance matrix.
% If method equals 'quadratic', it is a cell array containing the covariances per group.
% result.center : A vector in which the rows contain the estimated centers of the groups.
% result.rd : A vector of length n containing the robust distances of each observation from the training set
% to the center of its group.
% result.flagtrain : Observations from the training set whose robust distance exceeds a certain cut-off value
% can be considered as outliers and receive a flag equal to zero.
% The regular observations receive a flag 1. (See also mcdcov.m)
% result.flagvalid : Observations from the validation set whose robust distance (to the center of their group)
% exceeds a certain cut-off value can be considered as outliers and receive a
% flag equal to zero. The regular observations receive a flag 1.
% If there is no validation set, this field is equal to zero.
% result.grouppredict : If there is a prediction set, this vector contains the assigned group numbers
% for the observations of the prediction set.
% result.flagpredict : Observations from the new data set (predict) whose robust distance (to the center of their group)
% exceeds a certain cut-off value can be considered as overall outliers and receive a
% flag equal to zero. The regular observations receive a flag 1.
% If there is no prediction set, this field is equal to zero.
% result.membershipprob : A vector with the membership probabilities.
% result.misclassif : String containing the method used to estimate the misclassification probabilities
% (same as the input argument misclassif)
% result.groupmisclasprob : A vector containing the misclassification probabilities for each group.
% result.avemisclasprob : Overall probability of misclassification (weighted average of the misclassification
% probabilities over all groups).
% result.class : 'RDA'
% result.classic : If the input argument 'classic' is equal to one, this structure
% contains results of the classical discriminant analysis (see also cda.m).
% result.compare : If the input argument 'compare' is equal to one, this strucuture
% contains results for the classical discriminant analysis with the same weights
% and priors as in the robust analysis.
% result.x : The training data set (same as the input argument x).
% result.group : The group numbers of the training set (same as the input argument group).
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust.html
%
% Written by Nele Smets and Sabine Verboven on 01/03/2004
% Last Update: 01/07/2005
%
if nargin<2
error('There are too few input arguments.')
end
% assigning default-values
[n,p]=size(x);
if size(group,1)~=1
group=group';
end
if n ~= length(group)
error('The number of observations is not the same as the length of the group vector!')
end
g=group;
countsorig=tabulate(g); %contingency table (outputmatrix with 3 colums): value - number - percentage
[lev,levi,levj]=unique(g);
%Redefining the group number
if any(lev~= (1:length(lev)))
lev=1:length(lev);
g=lev(levj);
counts=tabulate(g);
else
counts=countsorig;
end
if ~all(counts(:,2)) %some groups have zero values, omit those groups
disp(['Warning: group(s) ', num2str(counts(counts(:,2)==0,1)'), 'are empty']);
empty=counts(counts(:,2)==0,:);
counts=counts(counts(:,2)~=0,:);
else
empty=[];
end
if any(counts(:,2)<5)%some groups have less than 5 observations
error(['Group(s) ', num2str(counts(counts(:,2)<5,1)'), ' have less than 5 observations.']);
end
proportions = zeros(size(counts,1),1);
y=0; %initial values of the validation data set and its groupsvector
groupy=0;
counter=1;
default=struct('alpha',0.75,'plots',1,'misclassif','training','method','linear','membershipprob',proportions,...
'valid',y,'groupvalid',groupy,'classic',0,'compare',0,'predictset',[]);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%reading the user's input
if nargin>2
%
%placing inputfields in array of strings
%
for j=1:nargin-2
if rem(j,2)~=0
chklist{i}=varargin{j};
i=i+1;
end
end
%
%Checking which default parameters have to be changed
% and keep them in the structure 'options'.
%
while counter<=IN
index=strmatch(list(counter,:),chklist,'exact');
if ~isempty(index) %in case of similarity
for j=1:nargin-2 %searching the index of the accompanying field
if rem(j,2)~=0 %fieldnames are placed on odd index
if strcmp(chklist{index},varargin{j})
I=j;
end
end
end
options=setfield(options,chklist{index},varargin{I+1});
index=[];
end
counter=counter+1;
end
end
%Checking prior (>0 )
prior=options.membershipprob;
if size(prior,1)~=1
prior=prior';
end
epsilon=10^-4;
if sum(prior) ~= 0
if (any(prior < 0) | (abs(sum(prior)-1)) > epsilon)
error('Invalid membership probabilities.')
end
end
ng=length(proportions);
if length(prior)~=ng
error('The number of membership probabilities is not the same as the number of groups.')
end
%%%%%%%%%%%%%%%%%%MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%
%Checking if a validation set is given
if strmatch(options.misclassif, 'valid','exact')
if options.valid==0
error(['The misclassification error will be estimated through a validation set',...
'but no validation set is given!'])
else
validx = options.valid;
validgrouping = options.groupvalid;
if size(validx,1)~=length(validgrouping)
error('The number of observations in the validation set is not the same as the length of its group vector!')
end
if size(validgrouping,1)~=1
validgrouping = validgrouping';
end
countsvalidorig=tabulate(validgrouping);
countsvalid=countsvalidorig(countsvalidorig(:,2)~=0,:);
if size(countsvalid,1)==1
error('The validation set must contain observations from more than one group!')
elseif any(ismember(empty,countsvalid(:,1)))
error(['Group(s) ' ,num2str(empty(ismember(empty,countsvalid(:,1)))), 'was/were empty in the original dataset.'])
end
end
elseif options.valid~=0
validx = options.valid;
validgrouping = options.groupvalid;
if size(validx,1) ~= length(validgrouping)
error('The number of observations in the validation set is not the same as the length of its group vector!')
end
if size(validgrouping,1)~=1
validgrouping = validgrouping';
end
options.misclassif='valid';
countsvalidorig=tabulate(validgrouping);
countsvalid=countsvalidorig(countsvalidorig(:,2)~=0);
if size(countsvalid,1)==1
error('The validation set must contain more than one group!')
elseif any(ismember(empty,countsvalid(:,1)))
error(['Group(s) ' , num2str(empty(ismember(empty,countsvalid(:,1)))), ' was/were empty in the original dataset.'])
end
end
%Discriminant rule based on the training set x
result1 = rawrule(x, g, prior, options.alpha, options.method);
%Apply discriminant rule on validation set
if strmatch(options.misclassif,'valid','exact')
result2 = rewrule(validx, result1);
finalgroup = result2.class;
else
result2 = rewrule(x, result1);
finalgroup = result2.class;
end
%Estimating the misclassification error
switch options.misclassif
case 'valid'
[v,vi,vj]=unique(validgrouping);
%Redefining the group number
if any(v~= (1:length(v)))
v=1:length(v);
validgrouping=v(vj);
end
if any(countsvalidorig(:,2)==0)
empty=setdiff(countsvalidorig(find(countsvalidorig(:,2)==0),1), countsorig(find(countsorig(:,2)==0)));
disp(['Warning: the test group(s) ' , num2str(empty), ' are empty']);
else
empty=[];
end
misclas=-ones(1,length(lev));
for i=1:size(validx,1)
if strmatch(options.method,'quadratic','exact')
dist(i) = libra_mahalanobis(validx(i,:), result1.center(vj(i),:),'invcov',result1.invcov{vj(i)});
else
dist(i) = libra_mahalanobis(validx(i, :), result1.center(vj(i),:),'invcov',result1.invcov);
end
end
weightsvalid=zeros(1,length(dist));
weightsvalid(dist <= chi2inv(0.975, p))=1;
for i=1:length(v)
if ~isempty(intersect(i,v))
misclas(i)=sum((validgrouping(weightsvalid==1)==finalgroup(weightsvalid==1)') & (validgrouping(weightsvalid==1)==repmat(lev(i),1,sum(weightsvalid))));
ingroup(i) = sum((validgrouping(weightsvalid == 1) == repmat(lev(i),1, sum(weightsvalid))));
misclas(i) = 1 - (misclas(i)./ingroup(i));
end
end
if any(misclas==-1)
misclas(misclas==-1)=0;
end
misclasprobpergroup=misclas;
misclas=misclas.*result1.prior;
misclasprob=sum(misclas);
case 'training'
for i=1:ng
misclas(i) = sum((g(result1.weights==1)==finalgroup(result1.weights==1)')&(g(result1.weights==1)==repmat(lev(i),1,sum(result1.weights))));
ingroup(i) = sum((g(result1.weights == 1) == repmat(lev(i),1,sum(result1.weights))));
end
misclas = (1 - (misclas./ingroup));
misclasprobpergroup = misclas;
misclas = misclas.*result1.prior;
misclasprob = sum(misclas);
weightsvalid=0;%only available with validation set
case 'cv'
finalgroup=[];
for i=1:length(x)
if (result1.weights(i) == 1)
xnew=removal(x,i,0);
groupnew=removal(group,0,i);
functie1res = rawrule(xnew, groupnew,prior,options.alpha, options.method);
functie2res = rewrule(x(i, :), functie1res);
finalgroup = [finalgroup; functie2res.class(1)];
end
end
for i=1:ng
misclas(i) = sum((g(result1.weights == 1) == finalgroup') & (g(result1.weights == 1) == repmat(lev(i),1,sum(result1.weights))));
ingroup(i) = sum(g(result1.weights == 1) == repmat(lev(i),1,sum(result1.weights)));
end
misclas = (1 - (misclas./ingroup));
misclasprobpergroup= misclas;
misclas = misclas.* result1.prior;
misclasprob = sum(misclas);
weightsvalid=0; %only available with validation set
end
%classify the new observations (predict)
if ~isempty(options.predictset)
resultpredict = rewrule(options.predictset, result1);
finalgrouppredict = resultpredict.class;
for i=1:size(options.predictset,1)
for j = 1:ng
if strmatch(options.method,'quadratic','exact')
distpredict(i,j) = libra_mahalanobis(options.predictset(i,:), result1.center(j,:),'invcov',result1.invcov{j});
else
distpredict(i,j) = libra_mahalanobis(options.predictset(i, :), result1.center(j,:),'invcov',result1.invcov);
end
end
end
weightspredict = zeros(1,size(distpredict,1));
weightspredict(min(distpredict,[],2) <= chi2inv(0.975, p))=1;
else
finalgrouppredict = 0;
weightspredict = 0;
end
if options.classic
classicout=cda(x,g,'method',result2.method,'misclassif',options.misclassif,'membershipprob',options.membershipprob,'valid',options.valid,...
'groupvalid',options.groupvalid,'plots',0,'predictset',options.predictset);
else
classicout=0;
end
if options.compare
compareout=cda(x,g,'method',result2.method,'misclassif',options.misclassif,'membershipprob',result1.prior,'valid',options.valid,...
'groupvalid',options.groupvalid,'plots',0,'predictset',options.predictset,'weightstrain',result1.weights,'weightsvalid',weightsvalid);
else
compareout=0;
end
%Output structure
result=struct('assignedgroup',{finalgroup'},'scores',{result2.scores'},'method',{result2.method},'cov',{result1.cov}, ...
'center',{result1.center},'rd',{result1.dist'},'flagtrain',{result1.weights},...
'flagvalid',weightsvalid,'grouppredict',finalgrouppredict,'flagpredict',weightspredict','membershipprob',{result1.prior},...
'misclassif',{options.misclassif},'groupmisclasprob',{misclasprobpergroup},'avemisclasprob',{misclasprob},...
'class',{'RDA'},'classic',{classicout},'compare',{compareout},'x',{x},'group',{group});
if size(x,2)~=2
result=rmfield(result,{'x','group'});
end
%Plotting the output
try
if options.plots & options.classic
makeplot(result,'classic',1)
elseif options.plots
makeplot(result)
end
catch %output must be given even if plots are interrupted
%> delete(gcf) to get rid of the menu
end
%--------------------------------------------------------------------------
function result=rawrule(x, g,prior, alfa, method)
%computes the discrimination rule based on the training set x.
[n,p]=size(x);
epsilon=10^-4;
counts=tabulate(g); %contingency table (outputmatrix with 3 colums): value - number - percentage
[lev,levi,levj]=unique(g);
if ~all(counts(:,2)) %some groups have zero values, omit those groups
empty=counts(counts(:,2)==0,1);
else
empty=[];
end
ng=size(counts,1);
switch method
case 'linear' %equal covariances supposed
[gun,gi,gj]=unique(g);
for j=1:length(gun)
group.mcd{j}=mcdcov(x(g==gun(j),:),'alpha',alfa,'plots',0); %covariance of group j
group.center(j,:)=group.mcd{j}.center; %center of all groups, matrix of ng x p
end
for i=1:n
zgeg(i,:)=x(i,:)-group.center(gj(i),:);
end
zmcd = mcdcov(zgeg,'alpha',alfa,'plots',0);
zmcdcenter = zmcd.center;
zmcdcov = zmcd.cov;
zgeg = zgeg - repmat(zmcdcenter,length(zgeg),1);
group.center = group.center + repmat(zmcdcenter,size(group.center,1),1);
dist=zeros(n,1);
for j=1:length(gun)
dist(g==gun(j))=libra_mahalanobis(x(g==gun(j),:),group.center(j,:),'invcov',inv(zmcd.cov));
end
weights=zeros(n,1);
weights(dist <= chi2inv(0.975,p))=1;
result.cov=zmcd.cov; %over all group
result.invcov=inv(zmcd.cov);
result.center=group.center; %all groups
result.weights=weights';
result.dist=dist;
result.method=method;
case 'quadratic'
[gun,gi,gj]=unique(g);
xmcdweights=zeros(n,1);
for j=1:length(gun)
[group.mcd{j} raw{j}]=mcdcov(x(g==gun(j),:),'alpha',alfa,'plots',0);
group.cov{j}=group.mcd{j}.cov; %covariance of group j
group.invcov{j}=inv(group.cov{j});
group.center(j,:)=group.mcd{j}.center; %center of all groups
xmcdweights(g==gun(j))=raw{j}.wt;
end
for i=1:n
xdist(i)=libra_mahalanobis(x(i,:), group.center(gj(i),:), 'invcov',group.invcov{gj(i)});
end
weights=xmcdweights;
result.cov=group.cov; %per group
result.invcov=group.invcov;
result.center = group.center; %all groups
result.weights = xmcdweights';
result.dist = xdist';
result.method = method;
end
%Define the prior
if sum(prior) ~= 0
result.prior = prior;
else
ngood=sum(weights);
%regular points are kept
ggood = g(weights==1);
countsgood=tabulate(ggood);
if empty
for i=1:length(empty)
countsgood(countsgood(:,1)==empty(i),:)= [];
end
end
if ~any(countsgood(:,2))
disp(['Warning: the group(s) ', num2str(countsgood(countsgood(:,2) == 0,1)'), 'contain only outliers']);
countsgood=countsgood(countsgood(:,2)~=0,:);
end
result.prior = (countsgood(:,3)/100)';
end
%--------------------------------------------------------------------------
function result=rewrule(x, rawobject)
epsilon=10^-4;
center=rawobject.center;
covar=rawobject.cov;
invcov=rawobject.invcov;
prior=rawobject.prior;
method=rawobject.method;
if (length(prior) == 0 | length(prior) ~= size(center,1))
error('invalid prior')
end
if sum(prior) ~= 0
if (any(prior < 0) | (abs(sum(prior)-1)) > epsilon)
error('invalid prior')
end
end
ngroup=length(prior);
[n,p]=size(x);
switch method
case 'linear'
for j=1:ngroup
for i=1:n
scores(i,j) = linclassification(x(i,:)', center(j,:)', invcov, prior(j));
end
end
[maxs,maxsI] = max(scores,[],2);
for i=1:n
maxscore(i,1) = scores(i,maxsI(i));
end
result.scores = maxscore;
result.class = maxsI;
result.method = method;
case 'quadratic'
for j=1:ngroup
for i=1:n
scores(i, j) = classification(x(i,:)', center(j,:)', covar{j}, invcov{j}, prior(j));
end
end
[maxs,maxsI] = max(scores,[],2);
for i=1:n
maxscore(i,1) = scores(i,maxsI(i));
end
result.scores = maxscore;
result.class = maxsI;
result.method = method;
end
%--------------make sure the input variables are column vectors!
function out=classification(x, center, covar,invcov, priorprob)
out=-0.5*log(abs(det(covar)))-0.5*(x - center)' * invcov *(x - center)+log(priorprob);
%-------------------
function out=linclassification(x, center, invcov, priorprob)
out=center'*invcov*x - 0.5*center'*invcov*center+log(priorprob);