abax 0.1.29

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
118

use crate::{betaincinv, norminv};
use std::f64::consts::PI;

/// Inverse of Student's T cumulative distribution function (CDF).
///
/// Given a probability `p` and degrees of freedom `v`, this function returns
/// the value `x` such that the probability of a Student's T random variable
/// being less than or equal to `x` is `p`.
///
/// # Mathematical Definition
/// - For `v = 1`, it is the standard Cauchy quantile function: `tan(pi * (p - 0.5))`.
/// - For other `v`, it uses the relationship with the Inverse Incomplete Beta function:
///   `x = sign(p - 0.5) * sqrt(v * (1 - z) / z)`, where `z` is the root of the
///   regularized incomplete beta function.
/// - For large `v`, it uses the Cornish-Fisher expansion (Abramowitz & Stegun 26.7.5).
///
/// # Domain
/// - `0 <= p <= 1`
/// - `v > 0`
/// - Returns `NaN` if `p` is out of range, `v` is non-positive, or any input is `NaN`.
///
/// # Examples
/// ```
/// use abax::tinv;
///
/// // Median of any T distribution is 0
/// assert!((tinv(0.5, 5.0) - 0.0).abs() < 1e-15);
///
/// // v=1 is the standard Cauchy distribution: p=0.75 -> x=1
/// assert!((tinv(0.75, 1.0) - 1.0).abs() < 1e-15);
/// ```
pub fn tinv(p: f64, v: f64) -> f64 {
    if p.is_nan() || v.is_nan() || !(0.0..=1.0).contains(&p) || v <= 0.0 {
        return f64::NAN;
    }

    if p == 0.0 {
        return f64::NEG_INFINITY;
    }
    if p == 1.0 {
        return f64::INFINITY;
    }

    if v == 1.0 {
        // Explicit inversion of Cauchy distribution
        return (PI * (p - 0.5)).tan();
    }

    if v < 1000.0 {
        let q = p - 0.5;
        let abs_2q = 2.0 * q.abs();

        // Use betaincinv to invert the T distribution
        // Relationship: P(|T| < x) = 1 - I_{v/(v+x^2)}(v/2, 0.5)
        let (z, oneminusz) = if q.abs() < 0.25 {
            // Central region: compute 1-z directly to avoid roundoff
            let omz = betaincinv(abs_2q, 0.5, v / 2.0, true);
            (1.0 - omz, omz)
        } else {
            // Tail region
            let z_val = betaincinv(abs_2q, v / 2.0, 0.5, false);
            (z_val, 1.0 - z_val)
        };

        return q.signum() * (v * (oneminusz / z)).sqrt();
    }

    // For large degrees of freedom, use Abramowitz & Stegun formula 26.7.5
    let xn = norminv(p, 0.0, 1.0);
    let xn2 = xn * xn;
    let xn3 = xn2 * xn;
    let xn5 = xn3 * xn2;
    let xn7 = xn5 * xn2;
    let xn9 = xn7 * xn2;

    xn + (xn3 + xn) / (4.0 * v)
        + (5.0 * xn5 + 16.0 * xn3 + 3.0 * xn) / (96.0 * v * v)
        + (3.0 * xn7 + 19.0 * xn5 + 17.0 * xn3 - 15.0 * xn) / (384.0 * v * v * v)
        + (79.0 * xn9 + 776.0 * xn7 + 1482.0 * xn5 - 1920.0 * xn3 - 945.0 * xn)
            / (92160.0 * v * v * v * v)
}

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

    #[test]
    fn test_tinv_basic() {
        let tol = 1e-14;
        assert!((tinv(0.5, 5.0) - 0.0).abs() < tol);
        // Reference value for v=5, p=0.975 (Standard 95% CI critical value)
        assert!((tinv(0.975, 5.0) - 2.570581835636315).abs() < tol);
    }

    #[test]
    fn test_tinv_cauchy() {
        let tol = 1e-15;
        assert!((tinv(0.75, 1.0) - 1.0).abs() < tol);
        assert!((tinv(0.25, 1.0) + 1.0).abs() < tol);
    }

    #[test]
    fn test_tinv_large_v() {
        // Should approach norminv
        let p = 0.95;
        let t = tinv(p, 1e6);
        let n = 1.6448551507220405;
        assert!((t - n).abs() < 1e-10);
    }

    #[test]
    fn test_tinv_limits() {
        assert_eq!(tinv(0.0, 10.0), f64::NEG_INFINITY);
        assert_eq!(tinv(1.0, 10.0), f64::INFINITY);
        assert!(tinv(0.5, 0.0).is_nan());
        assert!(tinv(-0.1, 5.0).is_nan());
    }
}