qualia_core_db/solvers/learning/resampling/
permutation.rs

1//! Permutation tests — assumption-free hypothesis testing by resampling.
2//!
3//! Under the null hypothesis that two groups are exchangeable, the labels carry no
4//! information, so the sampling distribution of any test statistic is obtained by
5//! **shuffling the pooled data** and recomputing it. The p-value is the fraction of
6//! shuffles whose statistic is at least as extreme as the observed one. No
7//! distributional assumption — the honest empirical twin of a parametric test.
8
9/// Result of a two-sample permutation test.
10#[derive(Debug, Clone, Copy, PartialEq)]
11pub struct PermutationResult {
12    /// The statistic on the original grouping.
13    pub observed: f64,
14    /// Two-sided p-value `(1 + #{|perm| ≥ |observed|}) / (n_perm + 1)`.
15    pub p_value: f64,
16    pub n_permutations: usize,
17}
18
19struct Lcg(u64);
20impl Lcg {
21    fn below(&mut self, bound: usize) -> usize {
22        self.0 = self
23            .0
24            .wrapping_mul(6364136223846793005)
25            .wrapping_add(1442695040888963407);
26        ((self.0 >> 33) as usize) % bound.max(1)
27    }
28}
29
30/// Two-sample permutation test of `statistic(a) − statistic(b)` (e.g. a difference
31/// of means). Pools the two samples, repeatedly shuffles and re-splits into the
32/// original sizes, and compares. `None` on an empty sample. The `statistic` closure
33/// maps a group slice to a scalar.
34pub fn two_sample_test(
35    a: &[f64],
36    b: &[f64],
37    n_perm: usize,
38    seed: u64,
39    statistic: impl Fn(&[f64]) -> f64,
40) -> Option<PermutationResult> {
41    let (na, nb) = (a.len(), b.len());
42    if na == 0 || nb == 0 || n_perm == 0 {
43        return None;
44    }
45    let observed = statistic(a) - statistic(b);
46    let obs_abs = observed.abs();
47
48    // Pool the samples.
49    let mut pool: Vec<f64> = Vec::with_capacity(na + nb);
50    pool.extend_from_slice(a);
51    pool.extend_from_slice(b);
52    let n = pool.len();
53
54    let mut rng = Lcg(seed ^ 0x9E3779B97F4A7C15);
55    let mut count_extreme = 0usize;
56    let mut ga = vec![0.0; na];
57    let mut gb = vec![0.0; nb];
58    for _ in 0..n_perm {
59        // Fisher–Yates shuffle of the pool, then split.
60        for i in (1..n).rev() {
61            let j = rng.below(i + 1);
62            pool.swap(i, j);
63        }
64        ga.copy_from_slice(&pool[..na]);
65        gb.copy_from_slice(&pool[na..]);
66        let stat = statistic(&ga) - statistic(&gb);
67        if stat.abs() >= obs_abs {
68            count_extreme += 1;
69        }
70    }
71    // Add-one smoothing so the p-value is never exactly 0.
72    let p_value = (1.0 + count_extreme as f64) / (n_perm as f64 + 1.0);
73    Some(PermutationResult {
74        observed,
75        p_value,
76        n_permutations: n_perm,
77    })
78}
79
80#[cfg(test)]
81mod tests {
82    use super::*;
83    use crate::solvers::statistics::descriptive::mean;
84
85    #[test]
86    fn detects_a_real_difference() {
87        // Group b clearly higher than a → small p-value.
88        let a = [1.0, 2.0, 1.5, 2.5, 1.8, 2.2];
89        let b = [8.0, 9.0, 8.5, 9.5, 8.2, 9.1];
90        let r = two_sample_test(&a, &b, 5000, 1, |s| mean(s).unwrap()).unwrap();
91        assert!(r.observed < 0.0); // mean(a) - mean(b) < 0
92        assert!(
93            r.p_value < 0.01,
94            "clear difference should be significant: p={}",
95            r.p_value
96        );
97    }
98
99    #[test]
100    fn no_difference_is_not_significant() {
101        // Two interleaved samples from the same distribution.
102        let a = [5.0, 6.0, 4.0, 5.5, 4.5, 6.5];
103        let b = [5.2, 5.8, 4.2, 5.6, 4.8, 6.2];
104        let r = two_sample_test(&a, &b, 5000, 7, |s| mean(s).unwrap()).unwrap();
105        assert!(
106            r.p_value > 0.2,
107            "similar groups should not be significant: p={}",
108            r.p_value
109        );
110    }
111
112    #[test]
113    fn works_for_a_difference_of_medians() {
114        use crate::solvers::statistics::descriptive::median_in_place;
115        let a = [1.0, 2.0, 3.0, 4.0, 100.0]; // outlier — medians are robust
116        let b = [10.0, 11.0, 12.0, 13.0, 14.0];
117        let r = two_sample_test(&a, &b, 3000, 2, |s| {
118            let mut v = s.to_vec();
119            median_in_place(&mut v).unwrap()
120        })
121        .unwrap();
122        // median(a)=3, median(b)=12 → observed difference negative.
123        assert!(r.observed < 0.0);
124    }
125
126    #[test]
127    fn guards() {
128        assert!(two_sample_test(&[], &[1.0], 100, 0, |_| 0.0).is_none());
129        assert!(two_sample_test(&[1.0], &[1.0], 0, 0, |_| 0.0).is_none());
130    }
131}
qualia_core_db/solvers/learning/resampling/permutation.rs

qualia_core_db/solvers/learning/resampling/
permutation.rs
