From 82bb33b585b3869094e1b7d47c62f59b0640d402 Mon Sep 17 00:00:00 2001 From: Jameson Quinn Date: Fri, 12 Feb 2021 12:39:44 -0500 Subject: [PATCH 1/2] 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. --- Project.toml | 1 + src/modeling_library/distributions/beta.jl | 10 +++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c8111b3a..ce0fdb2e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Marco Cusumano-Towner "] 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" diff --git a/src/modeling_library/distributions/beta.jl b/src/modeling_library/distributions/beta.jl index 012044fd..1503f6d7 100644 --- a/src/modeling_library/distributions/beta.jl +++ b/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) From dbcb12404c3041fb22a90aa0cd5f43bc63f9b71f Mon Sep 17 00:00:00 2001 From: Jameson Quinn Date: Fri, 12 Feb 2021 12:50:17 -0500 Subject: [PATCH 2/2] test cases --- test/modeling_library/distributions.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/modeling_library/distributions.jl b/test/modeling_library/distributions.jl index 9b894987..ae6d67b0 100644 --- a/test/modeling_library/distributions.jl +++ b/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)