From 22c06b7a4a394334e0d35f8880da2e9ce7d830e6 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Tue, 3 Nov 2026 23:49:52 +0100 Subject: [PATCH] BINOM.DIST, POISSON.DIST: Loader mass, recurrence sums, exact prefactors (#3324) The mass came from ln_gamma differences (relative error 1.1e-13 at n = 100, 1.5e-11 at n = 10,000) and the cumulative forms summed up to 10,000 of them. The mass now comes from Loader's saddle-point formula, the cumulative sum walks away from the mode by the ratio of neighbouring terms and stops when terms stop mattering, and the incomplete beta / gamma branches take the same mass as their prefactor. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/special.rs | 36 +++- crates/functions/src/stat.rs | 5 +- crates/functions/src/stat/distributions.rs | 191 +++++++++++++++-- crates/functions/src/stat/tests.rs | 230 +++++++++++++++++++++ 4 files changed, 436 insertions(+), 26 deletions(-) diff --git a/crates/functions/src/special.rs b/crates/functions/src/special.rs index 1541cd8a1..b173ff765 100644 --- a/crates/functions/src/special.rs +++ b/crates/functions/src/special.rs @@ -88,11 +88,17 @@ pub(crate) fn reg_gamma_p(a: f64, x: f64) -> f64 { if x.is_infinite() { return 1.0; } - if x < a + 1.0 { gamma_series(a, x) } else { 1.0 - gamma_cf(a, x) } + if x < a + 1.0 { gamma_series(a, x, front_of(a, x)) } else { 1.0 - gamma_cf(a, x, front_of(a, x)) } } /// Regularised upper incomplete gamma Q(a, x) = 1 − P(a, x). pub(crate) fn reg_gamma_q(a: f64, x: f64) -> f64 { + let front = if a > 0.0 && x > 0.0 && x.is_finite() { (-x + a * x.ln() - ln_gamma(a)).exp() } else { 0.0 }; + reg_gamma_q_front(a, x, front) +} + +/// Q(a, x) given the prefactor `e^-x x^a / Γ(a)` (the Poisson mass at a - 1, times x). +pub(crate) fn reg_gamma_q_front(a: f64, x: f64, front: f64) -> f64 { if a.is_nan() || a <= 0.0 || x.is_nan() || x < 0.0 || !a.is_finite() { return f64::NAN; } @@ -102,10 +108,15 @@ pub(crate) fn reg_gamma_q(a: f64, x: f64) -> f64 { if x.is_infinite() { return 0.0; } - if x < a + 1.0 { 1.0 - gamma_series(a, x) } else { gamma_cf(a, x) } + if x < a + 1.0 { 1.0 - gamma_series(a, x, front) } else { gamma_cf(a, x, front) } } -fn gamma_series(a: f64, x: f64) -> f64 { +/// `e^-x x^a / Γ(a)`. +fn front_of(a: f64, x: f64) -> f64 { + (-x + a * x.ln() - ln_gamma(a)).exp() +} + +fn gamma_series(a: f64, x: f64, front: f64) -> f64 { let mut ap = a; let mut sum = 1.0 / a; let mut del = sum; @@ -117,10 +128,10 @@ fn gamma_series(a: f64, x: f64) -> f64 { break; } } - sum * (-x + a * x.ln() - ln_gamma(a)).exp() + sum * front } -fn gamma_cf(a: f64, x: f64) -> f64 { +fn gamma_cf(a: f64, x: f64, front: f64) -> f64 { let tiny = 1e-300; let mut b = x + 1.0 - a; let mut c = 1.0 / tiny; @@ -144,7 +155,7 @@ fn gamma_cf(a: f64, x: f64) -> f64 { break; } } - (-x + a * x.ln() - ln_gamma(a)).exp() * h + front * h } /// Error function. @@ -187,6 +198,17 @@ pub(crate) fn erfc(x: f64) -> f64 { /// Regularised incomplete beta I_x(a, b). pub(crate) fn reg_inc_beta(x: f64, a: f64, b: f64) -> f64 { + if x > 0.0 && x < 1.0 && a > 0.0 && b > 0.0 { + let ln_bt = -ln_beta(a, b) + a * x.ln() + b * (1.0 - x).ln(); + return reg_inc_beta_front(x, a, b, ln_bt.exp()); + } + reg_inc_beta_front(x, a, b, 0.0) +} + +/// I_x(a, b) given the prefactor `bt = x^a (1-x)^b / B(a, b)` (symmetric in the two sides). A +/// caller that knows `bt` more accurately than the log-gamma route gives (the binomial mass, for +/// instance) passes it here. +pub(crate) fn reg_inc_beta_front(x: f64, a: f64, b: f64, bt: f64) -> f64 { if a.is_nan() || a <= 0.0 || b.is_nan() || b <= 0.0 || !(0.0..=1.0).contains(&x) { return f64::NAN; } @@ -196,8 +218,6 @@ pub(crate) fn reg_inc_beta(x: f64, a: f64, b: f64) -> f64 { if x == 1.0 { return 1.0; } - let ln_bt = -ln_beta(a, b) + a * x.ln() + b * (1.0 - x).ln(); - let bt = ln_bt.exp(); if x < (a + 1.0) / (a + b + 2.0) { bt * beta_cf(a, b, x) / a } else { 1.0 - bt * beta_cf(b, a, 1.0 - x) / b } } diff --git a/crates/functions/src/stat.rs b/crates/functions/src/stat.rs index 270bd2373..790915c95 100644 --- a/crates/functions/src/stat.rs +++ b/crates/functions/src/stat.rs @@ -3,7 +3,10 @@ use gridcraft_core::{CellError, Value}; -use crate::special::{erf, gamma, invert, ln_beta, ln_gamma, norm_cdf, norm_inv, norm_pdf, reg_gamma_p, reg_gamma_q, reg_inc_beta}; +use crate::special::{ + erf, gamma, invert, ln_beta, ln_gamma, norm_cdf, norm_inv, norm_pdf, reg_gamma_p, reg_gamma_q, reg_gamma_q_front, reg_inc_beta, + reg_inc_beta_front, +}; use crate::util::{A, R, S, arg, array_numbers, array_val, as_array, boolean, has, num, num_val, numbers, numbers_a, opt_num}; use crate::{Arg, Ctx, FnSpec, VAR}; diff --git a/crates/functions/src/stat/distributions.rs b/crates/functions/src/stat/distributions.rs index a3964ad80..7fa486e93 100644 --- a/crates/functions/src/stat/distributions.rs +++ b/crates/functions/src/stat/distributions.rs @@ -1,7 +1,7 @@ use super::descriptive::{arg_nums, arg_pairs, mean, var_s}; use super::{ Arg, C, CellError, R, Value, arg, as_array, boolean, erf, gamma, has, invert, ln_beta, ln_gamma, norm_cdf, norm_inv, norm_pdf, num, num_val, - reg_gamma_p, reg_gamma_q, reg_inc_beta, + reg_gamma_p, reg_gamma_q, reg_gamma_q_front, reg_inc_beta, reg_inc_beta_front, }; // --------------------------------------------------------------------------------------------- @@ -93,6 +93,95 @@ fn f_pdf(x: f64, d1: f64, d2: f64) -> f64 { const MAX_BINOM_TRIALS: f64 = 2_147_483_646.0; +/// Stirling-series error `ln n! - ((n + 1/2) ln n - n + ln(2π)/2)` for integer `n >= 1`: exact +/// table (computed with 60-digit Decimal arithmetic) up to 30, then the series, whose first +/// omitted term is below 1e-19 there. +fn stirlerr(n: f64) -> f64 { + const TABLE: [f64; 31] = [ + 0.0, + 0.08106146679532726, + 0.0413406959554093, + 0.02767792568499834, + 0.020790672103765093, + 0.016644691189821193, + 0.013876128823070748, + 0.01189670994589177, + 0.010411265261972096, + 0.009255462182712733, + 0.00833056343336287, + 0.007573675487951841, + 0.00694284010720953, + 0.006408994188004207, + 0.0059513701127588475, + 0.005554733551962801, + 0.0052076559196096404, + 0.004901395948434738, + 0.004629153749334028, + 0.004385560249232324, + 0.004166319691996922, + 0.00396795421864086, + 0.0037876180684444346, + 0.0036229602246830948, + 0.003472021382978767, + 0.003333155636728093, + 0.003204970228055038, + 0.0030862786826087773, + 0.002976063983550409, + 0.0028734493623524663, + 0.0027776749297526936, + ]; + if n <= 30.0 { + return TABLE.get(n as usize).copied().unwrap_or(0.0); + } + let r = 1.0 / n; + let r2 = r * r; + r * (1.0 / 12.0 - r2 * (1.0 / 360.0 - r2 * (1.0 / 1260.0 - r2 * (1.0 / 1680.0 - r2 * (1.0 / 1188.0))))) +} + +/// `x ln(x/np) + np - x` without the cancellation of the direct formula when `x` is near `np` +/// (Loader, 2000). +fn bd0(x: f64, np: f64) -> f64 { + if (x - np).abs() < 0.1 * (x + np) { + let v = (x - np) / (x + np); + let mut s = (x - np) * v; + if s.abs() < f64::MIN_POSITIVE { + return s; + } + let mut ej = 2.0 * x * v; + let v2 = v * v; + for j in 1..1000 { + ej *= v2; + let s1 = s + ej / f64::from(2 * j + 1); + if s1 == s { + return s1; + } + s = s1; + } + } + x * (x / np).ln() + np - x +} + +/// `(1 - p)^e` with the rounding error of `1 - p` (it is not exact for `p < 1/2`) folded back in, +/// which `e` would otherwise multiply. +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() +} + +/// `C(n, k)` for a small `k` as a product of `k` ratios. +fn choose_small(n: f64, k: f64) -> f64 { + let mut c = 1.0; + let mut j = 0.0; + while j < k { + c *= (n - j) / (j + 1.0); + j += 1.0; + } + c +} + +/// Binomial mass at integer `k`, `0 <= k <= n`, by Loader's saddle-point formula: no log-gamma +/// of large numbers, so the relative error stays near 1e-15 for any `n`. fn binom_pmf(k: f64, n: f64, p: f64) -> f64 { if k < 0.0 || k > n { return 0.0; @@ -103,7 +192,60 @@ fn binom_pmf(k: f64, n: f64, p: f64) -> f64 { if p == 1.0 { return if k == n { 1.0 } else { 0.0 }; } - (ln_comb(n, k) + k * p.ln() + (n - k) * (1.0 - p).ln()).exp() + if k == 0.0 { + return if n == 0.0 { 1.0 } else { pow_one_minus(p, n) }; + } + if k == n { + 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); + } + let q = 1.0 - p; + let lc = stirlerr(n) - stirlerr(k) - stirlerr(n - k) - bd0(k, n * p) - bd0(n - k, n * q); + let lf = std::f64::consts::TAU.ln() + k.ln() + (-k / n).ln_1p(); + (lc - 0.5 * lf).exp() +} + +/// Sum of the binomial masses 0..=k for `n <= 10_000`: starts at the term nearest the mode +/// (accurate by the formula above) and walks away from it by the ratio of neighbouring masses, +/// one multiply per term, stopping once the terms no longer matter. Above the mode it sums the +/// upper tail instead and subtracts it from 1. +fn binom_cdf_sum(k: f64, n: f64, p: f64) -> f64 { + let q = 1.0 - p; + if p == 0.0 { + return 1.0; + } + if p == 1.0 { + return 0.0; + } + if k <= n * p { + // mass(i - 1) = mass(i) * i / (n - i + 1) * q / p + let mut t = binom_pmf(k, n, p); + let mut s = t; + let mut i = k; + while i >= 1.0 && t > s * 1e-18 { + t *= i / (n - i + 1.0) * (q / p); + s += t; + i -= 1.0; + } + s.min(1.0) + } else { + // mass(i + 1) = mass(i) * (n - i) / (i + 1) * p / q + let mut t = binom_pmf(k + 1.0, n, p); + let mut s = t; + let mut i = k + 1.0; + while i < n && t > s * 1e-18 { + t *= (n - i) / (i + 1.0) * (p / q); + s += t; + i += 1.0; + } + (1.0 - s).max(0.0) + } } fn binom_cdf(k: f64, n: f64, p: f64) -> f64 { @@ -114,13 +256,7 @@ fn binom_cdf(k: f64, n: f64, p: f64) -> f64 { return 1.0; } if n <= 10_000.0 { - let mut s = 0.0; - let mut i = 0.0; - while i <= k { - s += binom_pmf(i, n, p); - i += 1.0; - } - return s.min(1.0); + return binom_cdf_sum(k, n, p); } if p == 0.0 { return 1.0; @@ -128,14 +264,19 @@ fn binom_cdf(k: f64, n: f64, p: f64) -> f64 { if p == 1.0 { return 0.0; } - reg_inc_beta(1.0 - p, n - k, k + 1.0) + // I_(1-p)(n-k, k+1) has prefactor (n-k) p mass(k). + reg_inc_beta_front(1.0 - p, n - k, k + 1.0, (n - k) * p * binom_pmf(k, n, p)) } +/// Poisson mass at integer `x >= 0` by Loader's formula. fn poisson_pmf(x: f64, m: f64) -> f64 { if m == 0.0 { return if x == 0.0 { 1.0 } else { 0.0 }; } - (-m + x * m.ln() - ln_gamma(x + 1.0)).exp() + if x == 0.0 { + return (-m).exp(); + } + (-stirlerr(x) - bd0(x, m)).exp() / (std::f64::consts::TAU * x).sqrt() } fn poisson_cdf(x: f64, m: f64) -> f64 { @@ -143,15 +284,31 @@ fn poisson_cdf(x: f64, m: f64) -> f64 { return 1.0; } if x <= 1000.0 { - let mut s = 0.0; - let mut i = 0.0; - while i <= x { - s += poisson_pmf(i, m); + if x <= m { + // mass(i - 1) = mass(i) * i / m, downwards from x while the terms matter + let mut t = poisson_pmf(x, m); + let mut s = t; + let mut i = x; + while i >= 1.0 && t > s * 1e-18 { + t *= i / m; + s += t; + i -= 1.0; + } + return s.min(1.0); + } + // Above the mean: 1 minus the upper tail, mass(i + 1) = mass(i) * m / (i + 1). + let mut t = poisson_pmf(x + 1.0, m); + let mut s = t; + let mut i = x + 1.0; + while t > s * 1e-18 { + t *= m / (i + 1.0); + s += t; i += 1.0; } - return s.min(1.0); + return (1.0 - s).max(0.0); } - reg_gamma_q(x + 1.0, m) + // Q(x+1, m) has prefactor m mass(x). + reg_gamma_q_front(x + 1.0, m, m * poisson_pmf(x, m)) } pub(super) fn norm_dist(a: &[Arg], _: C) -> R { diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index ba59116cd..f8a8c1890 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -617,3 +617,233 @@ fn aggregate_helper() { close(aggregate_values_k(9, &vals, &nv(2.0)), 9.0); is_err(aggregate_values_k(14, &vals, &nv(10.0)), CellError::Num); } + +// B1: BINOM.DIST and POISSON.DIST went through ln_gamma differences (relative error 1.1e-13 at +// n = 100, 1.5e-11 at n = 10,000). Reference values: 70-digit Decimal arithmetic on the exact +// double `p` (mass by the neighbour-ratio recurrence from q^n, cdf as the running sum), rounded +// once to f64. Columns: k, n, p, mass, cdf. +#[rustfmt::skip] +const BINOM_REF: &[(f64, f64, f64, f64, f64)] = &[ + (0.0, 10.0, 0.5, 0.0009765625, 0.0009765625), + (1.0, 10.0, 0.5, 0.009765625, 0.0107421875), + (3.0, 10.0, 0.5, 0.1171875, 0.171875), + (5.0, 10.0, 0.5, 0.24609375, 0.623046875), + (7.0, 10.0, 0.5, 0.1171875, 0.9453125), + (9.0, 10.0, 0.5, 0.009765625, 0.9990234375), + (0.0, 10.0, 0.3, 0.028247524900000005, 0.028247524900000005), + (1.0, 10.0, 0.3, 0.12106082100000001, 0.14930834590000003), + (2.0, 10.0, 0.3, 0.23347444050000002, 0.3827827864), + (3.0, 10.0, 0.3, 0.266827932, 0.6496107184000001), + (4.0, 10.0, 0.3, 0.200120949, 0.8497316674000001), + (5.0, 10.0, 0.3, 0.10291934519999998, 0.9526510126), + (7.0, 10.0, 0.3, 0.009001691999999999, 0.9984096136), + (9.0, 10.0, 0.3, 0.00013778099999999996, 0.9999940951), + (0.0, 10.0, 0.999, 1.0000000000000088e-30, 1.0000000000000088e-30), + (1.0, 10.0, 0.999, 9.99000000000008e-27, 9.99100000000008e-27), + (5.0, 10.0, 0.999, 2.507425174812609e-13, 2.50951798425701e-13), + (9.0, 10.0, 0.999, 0.00991035916125875, 0.009955119790251798), + (0.0, 100.0, 0.5, 7.888609052210118e-31, 7.888609052210118e-31), + (1.0, 100.0, 0.5, 7.888609052210118e-29, 7.967495142732219e-29), + (35.0, 100.0, 0.5, 0.0008638556657416528, 0.0017588208614850792), + (45.0, 100.0, 0.5, 0.048474296626430755, 0.18410080866334813), + (50.0, 100.0, 0.5, 0.07958923738717877, 0.5397946186935894), + (55.0, 100.0, 0.5, 0.048474296626430755, 0.8643734879630827), + (65.0, 100.0, 0.5, 0.0008638556657416528, 0.9991050348042566), + (99.0, 100.0, 0.5, 7.888609052210118e-29, 1.0), + (0.0, 100.0, 0.01, 0.3660323412732295, 0.3660323412732295), + (1.0, 100.0, 0.01, 0.3697296376497268, 0.7357619789229563), + (2.0, 100.0, 0.01, 0.1848648188248634, 0.9206267977478196), + (4.0, 100.0, 0.01, 0.014941714856895142, 0.9965676784122455), + (50.0, 100.0, 0.01, 6.103987557173001e-72, 1.0), + (99.0, 100.0, 0.01, 9.900000000000021e-197, 1.0), + (0.0, 100.0, 0.7, 5.15377520732019e-53, 5.15377520732019e-53), + (1.0, 100.0, 0.7, 1.2025475483747105e-50, 1.207701323582031e-50), + (50.0, 100.0, 0.7, 1.302622713144541e-05, 2.2060913327166086e-05), + (56.0, 100.0, 0.7, 0.001028871216633018, 0.0021146172813184595), + (65.0, 100.0, 0.7, 0.04677968235272982, 0.16285828837178737), + (70.0, 100.0, 0.7, 0.08678386475342809, 0.5376602639846381), + (75.0, 100.0, 0.7, 0.04955992276217012, 0.8864298182958186), + (84.0, 100.0, 0.7, 0.0005638655029474704, 0.9995950004580563), + (99.0, 100.0, 0.7, 1.3862042184106021e-14, 0.9999999999999997), + (0.0, 1000.0, 0.5, 9.332636185032189e-302, 9.332636185032189e-302), + (1.0, 1000.0, 0.5, 9.332636185032189e-299, 9.341968821217221e-299), + (453.0, 1000.0, 0.5, 0.00030353498700958767, 0.00162629990658169), + (484.0, 1000.0, 0.5, 0.015123665829158384, 0.1634683412157977), + (500.0, 1000.0, 0.5, 0.0252250181783608, 0.5126125090891804), + (516.0, 1000.0, 0.5, 0.015123665829158384, 0.8516553246133607), + (547.0, 1000.0, 0.5, 0.00030353498700958767, 0.9986772350804279), + (999.0, 1000.0, 0.5, 9.332636185032189e-299, 1.0), + (0.0, 1000.0, 0.3, 1.2532566399657381e-155, 1.2532566399657381e-155), + (1.0, 1000.0, 0.3, 5.371099885567449e-153, 5.383632451967106e-153), + (257.0, 1000.0, 0.3, 0.00030959977077928017, 0.0014835731103312242), + (286.0, 1000.0, 0.3, 0.0174263984921581, 0.1759288949498673), + (300.0, 1000.0, 0.3, 0.027521003821268385, 0.5155935198141203), + (314.0, 1000.0, 0.3, 0.017109129819617666, 0.8414952343056957), + (343.0, 1000.0, 0.3, 0.00036310112605077133, 0.9985016810663147), + (500.0, 1000.0, 0.3, 3.4791609391626546e-40, 1.0), + (999.0, 1000.0, 0.3, 0.0, 1.0), + (0.0, 1000.0, 0.001, 0.36769542477096406, 0.36769542477096406), + (1.0, 1000.0, 0.001, 0.36806348825922325, 0.7357589130301873), + (2.0, 1000.0, 0.001, 0.18403174412961162, 0.9197906571597989), + (4.0, 1000.0, 0.001, 0.01528995542083862, 0.9963631219704782), + (500.0, 1000.0, 0.001, 0.0, 1.0), + (999.0, 1000.0, 0.001, 0.0, 1.0), + (0.0, 10000.0, 0.5, 0.0, 0.0), + (1.0, 10000.0, 0.5, 0.0, 0.0), + (4850.0, 10000.0, 0.5, 8.861480551702583e-05, 0.0013939837588563594), + (4950.0, 10000.0, 0.5, 0.004839495142316416, 0.161087099897656), + (5000.0, 10000.0, 0.5, 0.007978646139382154, 0.5039893230696911), + (5050.0, 10000.0, 0.5, 0.004839495142316416, 0.8437523952446604), + (5150.0, 10000.0, 0.5, 8.861480551702583e-05, 0.9986946310466607), + (9999.0, 10000.0, 0.5, 0.0, 1.0), + (0.0, 10000.0, 0.3, 0.0, 0.0), + (1.0, 10000.0, 0.3, 0.0, 0.0), + (2863.0, 10000.0, 0.3, 9.71726742023389e-05, 0.0013927715709303387), + (2954.0, 10000.0, 0.3, 0.0052756073908582085, 0.16038574844431697), + (3000.0, 10000.0, 0.3, 0.008705361365067056, 0.504932983819298), + (3046.0, 10000.0, 0.3, 0.005244996814931252, 0.8448693061622872), + (3137.0, 10000.0, 0.3, 0.00010232641238302388, 0.9986013755801396), + (5000.0, 10000.0, 0.3, 0.0, 1.0), + (9999.0, 10000.0, 0.3, 0.0, 1.0), + (0.0, 10000.0, 0.0001, 0.3678610464329299, 0.3678610464329299), + (1.0, 10000.0, 0.0001, 0.3678978362165516, 0.7357588826494815), + (2.0, 10000.0, 0.0001, 0.1839489181082758, 0.9197078007577573), + (4.0, 10000.0, 0.0001, 0.015324477632794846, 0.996342452182816), + (5000.0, 10000.0, 0.0001, 0.0, 1.0), + (9999.0, 10000.0, 0.0001, 0.0, 1.0), + (0.0, 10000.0, 0.999, 0.0, 0.0), + (1.0, 10000.0, 0.999, 0.0, 0.0), + (5000.0, 10000.0, 0.999, 0.0, 0.0), + (9981.0, 10000.0, 0.999, 0.0037205945263045544, 0.007158161342006954), + (9987.0, 10000.0, 0.999, 0.07292252318084555, 0.20834867055478343), + (9990.0, 10000.0, 0.999, 0.12517263665023892, 0.5421328760209826), + (9993.0, 10000.0, 0.999, 0.09007019224850273, 0.8699847268232276), + (9999.0, 10000.0, 0.999, 0.00045218564541589875, 0.9999548266540229), + (0.0, 100000.0, 0.5, 0.0, 0.0), + (1.0, 100000.0, 0.5, 0.0, 0.0), + (49526.0, 100000.0, 0.5, 2.821098952769392e-05, 0.00137358872189177), + (49842.0, 100000.0, 0.5, 0.0015314620321580088, 0.15959653953914368), + (50000.0, 100000.0, 0.5, 0.00252312621419674, 0.5012615631070984), + (50158.0, 100000.0, 0.5, 0.0015314620321580088, 0.8419349224930143), + (50474.0, 100000.0, 0.5, 2.821098952769392e-05, 0.9986546222676359), + (99999.0, 100000.0, 0.5, 0.0, 1.0), + (0.0, 100000.0, 0.3, 0.0, 0.0), + (1.0, 100000.0, 0.3, 0.0, 0.0), + (29565.0, 100000.0, 0.3, 3.016598258904529e-05, 0.0013408650955754554), + (29855.0, 100000.0, 0.3, 0.0016703029966719045, 0.15934769308464344), + (30000.0, 100000.0, 0.3, 0.002752954648397428, 0.5015600059172672), + (30145.0, 100000.0, 0.3, 0.0016672321665529399, 0.8423208126454166), + (30435.0, 100000.0, 0.3, 3.067097885268188e-05, 0.9986570625478413), + (50000.0, 100000.0, 0.3, 0.0, 1.0), + (99999.0, 100000.0, 0.3, 0.0, 1.0), + (0.0, 100000.0, 0.99, 0.0, 0.0), + (1.0, 100000.0, 0.99, 0.0, 0.0), + (50000.0, 100000.0, 0.99, 0.0, 0.0), + (98906.0, 100000.0, 0.99, 0.00015962501729101418, 0.001680685850467731), + (98969.0, 100000.0, 0.99, 0.007723834924565287, 0.1660975153530963), + (99000.0, 100000.0, 0.99, 0.012678161323544589, 0.5042686879721229), + (99031.0, 100000.0, 0.99, 0.007885853372194001, 0.8416321140066325), + (99094.0, 100000.0, 0.99, 0.00013276363518518057, 0.9988439452168577), + (99999.0, 100000.0, 0.99, 0.0, 1.0), +]; + +// Columns: x, mean, mass, cdf. The last block has x > 1000, which takes the incomplete-gamma branch. +#[rustfmt::skip] +const POISSON_REF: &[(f64, f64, f64, f64)] = &[ + (0.0, 0.5, 0.6065306597126334, 0.6065306597126334), + (1.0, 0.5, 0.3032653298563167, 0.9097959895689501), + (3.0, 0.5, 0.012636055410679864, 0.9982483774437092), + (0.0, 5.0, 0.006737946999085467, 0.006737946999085467), + (2.0, 5.0, 0.08422433748856833, 0.12465201948308115), + (5.0, 5.0, 0.1754673697678507, 0.6159606548330632), + (10.0, 5.0, 0.018132788707821874, 0.986304731401617), + (25.0, 5.0, 1.294589320478486e-10, 0.9999999999695003), + (70.0, 100.0, 0.00031056089812001654, 0.0009714440282773882), + (100.0, 100.0, 0.039860996809147134, 0.5265621985299984), + (130.0, 100.0, 0.000575252683815463, 0.9982931596294985), + (900.0, 1000.0, 7.516954352125952e-05, 0.0006977673277963068), + (1000.0, 1000.0, 0.012614611348721499, 0.508409367168506), + (1100.0, 1000.0, 9.498944242299507e-05, 0.9991323590365564), + (1050.0, 1000.0, 0.003599490404803648, 0.9439711616363267), + (1200.0, 1000.0, 7.992642848843571e-11, 0.9999999996115061), + (1500.0, 2000.0, 1.878037542328086e-32, 7.468034373687017e-32), + (2000.0, 2000.0, 0.00892024889598624, 0.5059467004517837), + (2100.0, 2000.0, 0.0007442234118916615, 0.9872092799080792), + (2300.0, 2000.0, 4.011826753210622e-12, 0.999999999973895), + (4800.0, 5000.0, 9.987733294132511e-05, 0.002269403878343364), + (5000.0, 5000.0, 0.005641801804664022, 0.5037611677708467), +]; + +fn rel_err(got: Value, want: f64) -> f64 { + let g = got.as_f64().unwrap_or(f64::NAN); + if want == 0.0 { g.abs() } else { ((g - want) / want).abs() } +} + +/// Error bound the grid below must meet: 1e-13, plus 2e-16 per unit of |ln value| for far-tail +/// results (an exp of a sum of large logs cannot do better in double precision). Measured: at +/// most 9e-14 on the grid (the n = 100,000 cumulative cases), 2e-14 for n <= 10,000 away from the +/// far tails, 2e-15 on the Excel-sized central cases. The log-gamma code was off by 1.5e-11. +fn allowed(want: f64) -> f64 { + 1e-13 + 2e-16 * want.abs().ln().abs() +} + +#[test] +fn binom_and_poisson_match_exact_decimal_references() { + // Worst error relative to the allowance (above 1 is a failure) and plain worst error. + let (mut worst, mut worst_plain) = (0.0f64, 0.0f64); + let mut check = |what: &str, got: Value, want: f64| { + let e = rel_err(got, want); + let ratio = e / allowed(want); + if ratio > 1.0 { + println!("{what}: error {e:e}, allowed {:e}", allowed(want)); + } + worst = worst.max(ratio); + if want > 1e-3 { + worst_plain = worst_plain.max(e); + } + }; + for &(k, nn, p, pmf, cdf) in BINOM_REF { + check(&format!("BINOM.DIST({k},{nn},{p},FALSE)"), ev("BINOM.DIST", vec![n(k), n(nn), n(p), b(false)]), pmf); + check(&format!("BINOM.DIST({k},{nn},{p},TRUE)"), ev("BINOM.DIST", vec![n(k), n(nn), n(p), b(true)]), cdf); + } + for &(x, m, pmf, cdf) in POISSON_REF { + check(&format!("POISSON.DIST({x},{m},FALSE)"), ev("POISSON.DIST", vec![n(x), n(m), b(false)]), pmf); + check(&format!("POISSON.DIST({x},{m},TRUE)"), ev("POISSON.DIST", vec![n(x), n(m), b(true)]), cdf); + } + println!("worst error / allowance {worst:.3}; worst relative error where the result is above 1e-3: {worst_plain:e}"); + assert!(worst <= 1.0); + assert!(worst_plain < 1e-13); +} + +#[test] +fn binom_and_poisson_edge_cases() { + for cum in [false, true] { + let pm = |k: f64, nn: f64, p: f64| ev("BINOM.DIST", vec![n(k), n(nn), n(p), b(cum)]); + // p = 0 and p = 1 + close(pm(0.0, 50.0, 0.0), 1.0); + close(pm(3.0, 50.0, 1.0), 0.0); + close(pm(50.0, 50.0, 1.0), 1.0); + // n = 0 + close(pm(0.0, 0.0, 0.3), 1.0); + // k = n is certain when cumulative and p^n otherwise + close(pm(20.0, 20.0, 0.5), if cum { 1.0 } else { 0.5f64.powi(20) }); + close_tol(pm(10_000.0, 10_000.0, 0.999), if cum { 1.0 } else { 0.999f64.powi(10_000) }, 1e-12); + // k = 0 is q^n + close_tol(pm(0.0, 1000.0, 0.001), 0.999f64.powi(1000), 1e-12); + close_tol(pm(0.0, 100_000.0, 0.00001), (1.0f64 - 1e-5).powi(100_000), 1e-9); + } + close(ev("POISSON.DIST", vec![n(0.0), n(0.0), b(false)]), 1.0); + close(ev("POISSON.DIST", vec![n(4.0), n(0.0), b(false)]), 0.0); + close(ev("POISSON.DIST", vec![n(4.0), n(0.0), b(true)]), 1.0); + close_tol(ev("POISSON.DIST", vec![n(0.0), n(3.0), b(true)]), (-3.0f64).exp(), 1e-15); + // Fractions truncate as before. + close(ev("BINOM.DIST", vec![n(6.9), n(10.9), n(0.5), b(false)]), 0.205078125); + // The 10,000 trial switch of the method no longer shows: neighbours agree to 1e-14. + let lo = ev("BINOM.DIST", vec![n(5000.0), n(10_000.0), n(0.5), b(true)]).as_f64().unwrap(); + let hi = ev("BINOM.DIST", vec![n(5000.0), n(10_001.0), n(0.5), b(true)]).as_f64().unwrap(); + let pmf = ev("BINOM.DIST", vec![n(5000.0), n(10_000.0), n(0.5), b(false)]).as_f64().unwrap(); + // cdf(5000; 10,001) = cdf(5000; 10,000) - p * mass(5000; 10,000) by the recurrence on n + assert!(((hi - (lo - 0.5 * pmf)) / hi).abs() < 1e-13, "{lo} {hi}"); + // Inverse stays consistent: the smallest k with cdf >= alpha. + close(ev("BINOM.INV", vec![n(10_000.0), n(0.5), n(0.5)]), 5000.0); +}