sciforge_lib/maths/statistics/
distributions.rs

1use crate::constants::{
2    ERF_A1, ERF_A2, ERF_A3, ERF_A4, ERF_A5, ERF_P, STIRLING_LANCZOS_COEFFS, STIRLING_LANCZOS_G,
3};
4use std::f64::consts::PI;
5
6pub fn normal_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
7    let z = (x - mu) / sigma;
8    (-0.5 * z * z).exp() / (sigma * (2.0 * PI).sqrt())
9}
10
11pub fn normal_cdf(x: f64, mu: f64, sigma: f64) -> f64 {
12    0.5 * (1.0 + erf((x - mu) / (sigma * std::f64::consts::SQRT_2)))
13}
14
15pub fn erf(x: f64) -> f64 {
16    let t = 1.0 / (1.0 + ERF_P * x.abs());
17    let poly = t * (ERF_A1 + t * (ERF_A2 + t * (ERF_A3 + t * (ERF_A4 + t * ERF_A5))));
18    let result = 1.0 - poly * (-x * x).exp();
19    if x >= 0.0 { result } else { -result }
20}
21
22pub fn erfc(x: f64) -> f64 {
23    1.0 - erf(x)
24}
25
26pub fn poisson_pmf(k: u64, lambda: f64) -> f64 {
27    (-lambda).exp() * lambda.powi(k as i32) / factorial(k) as f64
28}
29
30pub fn binomial_pmf(n: u64, k: u64, p: f64) -> f64 {
31    binomial_coeff(n, k) as f64 * p.powi(k as i32) * (1.0 - p).powi((n - k) as i32)
32}
33
34pub fn exponential_pdf(x: f64, lambda: f64) -> f64 {
35    if x < 0.0 {
36        return 0.0;
37    }
38    lambda * (-lambda * x).exp()
39}
40
41pub fn exponential_cdf(x: f64, lambda: f64) -> f64 {
42    if x < 0.0 {
43        return 0.0;
44    }
45    1.0 - (-lambda * x).exp()
46}
47
48pub fn chi_squared_pdf(x: f64, k: u32) -> f64 {
49    if x <= 0.0 {
50        return 0.0;
51    }
52    let half_k = k as f64 / 2.0;
53    x.powf(half_k - 1.0) * (-x / 2.0).exp() / (2.0_f64.powf(half_k) * gamma(half_k))
54}
55
56pub fn student_t_pdf(x: f64, nu: f64) -> f64 {
57    let coeff = gamma((nu + 1.0) / 2.0) / ((nu * PI).sqrt() * gamma(nu / 2.0));
58    coeff * (1.0 + x * x / nu).powf(-(nu + 1.0) / 2.0)
59}
60
61pub fn cauchy_pdf(x: f64, x0: f64, gamma_val: f64) -> f64 {
62    1.0 / (PI * gamma_val * (1.0 + ((x - x0) / gamma_val).powi(2)))
63}
64
65pub fn uniform_pdf(x: f64, a: f64, b: f64) -> f64 {
66    if x >= a && x <= b { 1.0 / (b - a) } else { 0.0 }
67}
68
69pub fn beta_pdf(x: f64, alpha: f64, beta: f64) -> f64 {
70    if x <= 0.0 || x >= 1.0 {
71        return 0.0;
72    }
73    let b = gamma(alpha) * gamma(beta) / gamma(alpha + beta);
74    x.powf(alpha - 1.0) * (1.0 - x).powf(beta - 1.0) / b
75}
76
77pub fn gamma_pdf(x: f64, shape: f64, rate: f64) -> f64 {
78    if x <= 0.0 {
79        return 0.0;
80    }
81    rate.powf(shape) / gamma(shape) * x.powf(shape - 1.0) * (-rate * x).exp()
82}
83
84pub fn gamma(x: f64) -> f64 {
85    if x <= 0.5 {
86        PI / ((PI * x).sin() * gamma(1.0 - x))
87    } else {
88        let x = x - 1.0;
89        let mut sum = STIRLING_LANCZOS_COEFFS[0];
90        for (i, &ci) in STIRLING_LANCZOS_COEFFS.iter().enumerate().skip(1) {
91            sum += ci / (x + i as f64);
92        }
93        let t = x + STIRLING_LANCZOS_G + 0.5;
94        (2.0 * PI).sqrt() * t.powf(x + 0.5) * (-t).exp() * sum
95    }
96}
97
98pub fn ln_gamma(x: f64) -> f64 {
99    gamma(x).ln()
100}
101
102pub fn beta_function(a: f64, b: f64) -> f64 {
103    gamma(a) * gamma(b) / gamma(a + b)
104}
105
106pub fn incomplete_gamma_lower(a: f64, x: f64, terms: usize) -> f64 {
107    let mut sum = 0.0;
108    let mut term = 1.0 / a;
109    sum += term;
110    for n in 1..terms {
111        term *= x / (a + n as f64);
112        sum += term;
113    }
114    sum * x.powf(a) * (-x).exp()
115}
116
117fn factorial(n: u64) -> u64 {
118    (1..=n).product()
119}
120
121fn binomial_coeff(n: u64, k: u64) -> u64 {
122    if k > n {
123        return 0;
124    }
125    let k = k.min(n - k);
126    let mut result: u64 = 1;
127    for i in 0..k {
128        result = result * (n - i) / (i + 1);
129    }
130    result
131}
sciforge_lib/maths/statistics/distributions.rs

sciforge_lib/maths/statistics/
distributions.rs
