Fix underflow of logccdf of PGeneralizedGaussian (#1932)

Author: devmotion

src/univariate/continuous/pgeneralizedgaussian.jl

@@ -108,17 +108,37 @@ end
 # The CDF of the Gamma distribution provides this, with the necessary 1/Γ(a) normalization.
 function cdf(d::PGeneralizedGaussian, x::Real)
     μ, α, p = params(d)
-    v = cdf(Gamma(inv(p), 1), (abs(x - μ) / α)^p)
-    return (1 + copysign(v, x - μ)) / 2
+    return _normcdf_pgeneralizedgaussian(p, (x - μ) / α)
 end
+function ccdf(d::PGeneralizedGaussian, x::Real)
+    μ, α, p = params(d)
+    return _normcdf_pgeneralizedgaussian(p, (μ - x) / α)
+end
+function _normcdf_pgeneralizedgaussian(p::Real, z::Real)
+    r = abs(z)^p
+    d = Gamma(inv(p), 1)
+    if z < 0
+        return ccdf(d, r) / 2
+    else
+        return (1 + cdf(d, r)) / 2
+    end
+end
+
 function logcdf(d::PGeneralizedGaussian, x::Real)
     μ, α, p = params(d)
-    Δ = x - μ
-    logv = logcdf(Gamma(inv(p), 1), (abs(Δ) / α)^p)
-    if Δ < 0
-        return log1mexp(logv) - logtwo
+    return _normlogcdf_pgeneralizedgaussian(p, (x - μ) / α)
+end
+function logccdf(d::PGeneralizedGaussian, x::Real)
+    μ, α, p = params(d)
+    return _normlogcdf_pgeneralizedgaussian(p, (μ - x) / α)
+end
+function _normlogcdf_pgeneralizedgaussian(p::Real, z::Real)
+    r = abs(z)^p
+    d = Gamma(inv(p), 1)
+    if z < 0
+        return logccdf(d, r) - logtwo
     else
-        return log1pexp(logv) - logtwo
+        return log1pexp(logcdf(d, r)) - logtwo
     end
 end
 

test/univariate/continuous/pgeneralizedgaussian.jl

@@ -106,4 +106,20 @@ using Test
         # Additional tests, including sampling
         test_distr(d, 10^6)
     end
+
+    # issue #1916
+    @testset "Underflow in logccdf" begin
+        d = PGeneralizedGaussian(2.5)
+
+        # (c)cdf became more accurate
+        @test isfinite(log(cdf(d, -5.0)))
+        @test isfinite(log(ccdf(d, 5.0)))
+        @test ccdf(d, 5.0) ≈ cdf(d, -5.0)
+
+        # We have to use -20 to provoke the issue
+        @test log(cdf(d, -20.0)) === -Inf
+        @test log(ccdf(d, 20.0)) === -Inf
+        @test isfinite(logcdf(d, -20.0))
+        @test logccdf(d, 20.0) ≈ logcdf(d, -20.0)
+    end
 end