abax 0.1.39

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

use crate::erfc;
use std::f64::consts::SQRT_2;

/// Normal cumulative distribution function (CDF).
///
/// Given a value `x`, a mean `mu`, and a standard deviation `sigma`,
/// this function returns the probability that a normal random variable
/// is less than or equal to `x`.
///
/// # Mathematical Definition
/// For a normal distribution with mean <math><mi>μ</mi></math> and standard deviation <math><mi>σ</mi></math>:
/// - Lower tail (`upper = false`): <math><mi>Φ</mi><mo>(</mo><mi>x</mi><mo>;</mo><mi>μ</mi><mo>,</mo><mi>σ</mi><mo>)</mo><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>erfc</mi><mo>(</mo><mo>-</mo><mfrac><mrow><mi>x</mi><mo>-</mo><mi>μ</mi></mrow><mrow><mi>σ</mi><msqrt><mn>2</mn></msqrt></mrow></mfrac><mo>)</mo></math>
/// - Upper tail (`upper = true`): <math><mn>1</mn><mo>-</mo><mi>Φ</mi><mo>(</mo><mi>x</mi><mo>;</mo><mi>μ</mi><mo>,</mo><mi>σ</mi><mo>)</mo><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>erfc</mi><mo>(</mo><mfrac><mrow><mi>x</mi><mo>-</mo><mi>μ</mi></mrow><mrow><mi>σ</mi><msqrt><mn>2</mn></msqrt></mrow></mfrac><mo>)</mo></math>
///
/// # Domain
/// - `sigma >= 0`
/// - Returns `NaN` if `sigma` is negative or any input is `NaN`.
/// - If `sigma` is 0, the distribution is treated as a Dirac delta concentrated at `mu`.
///   - Returns 0 (lower) / 1 (upper) if `x < mu`.
///   - Returns 1 (lower) / 0 (upper) if `x >= mu`.
///
/// # Examples
/// ```
/// use abax::normcdf;
///
/// // Standard normal median is 0.5
/// assert!((normcdf(0.0, 0.0, 1.0, false) - 0.5).abs() < 1e-15);
/// // One sigma upper bound
/// assert!((normcdf(1.0, 0.0, 1.0, false) - 0.8413447460685429).abs() < 1e-15);
/// // Upper tail one sigma
/// assert!((normcdf(1.0, 0.0, 1.0, true) - 0.15865525393145705).abs() < 1e-15);
/// ```
pub fn normcdf(x: f64, mu: f64, sigma: f64, upper: bool) -> f64 {
    if x.is_nan() || mu.is_nan() || sigma.is_nan() || sigma < 0.0 {
        return f64::NAN;
    }
 
    if sigma == 0.0 {
        return if upper {
            if x < mu { 1.0 } else { 0.0 }
        } else {
            if x < mu { 0.0 } else { 1.0 }
        };
    }
 
    let z = (x - mu) / sigma;
    if upper {
        0.5 * erfc(z / SQRT_2)
    } else {
        0.5 * erfc(-z / SQRT_2)
    }
}

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

    #[test]
    fn test_normcdf_standard() {
        let tol = 1e-14;
        assert!((normcdf(0.0, 0.0, 1.0, false) - 0.5).abs() < tol);
        assert!((normcdf(1.959963984540054, 0.0, 1.0, false) - 0.975).abs() < tol);
        assert!((normcdf(-1.959963984540054, 0.0, 1.0, false) - 0.025).abs() < tol);
    }
 
    #[test]
    fn test_normcdf_upper() {
        let tol = 1e-14;
        assert!((normcdf(0.0, 0.0, 1.0, true) - 0.5).abs() < tol);
        assert!((normcdf(1.959963984540054, 0.0, 1.0, true) - 0.025).abs() < tol);
        assert!((normcdf(-1.959963984540054, 0.0, 1.0, true) - 0.975).abs() < tol);
    }

    #[test]
    fn test_normcdf_zero_sigma() {
        assert_eq!(normcdf(5.0, 5.0, 0.0, false), 1.0);
        assert_eq!(normcdf(4.99, 5.0, 0.0, false), 0.0);
        assert_eq!(normcdf(4.99, 5.0, 0.0, true), 1.0);
    }

    #[test]
    fn test_normcdf_invalid() {
        assert!(normcdf(0.5, 0.0, -1.0, false).is_nan());
    }
}