Stahel-Donoho Location Estimation

From Directional Outlyingness to a Location Estimate

A multivariate observation may look ordinary in one coordinate and extreme in another direction. Projection outlyingness addresses this by comparing an observation with the sample along every direction. The Stahel–Donoho location estimator then uses those scores to give less weight to observations that are extreme in at least one direction.

The construction has two distinct parts: first calculate how outlying each observation is, then turn that score into a weight and average the observations. The choice of weight function matters. A commonly used bounded inverse rule reduces an extreme observation’s weight without automatically discarding it.

Projection Outlyingness and the Weight Rule

Definition: Let $x_1,\ldots,x_n\in\mathbb{R}^p$ be a sample. For a unit direction $u\in\mathbb{R}^p$, project the observations to $u^\top x_1,\ldots,u^\top x_n$. Let $\operatorname{med}(u^\top X)$ denote the median of these projected values and let $\operatorname{MAD}(u^\top X)$ be the median of their absolute deviations from that median. When $\operatorname{MAD}(u^\top X)>0$ in every direction, the sample projection outlyingness of $x_i$ is $$ r_i=\sup_{\|u\|=1}\frac{|u^\top x_i-\operatorname{med}(u^\top X)|}{\operatorname{MAD}(u^\top X)}. $$

Here $X$ denotes the sample, not a random variable. The median absolute deviation is used without a normal-consistency multiplier; multiplying every directional scale by the same positive constant would simply rescale all the $r_i$ values and would require corresponding adjustment of the weight tuning constant. The definition requires positive projected MADs. If a sample has a direction with zero MAD, the displayed ratio is not defined in the ordinary way; a practical implementation must specify a convention, such as regularizing the scale, or use a sample for which the directional scales are positive.

The supremum searches all directions, not just the coordinate axes. In a computation, a finite set of directions usually approximates that search. Scores from such a direction grid are therefore approximate projection outlyingness scores; a grid can miss a direction in which an observation is more extreme.

Definition: Given a nonincreasing weight function $w(r)\geq 0$ and positive total weight, the Stahel–Donoho weighted location estimate is $$ \hat{\mu}_{\mathrm{SD}}=\frac{\sum_{i=1}^n w(r_i)x_i}{\sum_{i=1}^n w(r_i)}. $$ A common choice, for a tuning constant $c>0$, is $w(r)=\min(1,c/r)$ for $r>0$, with $w(0)=1$. Thus observations with $r\leq c$ receive weight one, while those with $r>c$ receive weight $c/r$.

The rule is bounded above by one and gives smaller weights to larger scores. The tuning constant $c$ sets where downweighting begins: a larger $c$ leaves more observations at full weight. Other nonincreasing weight functions can be used, so any reported estimate should be accompanied by its scoring and weighting conventions.

Worked Example: Calculating a Directional Score

Consider the seven observations $(0,2),(1,-1),(2,4),(3,0),(4,3),(5,-2),(20,1)$ in $\mathbb{R}^2$. Calculate the score of the last observation in the direction $u=(1,0)^\top$.

Project the observations. The projected values are $0,1,2,3,4,5,20$, so their median is $3$. The absolute deviations from $3$ are $3,2,1,0,1,2,17$. Ordered, these are $0,1,1,2,2,3,17$, with median $2$. Thus the projected MAD is $2$.

Standardize the last observation’s deviation. Its projection is $20$, so its absolute deviation from the projected median is $17$. Its directional score is therefore $17/2=8.5$. The full projection outlyingness is a supremum over directions, so the information from this one direction establishes $r_7\geq 8.5$; it does not by itself establish the exact value of $r_7$.

Forming Weights and Computing the Estimate

Worked Example: Turning Scores Into Weights

Suppose a projection search gives four observations the scores $0.5,1,2,6$, and choose $c=3$. Apply $w(r)=\min(1,3/r)$, with $w(0)=1$.

The first three scores are at most $3$, so their weights are all $1$. The fourth score exceeds $3$, so its weight is $3/6=0.5$. The weight vector is $(1,1,1,0.5)$. In particular, an outlyingness score above the cutoff does not mean that the observation is removed: it retains a positive, reduced weight under this rule.

Worked Example: A Weighted Location From a Direction Grid

Consider the six points $(0,0),(1,0),(2,0),(0,1),(0,2),(10,10)$. To make every calculation explicit, use only the two coordinate directions $e_1=(1,0)^\top$ and $e_2=(0,1)^\top$. This is a finite-grid approximation, not the full all-directions estimator. Use the larger of the two coordinate-direction scores for each point, and take $c=3$.

In each coordinate, the projected values are $0,0,0,1,2,10$, with median $(0+1)/2=0.5$. Their absolute deviations are $0.5,0.5,0.5,0.5,1.5,9.5$, so the MAD is $(0.5+0.5)/2=0.5$. Consequently, in either coordinate direction a coordinate value $a$ has score $|a-0.5|/0.5$.

The maximum of the two coordinate scores gives $r=(1,1,3,1,3,19)$, in the order listed. For example, $(2,0)$ has coordinate scores $3$ and $1$, hence maximum $3$; $(10,10)$ has score $19$ in each coordinate. The weights are therefore $(1,1,1,1,1,3/19)$.

The total weight is $5+3/19=98/19$. The weighted coordinate sums are equal, and each is $0+1+2+0+0+10(3/19)=87/19$. Thus the grid-based weighted location is $$ \left(\frac{87/19}{98/19},\frac{87/19}{98/19}\right) =\left(\frac{87}{98},\frac{87}{98}\right) \approx(0.888,0.888). $$ The outlying point is downweighted, but still contributes to the estimate. Adding more directions can change the scores, weights, and resulting location.

What the Construction Guarantees

Theorem: If the weights $w(r_i)$ are nonnegative and their sum is positive, the Stahel–Donoho weighted location estimate is a convex combination of the observations. In particular, it lies in their convex hull.

Proof. Write $W=\sum_{i=1}^n w(r_i)>0$ and $\alpha_i=w(r_i)/W$. Then $\alpha_i\geq0$ for every $i$, and $\sum_i\alpha_i=1$. The estimator is $\sum_i\alpha_i x_i$, which by the definition of a convex combination lies in the convex hull of $x_1,\ldots,x_n$. $\square$

This result is a geometric consequence of the final weighted average; it does not claim that the estimate is resistant to every possible sample configuration. Robustness depends on the outlyingness calculation and weight rule, as well as on the sample and the directions searched.

Theorem: Suppose $A$ is a nonsingular $p\times p$ matrix, and the projected medians and MADs are computed using conventions equivariant under translation and positive rescaling. If every directional MAD is positive, the all-directions Stahel–Donoho estimator with weights determined only by $r_i$ is affine equivariant: transforming every observation to $Ax_i+b$ transforms the estimate to $A\hat{\mu}_{\mathrm{SD}}+b$.

Proof. Let $y_i=Ax_i+b$. For a unit direction $v$ in the transformed space, set $a=A^\top v$. Since $A$ is nonsingular, $a\ne0$; write $u=a/\|a\|$, which is a unit direction in the original space. Then $v^\top y_i=v^\top b+\|a\|u^\top x_i$. Translation and positive-rescaling equivariance of the median give $\operatorname{med}(v^\top Y)=v^\top b+\|a\|\operatorname{med}(u^\top X)$. The corresponding MAD is $\operatorname{MAD}(v^\top Y)=\|a\|\operatorname{MAD}(u^\top X)$. Therefore the transformed standardized deviation for observation $i$ in direction $v$ equals its original standardized deviation in direction $u$. As $v$ ranges over the unit sphere, the normalized directions $u=A^\top v/\|A^\top v\|$ cover the unit sphere: for any unit $u$, one can choose $v$ proportional to $A^{-\top}u$. Taking suprema shows that each transformed score equals $r_i$. The weights are unchanged, so $\sum_i w(r_i)y_i/\sum_iw(r_i)=A\hat{\mu}_{\mathrm{SD}}+b$.

A finite direction grid need not preserve this affine equivariance: after a transformation, the grid may fail to correspond to the same set of original directions. This is one reason to distinguish the exact all-directions definition from a computational approximation.

Interpretation and Common Pitfalls

The estimator’s design is straightforward: measure extremeness relative to robust one-dimensional summaries in every direction, then limit the contribution of observations with large scores. The method matters when a small number of unusual observations could pull a conventional sample mean strongly away from the bulk of the data. The estimator is still a location summary; the weighting formula alone does not define a scatter estimate.

Using only coordinate-wise checks. An observation can be extreme in a direction that is not one of the coordinate axes. Coordinate directions alone define a grid approximation, not the supremum in projection outlyingness.
Treating a score as a weight. A score $r_i$ measures standardized extremeness; it is not itself the multiplier in the location average. Apply the stated function $w(r_i)$.
Assuming downweighting means deletion. With $w(r)=\min(1,c/r)$, every finite positive score receives a positive weight. Rejection requires a different weight rule.
Ignoring zero projected scales. A zero MAD makes the displayed directional ratio undefined or infinite unless an explicit convention is supplied. Check the scale condition or state how the implementation handles it.

Key takeaway: The Stahel–Donoho location estimator is constructed by maximizing robust standardized deviations over directions, converting the resulting observation-level scores to weights, and taking the weighted mean. State the direction search, projected center and scale, and weight rule to make the estimator unambiguous.

Check Your Understanding

Try these questions before moving on.

For one direction, the projected median is $4$, the projected MAD is $2$, and an observation projects to $10$. What is its directional score, and what does this tell you about its full projection outlyingness?
With $c=4$ and the bounded inverse weight rule, calculate the weights for scores $0,2,4,8$.
Why can a search over the coordinate axes underestimate full projection outlyingness?
Explain why nonnegative weights summing to one put the weighted estimate inside the convex hull of the observations.
What condition on the projected MADs is needed for the displayed all-directions score to be defined in the ordinary way?