From df74a2580a731ec465197d9f16d79d7239a3fe4c Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 22 Mar 2026 14:37:17 +0100 Subject: [PATCH] Implement binomial, Poisson, and negative binomial inverse functions 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. --- src/distrs/binom.jl | 148 ++++++++++++++++++++++++++++- src/distrs/nbinom.jl | 164 ++++++++++++++++++++++++++++++-- src/distrs/pois.jl | 153 +++++++++++++++++++++++++++++- test/rmath.jl | 219 ++++++++++++++++++++++++++++++++++++++++--- 4 files changed, 652 insertions(+), 32 deletions(-) diff --git a/src/distrs/binom.jl b/src/distrs/binom.jl index 7d16c48..78c12d7 100644 --- a/src/distrs/binom.jl +++ b/src/distrs/binom.jl @@ -30,20 +30,158 @@ for l in ("", "log"), compl in (false, true) end end -# Rmath implementations +# Inverse CDF: find smallest k such that binomcdf(n, p, k) >= cprob. +# Based on VBA critbinomial by Ian Smith. +# PMF ratios: PMF(k+1)/PMF(k) = (n-k)*p/((k+1)*(1-p)) +# PMF(k-1)/PMF(k) = k*(1-p)/((n-k+1)*p) +function _critbinomial(n::Float64, p::Float64, cprob::Float64) + q = 1.0 - p + + # Normal approximation for initial guess + σ = sqrt(n * p * q) + i = clamp(floor(n * p + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5), 0.0, n) + + # Compute CDF at the guess + pr = Float64(binomcdf(n, p, i)) + + if pr >= cprob + # Search left: find smallest k with CDF(k) >= cprob + while i > 0 + tpr = Float64(binompdf(n, p, i)) + if pr - tpr < cprob + return i # CDF(i-1) < cprob <= CDF(i) + end + pr -= tpr + i -= 1.0 + end + return 0.0 + else + # Search right + while i < n + i += 1.0 + tpr = Float64(binompdf(n, p, i)) + pr += tpr + if pr >= cprob + return i + end + end + return n + end +end + +# Inverse CCDF: find smallest k such that binomccdf(n, p, k) <= cprob. +# Based on VBA critcompbinomial by Ian Smith. +function _critcompbinomial(n::Float64, p::Float64, cprob::Float64) + q = 1.0 - p + + # Normal approximation + σ = sqrt(n * p * q) + i = clamp(floor(n * p - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5), 0.0, n) + + # Compute CCDF at the guess + pr = Float64(binomccdf(n, p, i)) + + if pr > cprob + # Search right + while i < n + i += 1.0 + tpr = Float64(binompdf(n, p, i)) + pr -= tpr + if pr <= cprob + return i + end + end + return n + else + # Search left: find smallest k with CCDF(k) <= cprob + while i > 0 + tpr = Float64(binompdf(n, p, i)) + if pr + tpr > cprob + return i # CCDF(i) <= cprob < CCDF(i-1) + end + pr += tpr + i -= 1.0 + end + return 0.0 + end +end + +# Wrapper with edge cases and post-correction (VBA crit_binomial) +function _binom_invcdf(n::Float64, p::Float64, q::Float64) + if q < 0 || q > 1 || n < 0 || p < 0 || p > 1 || isnan(q) || isnan(n) || isnan(p) + return NaN + elseif q == 0 || p == 0 + return 0.0 + elseif p == 1 + return n + end + + i = _critbinomial(n, p, q) + + # Post-correction + pr = Float64(binomcdf(n, p, i)) + if pr >= q + while i > 0 + pr2 = Float64(binomcdf(n, p, i - 1.0)) + if pr2 < q + return i + end + i -= 1.0 + pr = pr2 + end + return 0.0 + else + return i + 1.0 + end +end + +# Wrapper with edge cases and post-correction (VBA comp_crit_binomial) +function _binom_invccdf(n::Float64, p::Float64, q::Float64) + if q < 0 || q > 1 || n < 0 || p < 0 || p > 1 || isnan(q) || isnan(n) || isnan(p) + return NaN + elseif q == 1 || p == 0 + return 0.0 + elseif q == 0 || p == 1 + return n + end + + i = _critcompbinomial(n, p, q) + + # Post-correction + pr = Float64(binomccdf(n, p, i)) + if pr <= q + while i > 0 + pr2 = Float64(binomccdf(n, p, i - 1.0)) + if pr2 > q + return i + end + i -= 1.0 + pr = pr2 + end + return 0.0 + else + return i + 1.0 + end +end + +# Public API + function binominvcdf(n::Real, p::Real, q::Real) T = float(Base.promote_typeof(n, p, q)) - return convert(T, Rmath.qbinom(q, n, p, true, false)) + return convert(T, _binom_invcdf(Float64(n), Float64(p), Float64(q))) end + function binominvccdf(n::Real, p::Real, q::Real) T = float(Base.promote_typeof(n, p, q)) - return convert(T, Rmath.qbinom(q, n, p, false, false)) + return convert(T, _binom_invccdf(Float64(n), Float64(p), Float64(q))) end + function binominvlogcdf(n::Real, p::Real, lq::Real) T = float(Base.promote_typeof(n, p, lq)) - return convert(T, Rmath.qbinom(lq, n, p, true, true)) + return convert(T, _binom_invcdf(Float64(n), Float64(p), exp(Float64(lq)))) end + function binominvlogccdf(n::Real, p::Real, lq::Real) T = float(Base.promote_typeof(n, p, lq)) - return convert(T, Rmath.qbinom(lq, n, p, false, true)) + return convert(T, _binom_invccdf(Float64(n), Float64(p), exp(Float64(lq)))) end diff --git a/src/distrs/nbinom.jl b/src/distrs/nbinom.jl index 580e5db..343d6db 100644 --- a/src/distrs/nbinom.jl +++ b/src/distrs/nbinom.jl @@ -66,22 +66,172 @@ function nbinomlogccdf(r::Real, p::Real, k::Real) end end -# Rmath implementations -# TODO: implement https://arxiv.org/abs/2001.03953 -# for inverting the incomplete beta function wrt the 2nd argument. +# Inverse CDF: find smallest k such that nbinomcdf(r, p, k) >= cprob. +# Based on VBA critnegbinom by Ian Smith. +# PMF ratio: PMF(k+1)/PMF(k) = (k+r)/(k+1) * (1-p) +function _critnbinom(r::Float64, p::Float64, cprob::Float64) + q = 1.0 - p + + # Normal approximation: mean = r*q/p, var = r*q/p^2 + μ = r * q / p + σ = sqrt(μ / p) + i = max(0.0, floor(μ + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5)) + + # Compute CDF at the guess + pr = Float64(nbinomcdf(r, p, i)) + + if pr >= cprob + # Search left + while i > 0 + tpr = Float64(nbinompdf(r, p, i)) + if pr - tpr < cprob + return i + end + pr -= tpr + i -= 1.0 + end + return 0.0 + else + # Search right + for _ in 1:10_000 + i += 1.0 + tpr = Float64(nbinompdf(r, p, i)) + pr += tpr + if pr >= cprob + return i + end + end + return i + end +end + +# Inverse CCDF: find smallest k such that nbinomccdf(r, p, k) <= cprob. +# Based on VBA critcompnegbinom by Ian Smith. +function _critcompnbinom(r::Float64, p::Float64, cprob::Float64) + q = 1.0 - p + + # Normal approximation + μ = r * q / p + σ = sqrt(μ / p) + i = max(0.0, floor(μ - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5)) + + # Compute CCDF at the guess + pr = Float64(nbinomccdf(r, p, i)) + + if pr > cprob + # Search right + for _ in 1:10_000 + i += 1.0 + tpr = Float64(nbinompdf(r, p, i)) + pr -= tpr + if pr <= cprob + return i + end + end + return i + else + # Search left + while i > 0 + tpr = Float64(nbinompdf(r, p, i)) + if pr + tpr > cprob + return i + end + pr += tpr + i -= 1.0 + end + return 0.0 + end +end + +# Wrapper with edge cases +function _nbinom_invcdf(r::Float64, p::Float64, q::Float64) + if q < 0 || q > 1 || r <= 0 || p < 0 || p > 1 || isnan(q) || isnan(r) || isnan(p) + return NaN + elseif q == 0 || p == 1 + return 0.0 + elseif q == 1 + return Inf + elseif p == 0 + return Inf + end + + i = _critnbinom(r, p, q) + + # Post-correction + pr = Float64(nbinomcdf(r, p, i)) + if pr >= q + while i > 0 + pr2 = Float64(nbinomcdf(r, p, i - 1.0)) + if pr2 < q + return i + end + i -= 1.0 + end + return 0.0 + else + return i + 1.0 + end +end + +function _nbinom_invccdf(r::Float64, p::Float64, q::Float64) + if q < 0 || q > 1 || r <= 0 || p < 0 || p > 1 || isnan(q) || isnan(r) || isnan(p) + return NaN + elseif q == 0 || p == 0 + return Inf + elseif q == 1 + return 0.0 + elseif p == 1 + return 0.0 + end + + i = _critcompnbinom(r, p, q) + + # Post-correction + pr = Float64(nbinomccdf(r, p, i)) + if pr <= q + while i > 0 + pr2 = Float64(nbinomccdf(r, p, i - 1.0)) + if pr2 > q + return i + end + i -= 1.0 + end + return 0.0 + else + return i + 1.0 + end +end + +# Public API + function nbinominvcdf(r::Real, p::Real, q::Real) T = float(Base.promote_typeof(r, p, q)) - return convert(T, Rmath.qnbinom(q, r, p, true, false)) + return convert(T, _nbinom_invcdf(Float64(r), Float64(p), Float64(q))) end + function nbinominvccdf(r::Real, p::Real, q::Real) T = float(Base.promote_typeof(r, p, q)) - return convert(T, Rmath.qnbinom(q, r, p, false, false)) + return convert(T, _nbinom_invccdf(Float64(r), Float64(p), Float64(q))) end + function nbinominvlogcdf(r::Real, p::Real, lq::Real) T = float(Base.promote_typeof(r, p, lq)) - return convert(T, Rmath.qnbinom(lq, r, p, true, true)) + _lq = Float64(lq) + result = if _lq > -1 + _nbinom_invccdf(Float64(r), Float64(p), -expm1(_lq)) + else + _nbinom_invcdf(Float64(r), Float64(p), exp(_lq)) + end + return convert(T, result) end + function nbinominvlogccdf(r::Real, p::Real, lq::Real) T = float(Base.promote_typeof(r, p, lq)) - return convert(T, Rmath.qnbinom(lq, r, p, false, true)) + _lq = Float64(lq) + result = if _lq > -1 + _nbinom_invcdf(Float64(r), Float64(p), -expm1(_lq)) + else + _nbinom_invccdf(Float64(r), Float64(p), exp(_lq)) + end + return convert(T, result) end diff --git a/src/distrs/pois.jl b/src/distrs/pois.jl index 1c9b565..d65c187 100644 --- a/src/distrs/pois.jl +++ b/src/distrs/pois.jl @@ -18,20 +18,163 @@ poislogcdf(λ::Real, x::Real) = gammalogccdf(max(0, floor(x + 1)), 1, λ) poislogccdf(λ::Real, x::Real) = gammalogcdf(max(0, floor(x + 1)), 1, λ) -# Rmath implementations +# Inverse CDF: find smallest k such that poiscdf(λ, k) >= cprob. +# Based on VBA critpoiss by Ian Smith. +# PMF ratio: PMF(k+1)/PMF(k) = λ/(k+1), PMF(k-1)/PMF(k) = k/λ +function _critpoiss(λ::Float64, cprob::Float64) + # Normal approximation for initial guess + σ = sqrt(λ) + i = max(0.0, floor(λ + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5)) + + # Compute CDF at the guess + pr = Float64(poiscdf(λ, i)) + + if pr >= cprob + # Search left: find smallest k with CDF(k) >= cprob + while i > 0 + tpr = Float64(poispdf(λ, i)) + if pr - tpr < cprob + return i + end + pr -= tpr + i -= 1.0 + end + return 0.0 + else + # Search right (with guard against infinite loop) + for _ in 1:10_000 + i += 1.0 + tpr = Float64(poispdf(λ, i)) + pr += tpr + if pr >= cprob + return i + end + end + return i + end +end + +# Inverse CCDF: find smallest k such that poisccdf(λ, k) <= cprob. +# Based on VBA critcomppoiss by Ian Smith. +function _critcomppoiss(λ::Float64, cprob::Float64) + # Normal approximation + σ = sqrt(λ) + i = max(0.0, floor(λ - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5)) + + # Compute CCDF at the guess + pr = Float64(poisccdf(λ, i)) + + if pr > cprob + # Search right (with guard) + for _ in 1:10_000 + i += 1.0 + tpr = Float64(poispdf(λ, i)) + pr -= tpr + if pr <= cprob + return i + end + end + return i + else + # Search left + while i > 0 + tpr = Float64(poispdf(λ, i)) + if pr + tpr > cprob + return i + end + pr += tpr + i -= 1.0 + end + return 0.0 + end +end + +# Wrapper with edge cases +function _pois_invcdf(λ::Float64, q::Float64) + if q < 0 || q > 1 || λ < 0 || isnan(q) || isnan(λ) + return NaN + elseif q == 0 || λ == 0 + return 0.0 + elseif q == 1 + return Inf + end + + i = _critpoiss(λ, q) + + # Post-correction + pr = Float64(poiscdf(λ, i)) + if pr >= q + while i > 0 + pr2 = Float64(poiscdf(λ, i - 1.0)) + if pr2 < q + return i + end + i -= 1.0 + end + return 0.0 + else + return i + 1.0 + end +end + +function _pois_invccdf(λ::Float64, q::Float64) + if q < 0 || q > 1 || λ < 0 || isnan(q) || isnan(λ) + return NaN + elseif q == 0 + return Inf + elseif λ == 0 + return 0.0 + end + + i = _critcomppoiss(λ, q) + + # Post-correction + pr = Float64(poisccdf(λ, i)) + if pr <= q + while i > 0 + pr2 = Float64(poisccdf(λ, i - 1.0)) + if pr2 > q + return i + end + i -= 1.0 + end + return 0.0 + else + return i + 1.0 + end +end + +# Public API + function poisinvcdf(λ::Real, q::Real) T = float(Base.promote_typeof(λ, q)) - return convert(T, Rmath.qpois(q, λ, true, false)) + return convert(T, _pois_invcdf(Float64(λ), Float64(q))) end + function poisinvccdf(λ::Real, q::Real) T = float(Base.promote_typeof(λ, q)) - return convert(T, Rmath.qpois(q, λ, false, false)) + return convert(T, _pois_invccdf(Float64(λ), Float64(q))) end + function poisinvlogcdf(λ::Real, lq::Real) T = float(Base.promote_typeof(λ, lq)) - return convert(T, Rmath.qpois(lq, λ, true, true)) + _lq = Float64(lq) + # Use -expm1(lq) = 1-exp(lq) via invccdf for accuracy when lq ≈ 0 + result = if _lq > -1 + _pois_invccdf(Float64(λ), -expm1(_lq)) + else + _pois_invcdf(Float64(λ), exp(_lq)) + end + return convert(T, result) end + function poisinvlogccdf(λ::Real, lq::Real) T = float(Base.promote_typeof(λ, lq)) - return convert(T, Rmath.qpois(lq, λ, false, true)) + _lq = Float64(lq) + result = if _lq > -1 + _pois_invcdf(Float64(λ), -expm1(_lq)) + else + _pois_invccdf(Float64(λ), exp(_lq)) + end + return convert(T, result) end diff --git a/test/rmath.jl b/test/rmath.jl index 8d5bd60..9f8087e 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -219,8 +219,8 @@ end end end - rmathcomp_tests( - "binom", [ + @testset "binom" begin + configs = [ ((1, 0.5), 0.0:1.0), ((1, 0.7), 0.0:1.0), ((8, 0.6), 0.0:8.0), @@ -234,7 +234,52 @@ end ((10, 1 // 2), (0 // 1):(10 // 1)), ((1, 0.5), [-Inf, Inf]), ] - ) + @testset "params: $params" for (params, data) in configs + (n, p) = params + T = float(Base.promote_typeof(n, p, first(data))) + rtol = _default_rtol(T) + @testset "pdf with x=$x" for x in data + check_rmath(binompdf, (a, n, p) -> Rmath.dbinom(a, n, p, false), params, x, true, rtol) + end + @testset "logpdf with x=$x" for x in data + check_rmath(binomlogpdf, (a, n, p) -> Rmath.dbinom(a, n, p, true), params, x, false, rtol) + end + @testset "cdf with x=$x" for x in data + check_rmath(binomcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, false), params, x, true, rtol) + end + @testset "ccdf with x=$x" for x in data + check_rmath(binomccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, false), params, x, true, rtol) + end + @testset "logcdf with x=$x" for x in data + check_rmath(binomlogcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, true), params, x, false, rtol) + end + @testset "logccdf with x=$x" for x in data + check_rmath(binomlogccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, true), params, x, false, rtol) + end + if T === Float64 + @testset "invcdf round-trip with x=$x" for x in data + q = binomcdf(n, p, x) + (q == 0 || q == 1) && continue + @test binominvcdf(n, p, q) ≈ x rtol = rtol nans = true + end + @testset "invccdf round-trip with x=$x" for x in data + q = binomccdf(n, p, x) + (q == 0 || q == 1) && continue + @test binominvccdf(n, p, q) ≈ x rtol = rtol nans = true + end + @testset "invlogcdf round-trip with x=$x" for x in data + lq = binomlogcdf(n, p, x) + (lq == 0 || lq == -Inf) && continue + @test binominvlogcdf(n, p, lq) ≈ x atol = 1 rtol = rtol nans = true + end + @testset "invlogccdf round-trip with x=$x" for x in data + lq = binomlogccdf(n, p, x) + (lq == 0 || lq == -Inf) && continue + @test binominvlogccdf(n, p, lq) ≈ x atol = 1 rtol = rtol nans = true + end + end + end + end rmathcomp_tests( "chisq", [ @@ -307,8 +352,8 @@ end ] ) - rmathcomp_tests( - "nbinom", [ + @testset "nbinom" begin + configs = [ ((1, 0.5), 0.0:20.0), ((3, 0.5), 0.0:20.0), ((3, 0.2), 0.0:20.0), @@ -321,7 +366,60 @@ end ((1, 0.5), [-Inf, Inf, NaN]), ((1, 0.5f0), [-Inf32, Inf32, NaN32, -2, 1.1]), ] - ) + @testset "params: $params" for (params, data) in configs + (r, p) = params + T = float(Base.promote_typeof(r, p, first(data))) + rtol = _default_rtol(T) + @testset "pdf with x=$x" for x in data + check_rmath(nbinompdf, (a, r, p) -> Rmath.dnbinom(a, r, p, false), params, x, true, rtol) + end + @testset "logpdf with x=$x" for x in data + check_rmath(nbinomlogpdf, (a, r, p) -> Rmath.dnbinom(a, r, p, true), params, x, false, rtol) + end + @testset "cdf with x=$x" for x in data + check_rmath(nbinomcdf, (a, r, p) -> Rmath.pnbinom(a, r, p, true, false), params, x, true, rtol) + end + @testset "ccdf with x=$x" for x in data + check_rmath(nbinomccdf, (a, r, p) -> Rmath.pnbinom(a, r, p, false, false), params, x, true, rtol) + end + @testset "logcdf with x=$x" for x in data + check_rmath(nbinomlogcdf, (a, r, p) -> Rmath.pnbinom(a, r, p, true, true), params, x, false, rtol) + end + @testset "logccdf with x=$x" for x in data + check_rmath(nbinomlogccdf, (a, r, p) -> Rmath.pnbinom(a, r, p, false, true), params, x, false, rtol) + end + if T === Float64 + @testset "invcdf round-trip with x=$x" for x in data + (isnan(x) || isinf(x)) && continue + xf = floor(Float64(x)) + q = nbinomcdf(r, p, x) + (q == 0 || q == 1 || xf < 0) && continue + @test Float64(nbinominvcdf(r, p, q)) == xf + end + @testset "invccdf round-trip with x=$x" for x in data + (isnan(x) || isinf(x)) && continue + xf = floor(Float64(x)) + q = nbinomccdf(r, p, x) + (q == 0 || q == 1 || xf < 0) && continue + @test Float64(nbinominvccdf(r, p, q)) == xf + end + @testset "invlogcdf round-trip with x=$x" for x in data + (isnan(x) || isinf(x)) && continue + xf = floor(Float64(x)) + lq = nbinomlogcdf(r, p, x) + (!isfinite(lq) || lq == 0 || xf < 0) && continue + @test abs(Float64(nbinominvlogcdf(r, p, lq)) - xf) <= 1 + end + @testset "invlogccdf round-trip with x=$x" for x in data + (isnan(x) || isinf(x)) && continue + xf = floor(Float64(x)) + lq = nbinomlogccdf(r, p, x) + (!isfinite(lq) || lq == 0 || xf < 0) && continue + @test abs(Float64(nbinominvlogccdf(r, p, lq)) - xf) <= 1 + end + end + end + end rmathcomp_tests( "nchisq", [ @@ -386,8 +484,8 @@ end ] ) - rmathcomp_tests( - "pois", [ + @testset "pois" begin + configs = [ ((0.5,), 0.0:20.0), ((1.0,), 0.0:20.0), ((2.0,), 0.0:20.0), @@ -399,7 +497,52 @@ end ((0.5f0,), Float16(0):Float16(20)), ((1 // 2,), (0 // 1):(20 // 1)), ] - ) + @testset "params: $params" for (params, data) in configs + (λ,) = params + T = float(Base.promote_typeof(λ, first(data))) + rtol = _default_rtol(T) + @testset "pdf with x=$x" for x in data + check_rmath(poispdf, (a, λ) -> Rmath.dpois(a, λ, false), params, x, true, rtol) + end + @testset "logpdf with x=$x" for x in data + check_rmath(poislogpdf, (a, λ) -> Rmath.dpois(a, λ, true), params, x, false, rtol) + end + @testset "cdf with x=$x" for x in data + check_rmath(poiscdf, (a, λ) -> Rmath.ppois(a, λ, true, false), params, x, true, rtol) + end + @testset "ccdf with x=$x" for x in data + check_rmath(poisccdf, (a, λ) -> Rmath.ppois(a, λ, false, false), params, x, true, rtol) + end + @testset "logcdf with x=$x" for x in data + check_rmath(poislogcdf, (a, λ) -> Rmath.ppois(a, λ, true, true), params, x, false, rtol) + end + @testset "logccdf with x=$x" for x in data + check_rmath(poislogccdf, (a, λ) -> Rmath.ppois(a, λ, false, true), params, x, false, rtol) + end + if T === Float64 + @testset "invcdf round-trip with x=$x" for x in data + q = poiscdf(λ, x) + (q == 0 || q == 1) && continue + @test poisinvcdf(λ, q) ≈ x rtol = rtol nans = true + end + @testset "invccdf round-trip with x=$x" for x in data + q = poisccdf(λ, x) + (q == 0 || q == 1) && continue + @test poisinvccdf(λ, q) ≈ x rtol = rtol nans = true + end + @testset "invlogcdf round-trip with x=$x" for x in data + lq = poislogcdf(λ, x) + (lq == 0 || lq == -Inf) && continue + @test poisinvlogcdf(λ, lq) ≈ x atol = 1 rtol = rtol nans = true + end + @testset "invlogccdf round-trip with x=$x" for x in data + lq = poislogccdf(λ, x) + (lq == 0 || lq == -Inf) && continue + @test poisinvlogccdf(λ, lq) ≈ x atol = 1 rtol = rtol nans = true + end + end + end + end rmathcomp_tests( "tdist", [ @@ -483,13 +626,36 @@ end ((1, 1), [-10, -6, 2, 24]), ] ) - rmathcomp_tests( - "binom", [ + @testset "binom (edge cases)" begin + configs = [ ((5, 0.5), [-8, -2.3, 1.2, 5.4, 11.9]), ((5, 1 // 2), [-8, -23 // 10, 6 // 5, 27 // 5, 119 // 10]), ((5, 1 // 2), [-8, -2, 6, 12]), ] - ) + @testset "params: $params" for (params, data) in configs + (n, p) = params + T = float(Base.promote_typeof(n, p, first(data))) + rtol = _default_rtol(T) + @testset "pdf with x=$x" for x in data + check_rmath(binompdf, (a, n, p) -> Rmath.dbinom(a, n, p, false), params, x, true, rtol) + end + @testset "logpdf with x=$x" for x in data + check_rmath(binomlogpdf, (a, n, p) -> Rmath.dbinom(a, n, p, true), params, x, false, rtol) + end + @testset "cdf with x=$x" for x in data + check_rmath(binomcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, false), params, x, true, rtol) + end + @testset "ccdf with x=$x" for x in data + check_rmath(binomccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, false), params, x, true, rtol) + end + @testset "logcdf with x=$x" for x in data + check_rmath(binomlogcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, true), params, x, false, rtol) + end + @testset "logccdf with x=$x" for x in data + check_rmath(binomlogccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, true), params, x, false, rtol) + end + end + end rmathcomp_tests( "fdist", [ ((1.0, 1.0), [-10.0, -6.3]), @@ -504,11 +670,34 @@ end ((1, 1), [-10, -6]), ] ) - rmathcomp_tests( - "pois", [ + @testset "pois (edge cases)" begin + configs = [ ((0.5,), [-10, -2.5, 1.3, 8.7]), ((1 // 2,), [-10, -5 // 2, 13 // 10, 87 // 10]), ((1,), [-10, -3]), ] - ) + @testset "params: $params" for (params, data) in configs + (λ,) = params + T = float(Base.promote_typeof(λ, first(data))) + rtol = _default_rtol(T) + @testset "pdf with x=$x" for x in data + check_rmath(poispdf, (a, λ) -> Rmath.dpois(a, λ, false), params, x, true, rtol) + end + @testset "logpdf with x=$x" for x in data + check_rmath(poislogpdf, (a, λ) -> Rmath.dpois(a, λ, true), params, x, false, rtol) + end + @testset "cdf with x=$x" for x in data + check_rmath(poiscdf, (a, λ) -> Rmath.ppois(a, λ, true, false), params, x, true, rtol) + end + @testset "ccdf with x=$x" for x in data + check_rmath(poisccdf, (a, λ) -> Rmath.ppois(a, λ, false, false), params, x, true, rtol) + end + @testset "logcdf with x=$x" for x in data + check_rmath(poislogcdf, (a, λ) -> Rmath.ppois(a, λ, true, true), params, x, false, rtol) + end + @testset "logccdf with x=$x" for x in data + check_rmath(poislogccdf, (a, λ) -> Rmath.ppois(a, λ, false, true), params, x, false, rtol) + end + end + end end