qualia_core_db/solvers/statistics/distributions/
students_t.rs

1//! Student's t-distribution — pdf / cdf / quantile + the two-sided p-value the
2//! t-tests use. The CDF is exact via the regularized incomplete beta
3//! ([`super::special::betai`]); the quantile inverts it numerically.
4
5use super::special::{betai, ln_gamma};
6
7/// pdf of the t-distribution with `nu > 0` degrees of freedom.
8pub fn pdf(t: f64, nu: f64) -> f64 {
9    debug_assert!(nu > 0.0);
10    let c = (ln_gamma((nu + 1.0) / 2.0) - ln_gamma(nu / 2.0)).exp()
11        / (nu * std::f64::consts::PI).sqrt();
12    c * (1.0 + t * t / nu).powf(-(nu + 1.0) / 2.0)
13}
14
15/// cdf `P(T ≤ t)` with `nu` degrees of freedom. Uses
16/// `I_x(ν/2, ½)` with `x = ν/(ν+t²)` and the sign of `t`.
17pub fn cdf(t: f64, nu: f64) -> f64 {
18    debug_assert!(nu > 0.0);
19    let x = nu / (nu + t * t);
20    let ib = 0.5 * betai(nu / 2.0, 0.5, x);
21    if t >= 0.0 {
22        1.0 - ib
23    } else {
24        ib
25    }
26}
27
28/// Two-sided p-value for a t statistic: `2·(1 − P(T ≤ |t|)) = I_x(ν/2, ½)`.
29pub fn two_sided_p(t: f64, nu: f64) -> f64 {
30    let x = nu / (nu + t * t);
31    betai(nu / 2.0, 0.5, x)
32}
33
34/// One-sided upper-tail p-value `P(T ≥ t)`.
35pub fn upper_p(t: f64, nu: f64) -> f64 {
36    1.0 - cdf(t, nu)
37}
38
39/// Inverse cdf (quantile) for `0 < p < 1`.
40pub fn quantile(p: f64, nu: f64) -> f64 {
41    super::invert_cdf(p, None, |t| cdf(t, nu))
42}
43
44#[cfg(test)]
45mod tests {
46    use super::*;
47
48    #[test]
49    fn cdf_symmetry_and_center() {
50        assert!((cdf(0.0, 5.0) - 0.5).abs() < 1e-12);
51        // Symmetry: F(-t) = 1 - F(t).
52        for &(t, nu) in &[(1.3, 7.0), (2.5, 12.0), (0.8, 3.0)] {
53            assert!((cdf(-t, nu) - (1.0 - cdf(t, nu))).abs() < 1e-10);
54        }
55    }
56
57    #[test]
58    fn known_critical_values() {
59        // Two-sided 95% critical t: df=10 → 2.228; df=∞ → 1.96.
60        assert!((quantile(0.975, 10.0) - 2.228_138_851).abs() < 1e-4);
61        assert!((quantile(0.975, 1.0) - 12.706_204_736).abs() < 1e-3); // Cauchy
62        assert!((quantile(0.975, 1_000_000.0) - 1.959_963_98).abs() < 1e-3);
63    }
64
65    #[test]
66    fn two_sided_p_matches_tail() {
67        // p-value of t=2.228 at df=10 ≈ 0.05.
68        assert!((two_sided_p(2.228_138_851, 10.0) - 0.05).abs() < 1e-4);
69        assert!((two_sided_p(0.0, 5.0) - 1.0).abs() < 1e-12);
70    }
71
72    #[test]
73    fn quantile_inverts_cdf() {
74        for &p in &[0.01, 0.05, 0.5, 0.9, 0.99] {
75            let t = quantile(p, 8.0);
76            assert!((cdf(t, 8.0) - p).abs() < 1e-8, "p={p}");
77        }
78    }
79}
qualia_core_db/solvers/statistics/distributions/students_t.rs

qualia_core_db/solvers/statistics/distributions/
students_t.rs
