Skip to main content

robust_rs_core/scale/
qn.rs

1//! Rousseeuw & Croux's `Qn` robust scale.
2
3use crate::error::RobustError;
4use crate::scale::ScaleEstimator;
5use crate::types::Scale;
6
7/// `Qn = c · dₙ · {|rᵢ − rⱼ| : i < j}₍ₖ₎`, the `k`-th order statistic of the
8/// pairwise absolute differences with `k = C(⌊n/2⌋+1, 2)` (≈ the ¼-quantile).
9/// The constant `c = 2.2219 = 1/(√2·Φ⁻¹(⅝))` gives consistency for σ at the
10/// Gaussian; `dₙ` is a finite-sample correction (Croux & Rousseeuw 1992).
11///
12/// Uses no location step and stays smooth in the data, with ≈ 82% Gaussian
13/// efficiency and a 50% breakdown point (versus the MAD's ≈ 37%).
14///
15/// Evaluated directly over all pairs in `O(n²)` time / `O(n²)` space; Rousseeuw
16/// & Croux's `O(n log n)` algorithm is a possible follow-up (the datasets here
17/// are small).
18#[derive(Debug, Clone, Copy)]
19pub struct Qn {
20    /// Asymptotic consistency constant (default `2.2219`).
21    pub consistency: f64,
22    /// Apply the finite-sample correction `dₙ` (default `true`).
23    pub finite_sample_correction: bool,
24}
25
26impl Default for Qn {
27    fn default() -> Self {
28        Self {
29            consistency: 2.2219,
30            finite_sample_correction: true,
31        }
32    }
33}
34
35impl ScaleEstimator for Qn {
36    fn scale(&self, residuals: &[f64]) -> Result<Scale, RobustError> {
37        let n = residuals.len();
38        if n < 2 {
39            return Err(RobustError::InsufficientData { needed: 2, got: n });
40        }
41
42        let mut diffs = Vec::with_capacity(n * (n - 1) / 2);
43        for i in 0..n {
44            for j in (i + 1)..n {
45                diffs.push((residuals[i] - residuals[j]).abs());
46            }
47        }
48
49        // k-th smallest (1-indexed), k = C(h, 2) with h = ⌊n/2⌋ + 1.
50        let h = n / 2 + 1;
51        let k = h * (h - 1) / 2;
52        let (_, kth, _) = diffs.select_nth_unstable_by(k - 1, f64::total_cmp);
53        let kth = *kth;
54
55        let dn = if self.finite_sample_correction {
56            qn_correction(n)
57        } else {
58            1.0
59        };
60        let s = self.consistency * dn * kth;
61        if !s.is_finite() || s <= 0.0 {
62            return Err(RobustError::DegenerateScale);
63        }
64        Scale::new(s)
65    }
66}
67
68/// Finite-sample correction `dₙ` for `Qn` (Croux & Rousseeuw 1992): tabulated
69/// for `n ≤ 9`, else `n/(n+3.8)` (even) or `n/(n+1.4)` (odd); `→ 1` as `n → ∞`.
70fn qn_correction(n: usize) -> f64 {
71    const SMALL: [f64; 8] = [0.399, 0.994, 0.512, 0.844, 0.611, 0.857, 0.669, 0.872];
72    if n <= 9 {
73        SMALL[n - 2]
74    } else if n % 2 == 0 {
75        n as f64 / (n as f64 + 3.8)
76    } else {
77        n as f64 / (n as f64 + 1.4)
78    }
79}