function out=pcout(data,comp);

%PCOUT is the method proposed by Filzmoser et al. "Outlier detection in
%high dimensions",Computational Statistics&Data Analysis (2008). out.w  
% Method is claimed to be well for outlier detection in high dimensions.

%Required input:
%       data: The data matrix
%       comp: Determines the # of components:
%             less than 1 (percentage of explained variance)
%             greater than 1 (# of components)
%Output:
%    out.w  : Final weights
%  out.rd1  : distances based on kurtosis weighted components
%  out.rd2  : distances based on original components
%    out.V  : Coeffients for sphered data
%    out.T  : Scores for sphered data
%   out.med : Robust center
%   out.mad : Robust scatter


%Created by: Nedret Billor & Asuman Turkmen
%Updated on 04/01/2009

%MAIN PROGRAM:

[n p]=size(data);

%Step1: Robustly sphering the data and computing semi-robust covariance matrix


a1=medscale(data,1);
X=a1.sc;%Sphered data matrix, X*


%Step2: PCA analysis for X in step 1 and robustly sphered PC's (k
%components are chosen based on 99% default cutoff)

[Z,s,v1,V,pr]=pca(X);
I=(1:size(Z,2))';
e=cumsum(pr)/100;

% [U D V]=svd(X'*X,0);
% Z=X*V;
% I=(1:p)';
% e=cumsum(diag(D)/sum(diag(D)));

if comp<1
k=min(I(e>=comp));
else
k=comp;
end;

Told=Z(:,1:k);

a2=medscale(Told,1);
T=a2.sc;%Sphered principal components,T* 

%Step3: Computing robust kurtosis weights for each component and computing
%Euclidian norms for weighted components

T1=(T.^4);
we1=abs(mean(T1)-3);
We1=diag(we1./sum(we1));
Z=T*We1;

    if k==1;
        RD=abs(Z);
    else
        Z2=Z.^2;
    RD=(sqrt(sum(Z2')))'; %Equivalent to the Robust Mahalanobis distances
    end;%
    
d=(RD*sqrt((chi2inv(0.5,k))))/median(RD);%Location distance


%Step4: Determine the weights w1;

r=medscale(d,1);
c=r.med+2.5*r.mad;
M=quantile(d,.33);
cM=c-M;
w1=repmat(0,n,1);
for i=1:n;
    if (M<d(i) & d(i)<c)
        w1(i)=(1-((d(i)-M)/cM)^2)^2;
    elseif d(i)<=M
        w1(i)=1;
    end;
end;
w1=w1(:);

%Step 5: Compute robust distances based on un-weighted components

 if k==1;
        RD2=abs(T);
 else
        T2=T.^2;
   RD2=(sqrt(sum(T2')))'; %Equivalent to the Robust Mahalanobis distances
 end;
 
 d2=(RD2*sqrt((chi2inv(0.5,k))))/median(RD2);%Scatter weights

 %Step 6;Determine the weights w2;

c2=sqrt(chi2inv(0.99,k));
M2=sqrt(chi2inv(0.25,k));
cM2=c2-M2;
w2=repmat(0,n,1);
for i=1:n;
    if (M2<d2(i) & d2(i)<c2)
        w2(i)=(1-((d2(i)-M2)/cM2)^2)^2;
    elseif d2(i)<=M2
        w2(i)=1;
    end;
end;
w2=w2(:);

%Step 7; Combine w1 and w2

a=(1+0.25)^2;
svec=repmat(0.25,n,1);
w=((w1+svec).*(w2+svec))/a;
out.w=w;
out.rd1=d;
out.rd2=d2;
out.mad=a1.mad;
out.med=a1.med;
out.invmad=a1.invmad;
out.V=V(:,1:k);
out.T=Told;%pca based on semi-robust covariance marix


%other required programs:
%----------------------------------
function out=medscale(S,scale);

% Computes the robust standadized matrix 
% Center: using coordinate wise medians
% Scales: by MAD

%input: 
%S: Data matrix
%scale: 
%1: median
%2: l1-median

%Output:
%out.sc=Centered and scaled matrix (columnwise)
%out.med=median(S) %Robust center
%out.mad=MAD(S)    %Robust scatter
%out.invmad        %inverse of diag(scatter)

[n p]=size(S);

if scale==1;
mx=median(S);
elseif scale==2
mx=l1median(S,1e-4);
end;

Sdiff=S-repmat(mx,n,1);
diff=abs(Sdiff);
mad=1.4826*median(diff);
mad=mad(:);
w=[];
for i=1:p
    if (mad(i)~=0)
        a=1/mad(i);
    else
        a=1;
    end;
w=[w;a];
end;
result=Sdiff*diag(w);%Sphered data matrix 

out.sc=result;
out.med=mx;
out.mad=mad;
out.invmad=diag(w);%gives inverse of diag(mad)

%---------------------
function out=mad1(S)
%Computes the median absolute deviance for data set S

n=size(S,1);
mx=median(S);
Sdiff=S-repmat(mx,n,1);
diff=abs(Sdiff);
%out=1.4826*median(diff);
out=median(diff);

%---------------------------------------

function [pc,s,v,d,pr]=pca(x)

[m,n]=size(x);

if m<n                      	
   [s,v]=eig(x*x');
   v=abs(diag(v));
   [a,b]=sort(v);
   b=flipud(b);
   v=v(b);
   s=s(:,b);
   d=x'*(s*diag(sqrt(1./v)));
   pr=100*(v/sum(v));
   v=sqrt(v);
   pc=s*diag(v);
else                        
   [s,v,d]=svd(x,0);
   pc=s*v;
   v=diag(v);
   e=v.^2;
   pr=(e)*100/sum(e);
end