safer int for binomial loglik (#504)

* safer Int in Binomial log likelihood

* patch bump

* fallback method

* drop extra tests

* s/_safer/_safe/

Co-authored-by: Dave Kleinschmidt <[email protected]>
Co-authored-by: Alex Arslan <[email protected]>
Co-authored-by: Milan Bouchet-Valat <[email protected]>

Author: 4 people

Project.toml

@@ -1,6 +1,6 @@
 name = "GLM"
 uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
-version = "1.8.0"
+version = "1.8.1"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/glmtools.jl

@@ -499,6 +499,19 @@ false
 dispersion_parameter(D) = true
 dispersion_parameter(::Union{Bernoulli, Binomial, Poisson}) = false
 
+"""
+    _safe_int(x::T)
+    
+Convert to Int, when `x` is within 1 eps of an integer.
+"""
+function _safe_int(x::T) where {T<:AbstractFloat}
+    r = round(Int, x)
+    abs(x - r) <= eps(x) && return r
+    throw(InexactError(nameof(T), T, x))
+end
+
+_safe_int(x) = Int(x)
+
 """
     GLM.loglik_obs(D, y, μ, wt, ϕ)
 
@@ -512,7 +525,7 @@ The loglikelihood of a fitted model is the sum of these values over all the obse
 function loglik_obs end
 
 loglik_obs(::Bernoulli, y, μ, wt, ϕ) = wt*logpdf(Bernoulli(μ), y)
-loglik_obs(::Binomial, y, μ, wt, ϕ) = logpdf(Binomial(Int(wt), μ), Int(y*wt))
+loglik_obs(::Binomial, y, μ, wt, ϕ) = logpdf(Binomial(Int(wt), μ), _safe_int(y*wt))
 loglik_obs(::Gamma, y, μ, wt, ϕ) = wt*logpdf(Gamma(inv(ϕ), μ*ϕ), y)
 # In Distributions.jl, a Geometric distribution characterizes the number of failures before 
 # the first success in a sequence of independent Bernoulli trials with success rate p.

test/runtests.jl

@@ -1591,3 +1591,15 @@ end
     end
 end
 
+    
+@testset "Floating point error in Binomial loglik" begin
+    @test_throws InexactError GLM._safe_int(1.3)
+    @test GLM._safe_int(1) === 1
+    # see issue 503
+    y, μ, wt, ϕ = 0.6376811594202898, 0.8492925285671102, 69.0, NaN
+    # due to floating point:
+    # 1. y * wt == 43.99999999999999 
+    # 2. 44 / y == wt
+    # 3. 44 / wt == y
+    @test GLM.loglik_obs(Binomial(), y, μ, wt, ϕ) ≈ GLM.logpdf(Binomial(Int(wt), μ), 44)
+end