Merge pull request #359 from jamesonquinn/BetaDistInfIssue350

Avoid Inf from logpdf(::Beta,...) (issue #350)

Project.toml

@@ -4,6 +4,7 @@ authors = ["Marco Cusumano-Towner <[email protected]>"]
 version = "0.4.1"
 
 [deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/modeling_library/distributions/beta.jl

@@ -12,7 +12,15 @@ const beta = Beta()
 
 function logpdf(::Beta, x::Real, alpha::Real, beta::Real)
     (x < 0 || x > 1 ? -Inf :
-    (alpha - 1) * log(x) + (beta - 1) * log1p(-x) - logbeta(alpha, beta) )
+        ((x == 0 && alpha <= 1) ?      #Special case to avoid infinite PDF at 0
+            (alpha - 1) * log(eps(typeof(x))) -
+            alpha + 1 + #as if we were using x=(eps / e); even smaller than eps
+            (beta - 1) * log1p(-eps(typeof(x))) - logbeta(alpha, beta) :
+        ((x == 1 && beta <= 1) ?       #Special case to avoid infinite PDF at 1
+            (alpha - 1) * log1p(-eps(typeof(x))) -
+            beta + 1 + #as if we were using 1-x=(eps / e)
+            (beta - 1) * log(eps(typeof(x))) - logbeta(alpha, beta) :
+    (alpha - 1) * log(x) + (beta - 1) * log1p(-x) - logbeta(alpha, beta) )))
 end
 
 function logpdf_grad(::Beta, x::Real, alpha::Real, beta::Real)

test/modeling_library/distributions.jl

@@ -24,6 +24,12 @@ end
     # out of support
     @test logpdf(beta, -1, 0.5, 0.5) == -Inf
 
+    # avoid infinities, sanely
+    @test logpdf(beta, eps(typeof(0.)), 0.5, 0.5) < logpdf(beta, 0., 0.5, 0.5) < Inf
+    @test logpdf(beta, eps(typeof(0.)), 0.5, 1.5) < logpdf(beta, 0., 0.5, 1.5) < Inf
+    @test logpdf(beta, 1-eps(typeof(0.)), 0.5, 0.5) < logpdf(beta, 1., 0.5, 0.5) < Inf
+    @test logpdf(beta, 1-eps(typeof(0.)), 1.5, 0.5) < logpdf(beta, 1., 1.5, 0.5) < Inf
+
     # logpdf_grad
     f = (x, alpha, beta_param) -> logpdf(beta, x, alpha, beta_param)
     args = (0.4, 0.2, 0.3)