aleleba/gridcraft
1
0
Fork 0
mirror of https://github.com/storytold/gridcraft.git synced 2026-11-09 03:49:54 -06:00
Code Issues Packages Projects Releases Wiki Activity

BINOM.DIST: keep tiny masses and huge n out of the product form

Browse Source 
Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
vangel76
2026-11-03 23:54:53 +01:00
co-authored by Claude Sonnet 5.5
parent 22c06b7a4a
commit 7e8b94fc39
1 changed files with 12 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/functions/src/stat/distributions.rs
+12 -6
View File
@@ -166,7 +166,9 @@ fn bd0(x: f64, np: f64) -> f64 {
fn pow_one_minus(p: f64, e: f64) -> f64 {
let hi = 1.0 - p;
let lo = (1.0 - hi) - p;
hi.powf(e) * (e * (lo / hi).ln_1p()).exp()
let r = hi.powf(e) * (e * (lo / hi).ln_1p()).exp();
// An exponent so large that the two factors under- and overflow: the plain form is enough.
if r.is_nan() { (e * (-p).ln_1p()).exp() } else { r }
}
/// `C(n, k)` for a small `k` as a product of `k` ratios.
@@ -199,11 +201,15 @@ fn binom_pmf(k: f64, n: f64, p: f64) -> f64 {
return p.powf(n);
}
// A few successes or a few failures: the product form has no large cancelling terms.
if k < 16.0 {
return choose_small(n, k) * p.powi(k as i32) * pow_one_minus(p, n - k);
}
if n - k < 16.0 {
return choose_small(n, n - k) * p.powf(k) * pow_one_minus(p, n - k);
// (Not for astronomically large n or a power of p near underflow: the saddle-point form
// below handles those through logs.)
if n < 1e30 && p.powf(k) > 1e-280 {
if k < 16.0 {
return choose_small(n, k) * p.powi(k as i32) * pow_one_minus(p, n - k);
}
if n - k < 16.0 {
return choose_small(n, n - k) * p.powf(k) * pow_one_minus(p, n - k);
}
}
let q = 1.0 - p;
let lc = stirlerr(n) - stirlerr(k) - stirlerr(n - k) - bd0(k, n * p) - bd0(n - k, n * q);