solow_stats/
contingency.rs1use ndarray::Array2;
10use solow_distributions::chi2_sf;
11
12#[derive(Debug, Clone)]
14pub struct ContingencyResult {
15 pub statistic: f64,
17 pub df: usize,
19 pub pvalue: f64,
21 pub expected: Array2<f64>,
23}
24
25#[derive(Debug, Clone)]
27pub struct Table {
28 table: Array2<f64>,
29}
30
31impl Table {
32 pub fn new(table: Array2<f64>) -> Self {
34 Table { table }
35 }
36
37 pub fn marginal_probabilities(&self) -> (Vec<f64>, Vec<f64>) {
39 let n = self.table.sum();
40 let row: Vec<f64> = self
41 .table
42 .sum_axis(ndarray::Axis(1))
43 .iter()
44 .map(|&v| v / n)
45 .collect();
46 let col: Vec<f64> = self
47 .table
48 .sum_axis(ndarray::Axis(0))
49 .iter()
50 .map(|&v| v / n)
51 .collect();
52 (row, col)
53 }
54
55 pub fn fittedvalues(&self) -> Array2<f64> {
57 let (row, col) = self.marginal_probabilities();
58 let total = self.table.sum();
59 let (r, c) = self.table.dim();
60 let mut fit = Array2::<f64>::zeros((r, c));
61 for i in 0..r {
62 for j in 0..c {
63 fit[[i, j]] = total * row[i] * col[j];
64 }
65 }
66 fit
67 }
68
69 pub fn test_nominal_association(&self) -> ContingencyResult {
71 let expected = self.fittedvalues();
72 let (r, c) = self.table.dim();
73 let mut statistic = 0.0;
74 for i in 0..r {
75 for j in 0..c {
76 let e = expected[[i, j]];
77 let o = self.table[[i, j]];
78 let resid = (o - e) / e.sqrt();
79 statistic += resid * resid;
80 }
81 }
82 let df = (r - 1) * (c - 1);
83 let pvalue = chi2_sf(statistic, df as f64);
84 ContingencyResult {
85 statistic,
86 df,
87 pvalue,
88 expected,
89 }
90 }
91}
92
93#[cfg(test)]
94mod tests {
95 use super::*;
96 use ndarray::array;
97
98 #[test]
99 fn expected_row_sums_match_observed() {
100 let t = Table::new(array![[10.0, 20.0, 30.0], [6.0, 9.0, 17.0]]);
101 let exp = t.fittedvalues();
102 for i in 0..2 {
104 let so: f64 = (0..3).map(|j| t.table[[i, j]]).sum();
105 let se: f64 = (0..3).map(|j| exp[[i, j]]).sum();
106 assert!((so - se).abs() < 1e-9);
107 }
108 }
109
110 #[test]
111 fn statistic_nonnegative() {
112 let t = Table::new(array![[10.0, 20.0, 30.0], [6.0, 9.0, 17.0]]);
113 let res = t.test_nominal_association();
114 assert!(res.statistic >= 0.0);
115 assert_eq!(res.df, 2);
116 }
117}