Implement binomial, Poisson, and negative binomial inverse functions in pure Julia - #205
Draft
andreasnoack wants to merge 1 commit into
Draft
Implement binomial, Poisson, and negative binomial inverse functions in pure Julia#205andreasnoack wants to merge 1 commit into
andreasnoack wants to merge 1 commit into
Conversation
andreasnoack
force-pushed
the
an/normath-discrete
branch
from
March 22, 2026 13:22
cbc19d7 to
4b332a2
Compare
…in pure Julia Replace Rmath-dependent inverse functions for three discrete distributions with pure Julia implementations based on VBA code by Ian Smith: - Binomial: binominvcdf, binominvccdf, binominvlogcdf, binominvlogccdf - Poisson: poisinvcdf, poisinvccdf, poisinvlogcdf, poisinvlogccdf - Negative binomial: nbinominvcdf, nbinominvccdf, nbinominvlogcdf, nbinominvlogccdf All use normal approximation for initial guess followed by linear search with direct PMF evaluation. Each search function uses its own CDF variant (cdf or ccdf) directly, avoiding precision loss from 1-p subtraction. The invlog variants use -expm1(lq) when lq is close to 0. Tests use self-consistency round-trips for inverse functions, since our CDF implementation differs from Rmath's by ~1 ULP which can flip discrete answers.
andreasnoack
force-pushed
the
an/normath-discrete
branch
from
March 22, 2026 13:37
4b332a2 to
df74a25
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #205 +/- ##
==========================================
+ Coverage 72.39% 76.04% +3.65%
==========================================
Files 22 22
Lines 1007 1269 +262
==========================================
+ Hits 729 965 +236
- Misses 278 304 +26 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
andreasnoack
marked this pull request as draft
March 23, 2026 08:51
andreasnoack
added a commit
that referenced
this pull request
Jul 28, 2026
`inv_binomial` evaluated both `logabsbinomial` and `binomial` whenever the coefficient was small enough for the latter, and only used one of them. The exactly rounded reciprocal it bought is worth a few ulp in the cdf, which does not justify the branch, so the normaliser is now always evaluated as a beta function via `logabsbinomial`. The one thing that did depend on those last bits was the comparison of `wilcoxinvcdf` against `Rmath.qwilcox`. The inverses locate a discontinuity by searching for the first argument at which the cdf reaches `q`, and the `q` tested are themselves cdf values, so the outcome is decided by the last bit of the cdf; nudging the normaliser by a single ulp moves the answer in either direction. That is not a property to design around, and it is not what the package does elsewhere: as of #205 the natively implemented discrete distributions test their inverses by round tripping through their own cdf rather than against the Rmath quantile functions. `signrank` and `wilcox` are the only other native discrete implementations, they already have such round trips, so they are excluded from the Rmath inverse comparison outright. That also removes the `Sys.islinux()` special case, which was working around the same instability on one platform. Also drop the test that checked `exp(36.7) < 2^53` by brute force, which only restated the threshold it was given, and cover the asymmetric cases from #220, where the counts overflowed well before `binomial(nx + ny, nx)` did.
andreasnoack
added a commit
that referenced
this pull request
Sep 4, 2026
`inv_binomial` evaluated both `logabsbinomial` and `binomial` whenever the coefficient was small enough for the latter, and only used one of them. The exactly rounded reciprocal it bought is worth a few ulp in the cdf, which does not justify the branch, so the normaliser is now always evaluated as a beta function via `logabsbinomial`. The one thing that did depend on those last bits was the comparison of `wilcoxinvcdf` against `Rmath.qwilcox`. The inverses locate a discontinuity by searching for the first argument at which the cdf reaches `q`, and the `q` tested are themselves cdf values, so the outcome is decided by the last bit of the cdf; nudging the normaliser by a single ulp moves the answer in either direction. That is not a property to design around, and it is not what the package does elsewhere: as of #205 the natively implemented discrete distributions test their inverses by round tripping through their own cdf rather than against the Rmath quantile functions. `signrank` and `wilcox` are the only other native discrete implementations, they already have such round trips, so they are excluded from the Rmath inverse comparison outright. That also removes the `Sys.islinux()` special case, which was working around the same instability on one platform. Also drop the test that checked `exp(36.7) < 2^53` by brute force, which only restated the threshold it was given, and cover the asymmetric cases from #220, where the counts overflowed well before `binomial(nx + ny, nx)` did. (cherry picked from commit 7776c46)
andreasnoack
added a commit
that referenced
this pull request
Sep 4, 2026
…#230) * Avoid Int overflow in signrank and wilcox by normalising in place Both distributions accumulated combinatorial counts in `Int` and divided by the total count at the end. Those counts grow like `2^n` and `binomial(nx + ny, nx)`, so they wrapped around and made `signrankcdf` return silently invalid probabilities from `n = 72`, while `wilcoxcdf` threw an `OverflowError` from `nx + ny = 68`. Both recursions are linear and homogeneous, so the normaliser can be folded into the recursion instead of applied to its result, and the counts are never formed: - `signrankDP` halves at every step. Halving is exact in binary floating point and every count below `2^53` adds exactly, so the entries are bit-identical to the old ones divided by `2^n` for all `n <= 62`, i.e. across the whole range in which the old implementation was valid. - `wilcox_partitions` is seeded with `1 / binomial(nx + ny, nx)` rather than with `1`, and hence returns the probabilities themselves. It is renamed `wilcox_probabilities` accordingly. The seed uses `binomial` while the coefficient is below `2^53`, where it is exact, and `logabsbinomial`, which evaluates the coefficient as a beta function, above that. Agreement with Rmath across the support is now ~1e-15 up to at least `n = 200` and `nx + ny = 205`, and the pdfs sum to one again. Fixes #219 (cherry picked from commit 4a538b2) * Drop the exact `binomial` branch from the wilcox normaliser `inv_binomial` evaluated both `logabsbinomial` and `binomial` whenever the coefficient was small enough for the latter, and only used one of them. The exactly rounded reciprocal it bought is worth a few ulp in the cdf, which does not justify the branch, so the normaliser is now always evaluated as a beta function via `logabsbinomial`. The one thing that did depend on those last bits was the comparison of `wilcoxinvcdf` against `Rmath.qwilcox`. The inverses locate a discontinuity by searching for the first argument at which the cdf reaches `q`, and the `q` tested are themselves cdf values, so the outcome is decided by the last bit of the cdf; nudging the normaliser by a single ulp moves the answer in either direction. That is not a property to design around, and it is not what the package does elsewhere: as of #205 the natively implemented discrete distributions test their inverses by round tripping through their own cdf rather than against the Rmath quantile functions. `signrank` and `wilcox` are the only other native discrete implementations, they already have such round trips, so they are excluded from the Rmath inverse comparison outright. That also removes the `Sys.islinux()` special case, which was working around the same instability on one platform. Also drop the test that checked `exp(36.7) < 2^53` by brute force, which only restated the threshold it was given, and cover the asymmetric cases from #220, where the counts overflowed well before `binomial(nx + ny, nx)` did. (cherry picked from commit 7776c46) * Work in units of a power of two in the wilcox recurrence The seed `1 / binomial(nx + ny, nx)` is the probability of the least likely outcome, so it leaves the normal range long before the probabilities of interest do: it is subnormal from `nx + ny = 1030` and zero from `nx + ny = 1090`. Since the recurrence is homogeneous, a seed of zero zeroes every probability, so `wilcoxcdf` returned a plausible but wrong 0.776804 at `nx = ny = 540` and 0.0 from `nx = ny = 545`, where the answer is 0.5 either way by symmetry. The counts these replaced threw an `OverflowError` there, so this was a loud failure turning into a silent one. The recurrence is homogeneous, so the unit it works in cannot affect the result. It now seeds with `2^shift / binomial(nx + ny, nx)`, with `shift` chosen to place the seed in the normal range, and `ldexp` removes the unit from the scalar result. `shift` is zero whenever the seed is already normal, which leaves every input that worked before bit for bit unchanged. This moves the limit from `nx + ny = 1090` to about `nx + ny = 2050`, where the scaled values would overflow instead. The recurrence is O(U^2) in `U = nx * ny / 2`, so evaluating the cdf at the new limit takes minutes, and runtime rather than representability is now what stops you. Verified by shift invariance: recomputing with deliberately different units agrees to 1e-13, and the midpoints at `nx = ny = 540` and `560` now agree with the normal approximation to six digits. (cherry picked from commit 9910d20) * Address review comments on code comments Drop the claim that `logabsbinomial` evaluates the coefficient as a beta function — how it computes the logarithm is its own business — and the claim that `shift = 0` leaves results bit for bit unchanged, which compared against an intermediate state of this PR and means nothing in the merged code. Qualify the corresponding statement in the signrank comment: halving reproduces the counts divided by 2^n exactly only while every count stays below 2^53. (cherry picked from commit 64d880c) * Bump version from 1.5.2 to 1.5.3 * Run CI on release branches and PRs targeting them
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
2 participants
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.
Summary
Replace Rmath-dependent inverse functions for three discrete distributions with pure Julia implementations based on VBA code by Ian Smith:
binominvcdf,binominvccdf,binominvlogcdf,binominvlogccdfpoisinvcdf,poisinvccdf,poisinvlogcdf,poisinvlogccdfnbinominvcdf,nbinominvccdf,nbinominvlogcdf,nbinominvlogccdfAll use normal approximation for initial guess followed by linear search with direct PMF evaluation. Each search function uses its own CDF variant (cdf or ccdf) directly, avoiding precision loss from
1-psubtraction. Theinvlogvariants use-expm1(lq)whenlqis close to 0.Test plan
invcdf(cdf(k)) == kfor all support points with non-degenerate CDF valuesinvloground-trips allow ±1 tolerance due to inherentexp(log(cdf))precision loss🤖 Generated with Claude Code