qualia_core_db/solvers/statistics/distributions/
special.rs

1//! Special functions underpinning the probability distributions — the canonical,
2//! full-double-precision implementations the whole engine shares.
3//!
4//! These replace the scattered ad-hoc approximations (a coarse `normal_cdf` copied
5//! into `financial_modeling`, a `|t| > 1.96 ⇒ p = 0.05` placeholder in the t-test):
6//! every distribution CDF/quantile is built from `ln_gamma`, the regularized
7//! incomplete gamma `P(a,x)`, and the regularized incomplete beta `I_x(a,b)` here.
8//!
9//! Algorithms are the standard, well-conditioned ones (Lanczos for `ln_gamma`; a
10//! series + continued-fraction split for the incomplete gamma; a Lentz continued
11//! fraction for the incomplete beta) — they are ordinary numerical mathematics,
12//! accurate to ~1e-12 and verified against known closed forms in the tests.
13//!
14//! All scalar `f64` — these are pointwise special functions, **not** GPU-amenable
15//! (per-call scalar work below any dispatch crossover); the CPU path is the right
16//! and only path here (CLAUDE.md §13: "where GPU does not help, say so").
17
18const EPS: f64 = 1e-15;
19const ITMAX: usize = 300;
20/// Smallest positive normalized-ish guard for the continued-fraction Lentz start.
21const FPMIN: f64 = 1e-300;
22
23/// Natural log of the Gamma function, `ln Γ(x)`, via the Lanczos approximation
24/// (g = 7, 9 coefficients) with the reflection formula for `x < 0.5`. Accurate to
25/// ~15 significant figures for `x > 0`.
26pub fn ln_gamma(x: f64) -> f64 {
27    // Lanczos coefficients (g = 7).
28    const G: f64 = 7.0;
29    const C: [f64; 9] = [
30        0.999_999_999_999_809_93,
31        676.520_368_121_885_1,
32        -1_259.139_216_722_402_8,
33        771.323_428_777_653_13,
34        -176.615_029_162_140_59,
35        12.507_343_278_686_905,
36        -0.138_571_095_265_720_12,
37        9.984_369_578_019_571_6e-6,
38        1.505_632_735_149_311_6e-7,
39    ];
40    if x < 0.5 {
41        // Reflection: Γ(x)Γ(1-x) = π / sin(πx).
42        let pi = std::f64::consts::PI;
43        return (pi / (pi * x).sin()).ln() - ln_gamma(1.0 - x);
44    }
45    let x = x - 1.0;
46    let mut a = C[0];
47    let t = x + G + 0.5;
48    for (i, &ci) in C.iter().enumerate().skip(1) {
49        a += ci / (x + i as f64);
50    }
51    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + a.ln()
52}
53
54/// Γ(x) for convenience (small/moderate `x`).
55pub fn gamma(x: f64) -> f64 {
56    ln_gamma(x).exp()
57        * if x < 0.5 && (x.floor() == x) {
58            f64::NAN
59        } else {
60            1.0
61        }
62}
63
64/// Regularized lower incomplete gamma `P(a, x) = γ(a,x)/Γ(a)`, `a > 0`, `x ≥ 0`.
65/// Series for `x < a+1`, continued fraction (via `Q`) otherwise.
66pub fn gammp(a: f64, x: f64) -> f64 {
67    debug_assert!(a > 0.0);
68    if x <= 0.0 {
69        return 0.0;
70    }
71    if x < a + 1.0 {
72        gser(a, x)
73    } else {
74        1.0 - gcf(a, x)
75    }
76}
77
78/// Regularized upper incomplete gamma `Q(a, x) = 1 − P(a, x)`.
79pub fn gammq(a: f64, x: f64) -> f64 {
80    1.0 - gammp(a, x)
81}
82
83/// Series evaluation of `P(a, x)` (good for `x < a+1`).
84fn gser(a: f64, x: f64) -> f64 {
85    let mut ap = a;
86    let mut sum = 1.0 / a;
87    let mut del = sum;
88    for _ in 0..ITMAX {
89        ap += 1.0;
90        del *= x / ap;
91        sum += del;
92        if del.abs() < sum.abs() * EPS {
93            break;
94        }
95    }
96    sum * (-x + a * x.ln() - ln_gamma(a)).exp()
97}
98
99/// Continued-fraction evaluation of `Q(a, x)` (good for `x ≥ a+1`), Lentz's method.
100fn gcf(a: f64, x: f64) -> f64 {
101    let mut b = x + 1.0 - a;
102    let mut c = 1.0 / FPMIN;
103    let mut d = 1.0 / b;
104    let mut h = d;
105    for i in 1..=ITMAX {
106        let an = -(i as f64) * (i as f64 - a);
107        b += 2.0;
108        d = an * d + b;
109        if d.abs() < FPMIN {
110            d = FPMIN;
111        }
112        c = b + an / c;
113        if c.abs() < FPMIN {
114            c = FPMIN;
115        }
116        d = 1.0 / d;
117        let del = d * c;
118        h *= del;
119        if (del - 1.0).abs() < EPS {
120            break;
121        }
122    }
123    (-x + a * x.ln() - ln_gamma(a)).exp() * h
124}
125
126/// Regularized incomplete beta `I_x(a, b)`, `0 ≤ x ≤ 1`, `a,b > 0`. Continued
127/// fraction (Lentz) with the symmetry switch for fast convergence.
128pub fn betai(a: f64, b: f64, x: f64) -> f64 {
129    debug_assert!(a > 0.0 && b > 0.0);
130    if x <= 0.0 {
131        return 0.0;
132    }
133    if x >= 1.0 {
134        return 1.0;
135    }
136    let bt = (ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln()).exp();
137    if x < (a + 1.0) / (a + b + 2.0) {
138        bt * betacf(a, b, x) / a
139    } else {
140        1.0 - bt * betacf(b, a, 1.0 - x) / b
141    }
142}
143
144/// Lentz continued fraction for the incomplete beta.
145fn betacf(a: f64, b: f64, x: f64) -> f64 {
146    let qab = a + b;
147    let qap = a + 1.0;
148    let qam = a - 1.0;
149    let mut c = 1.0;
150    let mut d = 1.0 - qab * x / qap;
151    if d.abs() < FPMIN {
152        d = FPMIN;
153    }
154    d = 1.0 / d;
155    let mut h = d;
156    for m in 1..=ITMAX {
157        let m = m as f64;
158        let m2 = 2.0 * m;
159        // even step
160        let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
161        d = 1.0 + aa * d;
162        if d.abs() < FPMIN {
163            d = FPMIN;
164        }
165        c = 1.0 + aa / c;
166        if c.abs() < FPMIN {
167            c = FPMIN;
168        }
169        d = 1.0 / d;
170        h *= d * c;
171        // odd step
172        let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
173        d = 1.0 + aa * d;
174        if d.abs() < FPMIN {
175            d = FPMIN;
176        }
177        c = 1.0 + aa / c;
178        if c.abs() < FPMIN {
179            c = FPMIN;
180        }
181        d = 1.0 / d;
182        let del = d * c;
183        h *= del;
184        if (del - 1.0).abs() < EPS {
185            break;
186        }
187    }
188    h
189}
190
191/// Error function `erf(x)` via the incomplete gamma: `erf(x) = sign(x)·P(½, x²)`.
192pub fn erf(x: f64) -> f64 {
193    if x == 0.0 {
194        0.0
195    } else if x > 0.0 {
196        gammp(0.5, x * x)
197    } else {
198        -gammp(0.5, x * x)
199    }
200}
201
202/// Complementary error function `erfc(x) = 1 − erf(x)`.
203pub fn erfc(x: f64) -> f64 {
204    1.0 - erf(x)
205}
206
207#[cfg(test)]
208mod tests {
209    use super::*;
210    const TOL: f64 = 1e-9;
211
212    #[test]
213    fn ln_gamma_known_values() {
214        assert!((ln_gamma(1.0)).abs() < TOL); // Γ(1)=1
215        assert!((ln_gamma(2.0)).abs() < TOL); // Γ(2)=1
216        assert!((ln_gamma(5.0) - 24.0_f64.ln()).abs() < 1e-9); // Γ(5)=4!=24
217                                                               // Γ(1/2)=√π
218        assert!((ln_gamma(0.5) - std::f64::consts::PI.sqrt().ln()).abs() < 1e-9);
219    }
220
221    #[test]
222    fn erf_known_values() {
223        assert!((erf(0.0)).abs() < TOL);
224        assert!((erf(1.0) - 0.842_700_792_949_715).abs() < 1e-9);
225        assert!((erf(-1.0) + 0.842_700_792_949_715).abs() < 1e-9);
226        assert!((erf(2.0) - 0.995_322_265_018_953).abs() < 1e-9);
227        assert!((erfc(0.0) - 1.0).abs() < TOL);
228    }
229
230    #[test]
231    fn gammp_is_exponential_cdf_for_a_one() {
232        // P(1, x) = 1 - e^{-x} (the Exp(1) CDF).
233        for &x in &[0.5, 1.0, 2.5, 5.0] {
234            assert!((gammp(1.0, x) - (1.0 - (-x).exp())).abs() < 1e-10, "x={x}");
235        }
236        assert!((gammp(1.0, 0.0)).abs() < TOL);
237    }
238
239    #[test]
240    fn gammp_gammq_complementary() {
241        for &(a, x) in &[(0.5, 0.3), (2.0, 1.0), (3.5, 7.0), (10.0, 4.0)] {
242            assert!((gammp(a, x) + gammq(a, x) - 1.0).abs() < 1e-12);
243        }
244    }
245
246    #[test]
247    fn betai_endpoints_and_symmetry() {
248        // I_x(1,1) = x (uniform).
249        for &x in &[0.0, 0.25, 0.5, 0.9, 1.0] {
250            assert!((betai(1.0, 1.0, x) - x).abs() < 1e-10, "x={x}");
251        }
252        // Symmetry: I_x(a,b) = 1 - I_{1-x}(b,a).
253        for &(a, b, x) in &[(2.0, 3.0, 0.4), (0.5, 2.5, 0.7), (5.0, 1.5, 0.2)] {
254            assert!(
255                (betai(a, b, x) - (1.0 - betai(b, a, 1.0 - x))).abs() < 1e-10,
256                "{a},{b},{x}"
257            );
258        }
259    }
260}
qualia_core_db/solvers/statistics/distributions/special.rs

qualia_core_db/solvers/statistics/distributions/
special.rs
