Skip to content

fix: correct Geometric::inverse_cdf against its definition - #431

Open
teddytennant wants to merge 2 commits into
statrs-dev:mainfrom
teddytennant:fix/inverse-cdf-geometric-342
Open

fix: correct Geometric::inverse_cdf against its definition#431
teddytennant wants to merge 2 commits into
statrs-dev:mainfrom
teddytennant:fix/inverse-cdf-geometric-342

Conversation

@teddytennant

Copy link
Copy Markdown
Contributor

Summary

Geometric::inverse_cdf used the closed form ceil(ln1p(-x) / ln1p(-p)) alone. That quotient is only approximately integral at the step boundaries, so ceil often lands on the wrong side of the true quantile and breaks the discrete inverse definition (least k with cdf(k) >= x).

Over a sweep of 200 values of p and 400 of k, thousands of pairs failed inverse_cdf(cdf(k)) == k even though cdf separated adjacent integers.

Root cause

Floating-point evaluation of the closed form is not exact at bucket edges. Being an override, this path was also unreachable by later improvements to the default discrete inverse_cdf.

Fix

Keep the closed form as a fast candidate, then correct it against the definition:

  • check candidate, candidate - 1, and candidate + 1 with two or three cdf evaluations (the usual path);
  • when x is within a few ulp of 1, 1 - x loses significant bits and the closed form can be off without bound (cdf plateaus about 1/p wide), so fall back to exact integer bisection.

Test plan

  • cargo test --lib geometric::
  • cargo test --lib (775 passed, 2 ignored pre-existing pathological cases still marked #[ignore])
  • New regressions:
    • test_inverse_cdf_round_trips_and_matches_definition
    • test_inverse_cdf_is_least_k_with_cdf_at_least_x
    • test_inverse_cdf_long_plateau

Fixes #342

The closed form ceil(ln1p(-x)/ln1p(-p)) is only approximately integral at
step boundaries, so plain ceil often returns the wrong side of the true
quantile and breaks inverse_cdf(cdf(k)) == k (statrs-dev#342).

Correct the candidate with a short cdf check (and bisection near x ~ 1
where 1-x loses precision). Add round-trip and definition regression tests.
cdf(1) is mathematically p, but evaluating via expm1/ln1p can undershoot
by a ulp on MSVC, so the definition check in inverse_cdf returned 2 for
x=p (Windows CI: Expected 1, got 2 for p=0.2).

Special-case cdf(1)/sf(1) and short-circuit inverse_cdf for x <= p.
@teddytennant

Copy link
Copy Markdown
Contributor Author

The Coverage (nightly) failure here is unrelated to this change.

Recent nightly (1.99.0-nightly, 2026-08-06) started firing the deprecated lint on f64::NAN / INFINITY / NEG_INFINITY paths in the ~25 modules that do use core::f64; — the module import shadows the primitive in path resolution, so those resolve to the deprecated core::f64::* module constants rather than the associated constants. The Coverage job sets RUSTFLAGS: -D warnings, so it becomes a hard build failure.

It reproduces on a clean checkout of main (02fdab6) with nothing applied: 90 errors. This branch also produces exactly 90 — it adds no new instances. Without -D warnings it compiles fine (warnings only). Last green Coverage run on main was Jul 31.

cargo test --all-features passes on this branch: 780 passed, 0 failed.

Happy to send a separate PR swapping those paths over to the associated constants (or dropping the use core::f64; imports) if you want Coverage unblocked.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

issue in default inverse_cdf, found with simple case of Geometric

1 participant