qualia_core_db/solvers/learning/experiment/
power.rs1use crate::solvers::statistics::distributions::normal;
8
9pub 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
29pub 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 let ncp = d.abs() * (n as f64 / 2.0).sqrt();
38 Some(normal::standard_cdf(ncp - za))
39}
40
41pub 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 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 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 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}