Fix pdf/cdf domain handling in dist (StudentT negative cdf, Binomial overflow/infinite loop, off-support inputs)#112
Open
gaoflow wants to merge 1 commit into
Open
Conversation
…s in ops
- StudentT cdf returned cdf(-x) - 0.5 instead of 1 - cdf(-x) on the
negative half axis
- Binomial pdf computed C(n, m) in usize: silent overflow garbage from
n = 30 in release (pdf(15 | 30, 0.5) = 0.0131 instead of 0.1445),
panic in debug, and an infinite loop for x > n through C(n, n - r)
usize underflow; now evaluated in log-space via ln_gamma with a
support check
- Binomial cdf panicked in inc_beta for k < 0, k = n, k > n and
interpolated fractional k; now clamps the domain and floors k
(scipy/R semantics)
- Bernoulli pdf returned 1 - p for any x outside {0, 1}
- ops C / P return 0 for r > n and use checked multiplication (loud
panic instead of silent wraparound); C uses the multiplicative
formula, exact for all n <= 62 on 64-bit
Reference values in the new tests are cross-checked against scipy 1.17.1.
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
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.
While cross-checking
statistics::distagainst scipy 1.17.1 (all nine distributions, pdf + cdf, on and beyond their supports), I found the implementations agree to ~1e-7 on the happy path but four families break as soon as the input leaves it:StudentT cdf is wrong on the entire negative half. The negative branch used
self.cdf(-x) - 0.5where the reflection is1 - self.cdf(-x):Binomial pdf hangs the process for x > n and silently returns garbage from n ≈ 30. The pmf went through
ops::C(n, m)inusizearithmetic.P(30, 15)already overflows, so in releaseC(30, 15)came out 14052247 instead of 155117520 —Binomial(30, 0.5).pdf(15.0)returned 0.0131 instead of 0.1445 with no error (debug builds panic instead). Worse, for m > n the reductionC(n, n - r)underflowsusizeand the recursion cycles:Binomial(10, 0.3).pdf(12.0)spins at 100% CPU forever in a release build. The pmf is now computed in log-space via the existingln_gamma(exact-enough for any n, e.g. n = 100) behind a support check.Binomial cdf panics at the top of its own support.
cdfpassed out-of-domain arguments toinc_beta:Binomial(10, 0.3).cdf(10.0)panicked with "Bad a or b in routine betai", same for k > n and k < 0, and fractional k was interpolated instead of floored. Now clamped (k < 0 → 0, k ≥ n → 1) and floored, matching scipy/R semantics.Bernoulli pdf returned 1 − p for any x outside {0, 1} —
pdf(-1.0),pdf(0.5),pdf(2.0)all gave 0.7 for p = 0.3. Now 0 off support.At the root,
ops::C/Pnow return 0 for r > n (so the r > n/2 reduction can't underflow and recurse forever) and use checked multiplication, so genuine overflow is a clear panic instead of silent wraparound in release.Cuses the multiplicative formula, exact for everyC(n, r)with n ≤ 62 on 64-bit; the old path broke around n = 30.factorial/double_factorialget the same checked treatment.Tests: regression tests for every case above with scipy-cross-checked values, including the previously-hanging call, plus a new
tests/ops.rsforC/P. Re-running the full sweep (9 distributions × pdf/cdf, 101 points) against this branch gives zero divergences from scipy.cargo test(default and--no-default-features),cargo fmt,cargo clippy --all-targetsare clean — the remaining clippy warnings are the pre-existing doctestfn mainones.One behavior note to flag:
C/P/factorialpanicking on true overflow (instead of returning wrapped values in release) is a deliberate change; if you'd rather have a different convention there (saturating,Option), happy to adjust.