Tutorials › Multivariate Statistics (Graduate) › Tuning Constants For Robust M-Estimators

Home  /  Tutorials  /  Robust Location and Scatter

Robust Location and Scatter · Tutorial 712 of 975

Tuning Constants For Robust M-Estimators

Choose radial M-estimator tuning constants by setting a Gaussian efficiency target, then check the influence and weights that result.

Advanced 10 min read

What You'll Learn

  • Define radial score functions and their distance-dependent weights.
  • Derive the asymptotic Gaussian efficiency of a multivariate radial location M-estimator.
  • Calculate Huber and Tukey bisquare weights at representative Mahalanobis distances.
  • Calibrate a Huber tuning constant to a target efficiency in one dimension.
  • Distinguish bounded influence from redescending influence and identify limits of the tuning choice.

Why Tune a Radial Score?

A radial M-estimator of location uses the distance of an observation from the current location to determine how much it contributes. The tuning constant sets the distance at which that contribution is reduced. A large constant lets the estimator behave more like the sample mean under a Gaussian model; a smaller constant limits the effect of distant observations more aggressively. Choosing the constant is therefore a design decision, not a purely computational detail.

The calculations below use a Gaussian reference distribution to make that tradeoff explicit. They concern the radial score for location. Tuning a location score alone does not specify a joint location-scatter procedure: a scatter equation has its own score and calibration requirements.

Definition: Let \(X\in\mathbb{R}^p\) have location \(\mu\) and positive-definite scatter matrix \(\Sigma\), and write \(r(x;\mu,\Sigma)=\{(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu)\}^{1/2}\). A radial location score \(\psi\) gives the estimating equation $$ \sum_{i=1}^n w(r_i)(X_i-\mu)=0,\qquad w(r)=\frac{\psi(r)}{r}\quad(r>0), $$ where \(r_i=r(X_i;\mu,\Sigma)\). At \(r=0\), define \(w(0)\) by its continuous limit when that limit exists.

The scalar \(\psi(r)\) is the magnitude of the score contribution in the standardized radial direction; \(w(r)\) is the multiplier on the original residual vector. For example, when \(\psi(r)\) is bounded, the contribution \(\psi(r)(X-\mu)/r\) has bounded length. That is the basic mechanism behind bounded influence for location, subject to the usual nonsingularity and regularity conditions for the estimating equation.

Two Common Score Shapes

The Huber radial score is linear near zero and constant beyond its cutoff. Tukey’s bisquare score decreases back to zero at its cutoff. For a positive tuning constant \(c\), their definitions and weights are:

Definition: The Huber score is \(\psi_c(r)=\min(r,c)\), so \(w_c(r)=1\) for \(0\le r\le c\) and \(w_c(r)=c/r\) for \(r>c\). The Tukey bisquare score is \(\psi_c(r)=r\{1-(r/c)^2\}^2\) for \(0\le r<c\), and \(0\) for \(r\ge c\); its weight is \(w_c(r)=\{1-(r/c)^2\}^2\) for \(r<c\), and \(0\) otherwise.

Both scores are bounded, so their location contributions have bounded length. The Huber score does not discard an observation at any finite distance: its contribution keeps a constant magnitude in the outward radial direction. Tukey’s score is redescending, meaning that its contribution returns to zero for distances at or beyond \(c\). Redescending scores can create multiple solutions to the estimating equation, so boundedness alone does not guarantee an easy or unique computation.

Worked Example: Huber Weights at Three Distances

Take \(c=2.5\) and consider standardized Mahalanobis distances \(r=1,\ 2.5,\ 4\). For Huber’s score, \(w_c(r)=1\) up to the cutoff and \(w_c(r)=c/r\) above it.

At \(r=1\), \(r\le c\), so \(w_c(1)=1\). At \(r=2.5\), the observation is exactly at the cutoff, and \(w_c(2.5)=1\). At \(r=4\), \(r>c\), so \(w_c(4)=2.5/4=0.625\).

Thus the residual vectors at distances 1 and 2.5 retain their full multiplier, while the residual at distance 4 is multiplied by 0.625. Its score magnitude is \(r w_c(r)=4(0.625)=2.5=c\), the same magnitude as for any Huber residual beyond the cutoff.

Worked Example: Tukey Bisquare Weights

Set \(c=3\). For a point at distance \(r=1\), its weight is $$ w_3(1)=\left(1-\frac{1^2}{3^2}\right)^2 =\left(\frac{8}{9}\right)^2=\frac{64}{81}\approx0.7901. $$

At distance \(r=2\), $$ w_3(2)=\left(1-\frac{2^2}{3^2}\right)^2 =\left(\frac{5}{9}\right)^2=\frac{25}{81}\approx0.3086. $$ At \(r=3\), the definition gives \(w_3(3)=0\). The corresponding score magnitudes are \(1(64/81)=64/81\), \(2(25/81)=50/81\), and \(0\), respectively. Here the more distant point receives a smaller contribution and a point at the cutoff receives none.

Gaussian Efficiency and Its Derivation

To compare tuning constants, suppose the reference model is \(X\sim N_p(\mu,I_p)\), with known scatter \(I_p\), and let \(Z=X-\mu\) and \(R=\|Z\|\). Assume \(\psi(0)=0\), that \(\psi\) is bounded and piecewise continuously differentiable, and that the expectations below are finite. Also assume the population estimating equation is Fisher-consistent at \(\mu\), and that its derivative matrix is nonsingular. The Gaussian radial density assigns probability zero to the finitely many score kinks, so derivatives at those points need not be specified.

Theorem: For the radial estimating equation with vector score \(g(z)=\psi(\|z\|)z/\|z\|\), define $$ a=\frac{1}{p}E\left[\psi'(R)+(p-1)\frac{\psi(R)}{R}\right]. $$ If \(a\ne0\), the asymptotic covariance per observation of the location estimator is \(\{E[\psi(R)^2]/(p a^2)\}I_p\). Its asymptotic efficiency relative to the Gaussian sample mean is therefore $$ \operatorname{Eff}(\psi)=\frac{p a^2}{E[\psi(R)^2]}. $$

Proof. For \(z\ne0\), put \(r=\|z\|\) and \(u=z/r\). Differentiating \(g(z)=\psi(r)u\) gives a derivative with eigenvalue \(\psi'(r)\) in the radial direction \(u\), and eigenvalue \(\psi(r)/r\) in each of the \(p-1\) orthogonal directions. Thus the trace of its derivative is \(\psi'(r)+(p-1)\psi(r)/r\). Because \(Z\) has a spherical distribution, the expected derivative is a scalar multiple of \(I_p\); its trace shows that this matrix is \(aI_p\).

The derivative of the population estimating equation with respect to the location parameter is \(-aI_p\). Linearizing the sample equation therefore gives influence function \(a^{-1}g(Z)\). Rotational symmetry also implies \(E[g(Z)g(Z)^\mathsf{T}]=\{E[\psi(R)^2]/p\}I_p\): the trace of the left side is \(E[\|g(Z)\|^2]=E[\psi(R)^2]\), and all diagonal directions have the same expectation. The resulting asymptotic covariance is \(a^{-2}E[g(Z)g(Z)^\mathsf{T}]\), which is the stated expression. The sample mean has asymptotic covariance \(I_p\), so the ratio of its variance in any direction to the estimator’s variance in that direction is \(p a^2/E[\psi(R)^2]\). \(\square\)

This theorem gives a calibration rule: choose a reference dimension \(p\), compute the efficiency as a function of \(c\) under the corresponding Gaussian radial distribution, and select \(c\) to meet the desired efficiency. The reference dimension matters because \(R\) has a \(\chi_p\) distribution, not a common distribution for every \(p\). Efficiency is not generally the same across dimensions for a fixed numerical cutoff.

Worked Example: Calibrating a One-Dimensional Huber Score

For \(p=1\), \(R=|Z|\) with \(Z\sim N(0,1)\). Huber’s score has \(\psi_c'(R)=1\) when \(R<c\) and \(0\) when \(R>c\), so \(a=P(|Z|<c)=2\Phi(c)-1\). Also, $$ E[\psi_c(R)^2] =E[Z^2\,1\{|Z|\le c\}]+c^2P(|Z|>c). $$ Integration by parts gives \(E[Z^2\,1\{|Z|\le c\}]=2\Phi(c)-1-2c\phi(c)\), hence $$ \operatorname{Eff}(c)= \frac{\{2\Phi(c)-1\}^2} {2\Phi(c)-1-2c\phi(c)+2c^2\{1-\Phi(c)\}}. $$

At \(c=1.5\), use \(\Phi(1.5)\approx0.93319280\) and \(\phi(1.5)\approx0.12951760\). Then \(a=2(0.93319280)-1=0.86638560\), and $$ E[\psi_c(R)^2]\approx0.86638560-3(0.12951760)+4.5(1-0.93319280) =0.77846522. $$ Therefore \(\operatorname{Eff}(1.5)\approx(0.86638560)^2/0.77846522\approx0.9642\). This is the Gaussian efficiency for the one-dimensional location equation, not a dimension-free calibration.

Checking a Gaussian Tail Calculation

Efficiency calculations often contain Gaussian tail integrals. Keeping the tail probability and the unnormalized tail integral distinct prevents small numerical errors from propagating. For \(c>0\), integration by parts yields $$ \int_c^\infty r^2e^{-r^2/2}\,dr =ce^{-c^2/2}+\int_c^\infty e^{-r^2/2}\,dr. $$

Worked Example: Evaluating the Tail Terms at \(c=1.5\)

Let \(q=e^{-c^2/2}\) and \(T=\int_c^\infty e^{-r^2/2}\,dr=\sqrt{2\pi}\{1-\Phi(c)\}\). At \(c=1.5\), \(q\approx0.32465247\) and \(T\approx0.16746082\). Thus the tail integral is \(cq+T\approx1.5(0.32465247)+0.16746082\approx0.65443953\).

For example, in a decomposition whose interior contribution is \(J\approx0.62022701\), the expression \(J+c^2q+cT\) evaluates to $$ 0.62022701+2.25(0.32465247)+1.5(0.16746082) =0.62022701+0.73046806+0.25119123 \approx1.60188630. $$ The value of \(T\) must be computed as the unnormalized integral above; substituting the normal upper-tail probability alone would omit the factor \(\sqrt{2\pi}\).

Choosing a Constant in Practice

A tuning constant encodes how much Gaussian efficiency to exchange for resistance to large radial residuals. Under a fixed Gaussian reference model, increasing \(c\) makes more observations retain nearly full weight and generally moves the estimator closer to the mean. Decreasing \(c\) imposes stronger downweighting. For a redescending score, the tuning choice also determines where observations cease to contribute altogether.

A defensible choice states the reference dimension and model, the target efficiency or robustness goal, and the score family. It also checks that the resulting location equation is well behaved for the data and algorithm at hand. If the scatter matrix is estimated jointly, the distance itself depends on that estimate; the location efficiency formula above is not a complete joint location-scatter calibration.

  • Bounded is not the same as redescending. Huber contributions stay bounded but do not vanish at large distances; Tukey contributions do vanish beyond the cutoff.
  • Do not use a one-dimensional cutoff automatically in higher dimensions. The Gaussian radial distance changes with \(p\), so its calibration distribution changes too.
  • Do not confuse a weight with a score magnitude. The vector contribution has magnitude \(\psi(r)=r w(r)\), not \(w(r)\) alone.
  • Do not infer joint scatter performance from a location calculation. A scatter score and its calibration conditions are additional parts of a joint procedure.
Key takeaway: Choose a radial tuning constant against an explicit Gaussian reference and dimension, then inspect the resulting weight function and influence behavior. Efficiency quantifies the Gaussian cost; the score’s bounded or redescending shape describes how distant observations are controlled.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

  1. For Huber’s score with \(c=2\), what are the weights at distances \(r=1\) and \(r=5\), and what are the corresponding score magnitudes?
  2. For Tukey’s bisquare score with \(c=4\), calculate the weight at \(r=2\). What is the weight at \(r=4\)?
  3. Why does the Gaussian efficiency formula depend on the dimension \(p\), even if the score function and numerical tuning constant are held fixed?
  4. In one dimension, which two expectations determine the Huber efficiency, and how does the cutoff enter each?
  5. Explain why a bounded location score does not, by itself, establish bounded influence or calibrated efficiency for a jointly estimated scatter matrix.

Previous

Next

Back to Tutorials
Clinical Biostats