From a6b520beb6c9643a89f94be654abf9ad232038a3 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 22 Mar 2026 13:55:25 +0100 Subject: [PATCH] Replace Rmath hypergeometric with pure Julia implementation Port hypergeometric distribution functions from Ian Smith's VBA code, removing the dependency on Rmath for all 10 hyper functions (pdf, logpdf, cdf, ccdf, logcdf, logccdf, invcdf, invccdf, invlogcdf, invlogccdf). Add helper functions logfbit, lfbaccdif1, and ab_minus_cd to misc.jl for accurate Stirling correction and cross-product computation. --- src/distrs/hyper.jl | 541 +++++++++++++++++++++++++++++++++++++++++++- src/misc.jl | 58 +++++ test/rmath.jl | 146 ++++++++++++ 3 files changed, 733 insertions(+), 12 deletions(-) diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl index e7cf1c7..7993328 100644 --- a/src/distrs/hyper.jl +++ b/src/distrs/hyper.jl @@ -1,45 +1,562 @@ -# functions related to hyper-geometric distribution +# Functions related to hypergeometric distribution +# Pure Julia implementation based on VBA code by Ian Smith +# https://iandjmsmith.wordpress.com/ +# License: MIT + +# Constants +const _hyper_cfVSmall = 1.0e-15 +const _hyper_scalefactor = 1.1579208923731619542357098500869e+77 # 2^256 +const _hyper_scalefactor2 = 8.6361685550944446253863518628004e-78 # 2^-256 +const _hyper_minLog1Value = -0.79149064 +const _hyper_max_discrete = 9.007199254740991e15 # 2^53 +const _hyper_max_crit = 4.503599627370496e15 # 2^52 + +# Internal PMF computation +# ai = x, aji = n - x, aki = ms - x, amkji = mf - n + x +function _hypergeometric_term(ai::Float64, aji::Float64, aki::Float64, amkji::Float64) + ak = aki + ai # ms + amk = amkji + aji # mf + aj = aji + ai # n + am = amk + ak # ms + mf + amj = amkji + aki # ms + mf - n + + if am > _hyper_max_discrete + return NaN + end + + if ai == 0 && (aji <= 0 || aki <= 0 || amj < 0 || amk < 0) + return 1.0 + elseif ai > 0 && min(aki, aji) == 0 && max(amj, amk) == 0 + return 1.0 + elseif ai >= 0 && amkji > -1 && aki > -1 && aji >= 0 + c1 = lfbaccdif1(ak, amk) - lfbaccdif1(ai, aki) - lfbaccdif1(ai, aji) - lfbaccdif1(aki, amkji) - logfbit(ai) + + ai1 = ai + 1.0; aj1 = aj + 1.0; ak1 = ak + 1.0; am1 = am + 1.0 + aki1 = aki + 1.0; aji1 = aji + 1.0 + amk1 = amk + 1.0; amj1 = amj + 1.0; amkji1 = amkji + 1.0 + + cjkmi = ab_minus_cd(aji, aki, ai, amkji) + + c5 = (cjkmi - ai) / (amkji1 * am1) + c3 = if c5 < _hyper_minLog1Value + amkji * (log((amj1 * amk1) / (amkji1 * am1)) - c5) - c5 + else + amkji * log1pmx(c5) - c5 + end + + c5 = (-cjkmi - aji) / (aki1 * am1) + c4 = if c5 < _hyper_minLog1Value + aki * (log((ak1 * amj1) / (aki1 * am1)) - c5) - c5 + else + aki * log1pmx(c5) - c5 + end + c3 += c4 + + c5 = (-cjkmi - aki) / (aji1 * am1) + c4 = if c5 < _hyper_minLog1Value + aji * (log((aj1 * amk1) / (aji1 * am1)) - c5) - c5 + else + aji * log1pmx(c5) - c5 + end + c3 += c4 + + c5 = (cjkmi - amkji) / (ai1 * am1) + c4 = if c5 < _hyper_minLog1Value + ai * (log((aj1 * ak1) / (ai1 * am1)) - c5) - c5 + else + ai * log1pmx(c5) - c5 + end + c3 += c4 + + logterm = (c1 + 1.0 / am1) + c3 + sqrtterm = sqrt((amk1 * ak1) * (aj1 * amj1) / ((amkji1 * aki1 * aji1) * (am1 * ai1))) + return exp(logterm) * sqrtterm * Float64(invsqrt2π) + else + return 0.0 + end +end + +# Internal CDF computation +function _hypergeometric(ai::Float64, aji::Float64, aki::Float64, amkji::Float64, comp::Bool) + # Determine swap direction for numerical stability + if amkji > -1 && amkji < 0 + ip1 = -amkji + mkji = ip1 - 1.0 + allIntegral = false + else + ip1 = amkji + 1.0 + mkji = amkji + allIntegral = ai == floor(ai) && aji == floor(aji) && aki == floor(aki) && mkji == floor(mkji) + end + + if allIntegral + swapped = (ai + 0.5) * (mkji + 0.5) >= (aki - 0.5) * (aji - 0.5) + elseif (ai < 100 && ai == floor(ai)) || mkji < 0 + swapped = if comp + (ai + 0.5) * (mkji + 0.5) >= aki * aji + else + (ai + 0.5) * (mkji + 0.5) >= aki * aji + 1000 + end + elseif ai < 1 + swapped = (ai + 0.5) * (mkji + 0.5) >= aki * aji + elseif aji < 1 || aki < 1 || (ai < 1 && ai > 0) + swapped = false + else + swapped = (ai + 0.5) * (mkji + 0.5) >= (aki - 0.5) * (aji - 0.5) + end + + if !swapped + i = ai; ji = aji; ki = aki + else + i = aji - 1.0; ji = ai + 1.0; ki = ip1 + ip1 = aki; mkji = aki - 1.0 + end + + c2 = ji + i + c4_pop = mkji + ki + c2 # population size + + if c4_pop > _hyper_max_discrete + return NaN + end + + if (i >= 0 && (ji <= 0 || ki <= 0)) || (ip1 + ki <= 0) || (ip1 + ji <= 0) + exact = true + prob = i >= 0 ? 1.0 : 0.0 + elseif ip1 > 0 && ip1 < 1 + exact = false + prob = _hypergeometric_term(i, ji, ki, ip1) * (ip1 * (c4_pop + 1.0)) / ((ki + ip1) * (ji + ip1)) + else + exact = (i == 0 && (ji <= 0 || ki <= 0 || mkji + ki < 0 || mkji + ji < 0)) || + (i > 0 && min(ki, ji) == 0 && max(mkji + ki, mkji + ji) == 0) + prob = _hypergeometric_term(i, ji, ki, mkji) + end + + if exact || prob == 0.0 + return (swapped == comp) ? prob : 1.0 - prob + end + + a1 = 0.0 + c4 = c4_pop # working copy for CF + + if i < mkji + must_do_cf = i != floor(i) + maxSums = floor(i) + else + must_do_cf = mkji != floor(mkji) + maxSums = floor(max(mkji, 0.0)) + end + + if must_do_cf + sumAlways = 0 + sumFactor = 5 + else + sumAlways = 20 + sumFactor = 10 + end + + if maxSums > sumAlways || must_do_cf + numb = floor(sumFactor / c4 * exp(log((ki + i) * (ji + i) * (ip1 + ji) * (ip1 + ki)) / 3.0)) + numb = floor(i - (ki + i) * (ji + i) / c4 + numb) + numb = clamp(numb, 0.0, maxSums) + else + numb = maxSums + end + + if 2 * numb <= maxSums || must_do_cf + # Continued fraction evaluation + b1 = 1.0 + c1 = 0.0 + c2_cf = i - numb + c3 = mkji - numb + s = c3 + a2 = c2_cf + c3 -= 1.0 + b2 = ab_minus_cd(ki + numb + 1.0, ji + numb + 1.0, c2_cf - 1.0, c3) + bn = b2 + bnAdd = c3 + c4 + c2_cf - 2.0 + + while b2 > 0 && abs(a2 * b1 - a1 * b2) > abs(_hyper_cfVSmall * b1 * a2) + c1 += 1.0; c2_cf -= 1.0 + an = (c1 * c2_cf) * (c3 * c4) + c3 -= 1.0; c4 -= 1.0 + bn += bnAdd; bnAdd -= 4.0 + a1 = bn * a2 + an * a1 + b1 = bn * b2 + an * b1 + if b1 > _hyper_scalefactor + a1 *= _hyper_scalefactor2; b1 *= _hyper_scalefactor2 + a2 *= _hyper_scalefactor2; b2 *= _hyper_scalefactor2 + end + c1 += 1.0; c2_cf -= 1.0 + an = (c1 * c2_cf) * (c3 * c4) + c3 -= 1.0; c4 -= 1.0 + bn += bnAdd; bnAdd -= 4.0 + a2 = bn * a1 + an * a2 + b2 = bn * b1 + an * b2 + if b2 > _hyper_scalefactor + a1 *= _hyper_scalefactor2; b1 *= _hyper_scalefactor2 + a2 *= _hyper_scalefactor2; b2 *= _hyper_scalefactor2 + end + end + + if b1 < 0 || b2 < 0 + return NaN + else + a1 = a2 / b2 * s + end + else + numb = maxSums + end + + # Direct summation + c1_s = i - numb + 1.0 + c2_s = mkji - numb + 1.0 + c3_s = ki + numb + c4_s = ji + numb + for _ in 1:Int(numb) + a1 = (1.0 + a1) * ((c1_s * c2_s) / (c3_s * c4_s)) + c1_s += 1.0; c2_s += 1.0; c3_s -= 1.0; c4_s -= 1.0 + end + + a1 = (1.0 + a1) * prob + + if swapped == comp + return a1 + else + return a1 > 0.99 ? NaN : 1.0 - a1 + end +end + +# Inverse CDF search (lower tail) +# Faithful port of VBA crithyperg by Ian Smith +const _hyper_nearly_zero = 9.99999983659714e-317 +const _hyper_cSmall = 5.562684646268003457725581793331e-309 + +function _crithyperg(j::Float64, k::Float64, m::Float64, cprob::Float64) + if cprob > 0.5 + return _critcomphyperg(j, k, m, 1.0 - cprob) + end + + mx = min(j, k) + mn = max(0.0, j + k - m) + + # Normal approximation for initial guess + i = j * k / m + norminvcdf(cprob) * sqrt(j * k * (m - j) * (m - k) / (m * m * max(m - 1.0, 1.0))) + + while true + i = clamp(floor(i + 0.5), mn, mx) + if i >= _hyper_max_crit + return i + end + pr = _hypergeometric(i, j - i, k - i, m - k - j + i, false) + tpr = 0.0 + if pr >= cprob + if i == mn + return mn + end + tpr = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + if pr < 1.00001 * tpr + # PMF dominates: ratio stepping left + tpr *= ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + i -= 1.0 + while tpr > cprob + tpr *= ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + i -= 1.0 + end + # Falls through to top of while(true) for re-evaluation + else + pr -= tpr + if pr < cprob + return i + end + i -= 1.0 + if i == mn + return mn + end + temp = (pr - cprob) / tpr + if temp > 10 + # Large jump with parabolic refinement + temp = floor(temp + 0.5) + i -= temp + temp2 = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + i -= temp * (tpr - temp2) / (2.0 * temp2) + else + tpr *= ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + pr -= tpr + if pr < cprob + return i + end + i -= 1.0 + temp2 = (pr - cprob) / tpr + if temp2 < temp - 0.9 + # Linear stepping + while pr >= cprob + tpr *= ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + pr -= tpr + i -= 1.0 + end + return i + 1.0 + else + # Log-ratio jump + ratio = ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + temp = floor(log(cprob / pr) / log(ratio) + 0.5) + i -= temp + temp2 = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + if temp2 > _hyper_nearly_zero + temp = log((cprob / pr) * (tpr / temp2)) / log(ratio) + i -= temp + end + end + end + # Falls through to top of while(true) for re-evaluation + end + else + # Search right + while tpr < _hyper_cSmall && pr < cprob + i += 1.0 + tpr = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + pr += tpr + end + while pr < cprob + i += 1.0 + tpr *= ((k - i + 1.0) * (j - i + 1.0)) / (i * (m - j - k + i)) + pr += tpr + end + return i + end + end + return +end + +# Inverse CDF search (upper tail) +# Faithful port of VBA critcomphyperg by Ian Smith +function _critcomphyperg(j::Float64, k::Float64, m::Float64, cprob::Float64) + if cprob > 0.5 + return _crithyperg(j, k, m, 1.0 - cprob) + end + + mx = min(j, k) + mn = max(0.0, j + k - m) + + # Normal approximation + i = j * k / m - norminvcdf(cprob) * sqrt(j * k * (m - j) * (m - k) / (m * m * max(m - 1.0, 1.0))) + + while true + i = clamp(floor(i + 0.5), mn, mx) + if i >= _hyper_max_crit + return i + end + pr = _hypergeometric(i, j - i, k - i, m - k - j + i, true) + tpr = 0.0 + if pr > cprob + i += 1.0 + tpr = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + if pr < (1.0 + 0.00001) * tpr + # PMF dominates: ratio stepping right + while tpr > cprob + i += 1.0 + tpr *= ((k - i + 1.0) * (j - i + 1.0)) / (i * (m - j - k + i)) + end + # Falls through to top of while(true) for re-evaluation + else + pr -= tpr + if pr <= cprob + return i + end + temp = (pr - cprob) / tpr + if temp > 10 + # Large jump with parabolic refinement + temp = floor(temp + 0.5) + i += temp + temp2 = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + i += temp * (tpr - temp2) / (2.0 * temp2) + else + i += 1.0 + tpr *= ((k - i + 1.0) * (j - i + 1.0)) / (i * (m - j - k + i)) + pr -= tpr + if pr <= cprob + return i + end + temp2 = (pr - cprob) / tpr + if temp2 < temp - 0.9 + # Linear stepping + while pr > cprob + i += 1.0 + tpr *= ((k - i + 1.0) * (j - i + 1.0)) / (i * (m - j - k + i)) + pr -= tpr + end + return i + else + # Log-ratio jump + ratio = ((k - i + 1.0) * (j - i + 1.0)) / (i * (m - j - k + i)) + temp = floor(log(cprob / pr) / log(ratio) + 0.5) + i += temp + temp2 = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + if temp2 > _hyper_nearly_zero + temp = log((cprob / pr) * (tpr / temp2)) / log(ratio) + i += temp + end + end + end + # Falls through to top of while(true) for re-evaluation + end + else + # Search left + while tpr < _hyper_cSmall && pr <= cprob + tpr = _hypergeometric_term(i, j - i, k - i, m - k - j + i) + pr += tpr + i -= 1.0 + end + while pr <= cprob + tpr *= ((i + 1.0) * (m - j - k + i + 1.0)) / ((k - i) * (j - i)) + pr += tpr + i -= 1.0 + end + return i + 1.0 + end + end + return +end + +# Wrappers matching VBA crit_hypergeometric / comp_crit_hypergeometric +function _hyper_invcdf(ms::Float64, mf::Float64, n::Float64, p::Float64) + m = ms + mf + mn = max(0.0, n - mf) + mx = min(ms, n) + + if p < 0 || p > 1 || isnan(p) + return NaN + elseif p == 0 + return mn + elseif ms == 0 || n == 0 + return 0.0 + elseif mf == 0 || n == m + return mx + elseif p == 1 + return mx + end + + i = _crithyperg(n, ms, m, p) + + # Post-correction (from crit_hypergeometric) + pr = _hypergeometric(i, n - i, ms - i, mf - n + i, false) + if pr == p + return i + elseif pr > p + i2 = i - 1.0 + if i2 >= mn + pr2 = _hypergeometric(i2, n - i2, ms - i2, mf - n + i2, false) + if pr2 >= p + return i2 + end + end + return i + else + return i + 1.0 + end +end + +function _hyper_invccdf(ms::Float64, mf::Float64, n::Float64, q::Float64) + m = ms + mf + mn = max(0.0, n - mf) + mx = min(ms, n) + + if q < 0 || q > 1 || isnan(q) + return NaN + elseif q == 1 + return mn + elseif ms == 0 || n == 0 + return 0.0 + elseif mf == 0 || n == m + return mx + elseif q == 0 + return mx + end + + i = _critcomphyperg(n, ms, m, q) + + # Post-correction (from comp_crit_hypergeometric) + pr = _hypergeometric(i, n - i, ms - i, mf - n + i, true) + if pr == q + return i + elseif pr < q + i2 = i - 1.0 + if i2 >= mn + pr2 = _hypergeometric(i2, n - i2, ms - i2, mf - n + i2, true) + if pr2 <= q + return i2 + end + end + return i + else + return i + 1.0 + end +end + +# Public API -# Rmath implementations function hyperpdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.dhyper(x, ms, mf, n, false)) + _x = round(Float64(x)) + _ms = Float64(ms); _mf = Float64(mf); _n = Float64(n) + result = _hypergeometric_term(_x, _n - _x, _ms - _x, _mf - _n + _x) + return convert(T, result) end + function hyperlogpdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.dhyper(x, ms, mf, n, true)) + return convert(T, log(Float64(hyperpdf(ms, mf, n, x)))) end function hypercdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.phyper(x, ms, mf, n, true, false)) + _x = floor(Float64(x) + 1.0e-7) + _ms = Float64(ms); _mf = Float64(mf); _n = Float64(n) + result = _hypergeometric(_x, _n - _x, _ms - _x, _mf - _n + _x, false) + return convert(T, result) end + function hyperccdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.phyper(x, ms, mf, n, false, false)) + _x = floor(Float64(x) + 1.0e-7) + _ms = Float64(ms); _mf = Float64(mf); _n = Float64(n) + result = _hypergeometric(_x, _n - _x, _ms - _x, _mf - _n + _x, true) + return convert(T, result) end + function hyperlogcdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.phyper(x, ms, mf, n, true, true)) + return convert(T, log(Float64(hypercdf(ms, mf, n, x)))) end + function hyperlogccdf(ms::Real, mf::Real, n::Real, x::Real) T = float(Base.promote_typeof(ms, mf, n, x)) - return convert(T, Rmath.phyper(x, ms, mf, n, false, true)) + return convert(T, log(Float64(hyperccdf(ms, mf, n, x)))) end function hyperinvcdf(ms::Real, mf::Real, n::Real, q::Real) T = float(Base.promote_typeof(ms, mf, n, q)) - return convert(T, Rmath.qhyper(q, ms, mf, n, true, false)) + result = _hyper_invcdf(Float64(ms), Float64(mf), Float64(n), Float64(q)) + return convert(T, result) end + function hyperinvccdf(ms::Real, mf::Real, n::Real, q::Real) T = float(Base.promote_typeof(ms, mf, n, q)) - return convert(T, Rmath.qhyper(q, ms, mf, n, false, false)) + result = _hyper_invccdf(Float64(ms), Float64(mf), Float64(n), Float64(q)) + return convert(T, result) end + function hyperinvlogcdf(ms::Real, mf::Real, n::Real, lq::Real) T = float(Base.promote_typeof(ms, mf, n, lq)) - return convert(T, Rmath.qhyper(lq, ms, mf, n, true, true)) + _lq = Float64(lq) + isinf(_lq) && return convert(T, NaN) + result = _hyper_invcdf(Float64(ms), Float64(mf), Float64(n), exp(_lq)) + return convert(T, result) end + function hyperinvlogccdf(ms::Real, mf::Real, n::Real, lq::Real) T = float(Base.promote_typeof(ms, mf, n, lq)) - return convert(T, Rmath.qhyper(lq, ms, mf, n, false, true)) + _lq = Float64(lq) + isinf(_lq) && return convert(T, NaN) + result = _hyper_invccdf(Float64(ms), Float64(mf), Float64(n), exp(_lq)) + return convert(T, result) end diff --git a/src/misc.jl b/src/misc.jl index b7cbdff..bd4326f 100644 --- a/src/misc.jl +++ b/src/misc.jl @@ -97,3 +97,61 @@ function lstirling_asym(x::Float32) 8.417508417508f-4 ) / x # 1/1188 x^-9 end + +""" + logfbit(x) + +Stirling error term for log-factorial: + + logfbit(x) = log(x!) - log(√2π) + (x+1) - (x+0.5)*log(x+1) + +Equivalent to `lstirling_asym(x + 1)`. +""" +logfbit(x) = lstirling_asym(x + one(x)) + +""" + lfbaccdif1(a, b) + +Accurate computation of `logfbit(b) - logfbit(a + b)`. +Uses a polynomial expansion for `b ≥ 8` that avoids cancellation. +Based on VBA code by Ian Smith. +""" +function lfbaccdif1(a::Float64, b::Float64) + if a < 0 + return -lfbaccdif1(-a, b + a) + end + if b >= 8 + y1 = b + 1.0 + y2 = inv(y1 * y1) + x1 = a + b + 1.0 + x2 = inv(x1 * x1) + + # Initialize with innermost tuned coefficient (lfbc9) + x3 = x2 * 1.6769380337122674863 + y3 = y2 * 1.6769380337122674863 + acc = x2 * (a * (x1 + y1) * y3) + + # Unroll from lfbc8 down to lfbc2 + for c in (0.35068485511628418514, 1 / 13, 691 / 30030, 1 / 99, 1 / 140, 1 / 105, 1 / 30) + x3 = x2 * (c - x3) + y3 = y2 * (c - y3) + acc = x2 * (a * (x1 + y1) * y3 - acc) + end + + return (a * (1.0 - y3) - y1 * acc) / (12.0 * x1 * y1) + else + return logfbit(b) - logfbit(a + b) + end +end +lfbaccdif1(a::Real, b::Real) = lfbaccdif1(Float64(a), Float64(b)) + +""" + ab_minus_cd(a, b, c, d) + +Accurate computation of `a * b - c * d` using FMA. +""" +function ab_minus_cd(a::Float64, b::Float64, c::Float64, d::Float64) + w = c * d + return fma(a, b, -w) - fma(c, d, -w) +end +ab_minus_cd(a::Real, b::Real, c::Real, d::Real) = ab_minus_cd(Float64(a), Float64(b), Float64(c), Float64(d)) diff --git a/test/rmath.jl b/test/rmath.jl index 8d5bd60..f0ee9e7 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -291,6 +291,152 @@ end ] ) + # Additional hypergeometric tests with larger parameters for coverage. + # Uses self-consistency round-trips for inverse functions since our CDF + # differs from Rmath's by ~1 ULP which can flip discrete answers. + @testset "hyper (extended)" begin + for (ms, mf, n, xs) in [ + (10, 20, 15, 0.0:15.0), + (50, 100, 80, 0.0:50.0), + (0, 10, 5, 0.0:5.0), # ms=0 + (10, 0, 5, 0.0:5.0), # mf=0 + (5, 5, 10, 0.0:5.0), # n=ms+mf + ] + @testset "params: ($ms, $mf, $n)" begin + @testset "pdf/cdf match Rmath for x=$x" for x in xs + @test hyperpdf(ms, mf, n, x) ≈ + Rmath.dhyper(Float64(x), Float64(ms), Float64(mf), Float64(n), false) nans = true + @test hypercdf(ms, mf, n, x) ≈ + Rmath.phyper(Float64(x), Float64(ms), Float64(mf), Float64(n), true, false) nans = true + @test hyperccdf(ms, mf, n, x) ≈ + Rmath.phyper(Float64(x), Float64(ms), Float64(mf), Float64(n), false, false) nans = true + end + @testset "invcdf round-trip x=$x" for x in xs + xf = floor(Float64(x)) + q = hypercdf(ms, mf, n, x) + (0 < q < 1 && xf >= 0) || continue + @test Float64(hyperinvcdf(ms, mf, n, q)) == xf + end + @testset "invccdf round-trip x=$x" for x in xs + xf = floor(Float64(x)) + q = hyperccdf(ms, mf, n, x) + (0 < q < 1 && xf >= 0) || continue + @test Float64(hyperinvccdf(ms, mf, n, q)) == xf + end + end + end + end + + # Hypergeometric with non-integer parameters (exercises VBA continuous paths) + @testset "hyper (non-integer params)" begin + # Non-integer ms/mf exercises the non-allIntegral swap branches, + # fractional amkji paths, and non-integer CF summation limits + for (ms, mf, n) in [(2.5, 3.5, 4), (5.5, 10.5, 8), (0.5, 9.5, 5)] + mn = max(0, n - mf) + mx = min(ms, n) + @testset "params: ($ms, $mf, $n)" begin + for x in floor(mn):floor(mx) + pdf = @inferred hyperpdf(ms, mf, n, x) + cdf = @inferred hypercdf(ms, mf, n, x) + ccdf = @inferred hyperccdf(ms, mf, n, x) + @test pdf >= 0 + @test 0 <= cdf <= 1 + @test 0 <= ccdf <= 1 + @test cdf + ccdf ≈ 1 atol = 1.0e-3 + end + # CDF should be monotonically non-decreasing + xs = floor(mn):floor(mx) + cdfs = [hypercdf(ms, mf, n, x) for x in xs] + @test issorted(cdfs) + end + end + end + + # Hypergeometric edge cases + @testset "hyper (edge cases)" begin + # NaN for invalid params + @test isnan(hyperinvcdf(2, 3, 4, -0.1)) + @test isnan(hyperinvccdf(2, 3, 4, -0.1)) + # Boundary probabilities + @test hyperinvcdf(10, 20, 15, 0.0) == 0.0 + @test hyperinvcdf(10, 20, 15, 1.0) == 10.0 + @test hyperinvccdf(10, 20, 15, 0.0) == 10.0 + @test hyperinvccdf(10, 20, 15, 1.0) == 0.0 + # ms=0: only x=0 possible + @test hyperpdf(0, 10, 5, 0) == 1.0 + @test hypercdf(0, 10, 5, 0) == 1.0 + # mf=0: x must equal n (if n <= ms) — degenerate PMF=1 (line 30) + @test hyperpdf(10, 0, 5, 5) == 1.0 + @test hyperpdf(5, 0, 5, 5) == 1.0 + @test hypercdf(10, 0, 5, 5) == 1.0 + # Skewed params triggering minLog1Value fallback in PMF (lines 42, 57, 65) + @test hyperpdf(100, 5, 100, 99) > 0 # c5_1 < minLog1Value + @test hyperpdf(100, 5, 5, 4) > 0 # c5_3 < minLog1Value + @test hyperpdf(5, 100, 5, 1) > 0 # c5_4 < minLog1Value + @test hyperpdf(100, 5, 100, 100) >= 0 + # Skewed distribution where PMF dominates CDF (lines 258-282, 314-343) + @test hyperinvcdf(1, 1000, 500, 0.3) == 0.0 + @test hyperinvccdf(1, 1000, 500, 0.3) == 1.0 + @test hyperinvcdf(1, 1000, 500, 0.9) == 1.0 + @test hyperinvccdf(1, 1000, 500, 0.9) == 0.0 + # Inverse edge cases for ms=0, mf=0, n=0 + @test hyperinvcdf(0, 10, 5, 0.0) == 0.0 + @test hyperinvcdf(0, 10, 5, 0.5) == 0.0 # ms=0 → 0 (line 366) + @test hyperinvcdf(5, 10, 0, 0.5) == 0.0 # n=0 → 0 (line 366) + @test hyperinvcdf(10, 0, 5, 0.5) == 5.0 # mf=0 → ms (line 368) + @test hyperinvccdf(0, 10, 5, 0.5) == 0.0 # ms=0 → 0 (line 403) + @test hyperinvccdf(5, 10, 0, 0.5) == 0.0 # n=0 → 0 (line 403) + @test hyperinvccdf(10, 0, 5, 0.5) == 5.0 # mf=0 → ms (line 405) + @test hyperinvccdf(10, 20, 15, 0.0) == 10.0 + # Skewed distribution where initial guess overshoots (exercises search) + @test hyperinvcdf(100, 200, 150, 0.001) isa Float64 + @test hyperinvccdf(100, 200, 150, 0.001) isa Float64 + @test hyperinvcdf(100, 200, 150, 0.999) isa Float64 + @test hyperinvccdf(100, 200, 150, 0.999) isa Float64 + # Very large population triggers NaN (line 24, 119) + @test isnan(hyperpdf(2^53, 1, 1, 1)) + @test isnan(hypercdf(2^53, 1, 1, 1)) + # Exact boundary path in CDF (lines 123-124) + @test hypercdf(5, 5, 10, 10) == 1.0 + end + + # Tests calling internal functions directly with non-integer args to exercise + # code paths that are unreachable from the public integer-parameter API. + @testset "hyper (internal non-integer)" begin + using StatsFuns: _hypergeometric_term, _hypergeometric + + # Fractional amkji between -1 and 0 (lines 83-85) + @test _hypergeometric(1.0, 3.0, 2.0, -0.5, false) ≈ 0.09955360139478697 + @test _hypergeometric(1.0, 3.0, 2.0, -0.5, true) ≈ 0.900446398605213 + + # Fractional ai < 1 (lines 100-101) + @test 0 < _hypergeometric(0.5, 3.5, 2.0, 1.0, false) < 1 + + # Fractional aji < 1 (lines 102-103) + @test 0 < _hypergeometric(2.0, 0.5, 2.0, 1.0, false) < 1 + + # Fractional aki < 1 (lines 102-103) + @test 0 < _hypergeometric(2.0, 3.0, 0.5, 1.0, false) < 1 + + # Non-integer, ai >= 100 (line 105) + @test 0 < _hypergeometric(100.0, 50.0, 100.5, 50.5, false) < 1 + + # aji < 1 with ai >= 100, swapped=false (line 103) + @test 0 < _hypergeometric(100.0, 0.5, 50.0, 2.0, false) < 1 + + # Exact boundary with fractional params (lines 123-124) + @test _hypergeometric(1.0, 0.0, 2.0, -0.3, false) == 1.0 + + # minLog1Value fallback for c5_2 (line 49) + @test _hypergeometric_term(5.0, 30.0, 0.2, 0.1) > 0 + + # minLog1Value fallback for c5_3 (line 57) + @test _hypergeometric_term(5.0, 0.2, 30.0, 0.1) > 0 + + # minLog1Value fallback for c5_4 (line 65) + @test _hypergeometric_term(0.2, 0.2, 0.2, 30.0) > 0 + end + rmathcomp_tests( "nbeta", [ ((1.0, 1.0, 0.0), 0.01:0.01:0.99),