fix: correct Geometric::inverse_cdf against its definition#411
Open
agene0001 wants to merge 2 commits into
Open
fix: correct Geometric::inverse_cdf against its definition#411agene0001 wants to merge 2 commits into
agene0001 wants to merge 2 commits into
Conversation
`approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it
consults the ULPs distance. The crate's wrapper paired it with
`DEFAULT_EPS = 1e-9`, which made the ULPs bound unreachable, so every internal
`ulps_eq!(x, y)` was really a 1e-9 absolute comparison.
That matters because the macro is used to recognise *exact* parameter values, so
anything within 1e-9 of them silently took a degenerate branch:
Binomial(p = 1 - 2^-33, n = 100).pmf(99) 0 (true 1.16e-8)
Binomial(p = 1 - 2^-33, n = 100).ln_pmf(99) -inf
Geometric(p = 1 - 2^-33).max() 1 (true u64::MAX)
Geometric(p = 1 - 2^-33).skewness() inf (true 92681.9)
Beta(1 + 2^-33, 1 + 2^-33).pdf(0.0) 1 (true 0)
digamma(-1 + 2^-33) -inf (true -8.59e9)
Adds `DEFAULT_ULPS_EPS = f64::EPSILON` and uses it for `ulps_eq!` and
`assert_ulps_eq!`. The exact values still take the degenerate branches -
`Binomial(1.0, 5).pmf(5) == 1`, `Geometric(1.0).max() == 1`, `digamma` at the
poles is still -inf - all covered by existing tests.
Tests use `1 - 2^-33` rather than `1 - 1e-10` so that `1 - p` is exact and the
references are not limited by the representation of `p`.
Also fixes `Binomial::entropy`, which summed `-p * ln(p)` unguarded and so
returned `NaN` once any mass underflowed to zero (`0 * -inf`) - for `p = 0.5`
that is every `n` past about 1100. `Categorical::entropy` already filtered zero
terms; this brings `Binomial` in line.
`Geometric` overrides `inverse_cdf` with the closed form `ceil(ln1p(-x) / ln1p(-p))`. The quotient is only *approximately* integral at the step boundaries, so `ceil` lands on either side of the true one and the round-trip fails for most `p`: over a sweep of 200 values of `p` and 400 of `k`, 7255 of 16200 `(p, k)` pairs had `inverse_cdf(cdf(k)) != k`. Being an override, this was untouched by the work on the default `inverse_cdf` (statrs-dev#390, statrs-dev#391), which is why statrs-dev#342 stayed open after the explicit implementation it asks for had landed. The closed form is now corrected against the definition being inverted - the least `k` with `cdf(k) >= x`: * ordinary rounding puts it within one step, so the fast path is two or three `cdf` evaluations; * once `x` is within a few ulp of 1, `1 - x` has lost every significant bit and the closed form can be off without bound (`cdf` has a plateau roughly `1/p` wide there - at `p = 1e-9` the answer is 36331335444), so it falls back to a bisection that is exact by definition. After: 0 genuine round-trip failures (the 402 remaining are all cases where `cdf(k-1) == cdf(k)` in f64, so no integer is recoverable), and 0 violations of the definition over 30266 inputs including `x` within ulps of 0 and 1. This also un-ignores `test_inverse_cdf_small_p` and `test_inverse_cdf_extreme_tail`, marked `#[ignore = "Gaps in pathological corner cases and testing"]`. The former was masking two separate bugs: it needs both this fix and the `ulps_eq!` epsilon fix, because it constructs `p = 1 - 1e-14`, which the old 1e-9 epsilon snapped to the degenerate `p == 1` branch. Cost: 8 -> 23 ns/call for the verifying `cdf` evaluations. Closes statrs-dev#342
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Geometricoverridesinverse_cdfwith the closed formceil(ln1p(-x) / ln1p(-p)). The quotient is only approximately integral at the step boundaries, soceillands on either side of the true one: over a sweep of 200 values ofp× 400 ofk, 7255 of 16200 pairs hadinverse_cdf(cdf(k)) != k.Being an override, this was unreachable by the improvements to the default
inverse_cdf(#390, #391) — which is why #342 stayed open after the explicit implementation it asks for had landed.The closed form is now corrected against the definition being inverted — the least
kwithcdf(k) >= x:cdfevaluations (8 → 23 ns/call);xis within a few ulp of 1,1 − xhas lost every significant bit and the closed form can be off without bound (cdfhas a plateau ~1/p wide there; at p = 1e-9 the correct answer is 36331335444), so it falls back to a bisection that is exact by definition.After: 0 genuine round-trip failures (the 402 remaining are cases where
cdf(k-1) == cdf(k)in f64, so no integer is recoverable), and 0 violations of the definition over 30266 inputs includingxwithin ulps of 0 and 1.This also un-ignores
test_inverse_cdf_small_pandtest_inverse_cdf_extreme_tail(#[ignore = "Gaps in pathological corner cases and testing"]). The former was masking two independent bugs: it constructsp = 1 − 1e-14, which the oldulps_eq!epsilon (#410) snapped to the degeneratep == 1branch, returning 1 instead ofu64::MAX. It needs both fixes — hence the stacking.Closes #342
🤖 Generated with Claude Code