Fix quantiles of degenerate beta distributions (#158)

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/distrs/beta.jl

@@ -27,18 +27,22 @@ function betalogpdf(α::T, β::T, x::T) where {T<:Real}
 end
 
 function betacdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
+    # Handle degenerate cases
     if iszero(α) && β > 0
         return float(last(promote(α, β, x, x >= 0)))
+    elseif iszero(β) && α > 0
+        return float(last(promote(α, β, x, x >= 1)))
     end
 
     return first(beta_inc(α, β, clamp(x, 0, 1)))
 end
 
 function betaccdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
+    # Handle degenerate cases
     if iszero(α) && β > 0
         return float(last(promote(α, β, x, x < 0)))
+    elseif iszero(β) && α > 0
+        return float(last(promote(α, β, x, x < 1)))
     end
 
     last(beta_inc(α, β, clamp(x, 0, 1)))
@@ -47,9 +51,11 @@ end
 # The log version is currently based on non-log version. When the cdf is very small we shift
 # to an implementation based on the hypergeometric function ₂F₁ to avoid underflow.
 function betalogcdf(α::T, β::T, x::T) where {T<:Real}
-    # Handle a degenerate case
+    # Handle degenerate cases
     if iszero(α) && β > 0
         return log(last(promote(x, x >= 0)))
+    elseif iszero(β) && α > 0
+        return log(last(promote(x, x >= 1)))
     end
 
     _x = clamp(x, 0, 1)
@@ -67,9 +73,11 @@ end
 betalogcdf(α::Real, β::Real, x::Real) = betalogcdf(promote(α, β, x)...)
 
 function betalogccdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
+    # Handle degenerate cases
     if iszero(α) && β > 0
         return log(last(promote(α, β, x, x < 0)))
+    elseif iszero(β) && α > 0
+        return log(last(promote(α, β, x, x < 1)))
     end
 
     p, q = beta_inc(α, β, clamp(x, 0, 1))
@@ -80,6 +88,28 @@ function betalogccdf(α::Real, β::Real, x::Real)
     end
 end
 
-betainvcdf(α::Real, β::Real, p::Real) = first(beta_inc_inv(α, β, p))
+function betainvcdf(α::Real, β::Real, p::Real)
+    # Handle degenerate cases
+    if 0 ≤ p ≤ 1
+        if iszero(α) && β > 0
+            return last(promote(α, β, p, false))
+        elseif iszero(β) && α > 0
+            return last(promote(α, β, p, p > 0))
+        end
+    end
+
+    return first(beta_inc_inv(α, β, p))
+end
 
-betainvccdf(α::Real, β::Real, p::Real) = last(beta_inc_inv(β, α, p))
+function betainvccdf(α::Real, β::Real, p::Real)
+    # Handle degenerate cases
+    if 0 ≤ p ≤ 1
+        if iszero(α) && β > 0
+            return last(promote(α, β, p, p == 0))
+        elseif iszero(β) && α > 0
+            return last(promote(α, β, p, true))
+        end
+    end
+
+    return last(beta_inc_inv(β, α, p))
+end

test/misc.jl

@@ -96,3 +96,8 @@ end
 @testset "gammalogcdf: numerical issue" begin
     @test gammalogcdf(42648.50647826826, 2.2498007956420723e-5, 0.6991377135675367) ≈ -1933.2698959040617410
 end
+
+# https://github.com/JuliaStats/StatsFuns.jl/issues/154
+@testset "tvdistinvcdf: numerical issue" begin
+    @test isnan(@inferred(tdistinvcdf(0, 0.975)))
+end

test/rmath.jl

@@ -171,6 +171,41 @@ end
         x in (0.0, 0.5, 1.0)
         @test_throws DomainError f(0.0, 0.0, x)
     end
+    # Have to check them separately since Rmath is not completely consistent with our conventions:
+    # - betacdf(α, β, x) = P(X ≤ x) where X ~ Beta(α, β)
+    # - betaccdf(α, β, x) = P(X > x) where X ~ Beta(α, β)
+    # - betainvcdf(α, β, p) = inf { x : p ≤ P(X ≤ x) } where X ~ Beta(α, β)
+    # - betainvccdf(α, β, p) = sup { x : p ≤ P(X > x) } where X ~ Beta(α, β)
+    @testset "beta: degenerate cases" begin
+        # Check degenerate cases Beta(0, β) and Dirac(α, 0) with α, β > 0
+        # Beta(0, β) is a Dirac distribution at x=0
+        # Beta(α, 0) is a Dirac distribution at x=1
+        α = β =  1//2
+
+        for x in 0f0:0.01f0:1f0
+            # Check betacdf
+            @test @inferred(betacdf(0, β, x)) === 1f0
+            @test @inferred(betacdf(α, 0, x)) === (x < 1 ? 0f0 : 1f0)
+
+            # Check betaccdf, betalogcdf, and betalogccdf based on betacdf
+            @test @inferred(betaccdf(0, β, x)) === 1 - betacdf(0, β, x)
+            @test @inferred(betaccdf(α, 0, x)) === 1 - betacdf(α, 0, x)
+            @test @inferred(betalogcdf(0, β, x)) === log(betacdf(0, β, x))
+            @test @inferred(betalogcdf(α, 0, x)) === log(betacdf(α, 0, x))
+            @test @inferred(betalogccdf(0, β, x)) === log(betaccdf(0, β, x))
+            @test @inferred(betalogccdf(α, 0, x)) === log(betaccdf(α, 0, x))
+        end
+
+        for p in 0f0:0.01f0:1f0
+            # Check betainvcdf
+            @test @inferred(betainvcdf(0, β, p)) === 0f0
+            @test @inferred(betainvcdf(α, 0, p)) === (p > 0 ? 1f0 : 0f0)
+
+            # Check betainvccdf
+            @test @inferred(betainvccdf(0, β, p)) === (p > 0 ? 0f0 : 1f0)
+            @test @inferred(betainvccdf(α, 0, p)) === 1f0
+        end
+    end
 
     rmathcomp_tests("binom", [
         ((1, 0.5), 0.0:1.0),