sim_lib_numbers_stats/
parametric_distribution.rs

1//! Normal and Student-t distributions composed from shared numerical owners.
2
3use super::{StatsError, StatsResult};
4use sim_lib_numbers_special::{erfc, log_gamma, regularized_beta};
5
6fn valid_probability(metric: &'static str, value: f64) -> StatsResult<()> {
7    if value.is_finite() && (0.0..=1.0).contains(&value) {
8        Ok(())
9    } else {
10        Err(StatsError::ProbabilityOutOfRange {
11            metric,
12            index: None,
13            value,
14        })
15    }
16}
17fn quantile_root(mut lower: f64, mut upper: f64, p: f64, cdf: impl Fn(f64) -> f64) -> f64 {
18    for _ in 0..192 {
19        let middle = lower + (upper - lower) / 2.0;
20        if cdf(middle) < p {
21            lower = middle;
22        } else {
23            upper = middle;
24        }
25        if upper - lower <= 2.0 * f64::EPSILON * middle.abs().max(1.0) {
26            break;
27        }
28    }
29    lower + (upper - lower) / 2.0
30}
31
32/// Standard normal probability density.
33pub fn normal_density(x: f64) -> f64 {
34    (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt()
35}
36/// Standard normal cumulative probability.
37pub fn normal_cdf(x: f64) -> f64 {
38    0.5 * erfc(-x / std::f64::consts::SQRT_2).value
39}
40/// Standard normal survival probability, evaluated directly in the positive tail.
41pub fn normal_survival(x: f64) -> f64 {
42    0.5 * erfc(x / std::f64::consts::SQRT_2).value
43}
44/// Standard normal quantile, found by bounded monotone inversion.
45pub fn normal_quantile(p: f64) -> StatsResult<f64> {
46    valid_probability("normal_quantile", p)?;
47    if p == 0.0 {
48        return Ok(f64::NEG_INFINITY);
49    }
50    if p == 1.0 {
51        return Ok(f64::INFINITY);
52    }
53    Ok(quantile_root(-40.0, 40.0, p, normal_cdf))
54}
55
56/// Student-t probability density for positive degrees of freedom.
57pub fn student_t_density(x: f64, v: f64) -> StatsResult<f64> {
58    if !(v > 0.0 && v.is_finite()) {
59        return Err(StatsError::InvalidControl {
60            field: "degrees_of_freedom",
61            reason: "must be finite and positive",
62        });
63    }
64    let lg1 = log_gamma((v + 1.0) / 2.0)
65        .map_err(|_| StatsError::InvalidControl {
66            field: "degrees_of_freedom",
67            reason: "gamma evaluation failed",
68        })?
69        .value;
70    let lg0 = log_gamma(v / 2.0)
71        .map_err(|_| StatsError::InvalidControl {
72            field: "degrees_of_freedom",
73            reason: "gamma evaluation failed",
74        })?
75        .value;
76    Ok(
77        (lg1 - lg0 - 0.5 * (v * std::f64::consts::PI).ln() - 0.5 * (v + 1.0) * (x * x / v).ln_1p())
78            .exp(),
79    )
80}
81fn student_tail(x: f64, v: f64) -> StatsResult<f64> {
82    if !(v > 0.0 && v.is_finite()) {
83        return Err(StatsError::InvalidControl {
84            field: "degrees_of_freedom",
85            reason: "must be finite and positive",
86        });
87    }
88    regularized_beta(v / 2.0, 0.5, v / (v + x * x))
89        .map(|r| 0.5 * r.value)
90        .map_err(|_| StatsError::InvalidControl {
91            field: "degrees_of_freedom",
92            reason: "beta evaluation failed",
93        })
94}
95/// Student-t cumulative probability.
96pub fn student_t_cdf(x: f64, v: f64) -> StatsResult<f64> {
97    let tail = student_tail(x, v)?;
98    Ok(if x >= 0.0 { 1.0 - tail } else { tail })
99}
100/// Student-t survival probability, using the direct beta tail for positive values.
101pub fn student_t_survival(x: f64, v: f64) -> StatsResult<f64> {
102    let tail = student_tail(x, v)?;
103    Ok(if x >= 0.0 { tail } else { 1.0 - tail })
104}
105/// Student-t quantile, found by bounded monotone inversion.
106pub fn student_t_quantile(p: f64, v: f64) -> StatsResult<f64> {
107    valid_probability("student_t_quantile", p)?;
108    if p == 0.0 {
109        return Ok(f64::NEG_INFINITY);
110    }
111    if p == 1.0 {
112        return Ok(f64::INFINITY);
113    }
114    student_t_density(0.0, v)?;
115    let mut bound = 1.0;
116    while (student_t_cdf(-bound, v)? > p || student_t_cdf(bound, v)? < p) && bound < 1e16 {
117        bound *= 2.0;
118    }
119    Ok(quantile_root(-bound, bound, p, |x| {
120        student_t_cdf(x, v).unwrap_or(f64::NAN)
121    }))
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127    fn close(a: f64, b: f64, t: f64) {
128        assert!((a - b).abs() < t * b.abs().max(1.0), "{a} != {b}");
129    }
130    #[test]
131    fn normal_center_tails_and_inversion() {
132        close(normal_density(0.0), 0.3989422804014327, 1e-15);
133        close(normal_cdf(0.0), 0.5, 1e-15);
134        assert!(normal_survival(10.0) > 0.0);
135        for p in [1e-10, 0.01, 0.5, 0.99, 1.0 - 1e-10] {
136            close(normal_cdf(normal_quantile(p).unwrap()), p, 3e-8);
137        }
138    }
139    #[test]
140    fn student_symmetry_tails_and_inversion() {
141        for v in [1.0, 5.0, 30.0] {
142            close(student_t_cdf(0.0, v).unwrap(), 0.5, 1e-14);
143            assert!(student_t_survival(20.0, v).unwrap() > 0.0);
144            for p in [1e-5, 0.1, 0.9, 1.0 - 1e-5] {
145                let x = student_t_quantile(p, v).unwrap();
146                close(student_t_cdf(x, v).unwrap(), p, 2e-11);
147            }
148        }
149    }
150}
sim_lib_numbers_stats/parametric_distribution.rs

sim_lib_numbers_stats/
parametric_distribution.rs
