From f511a7c4020fcc18a5ebdd3ada6f59a951c23c6f Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 28 Jul 2026 10:36:52 +0200 Subject: [PATCH 1/6] 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 4a538b2426ad47acb71ebf8166422c49b4865e4f) --- src/StatsFuns.jl | 2 +- src/distrs/signrank.jl | 26 +++++++++++++++++------- src/distrs/wilcox.jl | 46 ++++++++++++++++++++++++++++-------------- test/misc.jl | 16 +++++++++++++++ test/rmath.jl | 41 +++++++++++++++++++++++++++++++++++++ 5 files changed, 108 insertions(+), 23 deletions(-) diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl index 22fdbc0..c57d411 100644 --- a/src/StatsFuns.jl +++ b/src/StatsFuns.jl @@ -3,7 +3,7 @@ module StatsFuns using Base: Math.@horner using Reexport: @reexport using SpecialFunctions: beta_inc, beta_inc_inv, digamma, - erfc, erfcinv, erfcx, gamma_inc, gamma_inc_inv, logbeta, loggamma + erfc, erfcinv, erfcx, gamma_inc, gamma_inc_inv, logabsbinomial, logbeta, loggamma # reexports @reexport using IrrationalConstants: diff --git a/src/distrs/signrank.jl b/src/distrs/signrank.jl index 96936db..3543aab 100644 --- a/src/distrs/signrank.jl +++ b/src/distrs/signrank.jl @@ -11,15 +11,27 @@ the number of subsets of {1,2,...,n-1} summing to W. This can be calculated bottom up using dynamic programming. The i'th element of DP in the j'th outer loop iteration represents: -the number of ways {1,2,...,j} can sum to W-i+1. +the probability that a subset of {1,2,...,j} sums to W-i+1, +i.e. the number of such subsets divided by the 2^j subsets in total. + +The counts themselves grow like 2^n and overflow `Int` well before the +distribution stops being of interest, so the recursion is normalised as it goes +rather than divided by 2^n at the end. Halving is exact in binary floating +point, so this is bit-for-bit what dividing the counts by 2^n would give, but +the entries stay in [0,1] and can never overflow. =# -@inline function signrankDP(n, W) - DP = zeros(Int, W + 1) - DP[W + 1] = 1 +function signrankDP(n, W) + DP = zeros(Float64, W + 1) + DP[W + 1] = 1.0 for j in 1:n for i in 1:(W + 1 - j) - DP[i] += DP[i + j] + DP[i] = (DP[i] + DP[i + j]) / 2 + end + # Sums smaller than j cannot contain j, so their counts are unchanged and + # the corresponding probabilities only pick up the factor 1/2 + for i in max(1, W + 2 - j):(W + 1) + DP[i] /= 2 end end return DP @@ -38,7 +50,7 @@ function signrankpdf(n::Int, W::Int) return signrankpdf(n, W2) end DP = signrankDP(n, W) - return ldexp(float(DP[1]), -n) + return DP[1] end function signranklogpdf(n::Int, W::Union{Float64, Int}) @@ -58,7 +70,7 @@ function signrankcdf(n::Int, W::Int) return 1.0 - signrankcdf(n, W2) end DP = signrankDP(n, W) - return sum(Base.Fix2(ldexp, -n) ∘ float, DP) + return sum(DP) end function signranklogcdf(n::Int, W::Union{Float64, Int}) diff --git a/src/distrs/wilcox.jl b/src/distrs/wilcox.jl index 64511cb..7bc898d 100644 --- a/src/distrs/wilcox.jl +++ b/src/distrs/wilcox.jl @@ -34,14 +34,30 @@ H. B. Mann, D. R. Whitney. "On a Test of Whether one of Two Random Variables is A. Löffler: "Über eine Partition der nat. Zahlen und ihre Anwendung beim U-Test." Wissenschaftliche Zeitschrift der Martin-Luther-Universität Halle-Wittenberg; Mathematisch-Naturwissenschaftliche Reihe, XXXII'83 M, Heft 5, 87–89; available as https://upload.wikimedia.org/wikipedia/commons/f/f5/LoefflerWilcoxonMannWhitneyTest.pdf =# -@inline function wilcox_partitions(nx::Int, ny::Int, U::Int) +# `1 / binomial(n, k)`, without ever forming a coefficient that overflows `Int`. +# `binomial(n, k)` is exact and converts to `Float64` exactly while it stays below `2^53`, +# so it is used whenever that is guaranteed (`exp(36.7) < 2^53`). Above that we fall back to +# `logabsbinomial`, which evaluates the coefficient as a beta function, +# `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`, and hence never overflows. +function inv_binomial(n::Int, k::Int) + logbinom = first(logabsbinomial(n, k)) + return logbinom < 36.7 ? inv(Float64(binomial(n, k))) : exp(-logbinom) +end + +#= +The recurrence below is linear and homogeneous in `pₘ,ₙ`, so seeding it with +`1 / binomial(nx + ny, nx)` instead of `1` yields the probabilities directly. +This matters because the counts themselves grow like `binomial(nx + ny, nx)` and overflow +`Int` from `nx + ny = 68` onwards, whereas the probabilities are bounded by 1. +=# +function wilcox_probabilities(nx::Int, ny::Int, U::Int) # This internal function expects 0 <= U <= nx * ny / 2 if !(0 <= U <= (nx * ny) / 2) - throw(ArgumentError("`wilcox_partitions(nx, ny, U)` is only implemented for 0 <= U <= (nx * ny) / 2")) + throw(ArgumentError("`wilcox_probabilities(nx, ny, U)` is only implemented for 0 <= U <= (nx * ny) / 2")) end - # Due to symmetry, `wilcox_partitions(nx, ny, U) = wilcox_partitions(ny, nx, U)` - # Hence for simplicity we only consider the case `wilcox_partitions(min(nx, ny), max(nx, ny), U)` here + # Due to symmetry, `wilcox_probabilities(nx, ny, U) = wilcox_probabilities(ny, nx, U)` + # Hence for simplicity we only consider the case `wilcox_probabilities(min(nx, ny), max(nx, ny), U)` here m, n = minmax(nx, ny) # Compute σ(k) = ∑_{d|k} ϵ(d) d where @@ -60,18 +76,18 @@ A. Löffler: "Über eine Partition der nat. Zahlen und ihre Anwendung beim U-Tes end end - # Recursively compute the number of partitions pₘ,ₙ(a) for 0 <= a <= U - partitions = Vector{Int}(undef, U + 1) - partitions[1] = 1 + # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U + probabilities = Vector{Float64}(undef, U + 1) + probabilities[1] = inv_binomial(nx + ny, nx) for a in 1:U - p = 0 + p = 0.0 for i in 1:a - p += partitions[i] * sigmas[a + 1 - i] + p += probabilities[i] * sigmas[a + 1 - i] end - partitions[a + 1] = p ÷ a + probabilities[a + 1] = p / a end - return partitions + return probabilities end function wilcoxpdf(nx::Int, ny::Int, U::Float64) @@ -83,8 +99,8 @@ function wilcoxpdf(nx::Int, ny::Int, U::Int) return 0.0 end U = min(U, max_U - U) - partitions = wilcox_partitions(nx, ny, U) - return partitions[end] / binomial(nx + ny, nx) + probabilities = wilcox_probabilities(nx, ny, U) + return probabilities[end] end function wilcoxlogpdf(nx::Int, ny::Int, U::Union{Float64, Int}) @@ -102,8 +118,8 @@ function wilcoxcdf(nx::Int, ny::Int, U::Int) return 1.0 else U2 = max_U - U - 1 - partitions = wilcox_partitions(nx, ny, min(U, U2)) - p = sum(float, partitions) / binomial(nx + ny, nx) + probabilities = wilcox_probabilities(nx, ny, min(U, U2)) + p = sum(probabilities) return U2 < U ? 1.0 - p : p end end diff --git a/test/misc.jl b/test/misc.jl index 40fa1c0..8261fb1 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -72,6 +72,22 @@ end end end +@testset "inv_binomial" begin + @testset "n = $n, k = $k" for (n, k) in ((1, 0), (4, 2), (22, 10), (52, 26), (56, 28), (58, 29), (68, 34), (200, 100), (1000, 3), (2000, 1000)) + @test StatsFuns.inv_binomial(n, k) ≈ Float64(1 / binomial(big(n), big(k))) + end + + # The fast branch computes `binomial(n, k)` in `Int` and converts it to `Float64`, + # which is only exact while the coefficient stays below 2^53 + @testset "the exact branch is only taken when it is exact" begin + @test all( + binomial(big(n), big(k)) <= 2^53 + for n in 1:200 for k in 0:n + if first(SpecialFunctions.logabsbinomial(n, k)) < 36.7 + ) + end +end + @testset "lstirling_asym" begin # can test for equality here because the lhs is the way the value is created @test Float32(lstirling_asym(1.0)) == @inferred lstirling_asym(1.0f0) diff --git a/test/rmath.jl b/test/rmath.jl index 373e34b..6b68dac 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -1,5 +1,6 @@ using StatsFuns using StatsFuns: RFunctions +using Rmath: Rmath using Test include("utils.jl") @@ -455,6 +456,26 @@ end @test isnan(signrankinvlogccdf.(50, signranklogccdf.(50, 1275))) @test isnan(signrankinvlogccdf.(50, signranklogccdf.(50, 1276))) + # The subset counts grow like 2^n, so accumulating them in `Int` silently wrapped + # around and returned invalid probabilities from n = 72 onwards (#219) + @testset "signrank: large n" begin + @testset "n = $n" for n in (71, 72, 73, 74, 80, 100, 120, 200) + # The distribution is symmetric about n * (n + 1) / 4 + W = (n * (n + 1)) ÷ 4 + @test signrankpdf(n, W) ≈ Rmath.dsignrank(W, n, false) + @test signrankcdf(n, W) ≈ Rmath.psignrank(W, n, true, false) + @test signrankccdf(n, W) ≈ Rmath.psignrank(W, n, false, false) + # Deep in the left tail, where the probabilities are tiny but still normal + W2 = W ÷ 3 + @test signrankpdf(n, W2) ≈ Rmath.dsignrank(W2, n, false) + @test signrankcdf(n, W2) ≈ Rmath.psignrank(W2, n, true, false) + end + + @testset "pdf sums to one, n = $n" for n in (70, 80, 100) + @test sum(signrankpdf(n, W) for W in 0:((n * (n + 1)) ÷ 2)) ≈ 1 + end + end + rmathcomp_tests( "srdist", [ ((1, 2), (0.0:0.2:5.0)), @@ -488,6 +509,26 @@ end @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 100))) @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 101))) + # The partition counts and their normaliser `binomial(nx + ny, nx)` exceed `typemax(Int)` + # from nx + ny = 68 onwards, which used to throw an `OverflowError` (#219) + @testset "wilcox: large nx and ny" begin + @testset "nx = $nx, ny = $ny" for (nx, ny) in ((33, 33), (34, 34), (40, 40), (50, 50), (5, 200), (150, 40)) + # The distribution is symmetric about nx * ny / 2 + U = (nx * ny) ÷ 2 + @test wilcoxpdf(nx, ny, U) ≈ Rmath.dwilcox(U, nx, ny, false) + @test wilcoxcdf(nx, ny, U) ≈ Rmath.pwilcox(U, nx, ny, true, false) + @test wilcoxccdf(nx, ny, U) ≈ Rmath.pwilcox(U, nx, ny, false, false) + # Deep in the left tail, where the probabilities are tiny but still normal + U2 = U ÷ 3 + @test wilcoxpdf(nx, ny, U2) ≈ Rmath.dwilcox(U2, nx, ny, false) + @test wilcoxcdf(nx, ny, U2) ≈ Rmath.pwilcox(U2, nx, ny, true, false) + end + + @testset "pdf sums to one, nx = $nx, ny = $ny" for (nx, ny) in ((34, 34), (50, 50), (5, 200)) + @test sum(wilcoxpdf(nx, ny, U) for U in 0:(nx * ny)) ≈ 1 + end + end + # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5. # Reduced df or k allows convergence. This test documents this behavior. x = 0.15 From d2d5b5ae86847a6a6b6de96a2e6078608b047a6a Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 28 Jul 2026 16:47:50 +0200 Subject: [PATCH 2/6] 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 7776c46fa0e52099301331d27f3e9bb99d1283ff) --- src/distrs/wilcox.jl | 18 ++++++------------ test/misc.jl | 16 ---------------- test/rmath.jl | 32 ++++++++++++++++++++++---------- 3 files changed, 28 insertions(+), 38 deletions(-) diff --git a/src/distrs/wilcox.jl b/src/distrs/wilcox.jl index 7bc898d..3f8e4da 100644 --- a/src/distrs/wilcox.jl +++ b/src/distrs/wilcox.jl @@ -34,16 +34,6 @@ H. B. Mann, D. R. Whitney. "On a Test of Whether one of Two Random Variables is A. Löffler: "Über eine Partition der nat. Zahlen und ihre Anwendung beim U-Test." Wissenschaftliche Zeitschrift der Martin-Luther-Universität Halle-Wittenberg; Mathematisch-Naturwissenschaftliche Reihe, XXXII'83 M, Heft 5, 87–89; available as https://upload.wikimedia.org/wikipedia/commons/f/f5/LoefflerWilcoxonMannWhitneyTest.pdf =# -# `1 / binomial(n, k)`, without ever forming a coefficient that overflows `Int`. -# `binomial(n, k)` is exact and converts to `Float64` exactly while it stays below `2^53`, -# so it is used whenever that is guaranteed (`exp(36.7) < 2^53`). Above that we fall back to -# `logabsbinomial`, which evaluates the coefficient as a beta function, -# `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`, and hence never overflows. -function inv_binomial(n::Int, k::Int) - logbinom = first(logabsbinomial(n, k)) - return logbinom < 36.7 ? inv(Float64(binomial(n, k))) : exp(-logbinom) -end - #= The recurrence below is linear and homogeneous in `pₘ,ₙ`, so seeding it with `1 / binomial(nx + ny, nx)` instead of `1` yields the probabilities directly. @@ -76,9 +66,13 @@ function wilcox_probabilities(nx::Int, ny::Int, U::Int) end end - # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U + # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U. + # The seed is evaluated with `logabsbinomial` since `binomial(nx + ny, nx)` itself + # overflows `Int` from `nx + ny = 68` onwards. `logabsbinomial` evaluates the coefficient + # as a beta function, `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`, + # and hence cannot overflow. probabilities = Vector{Float64}(undef, U + 1) - probabilities[1] = inv_binomial(nx + ny, nx) + probabilities[1] = exp(-first(logabsbinomial(nx + ny, nx))) for a in 1:U p = 0.0 for i in 1:a diff --git a/test/misc.jl b/test/misc.jl index 8261fb1..40fa1c0 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -72,22 +72,6 @@ end end end -@testset "inv_binomial" begin - @testset "n = $n, k = $k" for (n, k) in ((1, 0), (4, 2), (22, 10), (52, 26), (56, 28), (58, 29), (68, 34), (200, 100), (1000, 3), (2000, 1000)) - @test StatsFuns.inv_binomial(n, k) ≈ Float64(1 / binomial(big(n), big(k))) - end - - # The fast branch computes `binomial(n, k)` in `Int` and converts it to `Float64`, - # which is only exact while the coefficient stays below 2^53 - @testset "the exact branch is only taken when it is exact" begin - @test all( - binomial(big(n), big(k)) <= 2^53 - for n in 1:200 for k in 0:n - if first(SpecialFunctions.logabsbinomial(n, k)) < 36.7 - ) - end -end - @testset "lstirling_asym" begin # can test for equality here because the lhs is the way the value is created @test Float32(lstirling_asym(1.0)) == @inferred lstirling_asym(1.0f0) diff --git a/test/rmath.jl b/test/rmath.jl index 6b68dac..f28be21 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -28,17 +28,20 @@ function rmathcomp(basename, params, X::AbstractArray, rtol = _default_rtol(para end #= - R version of signrank has system variation + signrank and wilcox are implemented natively rather than by delegating to Rmath, so + their inverse functions are not expected to reproduce the Rmath quantile functions + bit for bit. The inverses locate a discontinuity by searching for the first argument + at which the cdf reaches `q`, and the `q` tested here are themselves cdf values, so + the comparison is decided by the last bit of the cdf and is not stable across + platforms. R's own signrank already varies this way: julia> psignrank(18,10,false,true) # windows -0.2076393647782445 julia> psignrank(18,10,false,true) # linux -0.20763936477824452 - This slight difference causes test failures for the inverse functions, - due to a slight shift in the location of the discontinuity. - - This also holds true for wilcox. + As for the other natively implemented discrete distributions, the inverse functions + are instead tested by round tripping through our own cdf further down. =# - test_inv = (basename != "signrank" && basename != "wilcox") || !Sys.islinux() + test_inv = basename != "signrank" && basename != "wilcox" if has_pdf pdf = string(basename, "pdf") @@ -509,10 +512,19 @@ end @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 100))) @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 101))) - # The partition counts and their normaliser `binomial(nx + ny, nx)` exceed `typemax(Int)` - # from nx + ny = 68 onwards, which used to throw an `OverflowError` (#219) + # Accumulating the partition counts in `Int` used to break down in two separate ways: + # the counts themselves wrapped around and returned silently invalid probabilities (#220), + # and their normaliser `binomial(nx + ny, nx)` threw an `OverflowError` from + # nx + ny = 68 onwards (#219). The silent limit is the lower of the two, and it is set by + # `binomial(nx + ny, nx)` rather than by `nx + ny`: the Löffler recurrence subtracts as + # well as adds, so its intermediates overshoot the counts and wrap first. Unbalanced + # groups therefore stayed correct much further out, which is why the cases below are not + # all balanced. @testset "wilcox: large nx and ny" begin - @testset "nx = $nx, ny = $ny" for (nx, ny) in ((33, 33), (34, 34), (40, 40), (50, 50), (5, 200), (150, 40)) + @testset "nx = $nx, ny = $ny" for (nx, ny) in ( + (32, 32), (33, 33), (34, 34), (40, 40), (50, 50), + (30, 36), (20, 48), (20, 60), (10, 325), (10, 326), (5, 200), (150, 40), + ) # The distribution is symmetric about nx * ny / 2 U = (nx * ny) ÷ 2 @test wilcoxpdf(nx, ny, U) ≈ Rmath.dwilcox(U, nx, ny, false) @@ -524,7 +536,7 @@ end @test wilcoxcdf(nx, ny, U2) ≈ Rmath.pwilcox(U2, nx, ny, true, false) end - @testset "pdf sums to one, nx = $nx, ny = $ny" for (nx, ny) in ((34, 34), (50, 50), (5, 200)) + @testset "pdf sums to one, nx = $nx, ny = $ny" for (nx, ny) in ((30, 30), (33, 33), (34, 34), (50, 50), (5, 200)) @test sum(wilcoxpdf(nx, ny, U) for U in 0:(nx * ny)) ≈ 1 end end From 431e67a57a8d3ceae4f10fba16b86585b156ed9e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 28 Jul 2026 21:41:27 +0200 Subject: [PATCH 3/6] 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 9910d20a9cd9c4658521c028d3913b2a8695bacc) --- src/distrs/wilcox.jl | 31 ++++++++++++++++++++----------- test/rmath.jl | 19 +++++++++++++++++++ 2 files changed, 39 insertions(+), 11 deletions(-) diff --git a/src/distrs/wilcox.jl b/src/distrs/wilcox.jl index 3f8e4da..85a7925 100644 --- a/src/distrs/wilcox.jl +++ b/src/distrs/wilcox.jl @@ -66,13 +66,22 @@ function wilcox_probabilities(nx::Int, ny::Int, U::Int) end end - # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U. - # The seed is evaluated with `logabsbinomial` since `binomial(nx + ny, nx)` itself - # overflows `Int` from `nx + ny = 68` onwards. `logabsbinomial` evaluates the coefficient - # as a beta function, `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`, - # and hence cannot overflow. + # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U, in units of + # 2^shift. The coefficient is evaluated with `logabsbinomial`, since `binomial(nx + ny, + # nx)` itself overflows `Int` from `nx + ny = 68` onwards, whereas `logabsbinomial` + # evaluates it as a beta function, `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`. + # + # The seed 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`, and since the recurrence is homogeneous a seed of zero + # would silently zero every probability. Working in units of 2^shift avoids that. The + # recurrence is homogeneous so the choice of unit cannot affect the result, `ldexp` + # removes it from the scalar result exactly, and `shift` is zero whenever the seed is + # normal, which leaves every input that already worked bit for bit unchanged. + logbinom = first(logabsbinomial(nx + ny, nx)) + shift = max(0, ceil(Int, (logbinom - 700) / logtwo)) probabilities = Vector{Float64}(undef, U + 1) - probabilities[1] = exp(-first(logabsbinomial(nx + ny, nx))) + probabilities[1] = exp(shift * logtwo - logbinom) for a in 1:U p = 0.0 for i in 1:a @@ -81,7 +90,7 @@ function wilcox_probabilities(nx::Int, ny::Int, U::Int) probabilities[a + 1] = p / a end - return probabilities + return probabilities, shift end function wilcoxpdf(nx::Int, ny::Int, U::Float64) @@ -93,8 +102,8 @@ function wilcoxpdf(nx::Int, ny::Int, U::Int) return 0.0 end U = min(U, max_U - U) - probabilities = wilcox_probabilities(nx, ny, U) - return probabilities[end] + probabilities, shift = wilcox_probabilities(nx, ny, U) + return ldexp(probabilities[end], -shift) end function wilcoxlogpdf(nx::Int, ny::Int, U::Union{Float64, Int}) @@ -112,8 +121,8 @@ function wilcoxcdf(nx::Int, ny::Int, U::Int) return 1.0 else U2 = max_U - U - 1 - probabilities = wilcox_probabilities(nx, ny, min(U, U2)) - p = sum(probabilities) + probabilities, shift = wilcox_probabilities(nx, ny, min(U, U2)) + p = ldexp(sum(probabilities), -shift) return U2 < U ? 1.0 - p : p end end diff --git a/test/rmath.jl b/test/rmath.jl index f28be21..4ac1c74 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -539,6 +539,25 @@ end @testset "pdf sums to one, nx = $nx, ny = $ny" for (nx, ny) in ((30, 30), (33, 33), (34, 34), (50, 50), (5, 200)) @test sum(wilcoxpdf(nx, ny, U) for U in 0:(nx * ny)) ≈ 1 end + + # `1 / binomial(nx + ny, nx)` is the probability of the least likely outcome and + # leaves the normal range from nx + ny = 1030. It seeds the recurrence, which is + # homogeneous, so once it reaches zero every probability follows it and the results + # are silently zero rather than merely imprecise. The recurrence therefore works in + # units of a power of two, which these check has not disturbed anything. + @testset "seed below the normal range, nx = ny = $m" for m in (515, 520, 530, 535) + # the pdf at U = 0 is exactly the seed, so this pins the scaling down directly + @test wilcoxpdf(m, m, 0) ≈ Float64(1 / BigFloat(binomial(big(2m), big(m)), precision = 512)) + end + + @testset "seed at zero without scaling" begin + # `1 / binomial(1080, 524)` is zero in `Float64`, so every one of these was + # exactly 0.0 before the recurrence was rescaled + @test 0 < wilcoxpdf(524, 556, 20000) < 1 + @test 0 < wilcoxcdf(524, 556, 20000) < 1 + @test wilcoxpdf(524, 556, 20000) <= wilcoxcdf(524, 556, 20000) + @test wilcoxcdf(524, 556, 20000) < wilcoxcdf(524, 556, 20001) + end end # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5. From 9f1bcf563392a349fc7517c84db4739a22973a89 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 29 Jul 2026 10:30:38 +0200 Subject: [PATCH 4/6] Address review comments on code comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 64d880c6e854ba33794051ec13068af854bd3bad) --- src/distrs/signrank.jl | 5 +++-- src/distrs/wilcox.jl | 7 +++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/distrs/signrank.jl b/src/distrs/signrank.jl index 3543aab..158c530 100644 --- a/src/distrs/signrank.jl +++ b/src/distrs/signrank.jl @@ -17,8 +17,9 @@ i.e. the number of such subsets divided by the 2^j subsets in total. The counts themselves grow like 2^n and overflow `Int` well before the distribution stops being of interest, so the recursion is normalised as it goes rather than divided by 2^n at the end. Halving is exact in binary floating -point, so this is bit-for-bit what dividing the counts by 2^n would give, but -the entries stay in [0,1] and can never overflow. +point, so while every count stays below 2^53 this gives exactly the counts +divided by 2^n; beyond that the additions round, but the entries stay in [0,1] +and can never overflow. =# function signrankDP(n, W) diff --git a/src/distrs/wilcox.jl b/src/distrs/wilcox.jl index 85a7925..9ebc033 100644 --- a/src/distrs/wilcox.jl +++ b/src/distrs/wilcox.jl @@ -67,9 +67,8 @@ function wilcox_probabilities(nx::Int, ny::Int, U::Int) end # Recursively compute pₘ,ₙ(a) / binomial(nx + ny, nx) for 0 <= a <= U, in units of - # 2^shift. The coefficient is evaluated with `logabsbinomial`, since `binomial(nx + ny, - # nx)` itself overflows `Int` from `nx + ny = 68` onwards, whereas `logabsbinomial` - # evaluates it as a beta function, `1 / binomial(n, k) = (n + 1) * beta(k + 1, n - k + 1)`. + # 2^shift. The logarithm of the coefficient is used since `binomial(nx + ny, nx)` itself + # overflows `Int` from `nx + ny = 68` onwards. # # The seed 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 @@ -77,7 +76,7 @@ function wilcox_probabilities(nx::Int, ny::Int, U::Int) # would silently zero every probability. Working in units of 2^shift avoids that. The # recurrence is homogeneous so the choice of unit cannot affect the result, `ldexp` # removes it from the scalar result exactly, and `shift` is zero whenever the seed is - # normal, which leaves every input that already worked bit for bit unchanged. + # normal. logbinom = first(logabsbinomial(nx + ny, nx)) shift = max(0, ceil(Int, (logbinom - 700) / logtwo)) probabilities = Vector{Float64}(undef, U + 1) From 85d3789515bb7f53aaba310f5d1b623247826c73 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Fri, 4 Sep 2026 09:26:06 +0200 Subject: [PATCH 5/6] Bump version from 1.5.2 to 1.5.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4e86690..865ee0f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "StatsFuns" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "1.5.2" +version = "1.5.3" [deps] HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" From 3b4cdf0e15b5e2bb953b3b0c33e278b061bfac8f Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Fri, 4 Sep 2026 09:29:27 +0200 Subject: [PATCH 6/6] Run CI on release branches and PRs targeting them --- .github/workflows/ci.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8b3fd62..f30f4e4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -4,9 +4,11 @@ on: pull_request: branches: - master + - 'release-*' push: branches: - master + - 'release-*' tags: '*' concurrency: