Robust statistics

Deep dive

How it works

In: a sample (or residuals), tuning constant k · Out: robust location and scale, or per-point weights · State: none (IRLS keeps mu between iterations) · Cost: median and MAD O(n) expected; Huber IRLS O(n) per iteration, 5-10 iterations

1. Median: the 50 % quantile; every point has bounded influence (only its rank matters).
2. MAD: median(|x_i - median|); multiply by 1.4826 (= 1 / Phi^-1(0.75)) to estimate sigma under normality.
3. Trimmed mean: sort, drop the lowest and highest fraction p, average the rest. Winsorizing clamps those points to the cut values instead of dropping them.
4. Huber M-estimator: minimize sum rho((x_i - mu) / s), where rho is quadratic for |r| <= k and linear beyond. Points near the centre count fully, far points get influence capped at k. Solve by iteratively reweighted least squares: w_i = min(1, k*s/|x_i - mu|), mu = sum(w_i x_i) / sum(w_i), repeat until mu moves less than a tolerance. Use a robust s (the MAD): with the plain standard deviation the outliers inflate s and disable the protection.
5. Redescending estimators (Tukey biweight, Hampel) give far points weight 0: more resistant to gross outliers, but the objective is non-convex and needs a good start (the median).

influence of a point at distance r:   mean: grows linearly   Huber: capped at k   biweight: back to 0

Code

/// Median in O(n) expected (introselect); reorders `v`. NaNs must be removed first.
pub fn median(v: &mut [f64]) -> f64 {
    let n = v.len();
    assert!(n > 0);
    let mid = n / 2;
    let (lo, hi, _) = v.select_nth_unstable_by(mid, |a, b| a.total_cmp(b));
    if n % 2 == 1 {
        *hi
    } else {
        // lower middle = max of the left partition
        let l = lo.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        0.5 * (l + *hi)
    }
}

/// Normal-consistent MAD scale: 1.4826 * median(|x - median|).
pub fn mad(xs: &[f64]) -> (f64, f64) {
    let mut tmp = xs.to_vec();
    let m = median(&mut tmp);
    for (t, x) in tmp.iter_mut().zip(xs) {
        *t = (x - m).abs();
    }
    (m, 1.482_602_218_505_602 * median(&mut tmp))
}

/// Mean after dropping the lowest and highest `frac` of the sorted sample.
pub fn trimmed_mean(xs: &[f64], frac: f64) -> f64 {
    let mut v = xs.to_vec();
    v.sort_unstable_by(|a, b| a.total_cmp(b));
    let g = (frac * v.len() as f64).floor() as usize;
    let kept = &v[g..v.len() - g];
    kept.iter().sum::<f64>() / kept.len() as f64
}

/// Winsorize in place: clamp to the sample's [frac, 1-frac] order statistics.
pub fn winsorize(xs: &mut [f64], frac: f64) {
    let mut v = xs.to_vec();
    v.sort_unstable_by(|a, b| a.total_cmp(b));
    let g = (frac * v.len() as f64).floor() as usize;
    let (lo, hi) = (v[g], v[v.len() - 1 - g]);
    for x in xs.iter_mut() {
        *x = x.clamp(lo, hi);
    }
}

/// Huber M-estimator of location with the scale fixed to the MAD;
/// IRLS with weights w = min(1, k*s/|r|).
/// k = 1.345 gives about 95 % efficiency at the normal.
pub fn huber_location(xs: &[f64], k: f64) -> f64 {
    let (mut mu, s) = mad(xs);
    if s == 0.0 {
        return mu; // more than half the points are identical: scale degenerate
    }
    for _ in 0..50 {
        let (mut num, mut den) = (0.0, 0.0);
        for &x in xs {
            let r = (x - mu).abs();
            let w = if r <= k * s { 1.0 } else { k * s / r };
            num += w * x;
            den += w;
        }
        let next = num / den;
        let done = (next - mu).abs() < 1e-10 * s;
        mu = next;
        if done {
            break;
        }
    }
    mu
}

Test on 20 points near 10 plus three outliers (1000, 1200, -900): mean 65.2, median 10.0, MAD scale 0.445, 20 % trimmed mean 10.01, Huber 10.02 (clean mean 9.99). Run and asserted.

Complexity and constants

• select_nth_unstable is introselect: O(n) expected, worst case O(n) in Rust's implementation; MAD needs two selections and one n-element scratch copy.
• Huber IRLS: one pass over n values per iteration (a subtraction, a compare, a division); converges in order of 5-10 iterations from the median start.
• Trimmed mean needs a sort or two selections at ranks g and n-g: O(n log n) or O(n).
• Measured here (1M random f64, Apple M-series, release build, includes the copy): median by selection about 2 ms, by full sort about 14 ms, MAD about 2.6 ms, Huber (includes the MAD) about 2.8 ms. Order of magnitude only.

Efficiency

• Selection beats sorting: about 7x here, and no need for sorted output. Never sort to get a median.
• Memory-bound: each pass streams 8 bytes per value; keep data as a contiguous f64/f32 slice and avoid per-element allocation. One scratch buffer can serve both medians of the MAD.
• Few passes, few iterations: start IRLS at the median with MAD scale and stop at a tolerance relative to s; there is no need to re-estimate s each iteration. The weight loop is branch-light and vectorizes (min form).
• Streaming and big data: don't materialize; per window keep a quantile sketch (quantile sketches) or two heaps for a running median. Compute the medians per shard and combine with a median of medians if exactness isn't required.
• Energy: the cost is dominated by memory traffic; process data where it lives, in one or two passes, rather than repeated sorts. Rolling windows: update incrementally (order-statistic tree or two heaps, O(log w) per step) instead of recomputing O(w).

Implementation notes

• Remove NaNs before selecting; total_cmp orders them but they are not meaningful. Use total_cmp, not partial_cmp().unwrap(), to avoid panics.
• Discrete data (integers, counts, rounded values): the MAD can be 0 when more than half the points tie. Guard s == 0 (code does) and fall back to the mean absolute deviation or an IQR-based scale.
• Even n: average the two middle values; for integer data decide whether a half-integer median is acceptable.
• Tuning constants (assuming normal data): Huber k = 1.345 for about 95 % efficiency; Tukey biweight c = 4.685 for about 95 %.
• For regression, the same weights give robust fits (statsmodels RLM, scikit-learn HuberRegressor); see linear and logistic regression. Huber loss is bounded in influence for y-outliers but not for leverage points.
• Two-sided anomaly rule of thumb: flag |x - median| > 3 * 1.4826 * MAD (a "modified z-score"); the threshold is a convention, tune against false positive cost.

Pitfalls

• Using the standard deviation for the scale in a Huber fit: the outliers inflate it and the "robust" estimator reverts to the mean.
• Applying 1.4826 to skewed data: the constant is only consistent at the normal; for log-normal latencies take logs first.
• MAD = 0 on discrete or heavily tied data, then dividing by it.
• Trimming or winsorizing with a fraction below the real contamination: the estimate still breaks.
• Estimating robustly, then testing with a formula that assumes the classical variance: standard errors need a matching robust or bootstrap form.
• Silently dropping outliers that are real events; log what was down-weighted.

Variants and alternatives

In production

Median and percentile summaries in monitoring systems, MAD- and median-based thresholds in anomaly detectors, trimmed means in judged scoring (e.g. sports judging), Huber loss in deep learning (PyTorch HuberLoss) and robust regression, MAD-based outlier screening in experimental science (Leys et al.), robust covariance (MCD) in outlier screening.

Implementations

Sources

• Peter J. Huber, Robust Estimation of a Location Parameter, Annals of Mathematical Statistics 1964. paper
• Frank R. Hampel, The Influence Curve and Its Role in Robust Estimation, JASA 1974. paper
• Peter J. Rousseeuw, Christophe Croux, Alternatives to the Median Absolute Deviation, JASA 1993. paper
• Christophe Leys, Christophe Ley, Olivier Klein, Philippe Bernard, Laurent Licata, Detecting outliers: do not use standard deviation around the mean, use absolute deviation around the median, J. Experimental Social Psychology 49(4), 2013. paper
• Peter J. Huber, Elvezio M. Ronchetti, Robust Statistics, 2nd ed., Wiley 2009.
• Frank R. Hampel, Elvezio M. Ronchetti, Peter J. Rousseeuw, Werner A. Stahel, Robust Statistics: The Approach Based on Influence Functions, Wiley 1986.