Skip to content

fix: compute Hypergeometric pmf in log space for large populations - #430

Open
teddytennant wants to merge 1 commit into
statrs-dev:mainfrom
teddytennant:fix/hypergeometric-pmf-large-426
Open

fix: compute Hypergeometric pmf in log space for large populations#430
teddytennant wants to merge 1 commit into
statrs-dev:mainfrom
teddytennant:fix/hypergeometric-pmf-large-426

Conversation

@teddytennant

Copy link
Copy Markdown
Contributor

Problem

Hypergeometric::pmf multiplies and divides three binomial coefficients held in f64. Once any coefficient exceeds f64::MAX, the result is wrong:

  • Hypergeometric::new(1030, 1, 515): pmf(0) and pmf(1) return 0.0 (exact value 0.5) because only the denominator overflows.
  • Hypergeometric::new(20000, 200, 300): pmf(0) and pmf(3) return NaN because numerator and denominator both overflow (inf/inf).

ln_pmf, cdf, and mean on the same objects were already correct.

Root cause

Direct evaluation of (K choose x) * (N-K choose n-x) / (N choose n) via factorial::binomial, which returns a rounded f64 that can be +inf for large arguments.

Fix

Compute the PMF in log space, matching the Multinomial fix style already in this repo:

  • pmf(x) delegates to self.ln_pmf(x).exp()
  • ln_pmf returns f64::NEG_INFINITY for out-of-support x (x > draws, x > successes, or x < min()), so pmf returns 0.0 consistently
  • cdf/sf already sum exp of log-binomial terms and are left unchanged

Test plan

  • Existing hypergeometric unit tests (pmf cases that used exact equality now use absolute tolerance where log-space rounding differs at ~1e-15)
  • Regression tests from Hypergeometric::pmf returns 0.0 or NaN for large populations #426 for populations 1030 and 20000
  • Out-of-support guards for pmf / ln_pmf
  • cargo test --lib (774 passed, 2 ignored)

Fixes #426

pmf multiplied and divided binomial coefficients in f64. Once a
coefficient exceeded f64::MAX the result became 0.0 (denominator-only
overflow) or NaN (inf/inf). cdf and ln_pmf already used ln_binomial.

Delegate pmf to ln_pmf(x).exp(), matching the Multinomial fix, and
guard ln_pmf for out-of-support x so it returns NEG_INFINITY instead of
relying on ln_binomial edge cases alone. Adds regression tests from statrs-dev#426.

Fixes statrs-dev#426
@teddytennant

Copy link
Copy Markdown
Contributor Author

The Coverage (nightly) failure here is essentially 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 gives 91 — the one extra is the f64::NEG_INFINITY added in hypergeometric.rs, which follows the same convention as the existing code in geometric.rs and the rest of the crate. Fixing just that line would not make the job green. Without -D warnings everything compiles (warnings only). Last green Coverage run on main was Jul 31.

cargo test --all-features passes on this branch: 778 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 — and I can rebase this on top of it.

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.

Hypergeometric::pmf returns 0.0 or NaN for large populations

1 participant