diff --git a/Project.toml b/Project.toml index d35a77b..29ce15b 100644 --- a/Project.toml +++ b/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" diff --git a/src/distrs/beta.jl b/src/distrs/beta.jl index 97a6a61..0276c57 100644 --- a/src/distrs/beta.jl +++ b/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 diff --git a/test/misc.jl b/test/misc.jl index d3b5d96..428c195 100644 --- a/test/misc.jl +++ b/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 diff --git a/test/rmath.jl b/test/rmath.jl index e650bdb..df74d1c 100644 --- a/test/rmath.jl +++ b/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),