mirror of
https://github.com/storytold/gridcraft.git
synced 2026-11-10 04:21:45 +00:00
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 <noreply@anthropic.com>
This commit is contained in:
vangel76andClaude Sonnet 5.5 committed
1 parent
a305229287
commit
22c06b7a4a
4 files changed
+436
-26
No files matched your search
@@ -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 }
|
||||
}
|
||||
|
||||
|
||||
@@ -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};
|
||||
|
||||
|
||||
@@ -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<Value> {
|
||||
|
||||
@@ -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);
|
||||
}
|
||||
Reference in new issue
Block a user