qualia_core_db/solvers/statistics/distributions/
chi_squared.rs

1//! χ² (chi-squared) distribution — pdf / cdf / quantile + upper-tail p-value, used
2//! by the goodness-of-fit and independence tests. CDF is exact via the regularized
3//! lower incomplete gamma `P(k/2, x/2)` ([`super::special::gammp`]).
4
5use super::special::{gammp, gammq, ln_gamma};
6
7/// pdf with `k > 0` degrees of freedom, `x ≥ 0`.
8pub fn pdf(x: f64, k: f64) -> f64 {
9    debug_assert!(k > 0.0);
10    if x < 0.0 {
11        return 0.0;
12    }
13    if x == 0.0 {
14        // Finite only for k = 2 (= 1/2); 0 for k > 2; +∞ for k < 2.
15        return if k < 2.0 {
16            f64::INFINITY
17        } else if (k - 2.0).abs() < 1e-12 {
18            0.5
19        } else {
20            0.0
21        };
22    }
23    let kh = k / 2.0;
24    ((kh - 1.0) * x.ln() - x / 2.0 - kh * std::f64::consts::LN_2 - ln_gamma(kh)).exp()
25}
26
27/// cdf `P(X ≤ x)` = `P(k/2, x/2)`.
28pub fn cdf(x: f64, k: f64) -> f64 {
29    debug_assert!(k > 0.0);
30    if x <= 0.0 {
31        return 0.0;
32    }
33    gammp(k / 2.0, x / 2.0)
34}
35
36/// Upper-tail p-value `P(X ≥ x)` — the usual χ² test p-value.
37pub fn upper_p(x: f64, k: f64) -> f64 {
38    if x <= 0.0 {
39        return 1.0;
40    }
41    gammq(k / 2.0, x / 2.0)
42}
43
44/// Inverse cdf (quantile) for `0 < p < 1`.
45pub fn quantile(p: f64, k: f64) -> f64 {
46    super::invert_cdf(p, Some(0.0), |x| cdf(x, k))
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52
53    #[test]
54    fn cdf_is_exponential_for_k_two() {
55        // χ²(2) = Exp with mean 2 → CDF(x) = 1 - e^{-x/2}.
56        for &x in &[0.5, 2.0, 5.0] {
57            assert!(
58                (cdf(x, 2.0) - (1.0 - (-x / 2.0).exp())).abs() < 1e-10,
59                "x={x}"
60            );
61        }
62    }
63
64    #[test]
65    fn known_critical_values() {
66        // 95th percentile: df=1 → 3.841; df=10 → 18.307.
67        assert!((quantile(0.95, 1.0) - 3.841_458_82).abs() < 1e-4);
68        assert!((quantile(0.95, 10.0) - 18.307_038_05).abs() < 1e-3);
69        // Upper-tail p of 3.841 at df=1 ≈ 0.05.
70        assert!((upper_p(3.841_458_82, 1.0) - 0.05).abs() < 1e-6);
71    }
72
73    #[test]
74    fn quantile_inverts_cdf() {
75        for &p in &[0.05, 0.5, 0.9, 0.99] {
76            let x = quantile(p, 7.0);
77            assert!((cdf(x, 7.0) - p).abs() < 1e-8, "p={p}");
78        }
79    }
80}
qualia_core_db/solvers/statistics/distributions/chi_squared.rs

qualia_core_db/solvers/statistics/distributions/
chi_squared.rs
