open_eeg_codec_standard/
stats.rs

1//! Statistical rigor for codec comparison (host-side; not on the grade path).
2//!
3//! Two tools the `bench` command uses to make a comparison *citable* rather
4//! than a single point estimate:
5//!
6//! - [`bootstrap_ci`] — a percentile bootstrap confidence interval for a
7//!   per-file metric (e.g. Pearson R across a corpus), so a codec's headline
8//!   number carries an uncertainty band. Seeded ([`rand::rngs::StdRng`]) so the
9//!   interval is reproducible run to run.
10//! - [`sign_test`] — a paired two-sided sign test (exact Binomial tail via
11//!   [`statrs`]) for "is codec A actually better than codec B on this corpus,
12//!   or is the gap noise?".
13//!
14//! These wrap mature crates (`rand`, `statrs`) rather than re-deriving
15//! distributions by hand.
16
17use rand::rngs::StdRng;
18use rand::{Rng, SeedableRng};
19use statrs::distribution::{Binomial, DiscreteCDF};
20
21/// A confidence interval: the point estimate plus its `[lo, hi]` band.
22#[derive(Clone, Copy, Debug, PartialEq)]
23pub struct Ci {
24    /// Point estimate (the sample mean).
25    pub mean: f64,
26    /// Lower bound of the interval.
27    pub lo: f64,
28    /// Upper bound of the interval.
29    pub hi: f64,
30}
31
32/// Percentile bootstrap confidence interval of the mean of `samples`.
33///
34/// Resamples `samples` with replacement `n_boot` times, takes the mean of each
35/// resample, and reports the `confidence` central interval of those means
36/// (e.g. `confidence = 0.95` → the 2.5th/97.5th percentiles). The RNG is
37/// seeded by `seed`, so the interval is deterministic. Degenerate inputs
38/// (0 or 1 sample, or `n_boot == 0`) collapse to a zero-width interval at the
39/// mean.
40pub fn bootstrap_ci(samples: &[f64], confidence: f64, n_boot: usize, seed: u64) -> Ci {
41    let n = samples.len();
42    if n == 0 {
43        return Ci { mean: 0.0, lo: 0.0, hi: 0.0 };
44    }
45    let mean = samples.iter().sum::<f64>() / n as f64;
46    if n == 1 || n_boot == 0 {
47        return Ci { mean, lo: mean, hi: mean };
48    }
49
50    let mut rng = StdRng::seed_from_u64(seed);
51    let mut means = Vec::with_capacity(n_boot);
52    for _ in 0..n_boot {
53        let mut acc = 0.0;
54        for _ in 0..n {
55            acc += samples[rng.gen_range(0..n)];
56        }
57        means.push(acc / n as f64);
58    }
59    means.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
60
61    let alpha = (1.0 - confidence) / 2.0;
62    Ci {
63        mean,
64        lo: percentile(&means, alpha),
65        hi: percentile(&means, 1.0 - alpha),
66    }
67}
68
69/// Linear-interpolated quantile of a sorted slice (`q` in `[0, 1]`).
70fn percentile(sorted: &[f64], q: f64) -> f64 {
71    if sorted.is_empty() {
72        return 0.0;
73    }
74    let q = q.clamp(0.0, 1.0);
75    let idx = q * (sorted.len() - 1) as f64;
76    let lo = idx.floor() as usize;
77    let hi = idx.ceil() as usize;
78    if lo == hi {
79        sorted[lo]
80    } else {
81        sorted[lo] + (idx - lo as f64) * (sorted[hi] - sorted[lo])
82    }
83}
84
85/// Paired two-sided sign-test p-value for "do `a` and `b` differ?".
86///
87/// For each paired index, the sign of `a[i] - b[i]` is a +/− vote (ties
88/// dropped). Under H0 (no systematic difference) the count of one sign is
89/// `Binomial(m, 0.5)` over the `m` non-tie pairs; the two-sided p-value is the
90/// exact Binomial tail (via [`statrs`]), clamped to `1.0`. Returns `1.0` when
91/// there are no non-tie pairs (no evidence either way). The slices are paired
92/// up to the shorter length.
93pub fn sign_test(a: &[f64], b: &[f64]) -> f64 {
94    let n = a.len().min(b.len());
95    let mut plus = 0u64;
96    let mut minus = 0u64;
97    for i in 0..n {
98        let d = a[i] - b[i];
99        if d > 0.0 {
100            plus += 1;
101        } else if d < 0.0 {
102            minus += 1;
103        }
104    }
105    let m = plus + minus;
106    if m == 0 {
107        return 1.0;
108    }
109    let k = plus.min(minus);
110    // P(X <= k) under Binomial(m, 0.5); doubled for the two-sided test.
111    let binom = Binomial::new(0.5, m).expect("0.5 is a valid Binomial p");
112    (2.0 * binom.cdf(k)).min(1.0)
113}
114
115#[cfg(test)]
116mod tests {
117    use super::*;
118
119    #[test]
120    fn ci_of_constant_is_zero_width() {
121        let ci = bootstrap_ci(&[5.0; 16], 0.95, 1000, 42);
122        assert_eq!(ci.mean, 5.0);
123        assert_eq!(ci.lo, 5.0);
124        assert_eq!(ci.hi, 5.0);
125    }
126
127    #[test]
128    fn ci_brackets_the_mean_and_is_reproducible() {
129        let xs: Vec<f64> = (0..50).map(|i| (i as f64) * 0.1).collect();
130        let a = bootstrap_ci(&xs, 0.95, 2000, 7);
131        let b = bootstrap_ci(&xs, 0.95, 2000, 7);
132        assert_eq!(a, b, "seeded bootstrap is reproducible");
133        assert!(a.lo <= a.mean && a.mean <= a.hi, "mean inside the interval");
134        assert!(a.lo < a.hi, "non-degenerate interval has positive width");
135    }
136
137    #[test]
138    fn edge_cases() {
139        assert_eq!(bootstrap_ci(&[], 0.95, 100, 1).mean, 0.0);
140        let one = bootstrap_ci(&[3.0], 0.95, 100, 1);
141        assert_eq!((one.lo, one.mean, one.hi), (3.0, 3.0, 3.0));
142    }
143
144    #[test]
145    fn sign_test_extremes() {
146        // 10 paired points, a always greater -> small p (~2 * 0.5^10).
147        let a = vec![1.0; 10];
148        let b = vec![0.0; 10];
149        let p = sign_test(&a, &b);
150        assert!(p < 0.01, "all-positive sign test should be significant, got {p}");
151        assert!((p - 2.0 * 0.5_f64.powi(10)).abs() < 1e-9);
152
153        // Identical -> all ties -> no evidence -> p = 1.0.
154        assert_eq!(sign_test(&a, &a), 1.0);
155
156        // Balanced -> p close to 1.0.
157        let x = vec![1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
158        let y = vec![0.0; 6];
159        assert!(sign_test(&x, &y) > 0.9);
160    }
161}
open_eeg_codec_standard/stats.rs

open_eeg_codec_standard/
stats.rs
