Repository navigation
Backport #221 to release-1: avoid Int overflow in signrank and wilcox - #230
Merged
Merged
Conversation
Member
That seems wrong? Neither is the version bumped on release-1 (AFAICT) nor does this PR target 1.5.2 (it bumps the version to 1.5.3). How did you create the backport? I think using |
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)
`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)
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)
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)
andreasnoack
force-pushed
the
backport-pr221-release-1
branch
from
September 4, 2026 07:26
116c810 to
85d3789
Compare
andreasnoack
changed the title
Member
Author
|
Done. I've updated the PR to use the |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## release-1 #230 +/- ##
=============================================
+ Coverage 69.71% 69.90% +0.18%
=============================================
Files 16 16
Lines 799 804 +5
=============================================
+ Hits 557 562 +5
Misses 242 242 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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.
Backports #221 ("Avoid Int overflow in signrank and wilcox by normalising in place") onto the
release-1branch, for av1.5.3release.Commits
The four commits of #221 are cherry-picked with
git cherry-pick -x, so each carries a(cherry picked from commit …)provenance line:f511a7c← 4a538b2 — Avoid Int overflow in signrank and wilcox by normalising in placed2d5b5a← 7776c46 — Drop the exactbinomialbranch from the wilcox normaliser431e67a← 9910d20 — Work in units of a power of two in the wilcox recurrence9f1bcf5← 64d880c — Address review comments on code commentsPlus one branch-specific commit:
85d3789— Bump version from 1.5.2 to 1.5.3Conflicts resolved during the cherry-pick
src/distrs/signrank.jl,src/distrs/wilcox.jl— applied cleanly; identical to Avoid Int overflow in signrank and wilcox by normalising in place #221's final state.src/StatsFuns.jl— the SpecialFunctions import line differs from master onrelease-1(noerf), sologabsbinomialwas merged into this branch's line by hand.test/rmath.jl—release-1still uses theRFunctions-based test harness (master'stest/rmath.jlwas later rewritten to callRmathdirectly). The newsignrank/wilcoxtestsets anchor on lines that exist on this branch and apply cleanly; ausing Rmath: Rmathimport was added for theRmath.*reference calls, and Avoid Int overflow in signrank and wilcox by normalising in place #221'stest_invchange (dropping theSys.islinux()special case) was reapplied to this branch's olderrmathcompstructure.test/misc.jl— the transientinv_binomialtestset is added by the first commit and removed by the second, netting to no change, as upstream.Fixes
signrankcdfreturns silently invalid probabilities from n = 72 (Intoverflow insignrankDP) #219 —signrankcdfreturned silently invalid probabilities fromn = 72wilcoxcdfreturns silently invalid probabilities for large samples (Intoverflow inwilcox_partitions) #220 —wilcoxcdfreturned silently invalid probabilities, and threwOverflowErrorfromnx + ny = 68Verified locally that the previously broken cases now agree with the Rmath reference: signrank at
n = 72/200, wilcox at(34,34)/(30,36), and the rescaled seed at(524,556). See #221 for the full analysis.