abax 0.1.33

A lightweight Rust library providing high-precision mathematical constants and special functions, including Bernoulli numbers, Riemann Zeta values, robust incomplete gamma functions, and probability distribution functions like normal and lognormal CDF.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

use crate::{betainc, normcdf};
use std::f64::consts::PI;

/// Student's T cumulative distribution function (CDF).
///
/// Given a value `x` and degrees of freedom `v`, this function returns
/// the probability that a Student's T random variable is less than or equal to `x`.
///
/// # Mathematical Definition
/// For a Student's T distribution with `v` degrees of freedom:
/// - If `x <= 0`: `F(x; v) = 0.5 * I_{v/(v+x^2)}(v/2, 0.5)`
/// - If `x > 0`: `F(x; v) = 1 - 0.5 * I_{v/(v+x^2)}(v/2, 0.5)`
///
/// where `I_z(a, b)` is the regularized incomplete beta function.
///
/// # Domain
/// - `v > 0`
/// - Returns `NaN` if `v <= 0` or any input is `NaN`.
/// - As `v` approaches infinity, the distribution approaches the standard normal distribution.
///
/// # Examples
/// ```
/// use abax::tcdf;
///
/// // Median of any T distribution is 0.5
/// assert!((tcdf(0.0, 5.0, false) - 0.5).abs() < 1e-15);
///
/// // v=1 is the standard Cauchy distribution
/// assert!((tcdf(1.0, 1.0, false) - 0.75).abs() < 1e-15);
/// ```
pub fn tcdf(x: f64, v: f64, upper: bool) -> f64 {
    if x.is_nan() || v.is_nan() || v <= 0.0 {
        return f64::NAN;
    }

    // Student's T is symmetric about 0.
    // P(T > x) = P(T < -x)
    let effective_x = if upper { -x } else { x };

    if v == 1.0 {
        // Special case: Standard Cauchy Distribution
        return 0.5 + effective_x.atan() / PI;
    }

    const NORM_CUTOFF: f64 = 1e7;
    if v > NORM_CUTOFF || v.is_infinite() {
        // Normal approximation for very large degrees of freedom
        return normcdf(effective_x, 0.0, 1.0, false);
    }

    if effective_x == 0.0 {
        return 0.5;
    }

    let abs_x = effective_x.abs();
    let xsq = abs_x * abs_x;

    // Use the relationship between Student's T and the Incomplete Beta function.
    // To maintain precision:
    // If v < x^2, v/(v+x^2) is small.
    // If v >= x^2, x^2/(v+x^2) is small.
    let p_low = if v < xsq {
        betainc(v / (v + xsq), v / 2.0, 0.5, true) / 2.0
    } else {
        // I_z(a, b) = 1 - I_{1-z}(b, a)
        // Here z = x^2/(v+x^2), 1-z = v/(v+x^2)
        // We calculate 0.5 * (1 - I_{x^2/(v+x^2)}(0.5, v/2))
        betainc(xsq / (v + xsq), 0.5, v / 2.0, false) / 2.0
    };

    if effective_x > 0.0 {
        1.0 - p_low
    } else {
        p_low
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_tcdf_cauchy() {
        let tol = 1e-15;
        assert!((tcdf(1.0, 1.0, false) - 0.75).abs() < tol);
        assert!((tcdf(-1.0, 1.0, false) - 0.25).abs() < tol);
        assert!((tcdf(0.0, 1.0, false) - 0.5).abs() < tol);
    }


    #[test]
    fn test_tcdf_symmetry_and_upper() {
        let x = 1.5;
        let v = 5.0;
        let p_lower = tcdf(x, v, false);
        let p_upper = tcdf(x, v, true);
        assert!((p_lower + p_upper - 1.0).abs() < 1e-15);
        assert_eq!(p_upper, tcdf(-x, v, false));
    }

    #[test]
    fn test_tcdf_limits() {
        assert_eq!(tcdf(f64::INFINITY, 10.0, false), 1.0);
        assert_eq!(tcdf(f64::NEG_INFINITY, 10.0, false), 0.0);
        assert!(tcdf(1.0, 0.0, false).is_nan());
        assert!(tcdf(1.0, -1.0, false).is_nan());
    }

    #[test]
    fn test_tcdf_large_v() {
        // Should be very close to Normal CDF
        let x = 1.96;
        let t = tcdf(x, 1e8, false);
        let n = normcdf(x, 0.0, 1.0, false);
        assert!((t - n).abs() < 1e-12);
    }
}