limma/
poolvar.rs

1//! Satterthwaite pooling of sample variances.
2//!
3//! Pure-Rust port of limma's `poolvar.R` ([`pool_var`]): combine sample
4//! variances with possibly unequal true variances into one pooled estimate with
5//! Satterthwaite-approximated degrees of freedom.
6
7/// Result of [`pool_var`]: the pooled `var`iance, the Satterthwaite effective
8/// `df`, and the summed `multiplier`.
9#[derive(Debug, Clone, Copy, PartialEq)]
10pub struct PoolVar {
11    pub var: f64,
12    pub df: f64,
13    pub multiplier: f64,
14}
15
16/// `poolVar(var, df, multiplier)`. Length-1 inputs are recycled against the
17/// longest. Panics on a negative multiplier; an all-zero multiplier yields
18/// `var = df = multiplier = 0`.
19pub fn pool_var(var: &[f64], df: &[f64], multiplier: &[f64]) -> PoolVar {
20    let len = var.len().max(df.len()).max(multiplier.len());
21    let at = |s: &[f64], i: usize| s[i % s.len()];
22
23    let mmin = (0..len)
24        .map(|i| at(multiplier, i))
25        .fold(f64::INFINITY, f64::min);
26    assert!(mmin >= 0.0, "Multipliers must be non-negative");
27    let mmax = (0..len)
28        .map(|i| at(multiplier, i))
29        .fold(f64::NEG_INFINITY, f64::max);
30    if mmax == 0.0 {
31        return PoolVar {
32            var: 0.0,
33            df: 0.0,
34            multiplier: 0.0,
35        };
36    }
37
38    let sm: f64 = (0..len).map(|i| at(multiplier, i)).sum();
39    let mut var_out = 0.0;
40    let mut denom = 0.0;
41    for i in 0..len {
42        let m = at(multiplier, i) / sm;
43        let s2 = at(var, i);
44        var_out += m * s2;
45        denom += m * m * s2 * s2 / at(df, i);
46    }
47    PoolVar {
48        var: var_out,
49        df: var_out * var_out / denom,
50        multiplier: sm,
51    }
52}
53
54#[cfg(test)]
55mod tests {
56    use super::*;
57
58    fn close(a: f64, b: f64, tol: f64) -> bool {
59        (a - b).abs() <= tol + tol * b.abs()
60    }
61
62    #[test]
63    fn equal_multipliers_match_r() {
64        // Reference: poolVar(c(1,2,3), df=c(9,9,9), multiplier=c(.1,.1,.1)).
65        let p = pool_var(&[1.0, 2.0, 3.0], &[9.0, 9.0, 9.0], &[0.1, 0.1, 0.1]);
66        assert!(close(p.var, 2.0, 1e-12));
67        assert!(close(p.df, 23.1428571428571, 1e-12));
68        assert!(close(p.multiplier, 0.3, 1e-12));
69    }
70
71    #[test]
72    fn unequal_inputs_match_r() {
73        // Reference: poolVar(c(2,4,8), df=c(5,10,20), multiplier=c(.2,.3,.5)).
74        let p = pool_var(&[2.0, 4.0, 8.0], &[5.0, 10.0, 20.0], &[0.2, 0.3, 0.5]);
75        assert!(close(p.var, 5.6, 1e-12));
76        assert!(close(p.df, 32.1311475409836, 1e-12));
77        assert!(close(p.multiplier, 1.0, 1e-12));
78    }
79
80    #[test]
81    fn all_zero_multiplier_is_zero() {
82        let p = pool_var(&[1.0, 2.0, 3.0], &[9.0, 9.0, 9.0], &[0.0, 0.0, 0.0]);
83        assert_eq!(
84            p,
85            PoolVar {
86                var: 0.0,
87                df: 0.0,
88                multiplier: 0.0
89            }
90        );
91    }
92}
limma/poolvar.rs

limma/
poolvar.rs
