subms_stats/
robust.rs

1//! Robust statistics: low sensitivity to outliers, useful when the
2//! raw distribution has a heavy tail (which latency typically does).
3//! Behind the `robust` Cargo feature (on by default).
4
5use crate::percentiles::{mean, percentile, stddev};
6
7/// Interquartile range (p75 - p25). The "middle 50%" spread. More
8/// robust to outliers than stddev.
9pub fn iqr(samples: &[u64]) -> u64 {
10    if samples.is_empty() {
11        return 0;
12    }
13    let mut sorted = samples.to_vec();
14    sorted.sort_unstable();
15    let p25 = percentile(&sorted, 0.25);
16    let p75 = percentile(&sorted, 0.75);
17    p75.saturating_sub(p25)
18}
19
20/// Median absolute deviation: the median of `|x - median(x)|`. Robust
21/// alternative to stddev. Returns 0 for empty input.
22pub fn median_absolute_deviation(samples: &[u64]) -> u64 {
23    if samples.is_empty() {
24        return 0;
25    }
26    let mut sorted = samples.to_vec();
27    sorted.sort_unstable();
28    let median = percentile(&sorted, 0.50);
29    let mut devs: Vec<u64> = sorted.iter().map(|&v| v.abs_diff(median)).collect();
30    devs.sort_unstable();
31    percentile(&devs, 0.50)
32}
33
34/// Coefficient of variation: stddev / mean. Unit-free measure of
35/// relative variability. `0.0` for empty / zero-mean input.
36pub fn coefficient_of_variation(samples: &[u64]) -> f64 {
37    let m = mean(samples) as f64;
38    if m <= 0.0 {
39        return 0.0;
40    }
41    stddev(samples) as f64 / m
42}
43
44/// Skewness (3rd standardised moment). Positive skew means a right
45/// tail (typical for latency distributions). Returns 0 with fewer
46/// than 3 samples.
47pub fn skewness(samples: &[u64]) -> f64 {
48    let n = samples.len();
49    if n < 3 {
50        return 0.0;
51    }
52    let m = mean(samples) as f64;
53    let mut s2 = 0.0f64;
54    let mut s3 = 0.0f64;
55    for &v in samples {
56        let d = v as f64 - m;
57        s2 += d * d;
58        s3 += d * d * d;
59    }
60    let variance = s2 / n as f64;
61    let std = variance.sqrt();
62    if std <= 0.0 {
63        return 0.0;
64    }
65    (s3 / n as f64) / std.powi(3)
66}
67
68/// Excess kurtosis (4th standardised moment minus 3). Positive
69/// excess kurtosis means heavier tails than a normal distribution
70/// (typical for latency). Returns 0 with fewer than 4 samples.
71pub fn kurtosis(samples: &[u64]) -> f64 {
72    let n = samples.len();
73    if n < 4 {
74        return 0.0;
75    }
76    let m = mean(samples) as f64;
77    let mut s2 = 0.0f64;
78    let mut s4 = 0.0f64;
79    for &v in samples {
80        let d = v as f64 - m;
81        s2 += d * d;
82        s4 += d * d * d * d;
83    }
84    let variance = s2 / n as f64;
85    if variance <= 0.0 {
86        return 0.0;
87    }
88    (s4 / n as f64) / variance.powi(2) - 3.0
89}
90
91#[cfg(test)]
92mod tests {
93    use super::*;
94
95    #[test]
96    fn iqr_known_distribution() {
97        let v: Vec<u64> = (0..100).collect();
98        assert_eq!(iqr(&v), 50);
99    }
100
101    #[test]
102    fn mad_basic() {
103        let v: Vec<u64> = (0..100).collect();
104        let mad = median_absolute_deviation(&v);
105        assert!((20..=30).contains(&mad), "MAD around 25: {}", mad);
106    }
107
108    #[test]
109    fn cov_constant_signal_is_zero() {
110        let v = vec![100u64; 100];
111        assert!(coefficient_of_variation(&v) < 0.001);
112    }
113
114    #[test]
115    fn skewness_right_tail_positive() {
116        let mut v: Vec<u64> = vec![100; 990];
117        v.resize(v.len() + 10, 10_000);
118        assert!(skewness(&v) > 0.0);
119    }
120
121    #[test]
122    fn kurtosis_heavy_tail_positive() {
123        let mut v: Vec<u64> = vec![100; 990];
124        v.resize(v.len() + 10, 10_000);
125        assert!(kurtosis(&v) > 0.0);
126    }
127
128    #[test]
129    fn iqr_empty_zero() {
130        assert_eq!(iqr(&[]), 0);
131    }
132}
subms_stats/robust.rs

subms_stats/
robust.rs
