From 98c77292bfca3ada572225505289bae78d7e560c Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Sun, 22 Mar 2026 14:31:23 +0100
Subject: [PATCH] Implement betainvlogcdf/ccdf and gammainvlogcdf/ccdf in pure
 Julia

Replace the Rmath-dependent invlogcdf/invlogccdf functions for beta and
gamma distributions with Newton iteration in log-space. Uses asymptotic
initial estimates for extreme tails and -expm1(lp) for accuracy near
lp=0. Leverages existing betainvcdf/gammainvcdf for moderate cases.
---
 src/distrs/beta.jl  |  99 ++++++++++++++++++++++++++++++++---
 src/distrs/gamma.jl | 123 +++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 208 insertions(+), 14 deletions(-)

diff --git a/src/distrs/beta.jl b/src/distrs/beta.jl
index 79c6045..7285422 100644
--- a/src/distrs/beta.jl
+++ b/src/distrs/beta.jl
@@ -101,12 +101,97 @@ function betainvccdf(α::Real, β::Real, p::Real)
     return last(beta_inc_inv(β, α, p))
 end
 
-# Rmath implementations
-function betainvlogcdf(α::Real, β::Real, lq::Real)
-    T = float(Base.promote_typeof(α, β, lq))
-    return convert(T, Rmath.qbeta(lq, α, β, true, true))
+# Newton iteration in log-space for inverting log-CDF.
+# Solves betalogcdf(α, β, x) = lp for x.
+# Uses the identity: d/dx log(F(x)) = f(x)/F(x) = exp(logpdf - logcdf)
+# Newton step: x_{n+1} = x_n - (logcdf(x_n) - lp) / exp(logpdf(x_n) - logcdf(x_n))
+#            = x_n - (logcdf(x_n) - lp) * exp(logcdf(x_n) - logpdf(x_n))
+function _betainvlogcdf(α::Float64, β::Float64, lp::Float64)
+    if lp > -1
+        # Moderate lp: use -expm1(lp) = 1-exp(lp) for accuracy near lp=0
+        return Float64(betainvccdf(α, β, -expm1(lp)))
+    end
+
+    # Use Newton iteration in log-space.
+    # Initial estimate from small-x asymptotic: I_x(α,β) ≈ x^α / (α·B(α,β))
+    # ⟹ log(I_x) ≈ α·log(x) - log(α) - logbeta(α,β)
+    # ⟹ x₀ ≈ exp((lp + log(α) + logbeta(α,β)) / α)
+    logx = (lp + log(α) + logbeta(α, β)) / α
+    if logx < log(nextfloat(0.0))
+        return 0.0  # answer underflows Float64
+    end
+    x = exp(logx)
+    x = min(x, 1.0 - eps(1.0))
+
+    for _ in 1:200
+        lcdf = Float64(betalogcdf(α, β, x))
+        residual = lcdf - lp
+        if abs(residual) <= 2 * eps(max(abs(lp), 1.0))
+            break
+        end
+        lpdf = Float64(betalogpdf(α, β, x))
+        # Newton step: Δx = residual * F(x)/f(x) = residual * exp(logcdf - logpdf)
+        dx = residual * exp(lcdf - lpdf)
+        # Guard: don't let x go negative or jump too far
+        x_new = x - clamp(dx, -x / 2, x / 2)
+        x_new = clamp(x_new, nextfloat(0.0), 1.0 - eps(1.0))
+        if x_new == x
+            break
+        end
+        x = x_new
+    end
+
+    return x
+end
+
+function betainvlogcdf(α::Real, β::Real, lp::Real)
+    T = float(Base.promote_typeof(α, β, lp))
+    _lp = Float64(lp)
+
+    if isnan(_lp) || _lp > 0
+        return convert(T, NaN)
+    elseif _lp == 0
+        return one(T)
+    elseif isinf(_lp)  # -Inf
+        return zero(T)
+    end
+
+    # Handle degenerate cases
+    _α = Float64(α); _β = Float64(β)
+    if iszero(_α) && _β > 0
+        return zero(T)
+    elseif iszero(_β) && _α > 0
+        return _lp > log1p(-1.0) ? one(T) : zero(T)
+    end
+
+    return convert(T, _betainvlogcdf(_α, _β, _lp))
 end
-function betainvlogccdf(α::Real, β::Real, lq::Real)
-    T = float(Base.promote_typeof(α, β, lq))
-    return convert(T, Rmath.qbeta(lq, α, β, false, true))
+
+function betainvlogccdf(α::Real, β::Real, lp::Real)
+    T = float(Base.promote_typeof(α, β, lp))
+    _lp = Float64(lp)
+
+    if isnan(_lp) || _lp > 0
+        return convert(T, NaN)
+    elseif _lp == 0
+        return zero(T)
+    elseif isinf(_lp)  # -Inf
+        return one(T)
+    end
+
+    # Handle degenerate cases
+    _α = Float64(α); _β = Float64(β)
+    if iszero(_α) && _β > 0
+        return _lp > log1p(-1.0) ? zero(T) : one(T)
+    elseif iszero(_β) && _α > 0
+        return one(T)
+    end
+
+    # For moderate lp (near 0): use -expm1(lp) = 1-exp(lp) for accuracy
+    if _lp > -1
+        return convert(T, Float64(betainvcdf(_α, _β, -expm1(_lp))))
+    end
+
+    # For very negative lp: betainvlogccdf(α, β, lp) = 1 - betainvlogcdf(β, α, lp)
+    return convert(T, 1.0 - _betainvlogcdf(_β, _α, _lp))
 end
diff --git a/src/distrs/gamma.jl b/src/distrs/gamma.jl
index cee2fca..5206a0c 100644
--- a/src/distrs/gamma.jl
+++ b/src/distrs/gamma.jl
@@ -95,12 +95,121 @@ function gammainvccdf(k::Real, θ::Real, p::Real)
     return _θ * gamma_inc_inv(_k, 1 - _p, _p)
 end
 
-# Rmath implementations
-function gammainvlogcdf(k::Real, θ::Real, lq::Real)
-    T = float(Base.promote_typeof(k, θ, lq))
-    return convert(T, Rmath.qgamma(lq, k, θ, true, true))
+# Newton iteration in log-space for inverting log-CDF.
+# Solves gammalogcdf(k, θ, x) = lp for x.
+function _gammainvlogcdf(k::Float64, θ::Float64, lp::Float64)
+    if lp > -1
+        # Use -expm1(lp) = 1-exp(lp) for accuracy near lp=0
+        return Float64(gammainvccdf(k, θ, -expm1(lp)))
+    end
+
+    # Use Newton iteration in log-space.
+    # Initial estimate from small-x asymptotic:
+    # P(k, x/θ) ≈ (x/θ)^k / (k·Γ(k)) for small x/θ
+    # ⟹ x₀/θ ≈ exp((lp + log(k) + loggamma(k)) / k)
+    logz = (lp + log(k) + loggamma(k)) / k
+    if logz + log(θ) < log(nextfloat(0.0))
+        return 0.0  # answer underflows Float64
+    end
+    z = exp(logz)
+    x = z * θ
+    x = max(x, nextfloat(0.0))
+
+    for _ in 1:200
+        lcdf = Float64(gammalogcdf(k, θ, x))
+        residual = lcdf - lp
+        if abs(residual) <= 2 * eps(max(abs(lp), 1.0))
+            break
+        end
+        lpdf = Float64(gammalogpdf(k, θ, x))
+        # Newton step: Δx = residual * F(x)/f(x) = residual * exp(logcdf - logpdf)
+        dx = residual * exp(lcdf - lpdf)
+        # Guard: don't let x go negative or jump too far
+        x_new = x - clamp(dx, -x / 2, x / 2)
+        x_new = max(x_new, nextfloat(0.0))
+        if x_new == x
+            break
+        end
+        x = x_new
+    end
+
+    return x
+end
+
+# Newton iteration in log-space for inverting log-CCDF.
+# Solves gammalogccdf(k, θ, x) = lp for x.
+function _gammainvlogccdf(k::Float64, θ::Float64, lp::Float64)
+    if lp > -1
+        # Use -expm1(lp) = 1-exp(lp) for accuracy near lp=0
+        return Float64(gammainvcdf(k, θ, -expm1(lp)))
+    end
+
+    # Extreme right tail: initial estimate
+    # Q(k, x/θ) ≈ (x/θ)^(k-1)·exp(-x/θ)/Γ(k) for large x/θ
+    # ⟹ log(Q) ≈ (k-1)·log(x/θ) - x/θ - loggamma(k)
+    # For large x/θ, the -x/θ term dominates: x₀/θ ≈ -lp - loggamma(k)
+    z = max(-lp - loggamma(k), 1.0)
+    x = z * θ
+
+    for _ in 1:200
+        lccdf = Float64(gammalogccdf(k, θ, x))
+        residual = lccdf - lp
+        if abs(residual) <= 2 * eps(max(abs(lp), 1.0))
+            break
+        end
+        lpdf = Float64(gammalogpdf(k, θ, x))
+        # d/dx log(ccdf(x)) = -f(x)/ccdf(x) = -exp(logpdf - logccdf)
+        # Newton step: Δx = (logccdf - lp) / (-exp(logpdf - logccdf))
+        #            = -(logccdf - lp) * exp(logccdf - logpdf)
+        dx = -residual * exp(lccdf - lpdf)
+        # Guard: don't overshoot
+        x_new = x - clamp(dx, -x, x / 2)
+        x_new = max(x_new, nextfloat(0.0))
+        if x_new == x
+            break
+        end
+        x = x_new
+    end
+
+    return x
 end
-function gammainvlogccdf(k::Real, θ::Real, lq::Real)
-    T = float(Base.promote_typeof(k, θ, lq))
-    return convert(T, Rmath.qgamma(lq, k, θ, false, true))
+
+function gammainvlogcdf(k::Real, θ::Real, lp::Real)
+    T = float(Base.promote_typeof(k, θ, lp))
+    _lp = Float64(lp)
+
+    if isnan(_lp) || _lp > 0
+        return convert(T, NaN)
+    elseif _lp == 0
+        return convert(T, Inf)
+    elseif isinf(_lp)  # -Inf
+        return zero(T)
+    end
+
+    _k = Float64(k); _θ = Float64(θ)
+    if iszero(_k)
+        return _lp >= 0 ? zero(T) : convert(T, NaN)
+    end
+
+    return convert(T, _gammainvlogcdf(_k, _θ, _lp))
+end
+
+function gammainvlogccdf(k::Real, θ::Real, lp::Real)
+    T = float(Base.promote_typeof(k, θ, lp))
+    _lp = Float64(lp)
+
+    if isnan(_lp) || _lp > 0
+        return convert(T, NaN)
+    elseif _lp == 0
+        return zero(T)
+    elseif isinf(_lp)  # -Inf
+        return convert(T, Inf)
+    end
+
+    _k = Float64(k); _θ = Float64(θ)
+    if iszero(_k)
+        return _lp >= 0 ? zero(T) : convert(T, NaN)
+    end
+
+    return convert(T, _gammainvlogccdf(_k, _θ, _lp))
 end