function result=robpca(x,varargin)
%ROBPCA is a 'ROBust method for Principal Components Analysis'.
% It is resistant to outliers in the data. The robust loadings are computed
% using projection-pursuit techniques and the MCD method.
% Therefore ROBPCA can be applied to both low and high-dimensional data sets.
% In low dimensions, the MCD method is applied (see mcdcov.m).
%
% The ROBPCA method is described in
% Hubert, M., Rousseeuw, P.J., Vanden Branden, K. (2005), ROBPCA: a
% new approach to robust principal components analysis, Technometrics, 47, 64-79.
%
%
% To select the number of principal components, a robust PRESS (predicted
% residual sum of squares) curve is drawn, based on a fast algorithm for
% cross-validation. This approach is described in:
%
% Hubert, M., Engelen, S. (2007),
% "Fast cross-validation of high-breakdown resampling algorithms for PCA",
% Computational Statistics and Data Analysis, 51, 5013-5024.
%
% ROBPCA is designed for normally distributed data. If the data are skewed,
% a modification of ROBPCA is available, based on the adjusted outlyingness
% (see adjustedoutlyingness.m). This method is described in:
%
% Hubert, M., Rousseeuw, P.J., Verdonck, T. (2008),
% "Robust PCA for skewed data and its outlier map", Computational Statistics
% and Data Analysis, in press.
%
% For the up-to-date reference, please consult the website:
% wis.kuleuven.be/stat/robust.html
%
% Required input arguments:
% x : Data matrix (observations in the rows, variables in the
% columns)
%
% Optional input arguments:
% k : Number of principal components to compute. If k is missing,
% or k = 0, a scree plot and a press curve are drawn which allows you to select
% the number of principal components.
% kmax : Maximal number of principal components to compute (default = 10).
% If k is provided, kmax does not need to be specified, unless k is larger
% than 10.
% alpha : (1-alpha) measures the fraction of outliers the algorithm should
% resist. Any value between 0.5 and 1 may be specified (default = 0.75).
% h : (n-h+1) measures the number of outliers the algorithm should
% resist. Any value between n/2 and n may be specified. (default = 0.75*n)
% Alpha and h may not both be specified.
% mcd : If equal to one: when the number of variables is sufficiently small,
% the loadings are computed as the eigenvectors of the MCD covariance matrix,
% hence the function 'mcdcov.m' is automatically called. The number of
% principal components is then taken as k = rank(x). (default)
% If equal to zero, the robpca algorithm is always applied.
% plots : If equal to one, a scree plot, a press curve and a robust score outlier map are
% drawn (default). If the input argument 'classic' is equal to one,
% the classical plots are drawn as well.
% If 'plots' is equal to zero, all plots are suppressed (unless k is missing,
% then the scree plot and press curve are still drawn).
% See also makeplot.m
% labsd : The 'labsd' observations with largest score distance are
% labeled on the outlier map. (default = 3)
% labod : The 'labod' observations with largest orthogonal distance are
% labeled on the outlier map. default = 3)
% classic : If equal to one, the classical PCA analysis will be performed
% (see also cpca.m). (default = 0)
% scree : If equal to one, a scree plot is drawn. If k is given as input, the default value is 0, else the default value is one.
% press : If equal to one, a plot of robust press-values is drawn.
% If k is given as input, the default value is 0, else the default value is one.
% If the input argument 'skew' is equal to one, no plot is
% drawn.
% robpcamcd : If equal to one (default), the whole robpca procedure is run (computation of outlyingness and
% MCD).
% If equal to zero, the program stops after the computation of the outlyingness. The
% robust eigenvectors then correspond with the eigenvectors of the covariance matrix
% of the h observations with smallest outlyingness. This yields the same
% PCA subspace as the full robpca, but not the same eigenvectors and eigenvalues.
% skew : If equal to zero the regular robpca is run. If equal to
% one, the adjusted robpca algorithm for skewed data is run.
%
% I/O: result=robpca(x,'k',k,'kmax',10,'alpha',0.75,'h',h,'mcd',1,'plots',1,'labsd',3,'labod',3,'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:
% result=robpca(x,'k',3,'alpha',0.65,'plots',0)
% result=robpca(x,'alpha',0.80,'kmax',15,'labsd',5)
%
% The output of ROBPCA is a structure containing
%
% result.P : Robust loadings (eigenvectors)
% result.L : Robust eigenvalues
% result.M : Robust center of the data
% result.T : Robust scores
% result.k : Number of (chosen) principal components
% result.kmax : Maximal number of principal components
% result.alpha : see interpretation in the list of input arguments
% result.h : The quantile h used throughout the algorithm
% result.Hsubsets : A structure that contains H0, H1 and Hfreq:
% H0 : The h-subset that contains the h points with the smallest outlyingness.
% H1 : The optimal h-subset of mcdcov.
% Hfreq : The subset of h points which are the most frequently selected during the mcdcov
% algorithm.
% result.sd : Robust score distances within the robust PCA subspace
% result.od : Orthogonal distances to the robust PCA subspace
% result.cutoff : Cutoff values for the robust score and orthogonal distances
% result.flag : The observations whose score distance is larger than result.cutoff.sd (==> result.flag.sd)
% or whose orthogonal distance is larger than result.cutoff.od (==> result.flag.od)
% can be considered as outliers and receive a flag equal to zero (result.flag.all).
% The regular observations receive a flag 1.
% result.class : 'ROBPCA'
% result.classic : If the input argument 'classic' is equal to one, this structure
% contains results of the classical PCA analysis (see also cpca.m).
%
% Short description of the method:
%
% Let n denote the number of observations, and p the number of original variables,
% then ROBPCA finds a robust center (p x 1) of the data M and a loading matrix P which
% is (p x k) dimensional. Its columns are orthogonal and define a new coordinate
% system. The scores (n x k) are the coordinates of the centered observations with
% respect to the loadings: T=(X-M)*P.
% Note that ROBPCA also yields a robust covariance matrix (often singular) which
% can be computed as
% cov=out.P*out.L*out.P'
%
% To select the number of principal components, it is useful to look at the scree plot which
% shows the eigenvalues, and the press curve which displays a weighted sum of the squared
% cross-validated orthogonal distances.
% The outlier map visualizes the observations by plotting their orthogonal
% distance to the robust PCA subspace versus their robust distances
% within the PCA subspace. This allows to classify the data points into 4 types:
% regular observations, good leverage points, bad leverage points and
% orthogonal outliers.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
% http://wis.kuleuven.be/stat/robust
%
% Written by Mia Hubert, Sabine Verboven, Karlien Vanden Branden, Sanne Engelen, Tim Verdonck
% Last Update: 17/06/2003, 03/07/2006, 31/07/2007
% Last Revision: 27/03/2008, 09/06/2008
%
% initialization with defaults
%
data=x;
[n,p]=size(data);
% First Step: classical PCA on data
if n < p
[P1,T1,L1,r,Xc,clm]=kernelEVD(data);
else
[P1,T1,L1,r,Xc,clm]=classSVD(data);
end
% dim(P1): p x r
if r==0
error('All data points collapse!')
end
niter=100;
counter=1;
kmax=min([10,floor(n/2),r]);
k=0;
alfa=0.75;
h=min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
labsd=3;
labod=3;
plots=1;
scree = 1;
press = 1;
mcd=1; % user wants the mcd approach (in case n>>p)
robpcamcd = 1;
cutoff = 0.975;
skew=0;
% default is a structure needed for input checking
default=struct('alpha',alfa,'h',h,'labsd',labsd,'labod',labod,...
'k',k,'plots',plots,'kmax',kmax,'mcd',mcd,'classic',0,'scree',scree,'press',press,'robpcamcd',robpcamcd,'cutoff',cutoff,'skew',0);
list=fieldnames(default);
options=default; %input by user
IN=length(list);
i=1;
%
if nargin==2
error('Incorrect number of input arguments!')
end
dummy = 0; %Assume we didn't get h or alpha, unless we find it below.
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
dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha')); %checking if h and alpha are provided both or not
switch dummy
case 0 % Take on default values
options.alpha=alfa;
if any(strcmp(chklist,'kmax'))
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('kmax',varargin{j})
I=j;
end
end
end
options=setfield(options,'kmax',varargin{I+1});
kmax = max(min([floor(options.kmax),floor(n/2),r]),1); %acceptable kmax
options.h=min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n); %depends on kmax, so if kmax is given by user, get it first!
else %kmax is not given by user
options.h=h;
end
case 3
error('Both inputarguments alpha and h are provided. Only one is required.')
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
options.h=floor(options.h);
options.kmax=floor(options.kmax);
options.k=floor(options.k);
kmax=max(min([options.kmax,floor(n/2),r]),1);
labod=max(0,min(floor(options.labod),n));
labsd=max(0,min(floor(options.labsd),n));
k=options.k;
if k<0
k=0;
elseif k > kmax
k=kmax;
mess=sprintf(['Attention (robpca.m): The number of principal components, k = ',num2str(options.k)...
,'\n is larger than kmax= ',num2str(kmax),'; k is set to ',num2str(kmax)]);
disp(mess)
end
if dummy==1 % checking input variable h
options.alpha=options.h/n;
if k==0
if options.h < floor((n+kmax+1)/2 )
options.h=floor((n+kmax+1)/2);
options.alpha=options.h/n;
mess=sprintf(['Attention (robpca.m): h should be larger than (n+kmax+1)/2.\n',...
'It is set to its minimum value ',num2str(options.h)]);
disp(mess)
end
else
if options.h < floor((n+k+1)/2)
options.h=floor((n+k+1)/2);
options.alpha=options.h/n;
mess=sprintf(['Attention (robpca.m): h should be larger than (n+k+1)/2.\n',...
'It is set to its minimum value ',num2str(options.h)]);
disp(mess)
end
end
if options.h > n
options.alpha=0.75;
if k==0
options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
else
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
end
mess=sprintf(['Attention (robpca.m): h should be smaller than n. \n',...
'It is set to its default value ',num2str(options.h)]);
disp(mess)
end
elseif dummy==2 %checking input variable alpha
if options.alpha < 0.5
options.alpha=0.5;
mess=sprintf(['Attention (robpca.m): Alpha should be larger than 0.5.\n',...
'It is set to 0.5.']);
disp(mess)
end
if options.alpha > 1
options.alpha=0.75;
mess=sprintf(['Attention (robpca.m): Alpha should be smaller than 1. \n',...
'It is set to 0.75.']);
disp(mess)
end
if k==0
options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
else
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
end
end
alfa=options.alpha;
dummyh = strcmp(chklist,'h');
dummykmax = strcmp(chklist,'kmax');
% if all(dummyh == 0) & any(dummykmax) & k==0
% h = min(floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*alfa),n);
% end
if all(dummyh == 0)&& any(dummykmax) %kmax was given by the user
if k==0
options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
else
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
end
elseif all(dummyh == 0) && ~any(dummykmax) %kmax is the default value
if k==0
options.h=floor(2*floor((n+kmax+1)/2)-n+2*(n-floor((n+kmax+1)/2))*options.alpha);
else
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
end
end
h=options.h;
dummyscree = strcmp(chklist,'scree');
dummypress = strcmp(chklist,'press');
if all(dummyscree == 0)
if k~=0
options.scree = 0;
end
end
if all(dummypress == 0)
if k~=0
options.press = 0;
end
end
scree = options.scree;
press = options.press;
labsd=floor(max(0,min(options.labsd,n)));
labod=floor(max(0,min(options.labod,n)));
plots=options.plots;
mcd=options.mcd;
robpcamcd = options.robpcamcd;
cutoff = options.cutoff;
skew = options.skew;
if skew==1
press=0; %press curve not available for skewed data
end
end
%
% MAIN PART
%
X=T1;
center=clm;
rot=P1;
% Depending on n and p, perform MCD or ROBPCA:
% p << n => MCD
p1=size(X,2);
if p1<=min(floor(n/5),kmax) && mcd && (skew==0)
options.h=h;
[res,raw]=mcdcov(X,'h',h,'plots',0);
[U,S,P]=svd(res.cov,0);
L=diag(S);
if k~=0
options.k=min(k,p1);
else
bdwidth=5;
topbdwidth=30;
set(0,'Units','pixels');
scnsize=get(0,'ScreenSize');
pos1=[bdwidth, 1/3*scnsize(4)+bdwidth, scnsize(3)/2-2*bdwidth, scnsize(4)/2-(topbdwidth+bdwidth)];
pos2=[pos1(1)+scnsize(3)/2, pos1(2), pos1(3), pos1(4)];
if press == 1
outcvMcd = cvMcd(X,p1,res,h);
figure('Position',pos1)
set(gcf,'Name', 'PRESS curve','NumberTitle', 'off');
plot(1:p1,outcvMcd.press,'o-')
title('MCD')
xlabel('number of LV')
ylabel('R-PRESS')
end
if scree == 1
figure('Position',pos2)
screeplot(L,'MCD');
end
if (scree == 1) || (press == 1)
cumperc = cumsum(L)./sum(L);
disp(['The cumulative percentage of variance explained by the first ',num2str(kmax),' components is:']);
disp(num2str(cumperc'));
disp(['How many principal components would you like to retain? Max = ',num2str(kmax),'. ']);
k=input('');
end
% to close the figures.
if scree == 1
close
end
if press == 1
close
end
end
options.k = k;
T=(X-repmat(res.center,size(X,1),1))*U;
out.M=center+res.center*rot';
out.L=L(1:options.k)';
out.P=rot*U(:,1:options.k);
out.T=T(:,1:options.k);
out.h=h;
out.k=options.k;
out.alpha=alfa;
out.Hsubsets.H0 = res.Hsubsets.Hopt;
out.Hsubsets.H1 = [];
out.Hsubsets.Hfreq = res.Hsubsets.Hfreq;
out.skew=skew;
else
% p > n => ROBPCA
niter=100;
seed=0;
if h~=n
if skew==0
B=twopoints(T1,250,seed); %n*ri
for i=1:size(B,1)
Bnorm(i)=norm(B(i,:),2);
end
Bnormr=Bnorm(Bnorm > 1.e-12); %ndirect*1
B=B(Bnorm > 1.e-12,:); %ndirect*n
A=diag(1./Bnormr)*B; %ndirect*n
%projected points in columns
Y=T1*A';%n*ndirect
m=length(Bnormr);
Z=zeros(n,m);
for i=1:m
[tmcdi,smcdi,weights]=unimcd(Y(:,i),h);
if smcdi<1.e-12
r2=rank(data(weights,:));
if r2==1
error(['At least ',num2str(sum(weights)),' obervations are identical.']);
end
else
Z(:,i)=abs(Y(:,i)-tmcdi)/smcdi;
end
end
d=max(Z,[],2);
else %adjusted robpca for skewed data
outAO=adjustedoutlyingness(T1,'ndir',min(250*p,2500));
d=outAO.adjout;
end
[ds,is]=sort(d);
Xh=T1(is(1:h),:); % Xh contains h (good) points out of Xcentr
[P2,T2,L2,r2,Xm,clmX]=classSVD(Xh);
out.Hsubsets.H0 = is(1:h);
Tn=(T1-repmat(clmX,n,1))*P2;
else
P2=eye(r);
Tn=T1;
L2=L1;
r2=r;
out.Hsubsets.H0=1:n;
Xm=T1;
clmX=zeros(1,size(T1,2));
end
%dim(P2) = r x r2
L=L2;
kmax=min(r2,kmax);
% choice of k:
%-------------
bdwidth=5;
topbdwidth=30;
set(0,'Units','pixels');
scnsize=get(0,'ScreenSize');
pos1=[bdwidth, 1/3*scnsize(4)+bdwidth, scnsize(3)/2-2*bdwidth, scnsize(4)/2-(topbdwidth+bdwidth)];
pos2=[pos1(1)+scnsize(3)/2, pos1(2), pos1(3), pos1(4)];
if press == 1
disp('The robust press curve based on cross-validation is now computed.')
outprMCDkmax = projectMCD(Tn,L,kmax,h,niter,rot,P1,P2,center,cutoff);
if size(out.Hsubsets.H0,2)==1
out.Hsubsets.H0=out.Hsubsets.H0';
end
outprMCDkmax.Hsubsets.H0 = out.Hsubsets.H0;
outpress = cvRobpca(data,kmax,outprMCDkmax,0,h);
figure('Position',pos1)
set(gcf,'Name', 'PRESS curve','NumberTitle', 'off');
plot(1:kmax,outpress.press,'o-')
title('ROBPCA')
xlabel('Number of LV')
ylabel('R-PRESS')
else
if size(out.Hsubsets.H0,2)==1
out.Hsubsets.H0=out.Hsubsets.H0';
end
end
if scree == 1
figure('Position',pos2)
screeplot(L(1:kmax),'ROBPCA')
end
if (scree == 1)||(press == 1)
cumperc = (cumsum(L(1:kmax))./sum(L))';
disp(['The cumulative percentage of variance explained by the first ',num2str(kmax),' components is:']);
disp(num2str(cumperc));
disp(['How many principal components would you like to retain? Max = ',num2str(kmax),'.']);
k=input('');
k=max(min(min(r2,k),kmax),1);
% we compute again the robpca results for a specific k value. alpha
% and h can change again, because until now they were based on the kmax
% value.
if dummy == 2
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*alfa);
%if dummy == 1 no changes needed
elseif dummy~=1
options.h=floor(2*floor((n+k+1)/2)-n+2*(n-floor((n+k+1)/2))*options.alpha);
end
h=options.h;
% to close the figures.
if scree == 1
close
end
if press == 1
close
end
else
k=min(min(r2,k),kmax);
end
if k~=r % extra reweighting step
XRc=T1-repmat(clmX,n,1);
Xtilde=XRc*P2(:,1:k)*P2(:,1:k)';
Rdiff = XRc-Xtilde;
for i=1:n
odh(i,1)=norm(Rdiff(i,:));
end
if skew==0
[m,s]=unimcd(odh.^(2/3),h);
cutoffodh = sqrt(norminv(cutoff,m,s).^3);
else %adjusted robpca for skewed data
mcodh=mc(odh);
if mcodh>0
cutoffodh = prctile(odh,75)+1.5*exp(3*mcodh)*iqr(odh);
else
cutoffodh = prctile(odh,75)+1.5*iqr(odh);
end
ttup = sort(-odh(odh<cutoffodh));
cutoffodh = -ttup(1);
end
indexset = find(odh<=cutoffodh)';
[P2,Th,Lh,rh,Xm,clmX]=classSVD(T1(indexset,:));
if k>rh
k = rh;
end
end
center=center+clmX*rot';
rot=rot*P2(:,1:k);
Tn=(T1-repmat(clmX,n,1))*P2;
% if only the subspace is important, not the PC themselves: do not
% perform MCD anymore.
if ~robpcamcd
out.P = rot; %=P1*P2(:,1:k);
out.T = Tn(:,1:k);
out.M = center;
out.L = Lh;
out.k = k;
out.h = h;
out.alpha = alfa;
out.kmax=kmax;
out.skew=skew;
end
% projection, mcd
%-----------------
if skew==0
outpr = projectMCD(Tn,L,k,h,niter,rot,P1,P2,center,cutoff);
else %adjusted robpca for skewed data
outpr = projectAO(Tn,k,h,rot,center);
end
out.T = outpr.T;
out.P = outpr.P;
out.M = outpr.M;
out.L = outpr.L;
out.k = k;
out.kmax=kmax;
out.h = h;
out.alpha = alfa;
out.skew = skew;
if skew==0
out.Hsubsets.H1 = outpr.Hsubsets.H1;
out.Hsubsets.Hfreq = outpr.Hsubsets.Hfreq;
else
out.AO=outpr.AO;
out.cutoff.AO=outpr.cutoff.AO;
end
end
% Classical analysis
if options.classic==1
out.classic.P=P1(:,1:out.k);
out.classic.L=L1(1:out.k)';
out.classic.M=clm;
out.classic.T=T1(:,1:out.k);
out.classic.k=out.k;
out.classic.Xc=Xc;
end
outpr = out;
% Calculation of the distances, flags
%-------------------------------------
if options.classic == 1
outDist = CompDist(data,r,outpr,cutoff,robpcamcd,out.classic);
else
outDist = CompDist(data,r,outpr,cutoff,robpcamcd);
end
out.sd = outDist.sd;
out.cutoff.sd = outDist.cutoff.sd;
out.od = outDist.od;
out.cutoff.od = outDist.cutoff.od;
out.flag = outDist.flag;
out.class = outDist.class;
out.classic = outDist.classic;
if options.classic == 1
out.classic.sd = outDist.classic.sd;
out.classic.od = outDist.classic.od;
out.classic.cutoff.sd = outDist.classic.cutoff.sd;
out.classic.cutoff.od = outDist.classic.cutoff.od;
out.classic.class = outDist.classic.class;
out.classic.flag = outDist.classic.flag;
end
result=struct('P',{out.P},'L',{out.L},'M',{out.M},'T',{out.T},'k',{out.k},'kmax',{kmax},'alpha',{out.alpha},...
'h',{out.h},'Hsubsets',{out.Hsubsets},'sd', {out.sd},'od',{out.od},'cutoff',{out.cutoff},'flag',out.flag',...
'class',{out.class},'classic',{out.classic});
% Plots
try
if plots && options.classic
makeplot(result,'classic',1,'labsd',labsd,'labod',labod)
elseif plots
makeplot(result,'labsd',labsd,'labod',labod)
end
catch %output must be given even if plots are interrupted
%> delete(gcf) to get rid of the menu
end
%--------------------------------------------------------------------------
function outprMCD = projectMCD(Tn,L,k,h,niter,rot,P1,P2,center,cutoff)
% this function performs the last part of ROBPCA when k is determined.
% input :
% Tn : the projected data
% L : the matrix of the eigenvalues
% k : the number of components
% h : lower bound for regular observations
% niter : the number of iterations
% rot : the rotation matrix
% P1, P2: the different eigenvector matrices after each transformation
% center : the classical center of the data
X2=Tn(:,1:k);
n = size(X2,1);
rot=rot(:,1:k);
% first apply c-step with h points from first step,i.e. those that
% determine the covariance matrix after the c-steps have converged.
mah=libra_mahalanobis(X2,zeros(size(X2,2),1),'cov',L(1:k));
oldobj=prod(L(1:k));
P4=eye(k);
korig=k;
for j=1:niter
[mahs,is]=sort(mah);
Xh=X2(is(1:h),:);
[P,T,L,r3,Xm,clmX]=classSVD(Xh);
obj=prod(L);
X2=(X2-repmat(clmX,n,1))*P;
center=center+clmX*rot';
rot=rot*P;
mah=libra_mahalanobis(X2,zeros(size(X2,2),1),'cov',diag(L));
P4=P4*P;
if ((r3==k) && (abs(oldobj-obj) < 1.e-12))
break;
else
oldobj=obj;
j=j+1;
if r3 < k
j=1;
k=r3;
end
end
end
% dim(P4): k x k0 with k0 <= k but denoted as k
% dim X2: n x k0
% perform mcdcov on X2
[zres,zraw]= mcdcov(X2,'plots',0,'ntrial',250,'h',h,'file',0);
out.resMCD = zres;
if zraw.objective < obj
z = zres;
out.Hsubsets.H1 = zres.Hsubsets.Hopt;
else
sortmah = sort(mah);
if h==n
factor=1;
else
factor = sortmah(h)/chi2inv(h/n,k);
end
mah = mah/factor;
weights = mah <= chi2inv(cutoff,k);
[center_noMCD,cov_noMCD] = weightmecov(X2,weights);
mah = libra_mahalanobis(X2,center_noMCD,'cov',cov_noMCD);
z.flag = (mah <= chi2inv(cutoff,k));
z.center = center_noMCD;
z.cov = cov_noMCD;
out.Hsubsets.H1 = is(1:h);
end
covf=z.cov;
centerf=z.center;
[P6,L]=eig(covf);
[L,I]=greatsort(diag(real(L)));
P6=P6(:,I);
out.T=(X2-repmat(centerf,n,1))*P6;
P=P1*P2;
out.P=P(:,1:korig)*P4*P6;
centerfp=center+centerf*rot';
out.M=centerfp;
out.L=L';
out.k=k;
out.h=h;
% creation of Hfreq
out.Hsubsets.Hfreq = zres.Hsubsets.Hfreq(1:h);
outprMCD = out;
%----------------------------------------------------------------------------
function outprAO = projectAO(Tn,k,h,rot,center)
% this function performs the last part of ROBPCA when k is determined.
% input :
% Tn : the projected data
% k : the number of components
% h : lower bound for regular observations
% rot : the rotation matrix
% center : the classical center of the data
seed=0;
X2=Tn(:,1:k);
[n,p]=size(X2);
outAO=adjustedoutlyingness(X2,'ndir',min(250*p,2500));
AO=outAO.adjout;
cutoffAO=outAO.cutoff;
indexset = find(AO<=cutoffAO)';
%SVD
[P6,T6,L6,r6,Xm6,clmX6]=classSVD(X2(indexset,:));
out.T = (X2-repmat(clmX6,n,1))*P6;
out.P = rot*P6;
out.M = center+clmX6*rot';
out.L = L6;
out.k = k;
out.h = h;
out.AO=AO;
out.cutoff.AO=cutoffAO;
outprAO=out;
%--------------------------------------------------------------------------------
function outDist = CompDist(data,r,out,cutoff,robpcamcd,classic)
% Calculates the distances.
% input: data : the original data
% r : the rank of the data
% out is a structure that contains the results of the PCA.
% classic: an optional structure:
% classic.P1
% classic.T1
% classic.L1
% classic.clm
% classic.Xc
if nargin < 6
options.classic = 0;
else
options.classic = 1;
end
n = size(data,1);
p = size(data,2);
k = out.k;
skew=out.skew;
% Computing distances
% Robust score distances in robust PCA subspace
if robpcamcd
if skew==0
out.sd=sqrt(libra_mahalanobis(out.T,zeros(size(out.T,2),1),'cov',out.L))';
out.cutoff.sd=sqrt(chi2inv(cutoff,out.k));
else
out.sd=out.AO;
out.cutoff.sd=out.cutoff.AO;
end
else
out.sd=zeros(n,1);
out.cutoff.sd=0;
end
% Orthogonal distances to robust PCA subspace
XRc=data-repmat(out.M,n,1);
Xtilde=out.T*out.P';
Rdiff=XRc-Xtilde;
for i=1:n
out.od(i,1)=norm(Rdiff(i,:));
end
% Robust cutoff-value for the orthogonal distance
if k~=r
if skew==0
[m,s]=unimcd(out.od.^(2/3),out.h);
out.cutoff.od = sqrt(norminv(cutoff,m,s).^3);
else
mcod=mc(out.od);
if mcod>0
out.cutoff.od = prctile(out.od,75)+1.5*exp(3*mcod)*iqr(out.od);
else
out.cutoff.od = prctile(out.od,75)+1.5*iqr(out.od);
end
ttup = sort(-out.od(out.od<out.cutoff.od));
out.cutoff.od = -ttup(1);
end
else
out.cutoff.od=0;
end
if options.classic==1
% Mahalanobis distance in classical PCA subspace
Tclas=classic.Xc*classic.P(:,1:out.k);
out.classic.sd=sqrt(libra_mahalanobis(Tclas,zeros(size(Tclas,2),1),'invcov',1./classic.L(1:out.k)))';
% Orthogonal distances to classical PCA subspace
Xtilde=Tclas*classic.P(:,1:out.k)';
Cdiff=classic.Xc-Xtilde;
for i=1:n
out.classic.od(i,1)=norm(Cdiff(i,:));
end
out.classic.cutoff.sd=sqrt(chi2inv(cutoff,out.k)); % should be defined after od to have the output in the correct order
% Classical cutoff-values
if k~=r
m=mean(out.classic.od.^(2/3));
s=sqrt(var(out.classic.od.^(2/3)));
out.classic.cutoff.od = sqrt(norminv(cutoff,m,s)^3);
else
out.classic.cutoff.od=0;
end
out.classic.cutoff.sd=sqrt(chi2inv(cutoff,out.k));
out.classic.flag.od=(out.classic.od<=out.classic.cutoff.od);
out.classic.flag.sd=(out.classic.sd<=out.classic.cutoff.sd);
out.classic.flag.all=(out.classic.flag.od)&(out.classic.flag.sd);
out.classic.class='CPCA';
else
out.classic=0;
end
out.class='ROBPCA';
if k~=r
out.flag.od=(out.od<=out.cutoff.od);
out.flag.sd=(out.sd<=out.cutoff.sd);
out.flag.all=(out.flag.od)&(out.flag.sd);
else
out.flag.od=(out.od<=out.cutoff.od);
out.flag.sd=(out.sd<=out.cutoff.sd);
out.flag.all=(out.sd<=out.cutoff.sd);
end
outDist = out;