robust_rs_core/scale/
mad.rs

1//! Median absolute deviation scale.
2
3use crate::error::RobustError;
4use crate::scale::ScaleEstimator;
5use crate::types::Scale;
6
7/// MAD: `s = 1.4826 · median(|rᵢ − median(r)|)`. The constant `1.4826 = 1/Φ⁻¹(¾)`
8/// makes it consistent for σ at the Gaussian model.
9#[derive(Debug, Clone, Copy)]
10pub struct Mad {
11    /// Consistency constant (default `1.4826`).
12    pub consistency: f64,
13}
14
15impl Default for Mad {
16    fn default() -> Self {
17        Self {
18            consistency: 1.482_602_218_505_602,
19        }
20    }
21}
22
23impl ScaleEstimator for Mad {
24    fn scale(&self, residuals: &[f64]) -> Result<Scale, RobustError> {
25        if residuals.is_empty() {
26            return Err(RobustError::InsufficientData { needed: 1, got: 0 });
27        }
28
29        // Copy so the caller's slice is untouched; we sort in place twice.
30        let mut buf = residuals.to_vec();
31
32        let center = median(&mut buf);
33        for x in buf.iter_mut() {
34            *x = (*x - center).abs(); // reuse buf; order no longer matters
35        }
36        let mad = median(&mut buf);
37
38        let s = mad * self.consistency;
39        // Broadened from the literal "if 0": also reject NaN/±inf. A majority of
40        // non-finite residuals reaches the median and yields a non-finite scale
41        // that `== 0.0` would miss, leaking out as Ok(Scale(NaN)).
42        if !s.is_finite() || s <= 0.0 {
43            return Err(RobustError::DegenerateScale);
44        }
45        Scale::new(s)
46    }
47}
48
49/// Median via in-place total-order sort. `total_cmp` never panics on NaN and
50/// for even `n` we average the two central order statistics, which is the
51/// definition the `1.4826` constant is derived for.
52fn median(v: &mut [f64]) -> f64 {
53    v.sort_unstable_by(f64::total_cmp);
54    let n = v.len();
55    let mid = n / 2;
56    if n % 2 == 1 {
57        v[mid]
58    } else {
59        0.5 * (v[mid - 1] + v[mid])
60    }
61}
robust_rs_core/scale/mad.rs

robust_rs_core/scale/
mad.rs
