Fix quantiles of degenerate beta distributions

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),