polydat_core/numeric/
special.rs

1// Copyright 2024-2026 Jonathan Shook
2// SPDX-License-Identifier: Apache-2.0
3
4//! Special functions behind the inverse-CDF samplers: the standard
5//! normal quantile and CDF, `ln Γ`, and the regularized incomplete
6//! beta and gamma functions with their inverses. The node library's
7//! distribution tables (`polydat-nodes`, `sampling::icd`) and the
8//! comprehension runtime's named-measure sampling
9//! (comprehension_forms.md §10.2 R2) compute from the same bodies,
10//! so a `normal` node and a `normal` comprehension source agree bit
11//! for bit on the same quantile.
12
13/// Rational approximation of the standard normal quantile function.
14/// Abramowitz & Stegun 26.2.23: the absolute error is below
15/// 4.5e-4 everywhere in (0, 1). `0` and `1` map to the infinities.
16pub fn probit(p: f64) -> f64 {
17    if p <= 0.0 {
18        return f64::NEG_INFINITY;
19    }
20    if p >= 1.0 {
21        return f64::INFINITY;
22    }
23
24    let t = if p < 0.5 {
25        (-2.0 * p.ln()).sqrt()
26    } else {
27        (-2.0 * (1.0 - p).ln()).sqrt()
28    };
29
30    let c0 = 2.515517;
31    let c1 = 0.802853;
32    let c2 = 0.010328;
33    let d1 = 1.432788;
34    let d2 = 0.189269;
35    let d3 = 0.001308;
36
37    let result = t - (c0 + c1 * t + c2 * t * t) / (1.0 + d1 * t + d2 * t * t + d3 * t * t * t);
38
39    if p < 0.5 { -result } else { result }
40}
41
42/// The standard normal CDF `Φ(x)`, through the regularized lower
43/// incomplete gamma function: `Φ(x) = ½(1 + sgn(x)·P(½, x²/2))`.
44pub fn normal_cdf(x: f64) -> f64 {
45    if x.is_nan() {
46        return f64::NAN;
47    }
48    if x == f64::NEG_INFINITY {
49        return 0.0;
50    }
51    if x == f64::INFINITY {
52        return 1.0;
53    }
54    let p = regularized_gamma_p(0.5, x * x / 2.0);
55    if x < 0.0 {
56        0.5 * (1.0 - p)
57    } else {
58        0.5 * (1.0 + p)
59    }
60}
61
62/// Lanczos approximation of ln(Gamma(x)) for x > 0.
63pub fn ln_gamma(x: f64) -> f64 {
64    let g = 7.0;
65    let c = [
66        0.999_999_999_999_809_9,
67        676.5203681218851,
68        -1259.1392167224028,
69        771.323_428_777_653_1,
70        -176.615_029_162_140_6,
71        12.507343278686905,
72        -0.13857109526572012,
73        9.984_369_578_019_572e-6,
74        1.5056327351493116e-7,
75    ];
76
77    if x < 0.5 {
78        let pi = std::f64::consts::PI;
79        return (pi / (pi * x).sin()).ln() - ln_gamma(1.0 - x);
80    }
81
82    let x = x - 1.0;
83    let mut sum = c[0];
84    for (i, &coeff) in c[1..].iter().enumerate() {
85        sum += coeff / (x + i as f64 + 1.0);
86    }
87
88    let t = x + g + 0.5;
89    0.5 * (2.0 * std::f64::consts::PI).ln() + (t.ln() * (x + 0.5)) - t + sum.ln()
90}
91
92/// Regularized incomplete beta function I_x(a, b) via series expansion.
93pub fn regularized_beta(x: f64, a: f64, b: f64) -> f64 {
94    if x <= 0.0 {
95        return 0.0;
96    }
97    if x >= 1.0 {
98        return 1.0;
99    }
100
101    // Use symmetry relation for better convergence when x > 0.5
102    if x > (a + 1.0) / (a + b + 2.0) {
103        return 1.0 - regularized_beta(1.0 - x, b, a);
104    }
105
106    let ln_prefix = ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln();
107    let prefix = ln_prefix.exp();
108
109    // Series expansion: I_x(a,b) = (x^a * (1-x)^b) / (a * B(a,b)) * sum
110    let mut sum = 0.0;
111    let mut term = 1.0;
112    for n in 0..300 {
113        sum += term;
114        term *= x * (a + b + n as f64) / (a + 1.0 + n as f64);
115        if term.abs() < 1e-15 * sum.abs() {
116            break;
117        }
118    }
119
120    (prefix * sum / a).clamp(0.0, 1.0)
121}
122
123/// Inverse regularized beta via bisection.
124pub fn inv_regularized_beta(p: f64, a: f64, b: f64) -> f64 {
125    if p <= 0.0 {
126        return 0.0;
127    }
128    if p >= 1.0 {
129        return 1.0;
130    }
131
132    let mut lo = 0.0_f64;
133    let mut hi = 1.0_f64;
134    for _ in 0..100 {
135        let mid = (lo + hi) / 2.0;
136        if regularized_beta(mid, a, b) < p {
137            lo = mid;
138        } else {
139            hi = mid;
140        }
141    }
142    (lo + hi) / 2.0
143}
144
145/// Regularized lower incomplete gamma function P(a, x) via series.
146pub fn regularized_gamma_p(a: f64, x: f64) -> f64 {
147    if x <= 0.0 {
148        return 0.0;
149    }
150    if x > a + 50.0 {
151        return 1.0;
152    } // far in the tail
153
154    let mut sum = 1.0 / a;
155    let mut term = 1.0 / a;
156    for n in 1..300 {
157        term *= x / (a + n as f64);
158        sum += term;
159        if term.abs() < 1e-14 * sum.abs() {
160            break;
161        }
162    }
163    (a * x.ln() - x - ln_gamma(a)).exp() * sum
164}
165
166/// Inverse regularized gamma P via bisection.
167pub fn inv_regularized_gamma_p(p: f64, a: f64) -> f64 {
168    if p <= 0.0 {
169        return 0.0;
170    }
171    if p >= 1.0 {
172        return f64::INFINITY;
173    }
174
175    // Bracket: upper bound heuristic
176    let mut hi = a.max(1.0);
177    while regularized_gamma_p(a, hi) < p {
178        hi *= 2.0;
179    }
180    let mut lo = 0.0_f64;
181
182    for _ in 0..100 {
183        let mid = (lo + hi) / 2.0;
184        if regularized_gamma_p(a, mid) < p {
185            lo = mid;
186        } else {
187            hi = mid;
188        }
189    }
190    (lo + hi) / 2.0
191}
192
193#[cfg(test)]
194mod tests {
195    use super::*;
196
197    #[test]
198    fn probit_is_odd_about_one_half_and_hits_the_infinities() {
199        assert!(probit(0.5).abs() < 1e-3);
200        assert!((probit(0.25) + probit(0.75)).abs() < 1e-3);
201        assert!((probit(0.975) - 1.96).abs() < 1e-3);
202        assert_eq!(probit(0.0), f64::NEG_INFINITY);
203        assert_eq!(probit(1.0), f64::INFINITY);
204    }
205
206    #[test]
207    fn normal_cdf_matches_the_table() {
208        assert!((normal_cdf(0.0) - 0.5).abs() < 1e-9);
209        assert!((normal_cdf(1.0) - 0.841_344_746).abs() < 1e-6);
210        assert!((normal_cdf(-1.96) - 0.024_997_895).abs() < 1e-6);
211        assert_eq!(normal_cdf(f64::NEG_INFINITY), 0.0);
212        assert_eq!(normal_cdf(f64::INFINITY), 1.0);
213    }
214
215    #[test]
216    fn normal_cdf_and_probit_are_inverses_to_the_probit_error() {
217        for p in [0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99] {
218            let back = normal_cdf(probit(p));
219            assert!((back - p).abs() < 5e-4, "p={p} back={back}");
220        }
221    }
222
223    #[test]
224    fn ln_gamma_matches_factorials() {
225        for n in 1..10u32 {
226            let fact: f64 = (1..n).map(f64::from).product();
227            assert!((ln_gamma(f64::from(n)) - fact.ln()).abs() < 1e-9, "n={n}");
228        }
229        assert!((ln_gamma(0.5) - std::f64::consts::PI.sqrt().ln()).abs() < 1e-9);
230    }
231
232    #[test]
233    fn regularized_beta_and_its_inverse_round_trip() {
234        assert!((regularized_beta(0.5, 2.0, 2.0) - 0.5).abs() < 1e-9);
235        for p in [0.05, 0.25, 0.5, 0.75, 0.95] {
236            let x = inv_regularized_beta(p, 2.0, 5.0);
237            assert!((regularized_beta(x, 2.0, 5.0) - p).abs() < 1e-9, "p={p}");
238        }
239    }
240
241    #[test]
242    fn regularized_gamma_and_its_inverse_round_trip() {
243        // P(1, x) = 1 - e^-x.
244        assert!((regularized_gamma_p(1.0, 1.0) - (1.0 - (-1.0f64).exp())).abs() < 1e-9);
245        for p in [0.05, 0.25, 0.5, 0.75, 0.95] {
246            let x = inv_regularized_gamma_p(p, 2.5);
247            assert!((regularized_gamma_p(2.5, x) - p).abs() < 1e-9, "p={p}");
248        }
249    }
250}
polydat_core/numeric/special.rs

polydat_core/numeric/
special.rs
