Summary: the only 'correct' definition is the original one [0][1], the other ones were designed to solve a problem with it that happens in a specific context and would probably have been best called 'pseudo-Stahel-Donoho' distances because they don't, in general, have the same interpretation.
Now, I will show you how to obtain the standard SD distances in R, explain why they are computed the way they are, then explain what these pseudo-SD 'distances' are and why the need for them arose.
Given an $n$ by $p$ data matrix $X=\{x_1,\ldots,x_n\}$ whose entries lie in general (linear) position in $\mathbb{R}^p$, the original definition is to obtain a $p$-vector $d$ as the vector of coefficients of the hyperplane through $p$ data points chosen randomly out of your sample of $n$ observations. In R this is done so:
moar_directions<-function(x){
n<-nrow(x)
p<-ncol(x)
P<-x[sample(n,p),,drop=FALSE]
G<-rep(NA,p)
E<-matrix(1,p,1)
if((qrP<-qr(P))$rank==p) G<-solve(qrP,E)
return(G)
}
n<-100
p<-5
x<-matrix(rnorm(n*p),nc=p)
moar_directions(x)
(see for example line 13 of robustbase::adjOutlyingness here, which you can check by tipping:
library(robustbase)
body(robustbase::adjOutlyingness)
and [2]).
Then, projecting the data unto the direction normal to $d$ is done so:
x%*%moar_directions(x)
The SD outlyingness of $x_i$ with respect to $X=\{x_1,\ldots,x_n\}$ on a single projection unto $d$ is ([0][1]):
$$O(x_i,d,X)=\frac{|x_i'd-\mbox{med}_i(x_i'd)|}{\mbox{mad}_i(x_i'd)}\;\;\;(0)$$
which in R is obtained as:
w<-x%*%moar_directions(x)
abs(w-median(w))/mad(w)
And computing the SD outlyingness as the maximum outlyingness over all the members of $B_p^n$, where
$B_p^n$ is the set of all ${n\choose p}$ such directions $d$ in a $n$ by $p$ data matrix whose entries lie in general (linear) position in $\mathbb{R}^p$. Often, $|B_p^n|={n\choose p}$ will be too large, and it will not be possible to consider all its members in which case one can sample randomly $K$ directions from it to compute an approximation to the SD outlyigness of $x_i$ w.r.t. $X$:
$$O(x_i,X)=\underset{\{d_k\}_{k=1}^K\in B_p^n}{\max}\frac{|x_i'd_k-\mbox{med}_i(x_i'd_k)|}{\mbox{mad}_i(x_i'd_k)}$$
which is R is obtained as:
library(matrixStats)
K<-100
AllThoseProj<-matrix(NA,n,K)
for(i in 1:K){
w<-moar_directions(x)
AllThoseProj[,i]<-abs(w-median(w))/mad(w)
}
rowMaxs(AllThoseProj)
(I will use for a couple of lines K=100 as my definition of 'many' then latter explain what is meant by many in this context)
To motivate the original definition, notice[3] that:
$$d(x_i,X)=\underset{d\in B_p}{\max}\frac{|x_i'd_k-\mbox{mean}_i(x_i'd_k)|}{\mbox{sd}_i(x_i'd_k)}\;\;(1)$$
where
$$d^2(x_i,X)=(x_i-\mbox{mean}_i(x_i))'\mbox{Cov}_i(x_i)^{-1}(x_i-\mbox{mean}_i(x_i))$$
is the vector of squared Mahalanobis distances. So the original SD outlyingness index was designed as a way to compute
a consistent estimator of $d^2(x_i)$ with 50% breakdown point (see [1]) and
equality (1) only holds when the directions $d$ are defined as above.
Now, I will delve a bit more on what is meant by 'many' in ' Often, $|B_p^n|$ will be too large, and ... one can sample randomly many directions'
The notion of many that should be used here depends on the objective being pursued. If you are using the SD to approximate the vector of $d(x_i,X)$'s, then $K$ in the high hundreds should already give a good approximation. If you are using the SD as a robust alternative to the vector of $d(x_i,X)$'s, then a much higher value of $K$ will be necessary (see page 13 of [0] for a discussion of this).
Now, I will delve a bit more on these pseudo-SD distances.
The problem with the way $d$ is defined above is that it is only uniquely defined in settings where $n>p$. In the last decade a lot of research started to be done on high dimensional robustness. The idea built up incrementally to search for outliers through projection pursuit (as was done with the SD distances) but using a type of projection that would also 'work' when $p>n$. There are many such notions, but, gradually, a consensus emerged around the the idea of using directions through two points (see [4] for a late implementation of this idea). Directions through two points are defined as (see [4], and line 52 of body(rrcov:::extradir) in package rrcov for example):
moar_directions_2points<-function(x){
n<-nrow(x)
p<-ncol(x)
P<-x[sample(n,2),,drop=FALSE]
G<-rep(NA,p)
E<-G
E<-P[1,]-P[2,]
N<-sqrt(crossprod(E))
if(N>1e-8) G<-E/N
return(G)
}
Now, if we substitute in equation (0) the original definition of the directions $d$ by the ones above, the new outlyingness index (let's call it pseudo-SD) is no longer (except when $p=2$) a consistent estimator of the vector of $d(x_i,X)$'s, so it is a bit hard to interpret what it is. On the other hand, it can still be computed when $p>n$ in which case the vector of $d(x_i)$'s is not even defined anyway. The fact that we lost consistency by using directions through 2 data points also means we no longer have a target (a quantity we are trying to estimate, like we did with the the vector of $d(x_i,X)$'s in the case of the SD) so the question of what is meant by 'many' here is not really addressable. rrcov for example uses 'many' (in the case of directions through two points) to mean 250.
References:
- [0] Stahel W. (1981).
Breakdown of Covariance Estimators.
Research Report 31, Fachgrupp fur Statistik, E.T.H. Zurich.
- [1] Donoho. D.L. (1982).
Breakdown properties of multivariate location estimators
Ph.D. Qualifying Paper Harvard University.
- [2] Hubert, M. and Van der Veeken, S. (2007). Outlier detection for skewed data. Journal of chemometrics vol:22 issue:3-4 pages:235-246.
- [3] Rousseeuw, P.J. and Leroy, A.M. (1987).
Robust Regression and Outlier Detection. Wiley, New York.
- [4] Hubert, M., Rousseeuw P. J. and Vanden Branden, K. (2005). ROBPCA: A New Approach to Robust Principal Component Analysis. Technometrics
Volume 47, Issue 1.
