qualia_core_db/solvers/statistics/hypothesis/
chi_square.rs1use super::super::distributions::chi_squared;
5
6#[derive(Debug, Clone, Copy, PartialEq)]
8pub struct ChiSquareResult {
9 pub statistic: f64,
10 pub p_value: f64,
11 pub dof: f64,
12}
13
14pub fn chi_square_gof(observed: &[f64], expected: &[f64]) -> Option<ChiSquareResult> {
17 if observed.len() != expected.len() || observed.len() < 2 {
18 return None;
19 }
20 if expected.iter().any(|&e| e <= 0.0) {
21 return None;
22 }
23 let stat: f64 = observed
24 .iter()
25 .zip(expected.iter())
26 .map(|(&o, &e)| (o - e).powi(2) / e)
27 .sum();
28 let dof = (observed.len() - 1) as f64;
29 Some(ChiSquareResult {
30 statistic: stat,
31 p_value: chi_squared::upper_p(stat, dof),
32 dof,
33 })
34}
35
36pub fn chi_square_independence(table: &[&[f64]]) -> Option<ChiSquareResult> {
40 let rows = table.len();
41 if rows < 2 {
42 return None;
43 }
44 let cols = table[0].len();
45 if cols < 2 || table.iter().any(|r| r.len() != cols) {
46 return None;
47 }
48 let row_sums: Vec<f64> = table.iter().map(|r| r.iter().sum()).collect();
49 let mut col_sums = vec![0.0; cols];
50 for r in table {
51 for (j, &v) in r.iter().enumerate() {
52 col_sums[j] += v;
53 }
54 }
55 let total: f64 = row_sums.iter().sum();
56 if total <= 0.0 {
57 return None;
58 }
59 let mut stat = 0.0;
60 for (i, r) in table.iter().enumerate() {
61 for (j, &o) in r.iter().enumerate() {
62 let e = row_sums[i] * col_sums[j] / total;
63 if e > 0.0 {
64 stat += (o - e).powi(2) / e;
65 }
66 }
67 }
68 let dof = ((rows - 1) * (cols - 1)) as f64;
69 Some(ChiSquareResult {
70 statistic: stat,
71 p_value: chi_squared::upper_p(stat, dof),
72 dof,
73 })
74}
75
76#[cfg(test)]
77mod tests {
78 use super::*;
79
80 #[test]
81 fn gof_fair_die_is_not_rejected() {
82 let observed = [9.0, 11.0, 10.0, 12.0, 8.0, 10.0];
84 let expected = [10.0; 6];
85 let r = chi_square_gof(&observed, &expected).unwrap();
86 assert_eq!(r.dof, 5.0);
87 assert!(
88 r.p_value > 0.5,
89 "fair die should not be rejected: p={}",
90 r.p_value
91 );
92 }
93
94 #[test]
95 fn gof_loaded_die_is_rejected() {
96 let observed = [5.0, 5.0, 5.0, 5.0, 5.0, 35.0]; let expected = [10.0; 6];
98 let r = chi_square_gof(&observed, &expected).unwrap();
99 assert!(r.statistic > 11.07, "statistic {}", r.statistic); assert!(r.p_value < 0.001);
101 }
102
103 #[test]
104 fn independence_known_example() {
105 let table: [&[f64]; 2] = [&[90.0, 60.0], &[30.0, 120.0]];
107 let r = chi_square_independence(&table).unwrap();
108 assert_eq!(r.dof, 1.0);
109 assert!(r.statistic > 30.0, "statistic {}", r.statistic);
111 assert!(r.p_value < 1e-6);
112 }
113
114 #[test]
115 fn independence_of_independent_table() {
116 let table: [&[f64]; 2] = [&[10.0, 20.0], &[20.0, 40.0]];
118 let r = chi_square_independence(&table).unwrap();
119 assert!(r.statistic < 1e-9);
120 assert!((r.p_value - 1.0).abs() < 1e-6);
121 }
122
123 #[test]
124 fn guards_bad_shapes() {
125 assert!(chi_square_gof(&[1.0], &[1.0]).is_none());
126 assert!(chi_square_gof(&[1.0, 2.0], &[1.0, 0.0]).is_none()); let ragged: [&[f64]; 2] = [&[1.0, 2.0], &[3.0]];
128 assert!(chi_square_independence(&ragged).is_none());
129 }
130}