qualia_core_db/solvers/statistics/
robust.rs

1//! Robust / exploratory estimators — location and spread measures that resist
2//! outliers, the EDA complement to the mean/variance in [`super::descriptive`].
3//! They reuse the descriptive median/quantile kernels (no re-implementation).
4
5use super::descriptive::{mean, median_in_place, quantile_in_place};
6
7/// Trimmed mean: drop a `proportion` (in `[0, 0.5)`) of the data from each end and
8/// average the rest. `proportion = 0` is the ordinary mean. `None` if empty or
9/// `proportion` is out of range.
10pub fn trimmed_mean(values: &[f64], proportion: f64) -> Option<f64> {
11    let n = values.len();
12    if n == 0 || !(0.0..0.5).contains(&proportion) {
13        return None;
14    }
15    let mut v = values.to_vec();
16    v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
17    let cut = (n as f64 * proportion).floor() as usize;
18    if 2 * cut >= n {
19        return median_in_place(&mut v); // everything trimmed → fall back to median
20    }
21    mean(&v[cut..n - cut])
22}
23
24/// Winsorized mean: clamp the lowest/highest `proportion` of the data to the
25/// boundary values (rather than dropping them), then average. `None` if empty or
26/// out of range.
27pub fn winsorized_mean(values: &[f64], proportion: f64) -> Option<f64> {
28    let n = values.len();
29    if n == 0 || !(0.0..0.5).contains(&proportion) {
30        return None;
31    }
32    let mut v = values.to_vec();
33    v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
34    let cut = (n as f64 * proportion).floor() as usize;
35    if 2 * cut >= n {
36        return median_in_place(&mut v);
37    }
38    let lo = v[cut];
39    let hi = v[n - 1 - cut];
40    for x in v.iter_mut() {
41        if *x < lo {
42            *x = lo;
43        } else if *x > hi {
44            *x = hi;
45        }
46    }
47    mean(&v)
48}
49
50/// Median absolute deviation `median(|xᵢ − median(x)|)`. With `scaled = true`,
51/// multiplied by 1.4826 so it is a consistent estimator of the standard deviation
52/// for normal data. `None` if empty.
53pub fn median_abs_deviation(values: &[f64], scaled: bool) -> Option<f64> {
54    if values.is_empty() {
55        return None;
56    }
57    let mut v = values.to_vec();
58    let med = median_in_place(&mut v)?;
59    let mut dev: Vec<f64> = values.iter().map(|&x| (x - med).abs()).collect();
60    let mad = median_in_place(&mut dev)?;
61    Some(if scaled {
62        1.482_602_218_505_602 * mad
63    } else {
64        mad
65    })
66}
67
68/// Interquartile range `Q3 − Q1` (the 0.75 and 0.25 quantiles). `None` if empty.
69pub fn iqr(values: &[f64]) -> Option<f64> {
70    if values.is_empty() {
71        return None;
72    }
73    let mut v = values.to_vec();
74    let q3 = quantile_in_place(&mut v, 0.75)?;
75    let q1 = quantile_in_place(&mut v, 0.25)?;
76    Some(q3 - q1)
77}
78
79#[cfg(test)]
80mod tests {
81    use super::*;
82    const EPS: f64 = 1e-9;
83
84    #[test]
85    fn trimmed_mean_ignores_outliers() {
86        // A wild outlier wrecks the mean but not the 20%-trimmed mean.
87        let data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 1000.0];
88        assert!(mean(&data).unwrap() > 100.0);
89        let tm = trimmed_mean(&data, 0.2).unwrap();
90        assert!((tm - 5.5).abs() < 1.0, "trimmed mean {tm}"); // ~ middle of 3..8
91                                                              // proportion 0 == ordinary mean.
92        assert!((trimmed_mean(&[1.0, 2.0, 3.0], 0.0).unwrap() - 2.0).abs() < EPS);
93    }
94
95    #[test]
96    fn winsorized_mean_pulls_in_tails() {
97        let data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 1000.0];
98        let wm = winsorized_mean(&data, 0.1).unwrap();
99        // The 1000 is clamped to 9 → mean far below the raw mean.
100        assert!(wm < 10.0, "winsorized mean {wm}");
101    }
102
103    #[test]
104    fn mad_and_iqr_measure_robust_spread() {
105        // 0..=8 symmetric: median 4, |x-4| = 4,3,2,1,0,1,2,3,4 → MAD median = 2.
106        let data: Vec<f64> = (0..=8).map(|i| i as f64).collect();
107        assert!((median_abs_deviation(&data, false).unwrap() - 2.0).abs() < EPS);
108        // Scaled MAD ≈ 1.4826·2.
109        assert!((median_abs_deviation(&data, true).unwrap() - 2.965204).abs() < 1e-5);
110        // IQR: Q1=2, Q3=6 → 4.
111        assert!((iqr(&data).unwrap() - 4.0).abs() < EPS);
112    }
113
114    #[test]
115    fn robust_resists_a_single_contaminant() {
116        let clean: Vec<f64> = (1..=20).map(|i| i as f64).collect();
117        let mut dirty = clean.clone();
118        dirty[0] = 1e6;
119        // MAD barely moves; std-dev explodes.
120        let m_clean = median_abs_deviation(&clean, false).unwrap();
121        let m_dirty = median_abs_deviation(&dirty, false).unwrap();
122        assert!((m_clean - m_dirty).abs() < 2.0);
123    }
124
125    #[test]
126    fn guards() {
127        assert_eq!(trimmed_mean(&[], 0.1), None);
128        assert_eq!(trimmed_mean(&[1.0], 0.6), None);
129        assert_eq!(iqr(&[]), None);
130        assert_eq!(median_abs_deviation(&[], false), None);
131    }
132}
qualia_core_db/solvers/statistics/robust.rs

qualia_core_db/solvers/statistics/
robust.rs
