From 0b6ca0962b6d8c879b4d131663ae82e28e69d559 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 3 Feb 2026 01:58:17 -0500 Subject: [PATCH 1/4] add hyper pdfs and cdfs --- src/distrs/hyper.jl | 75 +++++++++++++++++++++++++++++++++++++++++---- test/rmath.jl | 1 + 2 files changed, 70 insertions(+), 6 deletions(-) diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl index 6c2869f5..5fe0d452 100644 --- a/src/distrs/hyper.jl +++ b/src/distrs/hyper.jl @@ -2,13 +2,76 @@ # R implementations using .RFunctions: - hyperpdf, - hyperlogpdf, - hypercdf, - hyperccdf, - hyperlogcdf, - hyperlogccdf, + # hyperpdf, + # hyperlogpdf, + # hypercdf, + # hyperccdf, + # hyperlogcdf, + # hyperlogccdf, hyperinvcdf, hyperinvccdf, hyperinvlogcdf, hyperinvlogccdf + +function hyperpdf(ms, mf, n, x) + exp(hyperlogpdf(ms, mf, n, x)) +end + +hyperlogpdf(ms, mf, n, x) = hyperlogpdf(promote(ms,mf,n,x)...) + +function hyperlogpdf(ms::T, mf::T, n::T, x::T) where T<:Real + #loggamma(ms + 1) - loggamma(x + 1) - loggamma(ms - x + 1) + + #loggamma(mf + 1) - loggamma(n - x + 1) - loggamma(mf - n + x + 1) - + #loggamma(ms + mf + 1) + loggamma(n + 1) + loggamma(ms + mf - n + 1) + -log1p(ms) - logbeta(x + 1, ms - x + 1) - + log1p(mf) - logbeta(n - x + 1, mf - n + x + 1) + + log1p(ms+mf) + logbeta(n + 1, ms + mf - n + 1) +end + + +function _hypercdf(ms, mf, n, x, invert) + N = ms+mf + mode = fld((n + 1) * (ms + 1), N + 2) + local result + if x < mode + result = hyperpdf(ms, mf, n, x) + diff = result + lower_limit = max(0, n+ms-N) + while x != lower_limit && diff > (invert ? result : 1)*eps() + diff *= x * (N + x - n - ms) / ((1+n-x)*(1+ms -x)) + result += diff + x -= 1 + end + else + invert = !invert + x += 1 + result = hyperpdf(ms, mf, n, x) + diff = result + upper_limit = min(ms, n) + while x != upper_limit && diff > (invert ? 1 : result)*eps() + diff *= (n-x) * (ms-x) / ((1+x) * (N+x+1-n-ms)) + result += diff + x += 1 + end + end + result, invert +end + +function hypercdf(ms, mf, n, x) + result, invert = _hypercdf(ms, mf, n, x, false) + invert ? 1-result : result +end + +function hyperccdf(ms, mf, n, x) + result, invert = _hypercdf(ms, mf, n, x, true) + invert ? 1-result : result +end + +function hyperlogcdf(ms, mf, n, x) + result, invert = _hypercdf(ms, mf, n, x, false) + invert ? log1p(-result) : log(result) +end +function hyperlogccdf(ms, mf, n, x) + result, invert = _hypercdf(ms, mf, n, x, true) + invert ? log1p(-result) : log(result) +end diff --git a/test/rmath.jl b/test/rmath.jl index 373e34bc..f808b115 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -317,6 +317,7 @@ end ((2, 3, 4), 0.0f0:4.0f0), ((2, 3, 4), Float16(0):Float16(4)), ((2, 3, 4), (0 // 1):(4 // 1)), + ((40, 60, 40), 0.0:50.0), ] ) From 7b8e1df0adb99fa489fd017b92ed6543181f7173 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 13 Feb 2026 18:29:19 -0500 Subject: [PATCH 2/4] slow but accurate hyperpdf --- src/distrs/hyper.jl | 60 +++++++++++++++++++++++++++------------------ test/rmath.jl | 15 +++++++++--- 2 files changed, 48 insertions(+), 27 deletions(-) diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl index 5fe0d452..2766026c 100644 --- a/src/distrs/hyper.jl +++ b/src/distrs/hyper.jl @@ -2,42 +2,55 @@ # R implementations using .RFunctions: - # hyperpdf, - # hyperlogpdf, - # hypercdf, - # hyperccdf, - # hyperlogcdf, - # hyperlogccdf, hyperinvcdf, hyperinvccdf, hyperinvlogcdf, hyperinvlogccdf -function hyperpdf(ms, mf, n, x) - exp(hyperlogpdf(ms, mf, n, x)) -end -hyperlogpdf(ms, mf, n, x) = hyperlogpdf(promote(ms,mf,n,x)...) +hyperpdf(ms::Real, mf::Real, n::Real, x::Real) = hyperpdf(promote(ms,mf,n,x)...) + +function hyperpdf(ms::T, mf::T, n::T, x::T) where T<:Real + M, N = ms, ms+mf + fxnprev = fMNprev = one(T) + fxn = 1 - 2x/n + fMN = 1 - 2M/N + eiprev = ones(T, Int(n+1)) + k = 0:n + ei = 1 .- k .*(2/n) + square(x) = x^2 + res = fma(fxn, fMN/sum(square, ei), 1/(n+1)) + ei = collect(ei) + for i in 1:n-1 + fxnnext = ((2i + 1) * (n - 2x) * fxn - i * (n + i + 1) * fxnprev) / ((i + 1) * (n - i)) + fxn, fxnprev = fxnnext, fxn + fMNnext = ((2i + 1) * (N - 2M) * fMN - i * (N + i + 1) * fMNprev) / ((i + 1) * (N - i)) + fMN, fMNprev = fMNnext, fMN + @. eiprev = ((2i + 1) * (n - 2k) * ei - i * (n + i + 1) * eiprev) / ((i + 1) * (n - i)) + ei, eiprev = eiprev, ei + correction = fxn*fMN/sum(square, ei) + abs(correction) < res*eps(T) && break + res += correction + end + res +end +hyperlogpdf(ms::Real, mf::Real, n::Real, x::Real) = hyperlogpdf(promote(ms,mf,n,x)...) function hyperlogpdf(ms::T, mf::T, n::T, x::T) where T<:Real - #loggamma(ms + 1) - loggamma(x + 1) - loggamma(ms - x + 1) + - #loggamma(mf + 1) - loggamma(n - x + 1) - loggamma(mf - n + x + 1) - - #loggamma(ms + mf + 1) + loggamma(n + 1) + loggamma(ms + mf - n + 1) - -log1p(ms) - logbeta(x + 1, ms - x + 1) - - log1p(mf) - logbeta(n - x + 1, mf - n + x + 1) + - log1p(ms+mf) + logbeta(n + 1, ms + mf - n + 1) + return log1p(ms+mf) - log1p(ms) - log1p(mf) + + (logbeta(n + 1, ms + mf - n + 1) - logbeta(x + 1, ms - x + 1) - + logbeta(n - x + 1, mf - n + x + 1)) end - -function _hypercdf(ms, mf, n, x, invert) +_hypercdf(ms::Real, mf::Real, n::Real, x::Real, invert::Bool) = _hypercdf(promote(ms, mf, n, x)..., invert::Bool) +function _hypercdf(ms::T, mf::T, n::T, x::T, invert::Bool) where T<:Real N = ms+mf mode = fld((n + 1) * (ms + 1), N + 2) local result if x < mode - result = hyperpdf(ms, mf, n, x) - diff = result + result = diff = hyperpdf(ms, mf, n, x) lower_limit = max(0, n+ms-N) - while x != lower_limit && diff > (invert ? result : 1)*eps() + while x != lower_limit && diff > (invert ? result : 1)*eps(T) diff *= x * (N + x - n - ms) / ((1+n-x)*(1+ms -x)) result += diff x -= 1 @@ -45,10 +58,9 @@ function _hypercdf(ms, mf, n, x, invert) else invert = !invert x += 1 - result = hyperpdf(ms, mf, n, x) - diff = result + result = diff = hyperpdf(ms, mf, n, x) upper_limit = min(ms, n) - while x != upper_limit && diff > (invert ? 1 : result)*eps() + while x != upper_limit && diff > (invert ? 1 : result)*eps(T) diff *= (n-x) * (ms-x) / ((1+x) * (N+x+1-n-ms)) result += diff x += 1 diff --git a/test/rmath.jl b/test/rmath.jl index f808b115..ef5636e8 100644 --- a/test/rmath.jl +++ b/test/rmath.jl @@ -9,10 +9,19 @@ function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol) rv = @inferred(rmathfun(params..., a)) @test v isa float(Base.promote_typeof(params..., a)) @test rv isa float(Base.promote_typeof(params..., a)) - return if isprob - @test v ≈ rv rtol = rtol nans = true + # vbig == bigfloat version if available + vbig = try + statsfun(params..., big(a)) + catch + rv + end + # loose accuracy test to confirm shape is right + @test v ≈ rv nans = true + # tight accuracy test to confirm implimentation is correct numerically + if isprob + @test v ≈ vbig rtol = rtol nans = true else - @test v ≈ rv atol = rtol rtol = rtol nans = true + @test v ≈ vbig atol = rtol rtol = rtol nans = true end end From 492cd1f04b71732ee52d2095c2f88acfd20d90f8 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Thu, 2 Apr 2026 10:26:23 +0200 Subject: [PATCH 3/4] Fixup after conflict resolution --- src/distrs/hyper.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl index 2f6ae9a1..dbabedf7 100644 --- a/src/distrs/hyper.jl +++ b/src/distrs/hyper.jl @@ -78,6 +78,7 @@ end function hyperlogccdf(ms, mf, n, x) result, invert = _hypercdf(ms, mf, n, x, true) invert ? log1p(-result) : log(result) +end # Rmath implementations function hyperinvcdf(ms::Real, mf::Real, n::Real, q::Real) From 5566f0fa8c660c7b571006d53957ed0fff3cbc6e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 5 Apr 2026 11:28:04 +0200 Subject: [PATCH 4/4] Set timeout for CI jobs to 30 minutes --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 05f47886..b0d446fe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,6 +22,7 @@ jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} + timeout-minutes: 30 strategy: fail-fast: false matrix: