Skip to main content

qualia_core_db/solvers/learning/experiment/
power.rs

1//! Power analysis & sample-size (Practical Statistics ch 3) — how many observations
2//! an experiment needs to detect an effect, and the power it achieves at a given
3//! size. Honest experiment planning: state the effect you care about and the
4//! confidence you need, get the sample size you must collect. Reuses the Normal
5//! quantile/CDF from `statistics::distributions`.
6
7use crate::solvers::statistics::distributions::normal;
8
9/// Required sample size **per group** for a two-sample comparison of means, given a
10/// standardized effect size `d` (Cohen's d = mean difference / pooled SD), two-sided
11/// significance `alpha`, and desired `power` (e.g. 0.8). `None` if the inputs are
12/// out of range. Uses the standard normal approximation
13/// `n = 2·(z_{1−α/2} + z_power)² / d²`.
14pub fn required_sample_size_two_sample(d: f64, alpha: f64, power: f64) -> Option<usize> {
15    if d == 0.0
16        || !(0.0..1.0).contains(&alpha)
17        || alpha <= 0.0
18        || !(0.0..1.0).contains(&power)
19        || power <= 0.0
20    {
21        return None;
22    }
23    let za = normal::standard_quantile(1.0 - alpha / 2.0);
24    let zb = normal::standard_quantile(power);
25    let n = 2.0 * (za + zb).powi(2) / (d * d);
26    Some(n.ceil() as usize)
27}
28
29/// Achieved power of a two-sample mean test with `n` per group, effect size `d`,
30/// two-sided `alpha`. `None` on bad inputs.
31pub fn power_two_sample(n: usize, d: f64, alpha: f64) -> Option<f64> {
32    if n == 0 || !(0.0..1.0).contains(&alpha) || alpha <= 0.0 {
33        return None;
34    }
35    let za = normal::standard_quantile(1.0 - alpha / 2.0);
36    // Non-centrality on the standardized scale; one-direction normal approximation.
37    let ncp = d.abs() * (n as f64 / 2.0).sqrt();
38    Some(normal::standard_cdf(ncp - za))
39}
40
41/// Required sample size per group to detect a difference between two proportions
42/// `p1` vs `p2` at two-sided `alpha` and `power`. `None` on bad inputs / equal props.
43pub fn required_sample_size_two_proportion(
44    p1: f64,
45    p2: f64,
46    alpha: f64,
47    power: f64,
48) -> Option<usize> {
49    if !(0.0..=1.0).contains(&p1) || !(0.0..=1.0).contains(&p2) || (p1 - p2).abs() < 1e-12 {
50        return None;
51    }
52    if !(0.0..1.0).contains(&alpha) || alpha <= 0.0 || !(0.0..1.0).contains(&power) || power <= 0.0
53    {
54        return None;
55    }
56    let za = normal::standard_quantile(1.0 - alpha / 2.0);
57    let zb = normal::standard_quantile(power);
58    let pbar = 0.5 * (p1 + p2);
59    let num = (za * (2.0 * pbar * (1.0 - pbar)).sqrt()
60        + zb * (p1 * (1.0 - p1) + p2 * (1.0 - p2)).sqrt())
61    .powi(2);
62    let n = num / (p1 - p2).powi(2);
63    Some(n.ceil() as usize)
64}
65
66#[cfg(test)]
67mod tests {
68    use super::*;
69
70    #[test]
71    fn sample_size_grows_as_effect_shrinks() {
72        let big = required_sample_size_two_sample(0.8, 0.05, 0.8).unwrap();
73        let small = required_sample_size_two_sample(0.2, 0.05, 0.8).unwrap();
74        assert!(
75            small > big,
76            "smaller effect needs more data: {small} !> {big}"
77        );
78        // Known textbook value: d=0.5, α=0.05, power=0.8 → ~63–64 per group.
79        let mid = required_sample_size_two_sample(0.5, 0.05, 0.8).unwrap();
80        assert!((62..=65).contains(&mid), "n per group {mid}");
81    }
82
83    #[test]
84    fn power_and_sample_size_are_consistent() {
85        // Compute n for target power, then the achieved power should be ≥ target.
86        let d = 0.5;
87        let n = required_sample_size_two_sample(d, 0.05, 0.8).unwrap();
88        let achieved = power_two_sample(n, d, 0.05).unwrap();
89        assert!(achieved >= 0.8 - 1e-3, "achieved power {achieved}");
90    }
91
92    #[test]
93    fn power_rises_with_sample_size() {
94        let lo = power_two_sample(10, 0.4, 0.05).unwrap();
95        let hi = power_two_sample(200, 0.4, 0.05).unwrap();
96        assert!(hi > lo && hi <= 1.0);
97    }
98
99    #[test]
100    fn proportion_sample_size() {
101        // Detecting 0.10 vs 0.12 needs many samples; 0.10 vs 0.30 needs few.
102        let subtle = required_sample_size_two_proportion(0.10, 0.12, 0.05, 0.8).unwrap();
103        let obvious = required_sample_size_two_proportion(0.10, 0.30, 0.05, 0.8).unwrap();
104        assert!(subtle > obvious);
105    }
106
107    #[test]
108    fn guards() {
109        assert_eq!(required_sample_size_two_sample(0.0, 0.05, 0.8), None);
110        assert_eq!(power_two_sample(0, 0.5, 0.05), None);
111        assert_eq!(
112            required_sample_size_two_proportion(0.2, 0.2, 0.05, 0.8),
113            None
114        );
115    }
116}