diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..315c26ab 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -237,7 +237,11 @@ impl Distribution for Binomial { } else { (0..self.n + 1).fold(0.0, |acc, x| { let p = self.pmf(x); - acc - p * p.ln() + // A mass that underflows to zero contributes nothing, taking + // `0 * ln 0 == 0` as usual for entropy. Evaluating it would give + // `0 * -inf == NaN` and poison the whole sum, which made + // `entropy()` return `NaN` for any `n` past roughly 1100. + if p > 0.0 { acc - p * p.ln() } else { acc } }) }; Some(entr) @@ -412,6 +416,39 @@ mod tests { test_exact(0.3, 10, 10, max); } + /// `p` very close to but not equal to 1 must take the general path, not the + /// degenerate `p == 1` branch. The branch is guarded by `prec::ulps_eq!`, + /// whose default epsilon was `1e-9` *absolute*, so any `p` within `1e-9` of + /// 1 collapsed to a point mass at `n`. + /// + /// `1 - 2^-33` is used rather than `1 - 1e-10` so that `1 - p` is exact and + /// the reference values are not limited by the representation of `p`. + #[test] + fn test_pmf_p_near_one_is_not_degenerate() { + let n = Binomial::new(1.0 - f64::powi(2.0, -33), 100).unwrap(); + prec::assert_relative_eq!(n.pmf(99), 1.1641532048523463366e-8, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(n.pmf(100), 0.99999998835846788439, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(n.ln_pmf(99), -18.268686784015220704, epsilon = 0.0, max_relative = 5e-15); + // entropy of a non-degenerate distribution is strictly positive + assert!(n.entropy().unwrap() > 0.0); + } + + /// `entropy()` sums `-p * ln(p)` over the whole support, so it used to + /// return `NaN` as soon as one mass underflowed to zero (`0 * -inf`). For + /// `p = 0.5` that starts at around `n = 1100`. + #[test] + fn test_entropy_with_underflowing_masses() { + for n in [1_100, 2_000, 5_000, 20_000] { + let d = Binomial::new(0.5, n).unwrap(); + let h = d.entropy().unwrap(); + assert!(h.is_finite(), "entropy for n={n} was {h}"); + // 0.5 * ln(2 * pi * e * n * p * q) is the asymptotic form + let approx = 0.5 * (2.0 * f64::consts::PI * f64::consts::E * n as f64 * 0.25).ln(); + prec::assert_relative_eq!(h, approx, epsilon = 0.0, max_relative = 1e-3); + } + assert!(Binomial::new(0.1, 5_000).unwrap().entropy().unwrap().is_finite()); + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Binomial| x.pmf(arg); diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index f74f2a79..b52e1b19 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -203,7 +203,48 @@ impl DiscreteCDF for Geometric { self.p, x ); } - k as u64 + // `ln1p(-x) / ln1p(-p)` is only approximately integral at the step + // boundaries, so `ceil` lands on either side of the true one and + // `inverse_cdf(cdf(k)) != k` for most `p` (statrs-dev/statrs#342). + // Everything below corrects the closed form against the definition + // being inverted: the least `k` with `cdf(k) >= x`. + let candidate = (k.max(1.0) as u64).max(self.min()); + let is_answer = |k: u64| self.cdf(k) >= x && (k <= self.min() || self.cdf(k - 1) < x); + + // Ordinary rounding puts the closed form within one step, so this is + // the path essentially always taken - two or three `cdf` evaluations. + if is_answer(candidate) { + return candidate; + } + if candidate > self.min() && is_answer(candidate - 1) { + return candidate - 1; + } + if is_answer(candidate + 1) { + return candidate + 1; + } + + // Once `x` is within a few ulp of 1, `1 - x` has lost every significant + // bit and the closed form can be off by an unbounded number of steps + // (`cdf` has a long plateau there - about 1/p wide). Bisect instead, + // which is exact by definition regardless of how far off `candidate` is. + let mut hi = candidate; + while self.cdf(hi) < x { + let Some(doubled) = hi.checked_mul(2) else { + hi = u64::MAX; + break; + }; + hi = doubled; + } + let mut lo = self.min(); + while lo < hi { + let mid = lo + (hi - lo) / 2; + if self.cdf(mid) >= x { + hi = mid; + } else { + lo = mid + 1; + } + } + lo } } @@ -572,6 +613,74 @@ mod tests { density_util::check_discrete_distribution(&create_ok(1.0), 1); } + /// `inverse_cdf` must return the least `k` with `cdf(k) >= x`, and so must + /// round-trip `cdf` wherever `cdf` is injective. + /// + /// The closed form `ceil(ln1p(-x) / ln1p(-p))` is only approximately + /// integral at the step boundaries, so on its own it landed on the wrong + /// side for 7255 of 16200 (p, k) pairs (statrs-dev/statrs#342). + #[test] + fn test_inverse_cdf_round_trips_and_matches_definition() { + let mut mismatches = 0; + for i in 1..200 { + let p = i as f64 / 200.0; + let g = create_ok(p); + for k in 1..=400u64 { + let c = g.cdf(k); + if c >= 1.0 { + break; + } + // only meaningful where `cdf` actually separates k-1 from k + if g.cdf(k - 1) < c && g.inverse_cdf(c) != k { + mismatches += 1; + } + } + } + assert_eq!(mismatches, 0, "cdf/inverse_cdf round-trip failures"); + } + + #[test] + fn test_inverse_cdf_is_least_k_with_cdf_at_least_x() { + for i in 1..60 { + let p = i as f64 / 60.0; + let g = create_ok(p); + // Chained rather than collected, so the test also builds under + // `no_std`, where there is no `Vec`. The second and third groups + // put x within a few ulp of 1, where `1 - x` has lost every + // significant bit and the closed form needs the bisection fallback. + let xs = (1..400) + .map(|j| j as f64 / 400.0) + .chain((1..6).map(|u| 1.0 - u as f64 * f64::EPSILON / 2.0)) + .chain((1..6).map(|u| u as f64 * f64::MIN_POSITIVE)) + .chain((1..400).map(|j| g.cdf(j))); + for x in xs { + // `x == 1` saturates to `u64::MAX` by design, iCDF(1) = +inf + if !(0.0..=1.0).contains(&x) || x == 1.0 { + continue; + } + let k = g.inverse_cdf(x); + assert!(g.cdf(k) >= x, "p={p} x={x:e}: cdf({k}) < x"); + assert!( + k <= g.min() || g.cdf(k - 1) < x, + "p={p} x={x:e}: {k} is not the least such k" + ); + } + } + } + + /// Tiny `p` gives `cdf` a plateau roughly `1/p` wide near 1, so the + /// bisection fallback has to cover an unbounded number of steps. + #[test] + fn test_inverse_cdf_long_plateau() { + for p in [1e-3f64, 1e-6, 1e-9] { + let g = create_ok(p); + let x = 1.0 - f64::EPSILON / 2.0; + let k = g.inverse_cdf(x); + assert!(g.cdf(k) >= x, "p={p:e}: cdf({k}) < x"); + assert!(k <= g.min() || g.cdf(k - 1) < x, "p={p:e}: {k} not least"); + } + } + #[test] fn test_inverse_cdf() { let invcdf = |arg: f64| move |x: Geometric| x.inverse_cdf(arg); @@ -605,8 +714,9 @@ mod tests { distribution.inverse_cdf(0.5); } + /// Un-ignored: this failed until `inverse_cdf` was corrected against the + /// definition (statrs-dev/statrs#342). #[test] - #[ignore = "Gaps in pathological corner cases and testing"] fn test_inverse_cdf_small_p() { let ident_prob = |arg: f64| move |x: Geometric| x.cdf(x.inverse_cdf(arg)); let ident_count = |count: u64| move |x: Geometric| x.inverse_cdf(x.cdf(count)); @@ -619,9 +729,8 @@ mod tests { test_relative(0.5, q, ident_prob(q)); } - #[test] - #[ignore = "Gaps in pathological corner cases and testing"] /// large k, small p regime + #[test] fn test_inverse_cdf_extreme_tail() { // k_below(k): probability midpoint of the k-th bucket — strictly between cdf(k-1) and // cdf(k) — so inverse_cdf must return k. diff --git a/src/function/gamma.rs b/src/function/gamma.rs index bff5c45a..771a62d1 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -1195,6 +1195,25 @@ mod tests { assert!(super::checked_gamma_ui(1.0, f64::INFINITY).is_err()); } + /// `digamma` has poles at the non-positive integers and returns -inf there. + /// The pole test is `prec::ulps_eq!(x.floor(), x)`, whose default epsilon + /// was `1e-9` absolute, so a whole neighbourhood around each pole returned + /// -inf instead of a large finite value. + /// + /// The tolerance is loose because `digamma(x < 0)` goes through the + /// reflection formula and `(PI * x).tan()` loses roughly + /// `ulp(PI) / dist_to_pole` of relative accuracy; what is pinned here is + /// that the values are finite and of the right magnitude. + #[test] + fn test_digamma_near_negative_integer_poles_is_finite() { + prec::assert_relative_eq!(digamma(-1.0 + f64::powi(2.0, -33)), -8589934591.5772156646, epsilon = 0.0, max_relative = 1e-5); + prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)), 8589934592.4227843348, epsilon = 0.0, max_relative = 1e-5); + // the poles themselves are unchanged + assert_eq!(digamma(-1.0), f64::NEG_INFINITY); + assert_eq!(digamma(0.0), f64::NEG_INFINITY); + assert_eq!(digamma(-5.0), f64::NEG_INFINITY); + } + // TODO: precision testing could be more accurate #[test] fn test_digamma() { diff --git a/src/prec.rs b/src/prec.rs index ebf1b6d4..5e23b05b 100644 --- a/src/prec.rs +++ b/src/prec.rs @@ -29,6 +29,7 @@ //! - `DEFAULT_RELATIVE_ACC`: 1e-14 for relative comparisons //! - `DEFAULT_EPS`: 1e-9 for absolute comparisons //! - `DEFAULT_ULPS`: 5 for ULPs comparisons +//! - `DEFAULT_ULPS_EPS`: `f64::EPSILON`, the absolute floor paired with `DEFAULT_ULPS` //! //! These defaults should be used unless there is a specific reason to use different //! precision levels. @@ -62,6 +63,16 @@ pub const DEFAULT_EPS: f64 = 1e-9; /// Default and target ULPs accuracy for f64 operations pub const DEFAULT_ULPS: u32 = 5; +/// Default absolute epsilon for ULPs comparisons. +/// +/// `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it looks +/// at the ULPs distance, so pairing it with [`DEFAULT_EPS`] (`1e-9`) would make +/// the ULPs bound unreachable and turn `ulps_eq!(x, y)` into a `1e-9` absolute +/// comparison. `ulps_eq!` is used inside the crate to recognise exact parameter +/// values (`p == 1.0`, `x == x.floor()`), so its epsilon has to stay at the +/// scale of a genuine rounding error. +pub const DEFAULT_ULPS_EPS: f64 = f64::EPSILON; + /// Compares if two floats are close via `approx::abs_diff_eq` /// using a maximum absolute difference (epsilon) of `acc`. #[deprecated(since = "0.19.0", note = "Use abs_diff_eq! macro instead")] @@ -144,7 +155,7 @@ mod macros { ); redefine_two_opt_approx_macro!( ulps_eq, - { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS } + { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS } ); pub(crate) use abs_diff_eq; @@ -162,7 +173,7 @@ mod macros { ); redefine_two_opt_approx_macro!( assert_ulps_eq, - { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS } + { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS } ); pub(crate) use assert_abs_diff_eq;