issue 350: avoid logpdf=inf for beta

Adds two special cases for logpdf(::Beta...), when x is 0 or 1. Effectively, uses (machine epsilon / e) or one minus that instead of 0 or 1.

Author: jamesonquinn

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)