nereids_core/
stats.rs

1//! Small robust-statistics helpers shared across NEREIDS crates.
2//!
3//! These live in `nereids-core` (the dependency-free foundation crate) so that
4//! the same estimator is computed identically everywhere instead of being
5//! re-implemented per crate.  Each helper returns a plain `Option`/value
6//! rather than a formatted error, so the calling crate can map degenerate
7//! inputs onto its own error type / message wording without `nereids-core`
8//! having to know about `IoError`, `FittingError`, etc.
9//!
10//! Precondition: callers pass **finite** inputs (they enforce this via
11//! [`crate::validation`] before calling in).  The helpers nevertheless sort
12//! with [`f64::total_cmp`], which is a total order even over NaN (NaN sorts
13//! last), so a violated precondition degrades to a deterministic — if
14//! meaningless — result instead of a `partial_cmp` panic path.
15
16/// Consistency factor converting a median absolute deviation into a Gaussian
17/// standard-deviation estimate: `sigma ≈ MAD_TO_SIGMA * MAD`.
18///
19/// Derivation: for X ~ N(μ, σ²), the MAD about the median satisfies
20/// P(|X − μ| ≤ MAD) = 1/2, i.e. MAD = σ·Φ⁻¹(3/4) where Φ is the standard
21/// normal CDF.  The consistency factor is therefore
22///
23///   1 / Φ⁻¹(3/4) = 1 / 0.674489750196081… = 1.482602218505601…
24///
25/// so that `MAD_TO_SIGMA * MAD` is an unbiased estimate of σ for Gaussian
26/// data while staying robust (50% breakdown point) against outliers.
27pub const MAD_TO_SIGMA: f64 = 1.482_602_218_505_601_8;
28
29/// Median of a slice of values.
30///
31/// Copies the input and sorts with [`f64::total_cmp`] (total, deterministic —
32/// see the module-level precondition note on NaN).  For even `n` the median
33/// is the midpoint mean of the two central order statistics.
34///
35/// # Arguments
36/// * `values` — Sample values.  Precondition: finite (callers enforce via
37///   [`crate::validation`]).
38///
39/// # Returns
40/// `Some(median)`, or `None` if `values` is empty.
41pub fn median(values: &[f64]) -> Option<f64> {
42    if values.is_empty() {
43        return None;
44    }
45    let mut sorted = values.to_vec();
46    sorted.sort_unstable_by(f64::total_cmp);
47    let n = sorted.len();
48    if n % 2 == 1 {
49        Some(sorted[n / 2])
50    } else {
51        // Even n: midpoint mean of the two central order statistics.
52        Some((sorted[n / 2 - 1] + sorted[n / 2]) / 2.0)
53    }
54}
55
56/// Median absolute deviation of a slice about a given center.
57///
58/// Computes `median(|v − center|)` with the same conventions as [`median`]
59/// (total-order sort, even-`n` midpoint mean).  Note this returns the *raw*
60/// MAD — multiply by [`MAD_TO_SIGMA`] to obtain a Gaussian-consistent σ
61/// estimate.
62///
63/// # Arguments
64/// * `values` — Sample values.  Precondition: finite (callers enforce via
65///   [`crate::validation`]).
66/// * `center` — Center to take deviations about (typically the sample
67///   median).
68///
69/// # Returns
70/// `Some(mad)`, or `None` if `values` is empty.
71pub fn median_abs_deviation(values: &[f64], center: f64) -> Option<f64> {
72    if values.is_empty() {
73        return None;
74    }
75    let deviations: Vec<f64> = values.iter().map(|v| (v - center).abs()).collect();
76    median(&deviations)
77}
78
79#[cfg(test)]
80mod tests {
81    use super::*;
82
83    #[test]
84    fn median_empty_is_none() {
85        assert_eq!(median(&[]), None);
86    }
87
88    #[test]
89    fn median_single_element() {
90        assert_eq!(median(&[3.5]), Some(3.5));
91    }
92
93    #[test]
94    fn median_odd_n() {
95        // Unsorted input; median of {1, 2, 9} is 2.
96        assert_eq!(median(&[9.0, 1.0, 2.0]), Some(2.0));
97    }
98
99    #[test]
100    fn median_even_n_midpoint_mean() {
101        // Sorted {1, 2, 3, 10} → midpoint mean of 2 and 3 = 2.5.
102        assert_eq!(median(&[10.0, 3.0, 1.0, 2.0]), Some(2.5));
103    }
104
105    #[test]
106    fn mad_empty_is_none() {
107        assert_eq!(median_abs_deviation(&[], 0.0), None);
108    }
109
110    #[test]
111    fn mad_basic() {
112        // Values {1, 2, 4, 8}, center 3 → |dev| = {2, 1, 1, 5} → sorted
113        // {1, 1, 2, 5} → midpoint mean of 1 and 2 = 1.5.
114        assert_eq!(median_abs_deviation(&[1.0, 2.0, 4.0, 8.0], 3.0), Some(1.5));
115    }
116
117    #[test]
118    fn mad_about_own_median() {
119        // Values {1, 2, 3, 4, 100}, median 3 → |dev| = {2, 1, 0, 1, 97}
120        // → median 1: the outlier does not blow up the scale estimate.
121        let vals = [1.0, 2.0, 3.0, 4.0, 100.0];
122        let med = median(&vals).unwrap();
123        assert_eq!(med, 3.0);
124        assert_eq!(median_abs_deviation(&vals, med), Some(1.0));
125    }
126
127    #[test]
128    fn mad_of_constant_slice_is_zero() {
129        let vals = [7.0; 5];
130        let med = median(&vals).unwrap();
131        assert_eq!(median_abs_deviation(&vals, med), Some(0.0));
132    }
133
134    #[test]
135    fn mad_to_sigma_matches_inverse_normal_quantile() {
136        // Φ(1/MAD_TO_SIGMA) must equal 3/4: check via the error function
137        // relation Φ(x) = (1 + erf(x/√2))/2 evaluated numerically with a
138        // rational-approximation-free identity — instead verify the inverse
139        // direction: for the standard normal, P(|X| <= q) = 1/2 at
140        // q = 1/MAD_TO_SIGMA ≈ 0.67449.  We hard-check the published value
141        // of Φ⁻¹(3/4) to 15 significant digits.
142        let q = 1.0 / MAD_TO_SIGMA;
143        assert!((q - 0.674_489_750_196_081_7).abs() < 1e-15);
144    }
145}
nereids_core/stats.rs

nereids_core/
stats.rs
