function result = adjustedoutlyingness(x,varargin)
%ADJUSTEDOUTLYINGNESS calculates the 'Skewness-Adjusted Outlyingness'.
% The method searches for outliers in multivariate skewed data
% (thus without assuming elliptical symmetric).
% It is based on the outlyingness measure of Stahel and Donoho (1981, 1982)
% and on the skewness-adjusted boxplot of Hubert and Vandervieren (2007).
%
% The Skewness-Adjusted Outlyingness is described in:
% Hubert, M., and Van der Veeken, S. (2008),
% "Outlier detection for skewed data",
% Journal of Chemometrics, 22, 235-246.
%
% The method can be useful as preprocessing in
% FASTICA (www.cis.hut.fi/projects/ica/fastica/), see
%
% Brys, G., Hubert, M., and Rousseeuw, P.J. (2005),
% "A Robustification of Independent Component Analysis",
% Journal of Chemometrics, 19, 364-375.
%
%
% Required input arguments:
% x : Data matrix (rows=observations, columns=variables)
%
% Optional input arguments:
% ndir : Number of directions in which de outlyingness will be computed
% (default = 250*number of variables)
% classic : If equal to one, the classical Stahel Donoho outlyingness is
% calculated as well (default = 0)
%
% When the number of observations n is sufficiently large compared to the
% dimension p, i.e. when n > 5*p, directions are generated orthogonal to the
% hyperplane through p random observations. The adjusted outlyingness is
% then affine invariant. When n is small compared to p, we take random
% directions through two data points. This procedure is orthogonal
% invariant.
%
% I/O:
% result=adjustedoutlyingness(x,'ndir',250);
%
% Example:
% x = [chi2rnd(2,1000,1) trnd(3,1000,1)];
% x(1:10,:) = mvnrnd([-2 -3],eye(2)/10,10);
% result = adjustedoutlyingness(x);
%
% The output of ADJUSTEDOUTLYINGNESS is a structure containing:
% result.adjout : skewness-adjusted outlyingness values for all observations
% result.cutoff : cutoff value for the AO-values
% result.flag : The observations whose AO-value exceeds the cutoff value can be
% considered as outliers and receive a flag equal to zero. The regular observations
% receive a flag 1.
% result.classic : If the input argument 'classic' is equal to one, this structure
% contains results of the classical analysis:
% .outl: Stahel-Donoho outlyingness for all observations
% .cutoff: cutoff value for the outlyingness values.
% .flag: The observations whose outlyingness exceeds the cutoff receive flag zero, the others receive flag 1.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust
%
% Written by Guy Brys, Mia Hubert, Stephan Van der Veeken, Tim Verdonck
% Last Update: 09/06/2008
if (nargin<1)
error('Input matrix x is required');
end
[n,p]=size(x);
counter=1;
default=struct('ndir',250*p,'a',-4,'b',3,'classic',0);
list=fieldnames(default);
result=default;
IN=length(list);
i=1;
%reading the user's input
if nargin>1
%
%placing inputfields in array of strings
%
for j=1:nargin-1
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 'result'.
%
while counter<=IN
index=strmatch(list(counter,:),chklist,'exact');
if ~isempty(index) %in case of similarity
for j=1:nargin-1 %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
result=setfield(result,chklist{index},varargin{I+1});
index=[];
end
counter=counter+1;
end
end
ndir=result.ndir;
a=result.a;
b=result.b;
classic=result.classic;
if (ndir<=0)
error('The number of directions should be positive.');
end
if isempty(ndir)
ndir=250*p;
end
% a and b must be numeric scalars, classic should be a binary variable
if isempty(a)
a = -4;
elseif ~isscalar(a) || ~isnumeric(a)
error('The ''a'' parameter value must be a numeric scalar.');
end
if isempty(b)
b = 3;
elseif ~isscalar(b) || ~isnumeric(b)
error('The ''b'' parameter value must be a numeric scalar.');
end
if isempty(classic)
classic = 0;
elseif ~isscalar(classic) || ~ismember(classic,0:1)
error('Invalid value for ''classic'' parameter.');
end
if n>5*p
B=[];
for i =1:ndir
xx = randperm(n);
P = x(xx(1:p),:);
if (rank(P)==p)
B = [B ; (P\ones(p,1))'];
end
end
else
B=twopoints(x,ndir,0);
end
for i=1:size(B,1)
Bnorm(i)=norm(B(i,:),2);
end
Bnormr=Bnorm(Bnorm > 1.e-12);
B=B(Bnorm > 1.e-12,:);
A=diag(1./Bnormr)*B;
%Looking in ndir directions for skewness-adjusted outlyingness
Y=x*A';
[n,p]=size(Y);
YZ=zeros(n,p);
tmc = mc(Y);
h=find(abs(tmc)==1);
if(sum(h)>1)
error('There are too many ties in the data. The adjusted outlyingness can not be computed.')
end
tme = median(Y);
tiq = iqr(Y);
s=find(tiq==0);
if(sum(s)>1)
error('There are too many ties in the data. The adjusted outlyingness can not be computed.')
end
tp1 = prctile(Y,25);
tp3 = prctile(Y,75);
if (tmc>=0)
tup = (tp3+1.5*exp(b*tmc).*tiq)-tme;%3
tlo = tme-(tp1-1.5*exp(a*tmc).*tiq);%-4
else
tup = (tp3+1.5*exp(-a*tmc).*tiq)-tme;%4
tlo = tme-(tp1-1.5*exp(-b*tmc).*tiq);%-3
end
for k = 1:p
ttup = sort(-Y(Y(:,k)<(tup(k)+tme(k)),k));
tup(k) = -ttup(1)-tme(k);
ttlo = sort(Y(Y(:,k)>(tme(k)-tlo(k)),k));
tlo(k) = tme(k)-ttlo(1);
end
tmp = ((Y>=repmat(tme,n,1)).*(repmat(tup,n,1)) + (Y<repmat(tme,n,1)).*(repmat(tlo,n,1)));
YZ = abs((Y-repmat(tme,n,1)))./tmp;
adjout = max(YZ,[],2);
mcadjout=mc(adjout);
if mcadjout>0
cutoff = prctile(adjout,75)+1.5*exp(b*mcadjout)*iqr(adjout);
else
cutoff = prctile(adjout,75)+1.5*iqr(adjout);
end
ttup=sort(-adjout(adjout<cutoff));
cutoff=-ttup(1);
flag=(adjout<=cutoff);
if classic==1
% The classical estimates are computed.
umad=madc(Y);
g=find(umad==0);
if(sum(g)>1)
error('There are too many ties in the data. The Stahel-Donoho outlyingness can not be computed.')
end
clYZ=abs((Y-repmat(tme,n,1)))./repmat(umad,n,1);
classout=max(clYZ,[],2);
%classical chi-square cutoff and flag
quant=chi2inv(0.99,size(x,2));
classcutoff=sqrt(quant);
classflag=(classout<=classcutoff);
classic=struct('outl',classout,'cutoff',classcutoff,'flag',classflag);
end
%Putting things together
result = struct('adjout',adjout,'cutoff',cutoff,'flag',flag,'classic',classic);