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
86
87
88
89

use crate::gampdf;

/// Chi-squared probability density function (PDF).
///
/// Returns the probability density at `x` for the Chi-squared distribution
/// with <math><mi>v</mi></math> degrees of freedom.
///
/// # Mathematical Definition
/// The Chi-squared distribution with <math><mi>v</mi></math> degrees of freedom is a special case
/// of the Gamma distribution with shape <math><mi>a</mi><mo>=</mo><mi>v</mi><mo>/</mo><mn>2</mn></math> and scale <math><mi>b</mi><mo>=</mo><mn>2</mn></math>:
/// <math display="block">
///   <mi>f</mi><mo>(</mo><mi>x</mi><mo>;</mo><mi>v</mi><mo>)</mo>
///   <mo>=</mo>
///   <mfrac>
///     <mrow>
///       <msup><mi>x</mi><mrow><mi>v</mi><mo>/</mo><mn>2</mn><mo>-</mo><mn>1</mn></mrow></msup>
///       <msup><mi>e</mi><mrow><mo>-</mo><mi>x</mi><mo>/</mo><mn>2</mn></mrow></msup>
///     </mrow>
///     <mrow>
///       <msup><mn>2</mn><mrow><mi>v</mi><mo>/</mo><mn>2</mn></mrow></msup>
///       <mi>Γ</mi><mo>(</mo><mi>v</mi><mo>/</mo><mn>2</mn><mo>)</mo>
///     </mrow>
///   </mfrac>
/// </math>
/// for <math><mi>x</mi><mo>&gt;</mo><mn>0</mn></math>.
///
/// # Examples
/// ```
/// use abax::chi2pdf;
///
/// // For v=2, the distribution is Exponential with mean 2: f(x) = 0.5 * exp(-x/2)
/// let x = 2.0;
/// let pdf = chi2pdf(x, 2.0);
/// assert!((pdf - 0.5 * (-1.0f64).exp()).abs() < 1e-15);
/// ```
pub fn chi2pdf(x: f64, v: f64) -> f64 {
    gampdf(x, v / 2.0, 2.0)
}

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

    #[test]
    fn test_chi2pdf_standard_cases() {
        let tol = 1e-14;

        // v = 2: Exponential(0.5) => 0.5 * exp(-x/2)
        // at x = 0, pdf = 0.5
        assert!((chi2pdf(0.0, 2.0) - 0.5).abs() < tol);
        // at x = 2, pdf = 0.5 * exp(-1)
        assert!((chi2pdf(2.0, 2.0) - 0.5 * (-1.0f64).exp()).abs() < tol);

        // v = 4: f(x) = 0.25 * x * exp(-x/2)
        // at x = 2, pdf = 0.25 * 2 * exp(-1) = 0.5 * exp(-1)
        assert!((chi2pdf(2.0, 4.0) - 0.5 * (-1.0f64).exp()).abs() < tol);
    }

    #[test]
    fn test_chi2pdf_v1_relation() {
        // v = 1: f(x) = (1/sqrt(2*pi)) * x^(-1/2) * exp(-x/2)
        // at x = 1, pdf = (1/sqrt(2*pi)) * exp(-0.5)
        let expected = (1.0 / (2.0 * std::f64::consts::PI).sqrt()) * (-0.5f64).exp();
        assert!((chi2pdf(1.0, 1.0) - expected).abs() < 1e-15);
    }

    #[test]
    fn test_chi2pdf_poles_at_zero() {
        // if v < 2, pdf(0) is infinity
        assert_eq!(chi2pdf(0.0, 1.0), f64::INFINITY);
        // if v = 2, pdf(0) is 0.5
        assert_eq!(chi2pdf(0.0, 2.0), 0.5);
        // if v > 2, pdf(0) is 0
        assert_eq!(chi2pdf(0.0, 3.0), 0.0);
    }

    #[test]
    fn test_chi2pdf_invalid_params() {
        assert!(chi2pdf(1.0, 0.0) == 0.0);
        assert!(chi2pdf(1.0, -1.0).is_nan());
        assert!(chi2pdf(f64::NAN, 2.0).is_nan());
        assert!(chi2pdf(1.0, f64::NAN).is_nan());
    }

    #[test]
    fn test_chi2pdf_negative_x() {
        assert_eq!(chi2pdf(-1.0, 2.0), 0.0);
    }
}