robust_rs_core/scale/
sn.rs

1//! Rousseeuw & Croux's `Sn` robust scale.
2
3use crate::error::RobustError;
4use crate::scale::ScaleEstimator;
5use crate::types::Scale;
6
7/// `Sn = c · dₙ · med_i { med_j |rᵢ − rⱼ| }`: a low-median over points of the
8/// per-point high-median of the pairwise absolute differences. The constant
9/// `c = 1.1926` gives consistency for σ at the Gaussian; `dₙ` is a finite-sample
10/// correction (Croux & Rousseeuw 1992).
11///
12/// Like `Qn` it needs no location reference and has a 50% breakdown point, with
13/// ≈ 58% Gaussian efficiency; it is often preferred for its simpler influence
14/// function. Evaluated in `O(n²)` time / `O(n)` space here (the `O(n log n)`
15/// algorithm is a possible follow-up).
16#[derive(Debug, Clone, Copy)]
17pub struct Sn {
18    /// Asymptotic consistency constant (default `1.1926`).
19    pub consistency: f64,
20    /// Apply the finite-sample correction `dₙ` (default `true`).
21    pub finite_sample_correction: bool,
22}
23
24impl Default for Sn {
25    fn default() -> Self {
26        Self {
27            consistency: 1.1926,
28            finite_sample_correction: true,
29        }
30    }
31}
32
33impl ScaleEstimator for Sn {
34    fn scale(&self, residuals: &[f64]) -> Result<Scale, RobustError> {
35        let n = residuals.len();
36        if n < 2 {
37            return Err(RobustError::InsufficientData { needed: 2, got: n });
38        }
39
40        let mut inner = vec![0.0_f64; n];
41        let mut a = vec![0.0_f64; n];
42        for (i, &ri) in residuals.iter().enumerate() {
43            for (j, &rj) in residuals.iter().enumerate() {
44                inner[j] = (ri - rj).abs();
45            }
46            // high median of the n differences (includes the 0 at j = i).
47            let (_, hm, _) = inner.select_nth_unstable_by(n / 2, f64::total_cmp);
48            a[i] = *hm;
49        }
50        // low median of the per-point high-medians.
51        let (_, lm, _) = a.select_nth_unstable_by((n - 1) / 2, f64::total_cmp);
52        let med = *lm;
53
54        let dn = if self.finite_sample_correction {
55            sn_correction(n)
56        } else {
57            1.0
58        };
59        let s = self.consistency * dn * med;
60        if !s.is_finite() || s <= 0.0 {
61            return Err(RobustError::DegenerateScale);
62        }
63        Scale::new(s)
64    }
65}
66
67/// Finite-sample correction `dₙ` for `Sn` (Croux & Rousseeuw 1992): tabulated
68/// for `n ≤ 9`, else `1` (odd) or `n/(n−0.9)` (even); `→ 1` as `n → ∞`.
69fn sn_correction(n: usize) -> f64 {
70    const SMALL: [f64; 8] = [0.743, 1.851, 0.954, 1.351, 0.993, 1.198, 1.005, 1.131];
71    if n <= 9 {
72        SMALL[n - 2]
73    } else if n % 2 == 1 {
74        1.0
75    } else {
76        n as f64 / (n as f64 - 0.9)
77    }
78}
robust_rs_core/scale/sn.rs

robust_rs_core/scale/
sn.rs
