qualia_core_db/solvers/statistics/distributions/
normal.rs

1//! Normal (Gaussian) distribution — pdf / cdf / quantile, the canonical engine-wide
2//! implementation. The CDF is `½·erfc(−z/√2)` over the shared [`super::special::erfc`]
3//! (full precision), and the quantile is Acklam's rational inverse refined by one
4//! Halley step (≈ machine precision). Domain libraries (`financial_modeling`'s
5//! Black–Scholes, etc.) should call these instead of re-deriving a local `normal_cdf`.
6
7use super::special::erfc;
8use std::f64::consts::{PI, SQRT_2};
9
10/// Standard-normal pdf `φ(z)`.
11pub fn standard_pdf(z: f64) -> f64 {
12    (-(z * z) / 2.0).exp() / (2.0 * PI).sqrt()
13}
14
15/// Standard-normal cdf `Φ(z) = ½·erfc(−z/√2)`.
16pub fn standard_cdf(z: f64) -> f64 {
17    0.5 * erfc(-z / SQRT_2)
18}
19
20/// Normal pdf with mean `mu`, std-dev `sigma` (> 0).
21pub fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
22    debug_assert!(sigma > 0.0);
23    standard_pdf((x - mu) / sigma) / sigma
24}
25
26/// Normal cdf with mean `mu`, std-dev `sigma` (> 0).
27pub fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
28    standard_cdf((x - mu) / sigma)
29}
30
31/// Inverse standard-normal cdf `Φ⁻¹(p)`, `0 < p < 1` (Acklam + one Halley refinement).
32/// `±∞` at the endpoints.
33pub fn standard_quantile(p: f64) -> f64 {
34    if p <= 0.0 {
35        return f64::NEG_INFINITY;
36    }
37    if p >= 1.0 {
38        return f64::INFINITY;
39    }
40    const A: [f64; 6] = [
41        -3.969_683_028_665_376e1,
42        2.209_460_984_245_205e2,
43        -2.759_285_104_469_687e2,
44        1.383_577_518_672_690e2,
45        -3.066_479_806_614_716e1,
46        2.506_628_277_459_239e0,
47    ];
48    const B: [f64; 5] = [
49        -5.447_609_879_822_406e1,
50        1.615_858_368_580_409e2,
51        -1.556_989_798_598_866e2,
52        6.680_131_188_771_972e1,
53        -1.328_068_155_288_572e1,
54    ];
55    const C: [f64; 6] = [
56        -7.784_894_002_430_293e-3,
57        -3.223_964_580_411_365e-1,
58        -2.400_758_277_161_838e0,
59        -2.549_732_539_343_734e0,
60        4.374_664_141_464_968e0,
61        2.938_163_982_698_783e0,
62    ];
63    const D: [f64; 4] = [
64        7.784_695_709_041_462e-3,
65        3.224_671_290_700_398e-1,
66        2.445_134_137_142_996e0,
67        3.754_408_661_907_416e0,
68    ];
69    const PLOW: f64 = 0.02425;
70    const PHIGH: f64 = 1.0 - PLOW;
71
72    let mut x = if p < PLOW {
73        let q = (-2.0 * p.ln()).sqrt();
74        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
75            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
76    } else if p <= PHIGH {
77        let q = p - 0.5;
78        let r = q * q;
79        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
80            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
81    } else {
82        let q = (-2.0 * (1.0 - p).ln()).sqrt();
83        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
84            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
85    };
86    // One Halley refinement step (uses the full-precision cdf/pdf).
87    let e = standard_cdf(x) - p;
88    let u = e * (2.0 * PI).sqrt() * (x * x / 2.0).exp();
89    x -= u / (1.0 + x * u / 2.0);
90    x
91}
92
93/// Inverse normal cdf with mean `mu`, std-dev `sigma`.
94pub fn quantile(p: f64, mu: f64, sigma: f64) -> f64 {
95    mu + sigma * standard_quantile(p)
96}
97
98/// Two-sided p-value for a standard-normal (z) statistic: `2·(1 − Φ(|z|))`.
99pub fn two_sided_p(z: f64) -> f64 {
100    2.0 * (1.0 - standard_cdf(z.abs()))
101}
102
103#[cfg(test)]
104mod tests {
105    use super::*;
106    const TOL: f64 = 1e-9;
107
108    #[test]
109    fn cdf_known_quantiles() {
110        assert!((standard_cdf(0.0) - 0.5).abs() < TOL);
111        assert!((standard_cdf(1.0) - 0.841_344_746_068_543).abs() < 1e-9);
112        assert!((standard_cdf(1.959_963_984_540_054) - 0.975).abs() < 1e-9);
113        assert!((standard_cdf(-2.0) - 0.022_750_131_948_179).abs() < 1e-9);
114    }
115
116    #[test]
117    fn pdf_peak_and_symmetry() {
118        assert!((standard_pdf(0.0) - 1.0 / (2.0 * PI).sqrt()).abs() < TOL);
119        assert!((standard_pdf(1.5) - standard_pdf(-1.5)).abs() < TOL);
120    }
121
122    #[test]
123    fn quantile_inverts_cdf() {
124        for &p in &[0.001, 0.025, 0.1, 0.5, 0.84, 0.975, 0.999] {
125            let z = standard_quantile(p);
126            assert!((standard_cdf(z) - p).abs() < 1e-10, "p={p} z={z}");
127        }
128        // The canonical 1.96 ≈ Φ⁻¹(0.975).
129        assert!((standard_quantile(0.975) - 1.959_963_984_540_054).abs() < 1e-7);
130    }
131
132    #[test]
133    fn general_params_shift_and_scale() {
134        // X ~ N(10, 2): cdf(10)=0.5, quantile(0.5)=10.
135        assert!((cdf(10.0, 10.0, 2.0) - 0.5).abs() < TOL);
136        assert!((quantile(0.5, 10.0, 2.0) - 10.0).abs() < 1e-9);
137        assert!((cdf(12.0, 10.0, 2.0) - standard_cdf(1.0)).abs() < TOL);
138    }
139
140    #[test]
141    fn two_sided_p_value() {
142        // |z| = 1.96 → p ≈ 0.05.
143        assert!((two_sided_p(1.959_963_984_540_054) - 0.05).abs() < 1e-9);
144        assert!((two_sided_p(0.0) - 1.0).abs() < TOL);
145    }
146}
qualia_core_db/solvers/statistics/distributions/normal.rs

qualia_core_db/solvers/statistics/distributions/
normal.rs
