Fix InverseGaussian cdf overflows (#1882)

* Fix inversegaussian cdf overflows * Update src/univariate/continuous/inversegaussian.jl Co-authored-by: David Widmann <[email protected]> * Update src/univariate/continuous/inversegaussian.jl Co-authored-by: David Widmann <[email protected]> * Fix missing parentheses * Move parentheses * Move new tests to dedicated inversegaussian file * Add ccdf tests --------- Co-authored-by: David Widmann <[email protected]>
JuliaStats · Aug 9, 2024 · 6e0f1dc · 6e0f1dc
1 parent 47c040b
commit 6e0f1dc
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 3 deletions.
diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl
@@ -99,7 +99,10 @@ function cdf(d::InverseGaussian, x::Real)
     y = max(x, 0)
     u = sqrt(λ / y)
     v = y / μ
-    z = normcdf(u * (v - 1)) + exp(2λ / μ) * normcdf(-u * (v + 1))
+    # 2λ/μ and normlogcdf(-u*(v+1)) are similar magnitude, opp. sign
+    # truncating to [0, 1] as an additional precaution
+    # Ref https://github.com/JuliaStats/Distributions.jl/issues/1873
+    z = clamp(normcdf(u * (v - 1)) + exp(2λ / μ + normlogcdf(-u * (v + 1))), 0, 1)
 
     # otherwise `NaN` is returned for `+Inf`
     return isinf(x) && x > 0 ? one(z) : z
@@ -110,7 +113,10 @@ function ccdf(d::InverseGaussian, x::Real)
     y = max(x, 0)
     u = sqrt(λ / y)
     v = y / μ
-    z = normccdf(u * (v - 1)) - exp(2λ / μ) * normcdf(-u * (v + 1))
+    # 2λ/μ and normlogcdf(-u*(v+1)) are similar magnitude, opp. sign
+    # truncating to [0, 1] as an additional precaution
+    # Ref https://github.com/JuliaStats/Distributions.jl/issues/1873
+    z = clamp(normccdf(u * (v - 1)) - exp(2λ / μ + normlogcdf(-u * (v + 1))), 0, 1)
 
     # otherwise `NaN` is returned for `+Inf`
     return isinf(x) && x > 0 ? zero(z) : z

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -96,6 +96,7 @@ const tests = [
     "eachvariate",
     "univariate/continuous/triangular",
     "statsapi",
+    "univariate/continuous/inversegaussian",
 
     ### missing files compared to /src:
     # "common",
@@ -137,7 +138,6 @@ const tests = [
     # "univariate/continuous/generalizedextremevalue",
     # "univariate/continuous/generalizedpareto",
     # "univariate/continuous/inversegamma",
-    # "univariate/continuous/inversegaussian",
     # "univariate/continuous/ksdist",
     # "univariate/continuous/ksonesided",
     # "univariate/continuous/levy",

diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl
@@ -0,0 +1,18 @@
+@testset "InverseGaussian cdf outside of [0, 1] (#1873)" begin
+    for d in [
+        InverseGaussian(1.65, 590),
+        InverseGaussian(0.5, 1000)
+    ]
+        for x in [0.02, 1.0, 20.0, 300.0]
+            p = cdf(d, x)
+            @test 0.0 <= p <= 1.0
+            @test p ≈ exp(logcdf(d, x))
+
+            q = ccdf(d, x)
+            @test 0.0 <= q <= 1.0
+            @test q ≈ exp(logccdf(d, x))
+
+            @test (p + q) ≈ 1
+        end
+    end
+end