Handle tiny θ in gammalogpdf

JuliaStats · Feb 13, 2024 · f4219ac · f4219ac
1 parent 0079530
commit f4219ac
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 1 deletion.
diff --git a/src/distrs/gamma.jl b/src/distrs/gamma.jl
@@ -22,7 +22,12 @@ gammalogpdf(k::Real, θ::Real, x::Real) = gammalogpdf(promote(k, θ, x)...)
 function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
     # we ensure that `log(x)` does not error if `x < 0`
     xθ = max(x, 0) / θ
-    val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
+    val = -loggamma(k) - log(θ) - xθ
+    # xlogy(k - 1, xθ) - xθ -> -∞ for xθ -> ∞ so we only add the first term
+    # when it's safe
+    if isfinite(xθ)
+        val += xlogy(k - 1, xθ)
+    end
     return x < 0 ? oftype(val, -Inf) : val
 end
 

diff --git a/test/rmath.jl b/test/rmath.jl
@@ -255,6 +255,7 @@ end
         ((Float16(1), Float16(1)), (Float16(0):Float16(0.05):Float16(12))),
         ((1f0, 1f0), (Float16(0):Float16(0.05):Float16(12))),
         ((2, 3), (0//1:12//1)),
+        ((3.0, 1e-310), (0.0:0.05:12.0))
     ])
 
     rmathcomp_tests("hyper", [