limma-rust 0.1.0

//! Satterthwaite pooling of sample variances.
//!
//! Pure-Rust port of limma's `poolvar.R` ([`pool_var`]): combine sample
//! variances with possibly unequal true variances into one pooled estimate with
//! Satterthwaite-approximated degrees of freedom.

/// Result of [`pool_var`]: the pooled `var`iance, the Satterthwaite effective
/// `df`, and the summed `multiplier`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PoolVar {
    pub var: f64,
    pub df: f64,
    pub multiplier: f64,
}

/// `poolVar(var, df, multiplier)`. Length-1 inputs are recycled against the
/// longest. Panics on a negative multiplier; an all-zero multiplier yields
/// `var = df = multiplier = 0`.
pub fn pool_var(var: &[f64], df: &[f64], multiplier: &[f64]) -> PoolVar {
    let len = var.len().max(df.len()).max(multiplier.len());
    let at = |s: &[f64], i: usize| s[i % s.len()];

    let mmin = (0..len)
        .map(|i| at(multiplier, i))
        .fold(f64::INFINITY, f64::min);
    assert!(mmin >= 0.0, "Multipliers must be non-negative");
    let mmax = (0..len)
        .map(|i| at(multiplier, i))
        .fold(f64::NEG_INFINITY, f64::max);
    if mmax == 0.0 {
        return PoolVar {
            var: 0.0,
            df: 0.0,
            multiplier: 0.0,
        };
    }

    let sm: f64 = (0..len).map(|i| at(multiplier, i)).sum();
    let mut var_out = 0.0;
    let mut denom = 0.0;
    for i in 0..len {
        let m = at(multiplier, i) / sm;
        let s2 = at(var, i);
        var_out += m * s2;
        denom += m * m * s2 * s2 / at(df, i);
    }
    PoolVar {
        var: var_out,
        df: var_out * var_out / denom,
        multiplier: sm,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn close(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol + tol * b.abs()
    }

    #[test]
    fn equal_multipliers_match_r() {
        // Reference: poolVar(c(1,2,3), df=c(9,9,9), multiplier=c(.1,.1,.1)).
        let p = pool_var(&[1.0, 2.0, 3.0], &[9.0, 9.0, 9.0], &[0.1, 0.1, 0.1]);
        assert!(close(p.var, 2.0, 1e-12));
        assert!(close(p.df, 23.1428571428571, 1e-12));
        assert!(close(p.multiplier, 0.3, 1e-12));
    }

    #[test]
    fn unequal_inputs_match_r() {
        // Reference: poolVar(c(2,4,8), df=c(5,10,20), multiplier=c(.2,.3,.5)).
        let p = pool_var(&[2.0, 4.0, 8.0], &[5.0, 10.0, 20.0], &[0.2, 0.3, 0.5]);
        assert!(close(p.var, 5.6, 1e-12));
        assert!(close(p.df, 32.1311475409836, 1e-12));
        assert!(close(p.multiplier, 1.0, 1e-12));
    }

    #[test]
    fn all_zero_multiplier_is_zero() {
        let p = pool_var(&[1.0, 2.0, 3.0], &[9.0, 9.0, 9.0], &[0.0, 0.0, 0.0]);
        assert_eq!(
            p,
            PoolVar {
                var: 0.0,
                df: 0.0,
                multiplier: 0.0
            }
        );
    }
}