modewise_stats/
dist.rs

1//! Student's t distribution through the regularised incomplete beta function.
2
3/// Upper tail probability P(T > t) for Student's t with `df` degrees of freedom.
4pub fn t_upper(t: f64, df: f64) -> f64 {
5    if t.is_nan() || df.is_nan() || df <= 0.0 {
6        return f64::NAN;
7    }
8    if t.is_infinite() {
9        return if t > 0.0 { 0.0 } else { 1.0 };
10    }
11    let x = df / (df + t * t);
12    let half = 0.5 * incomplete_beta(0.5 * df, 0.5, x);
13    if t > 0.0 { half } else { 1.0 - half }
14}
15
16/// Regularised incomplete beta I_x(a, b) by its continued fraction (Lentz's method), using the
17/// symmetry I_x(a, b) = 1 − I_{1−x}(b, a) where the fraction converges faster.
18pub fn incomplete_beta(a: f64, b: f64, x: f64) -> f64 {
19    if x <= 0.0 {
20        return 0.0;
21    }
22    if x >= 1.0 {
23        return 1.0;
24    }
25    let front =
26        (ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln()).exp();
27    if x < (a + 1.0) / (a + b + 2.0) {
28        front * continued_fraction(a, b, x) / a
29    } else {
30        1.0 - front * continued_fraction(b, a, 1.0 - x) / b
31    }
32}
33
34fn continued_fraction(a: f64, b: f64, x: f64) -> f64 {
35    const TINY: f64 = 1e-300;
36    let mut c = 1.0;
37    let mut d = 1.0 - (a + b) * x / (a + 1.0);
38    if d.abs() < TINY {
39        d = TINY;
40    }
41    d = 1.0 / d;
42    let mut h = d;
43    for m in 1..=1000 {
44        let m = m as f64;
45        for numerator in [
46            m * (b - m) * x / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)),
47            -(a + m) * (a + b + m) * x / ((a + 2.0 * m) * (a + 2.0 * m + 1.0)),
48        ] {
49            d = 1.0 + numerator * d;
50            if d.abs() < TINY {
51                d = TINY;
52            }
53            c = 1.0 + numerator / c;
54            if c.abs() < TINY {
55                c = TINY;
56            }
57            d = 1.0 / d;
58            h *= d * c;
59        }
60        if (d * c - 1.0).abs() < 1e-15 {
61            break;
62        }
63    }
64    h
65}
66
67/// ln Γ(x) for x > 0, by the Lanczos approximation (g = 7, nine coefficients).
68fn ln_gamma(x: f64) -> f64 {
69    const G: [f64; 9] = [
70        0.999_999_999_999_809_9,
71        676.520_368_121_885_1,
72        -1_259.139_216_722_402_8,
73        771.323_428_777_653_1,
74        -176.615_029_162_140_6,
75        12.507_343_278_686_905,
76        -0.138_571_095_265_720_12,
77        9.984_369_578_019_572e-6,
78        1.505_632_735_149_311_6e-7,
79    ];
80    if x < 0.5 {
81        let pi = std::f64::consts::PI;
82        return (pi / (pi * x).sin()).ln() - ln_gamma(1.0 - x);
83    }
84    let x = x - 1.0;
85    let t = x + 7.5;
86    let series = G[1..]
87        .iter()
88        .enumerate()
89        .fold(G[0], |s, (i, g)| s + g / (x + i as f64 + 1.0));
90    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + series.ln()
91}
92
93#[cfg(test)]
94mod tests {
95    use super::*;
96
97    #[test]
98    fn known_values() {
99        assert!((ln_gamma(5.0) - 24f64.ln()).abs() < 1e-12);
100        assert!((ln_gamma(0.5) - std::f64::consts::PI.sqrt().ln()).abs() < 1e-12);
101        assert!(
102            (ln_gamma(0.25) - 3.625_609_908_221_908f64.ln()).abs() < 1e-12,
103            "reflection below 0.5"
104        );
105        assert_eq!(
106            (
107                incomplete_beta(2.0, 3.0, 0.0),
108                incomplete_beta(2.0, 3.0, 1.0)
109            ),
110            (0.0, 1.0)
111        );
112        // I_x(1, 1) = x and I_x(a, 1) = x^a.
113        assert!((incomplete_beta(1.0, 1.0, 0.3) - 0.3).abs() < 1e-14);
114        assert!((incomplete_beta(3.0, 1.0, 0.6) - 0.216).abs() < 1e-14);
115        assert!((t_upper(0.0, 7.0) - 0.5).abs() < 1e-15);
116        assert_eq!(
117            (t_upper(f64::INFINITY, 3.0), t_upper(f64::NEG_INFINITY, 3.0)),
118            (0.0, 1.0)
119        );
120        assert!(t_upper(f64::NAN, 3.0).is_nan() && t_upper(1.0, 0.0).is_nan());
121        // Cauchy (df = 1): P(T > 1) = 1/4.
122        assert!((t_upper(1.0, 1.0) - 0.25).abs() < 1e-14);
123        assert!((t_upper(-1.0, 1.0) - 0.75).abs() < 1e-14);
124    }
125}
modewise_stats/dist.rs

modewise_stats/
dist.rs
