qualia_core_db/solvers/statistics/distributions/
mod.rs

1//! Probability distributions — the canonical, full-precision pdf / cdf / quantile
2//! for the whole engine, built on the shared special functions ([`special`]).
3//!
4//! This is what makes p-values and confidence intervals **honest**: hypothesis
5//! tests ([`super::hypothesis`]) get their tail probabilities from a real
6//! Student-t / χ² / F CDF here, not a `|t| > 1.96 ⇒ p = 0.05` placeholder, and
7//! domain libraries reuse [`normal`] instead of copying a local `normal_cdf`.
8//!
9//! Submodules (one distribution each, PROJECT RULE §11): [`special`] (erf / incomplete
10//! gamma & beta), [`normal`], [`students_t`], [`chi_squared`], [`fisher_f`].
11//!
12//! Everything here is scalar `f64` special-function evaluation — pointwise, not
13//! GPU-amenable (CLAUDE.md §13: the CPU path is the right one; the *data-aggregate*
14//! kernels that feed these, e.g. mean/variance, are the `Reduction`-class work that
15//! routes through `ComputePolicy`).
16
17pub mod chi_squared;
18pub mod fisher_f;
19pub mod multivariate_normal;
20pub mod normal;
21pub mod special;
22pub mod students_t;
23
24// Lightweight additional distributions for computational economics (5.1-A).
25// These are scalar, no allocation. Full families can grow into dedicated files.
26
27/// Binomial PMF: P(K = k | n, p).
28#[inline]
29pub fn binomial_pmf(k: u32, n: u32, p: f64) -> f64 {
30    if p < 0.0 || p > 1.0 || k > n {
31        return f64::NAN;
32    }
33    if n == 0 {
34        return if k == 0 { 1.0 } else { 0.0 };
35    }
36    let ln_c = ln_binom(n, k);
37    let ln_p = (k as f64) * p.ln() + ((n - k) as f64) * (1.0 - p).ln();
38    (ln_c + ln_p).exp()
39}
40
41/// Binomial CDF via direct sum (small n only; for large use normal approx in caller).
42pub fn binomial_cdf(k: u32, n: u32, p: f64) -> f64 {
43    if p < 0.0 || p > 1.0 {
44        return f64::NAN;
45    }
46    let kk = k.min(n);
47    let mut s = 0.0;
48    for i in 0..=kk {
49        s += binomial_pmf(i, n, p);
50        if !s.is_finite() {
51            return f64::NAN;
52        }
53    }
54    s
55}
56
57fn ln_binom(n: u32, k: u32) -> f64 {
58    // ln(n! / (k!(n-k)! )) using sum of logs
59    if k > n {
60        return f64::NEG_INFINITY;
61    }
62    let mut s = 0.0;
63    for i in 0..k {
64        s += ((n - i) as f64).ln() - ((i + 1) as f64).ln();
65    }
66    s
67}
68
69/// Poisson PMF.
70#[inline]
71pub fn poisson_pmf(k: u32, lambda: f64) -> f64 {
72    if lambda < 0.0 {
73        return f64::NAN;
74    }
75    if k == 0 {
76        return (-lambda).exp();
77    }
78    let mut pmf = (-lambda).exp();
79    for i in 1..=k {
80        pmf *= lambda / (i as f64);
81    }
82    pmf
83}
84
85/// Poisson CDF.
86pub fn poisson_cdf(k: u32, lambda: f64) -> f64 {
87    if lambda < 0.0 {
88        return f64::NAN;
89    }
90    let mut s = 0.0;
91    for i in 0..=k {
92        s += poisson_pmf(i, lambda);
93    }
94    s
95}
96
97/// Lognormal PDF (mu, sigma>0).
98#[inline]
99pub fn lognormal_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
100    if x <= 0.0 || sigma <= 0.0 || !x.is_finite() || !mu.is_finite() || !sigma.is_finite() {
101        return f64::NAN;
102    }
103    let z = (x.ln() - mu) / sigma;
104    (1.0 / (x * sigma * core::f64::consts::TAU.sqrt())) * (-0.5 * z * z).exp()
105}
106
107/// Lognormal CDF via normal cdf of ln(x).
108pub fn lognormal_cdf(x: f64, mu: f64, sigma: f64) -> f64 {
109    if x <= 0.0 || sigma <= 0.0 {
110        return if x <= 0.0 { 0.0 } else { f64::NAN };
111    }
112    normal::cdf((x.ln() - mu) / sigma, 0.0, 1.0)
113}
114
115/// Exponential PDF (rate > 0).
116#[inline]
117pub fn exponential_pdf(x: f64, rate: f64) -> f64 {
118    if x < 0.0 || rate <= 0.0 || !x.is_finite() || !rate.is_finite() {
119        return f64::NAN;
120    }
121    rate * (-rate * x).exp()
122}
123
124/// Exponential CDF.
125#[inline]
126pub fn exponential_cdf(x: f64, rate: f64) -> f64 {
127    if x < 0.0 || rate <= 0.0 {
128        return if x < 0.0 { 0.0 } else { f64::NAN };
129    }
130    1.0 - (-rate * x).exp()
131}
132
133/// Uniform PDF on [a, b].
134#[inline]
135pub fn uniform_pdf(x: f64, a: f64, b: f64) -> f64 {
136    if a >= b || !x.is_finite() || !a.is_finite() || !b.is_finite() {
137        return f64::NAN;
138    }
139    if x >= a && x <= b {
140        1.0 / (b - a)
141    } else {
142        0.0
143    }
144}
145
146/// Uniform CDF.
147#[inline]
148pub fn uniform_cdf(x: f64, a: f64, b: f64) -> f64 {
149    if a >= b {
150        return f64::NAN;
151    }
152    if x < a {
153        0.0
154    } else if x > b {
155        1.0
156    } else {
157        (x - a) / (b - a)
158    }
159}
160
161/// Laplace (double exponential) PDF.
162#[inline]
163pub fn laplace_pdf(x: f64, mu: f64, b: f64) -> f64 {
164    if b <= 0.0 || !x.is_finite() || !mu.is_finite() || !b.is_finite() {
165        return f64::NAN;
166    }
167    (1.0 / (2.0 * b)) * (-((x - mu).abs() / b)).exp()
168}
169
170/// Laplace CDF.
171pub fn laplace_cdf(x: f64, mu: f64, b: f64) -> f64 {
172    if b <= 0.0 {
173        return f64::NAN;
174    }
175    let z = (x - mu) / b;
176    if z < 0.0 {
177        0.5 * (z).exp()
178    } else {
179        1.0 - 0.5 * (-z).exp()
180    }
181}
182
183/// Gamma PDF (shape k>0, scale theta>0).
184pub fn gamma_pdf(x: f64, shape: f64, scale: f64) -> f64 {
185    if x <= 0.0 || shape <= 0.0 || scale <= 0.0 {
186        return if x <= 0.0 { 0.0 } else { f64::NAN };
187    }
188    let log_pdf =
189        (shape - 1.0) * x.ln() - x / scale - shape * scale.ln() - special::ln_gamma(shape);
190    log_pdf.exp()
191}
192
193/// Basic Beta PDF (alpha, beta >0) on (0,1).
194pub fn beta_pdf(x: f64, alpha: f64, beta: f64) -> f64 {
195    if x <= 0.0 || x >= 1.0 || alpha <= 0.0 || beta <= 0.0 {
196        return 0.0;
197    }
198    let log_b =
199        special::ln_gamma(alpha) + special::ln_gamma(beta) - special::ln_gamma(alpha + beta);
200    ((alpha - 1.0) * x.ln() + (beta - 1.0) * (1.0 - x).ln() - log_b).exp()
201}
202
203/// Weibull PDF (shape k>0, scale lambda>0).
204#[inline]
205pub fn weibull_pdf(x: f64, shape: f64, scale: f64) -> f64 {
206    if x < 0.0 || shape <= 0.0 || scale <= 0.0 {
207        return 0.0;
208    }
209    (shape / scale) * (x / scale).powf(shape - 1.0) * (-(x / scale).powf(shape)).exp()
210}
211
212/// Empirical CDF from sorted samples (for caller-sorted data).
213pub fn empirical_cdf(sorted_samples: &[f64], x: f64) -> f64 {
214    if sorted_samples.is_empty() {
215        return f64::NAN;
216    }
217    let mut count = 0usize;
218    for &s in sorted_samples {
219        if s <= x {
220            count += 1;
221        }
222    }
223    count as f64 / sorted_samples.len() as f64
224}
225
226/// Invert a monotone-increasing CDF: find `x` with `cdf(x) ≈ p` by adaptive
227/// bracketing + bisection. `lower` bounds the support below (`Some(0.0)` for χ²/F,
228/// `None` for a doubly-infinite support like Student-t). Used by the distribution
229/// quantiles that have no closed-form inverse. Robust and ~machine-precision.
230pub(crate) fn invert_cdf(p: f64, lower: Option<f64>, cdf: impl Fn(f64) -> f64) -> f64 {
231    if p <= 0.0 {
232        return lower.unwrap_or(f64::NEG_INFINITY);
233    }
234    if p >= 1.0 {
235        return f64::INFINITY;
236    }
237    // Establish a bracket [lo, hi] with cdf(lo) ≤ p ≤ cdf(hi).
238    let mut lo = lower.unwrap_or(-1.0);
239    let mut hi = lower.map(|l| l + 1.0).unwrap_or(1.0);
240    if lower.is_none() {
241        let mut guard = 0;
242        while cdf(lo) > p && guard < 80 {
243            lo *= 2.0;
244            guard += 1;
245        }
246    }
247    let mut guard = 0;
248    while cdf(hi) < p && guard < 80 {
249        hi *= 2.0;
250        guard += 1;
251    }
252    // Bisection.
253    for _ in 0..200 {
254        let mid = 0.5 * (lo + hi);
255        if cdf(mid) < p {
256            lo = mid;
257        } else {
258            hi = mid;
259        }
260    }
261    0.5 * (lo + hi)
262}
263
264#[cfg(test)]
265mod tests {
266    use super::*;
267
268    #[test]
269    fn invert_cdf_recovers_a_linear_cdf() {
270        // CDF(x) = x on [0,1] → quantile(p) = p.
271        let q = invert_cdf(0.37, Some(0.0), |x| x.clamp(0.0, 1.0));
272        assert!((q - 0.37).abs() < 1e-9);
273    }
274}
qualia_core_db/solvers/statistics/distributions/mod.rs

qualia_core_db/solvers/statistics/distributions/
mod.rs
