Merge branch 'master' into dw/version

.github/workflows/ci.yml

@@ -29,6 +29,7 @@ jobs:
           - '1.3' # oldest supported version
           - '1.6' # LTS
           - '1' # latest release
+          - 'nightly'
         os:
           - ubuntu-latest
           - macos-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.1.1"
+version = "1.3.1"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -12,22 +12,32 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
+[weakdeps]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+
+[extensions]
+StatsFunsChainRulesCoreExt = "ChainRulesCore"
+StatsFunsInverseFunctionsExt = "InverseFunctions"
+
 [compat]
 ChainRulesCore = "1"
 HypergeometricFunctions = "0.3.10"
 InverseFunctions = "0.1"
-IrrationalConstants = "0.1"
+IrrationalConstants = "0.1, 0.2"
 LogExpFunctions = "0.3.26"
 Reexport = "1"
 Rmath = "0.6.1, 0.7"
 SpecialFunctions = "1.8.4, 2.1.4"
 julia = "1.3"
 
 [extras]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ChainRulesTestUtils", "ForwardDiff", "Random", "Test"]
+test = ["ChainRulesCore", "ChainRulesTestUtils", "ForwardDiff", "InverseFunctions", "Random", "Test"]

src/chainrules.jl → ext/StatsFunsChainRulesCoreExt.jl

@@ -1,3 +1,9 @@
+module StatsFunsChainRulesCoreExt
+
+using StatsFuns
+using StatsFuns: digamma
+import ChainRulesCore
+
 ChainRulesCore.@scalar_rule(
     betalogpdf(α::Real, β::Real, x::Number),
     @setup(z = digamma(α + β)),
@@ -76,3 +82,5 @@ ChainRulesCore.@scalar_rule(
         - x * b,
     ),
 )
+
+end # module

src/inverse.jl → ext/StatsFunsInverseFunctionsExt.jl

@@ -1,3 +1,8 @@
+module StatsFunsInverseFunctionsExt
+
+using StatsFuns
+import InverseFunctions
+
 InverseFunctions.inverse(::typeof(normcdf)) = norminvcdf
 InverseFunctions.inverse(::typeof(norminvcdf)) = normcdf
 
@@ -9,3 +14,5 @@ InverseFunctions.inverse(::typeof(norminvlogcdf)) = normlogcdf
 
 InverseFunctions.inverse(::typeof(normlogccdf)) = norminvlogccdf
 InverseFunctions.inverse(::typeof(norminvlogccdf)) = normlogccdf
+
+end # module

src/StatsFuns.jl

@@ -3,8 +3,6 @@ module StatsFuns
 using Base: Math.@horner
 using Reexport
 using SpecialFunctions
-import ChainRulesCore
-import InverseFunctions
 
 # reexports
 @reexport using IrrationalConstants:
@@ -262,7 +260,9 @@ include(joinpath("distrs", "pois.jl"))
 include(joinpath("distrs", "tdist.jl"))
 include(joinpath("distrs", "srdist.jl"))
 
-include("chainrules.jl")
-include("inverse.jl")
+if !isdefined(Base, :get_extension)
+    include("../ext/StatsFunsChainRulesCoreExt.jl")
+    include("../ext/StatsFunsInverseFunctionsExt.jl")
+end
 
 end # module

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

src/distrs/gamma.jl

@@ -22,7 +22,12 @@ gammalogpdf(k::Real, θ::Real, x::Real) = gammalogpdf(promote(k, θ, x)...)
 function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
     # we ensure that `log(x)` does not error if `x < 0`
     xθ = max(x, 0) / θ
-    val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
+    val = -loggamma(k) - log(θ) - xθ
+    # xlogy(k - 1, xθ) - xθ -> -∞ for xθ -> ∞ so we only add the first term
+    # when it's safe
+    if isfinite(xθ)
+        val += xlogy(k - 1, xθ)
+    end
     return x < 0 ? oftype(val, -Inf) : val
 end
 

src/tvpack.jl

@@ -427,7 +427,7 @@ function bvnuppercdf(dh::T, dk::T, r::T)::T where {T<:Float64}
         elseif dk == typemin(T)
             return normcdf(-dh)
         else
-            thow(error())
+            throw(error())
         end
     end
 end

test/generic.jl

@@ -1,38 +1,35 @@
 using StatsFuns
 using StatsFuns: RFunctions
 using ForwardDiff: Dual
+using Test
+
+include("utils.jl")
 
 function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
     v = @inferred(rmathfun(params..., a))
     rv = @inferred(statsfun(params..., Dual(a))).value
+    @test v isa float(Base.promote_typeof(params..., a))
+    @test rv isa float(Base.promote_typeof(params..., a))
     if isprob
-        rd = abs(v / rv - 1.0)
-        if rd > rtol
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v/rv - 1| = $rd > $rtol.")
-        end
+        @test v ≈ rv rtol=rtol nans=true
     else
-        τ = (1.0 + abs(rv)) * rtol
-        ad = abs(v - rv)
-        if ad > τ
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v - rv| = $ad > $τ.")
-        end
+        @test v ≈ rv atol=rtol rtol=rtol nans=true
     end
 end
 
-function genericcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
+function genericcomp(basename, params, X::AbstractArray, rtol=_default_rtol(params, X))
   pdf = string(basename, "pdf")
   logpdf = string(basename, "logpdf")
   stats_pdf = eval(Symbol(pdf))
   stats_logpdf = eval(Symbol(logpdf))
   rmath_pdf = eval(Meta.parse(string("RFunctions.", pdf)))
   rmath_logpdf = eval(Meta.parse(string("RFunctions.", logpdf)))
-  for i = 1:length(X)
-    x = X[i]
+
+  @testset "pdf with params=$params, x=$x" for x in X
     check_rmath(pdf, stats_pdf, rmath_pdf, params, "x", x, true, rtol)
+  end
+
+  @testset "logpdf with params=$params, x=$x" for x in X
     check_rmath(logpdf, stats_logpdf, rmath_logpdf, params, "x", x, false, rtol)
   end
 end
@@ -121,8 +118,11 @@ end
 
     genericcomp_tests("pois", [
         ((1.0,), 0:30),
-        ((10.0,), 0:42)
+        ((10.0,), 0:37),
+        ((10.0,), 39:42),
     ])
+    # Requires slightly larger tolerance: #157
+    genericcomp("pois", (10.0,), 38:38, 2.5e-14)
 
     genericcomp_tests("tdist", [
         ((1,), -5.0:0.1:5.0),

test/misc.jl

@@ -1,4 +1,5 @@
 using SpecialFunctions, StatsFuns
+using Test
 
 @testset "logmvgamma" begin
     @testset "type behavior" for eltya in (Float32, Float64)
@@ -50,7 +51,7 @@ end
                                       for eltyb in (Float32, Float64)
             a = rand(eltya)
             b = rand(eltyb)
-            T = Base.promote_eltype(eltya, eltyb)
+            T = Base.promote_type(eltya, eltyb)
             @test typeof(logmvbeta(1, a, b)) == T
         end
     end
@@ -90,3 +91,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

@@ -2,9 +2,20 @@ using StatsFuns
 using StatsFuns: RFunctions
 using Test
 
+include("utils.jl")
+
 function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
-    v = @inferred(statsfun(params..., a))
-    rv = @inferred(rmathfun(params..., a))
+    # HypergeometricFunctions@0.3.18 made some changes that trips inference
+    # on older versions of Julia. Not worth the the time trying to debug/fix
+    # as we might drop support for these version soon so just putting the
+    # inference check behind a VERSION branch
+    if VERSION >= v"1.7"
+        v = @inferred(statsfun(params..., a))
+        rv = @inferred(rmathfun(params..., a))
+    else
+        v = statsfun(params..., a)
+        rv = rmathfun(params..., a)
+    end
     @test v isa float(Base.promote_typeof(params..., a))
     @test rv isa float(Base.promote_typeof(params..., a))
     if isprob
@@ -17,15 +28,7 @@ end
 get_statsfun(fname) = eval(Symbol(fname))
 get_rmathfun(fname) = eval(Meta.parse(string("RFunctions.", fname)))
 
-function rmathcomp(basename, params, X::AbstractArray)
-    # compute default tolerance:
-    # has to take into account `params` as well since otherwise e.g. `X::Array{<:Rational}`
-    # always uses a tolerance based on `eps(one(Float64))` even when parameters are of type
-    # Float32
-    rtol = 100 * eps(float(one(promote_type(Base.promote_typeof(params...), eltype(X)))))
-    rmathcomp(basename, params, X, rtol)
-end
-function rmathcomp(basename, params, X::AbstractArray, rtol)
+function rmathcomp(basename, params, X::AbstractArray, rtol=_default_rtol(params, X))
     # tackle pdf specially
     has_pdf = true
     if basename == "srdist"
@@ -177,6 +180,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),
@@ -226,6 +264,7 @@ end
         ((Float16(1), Float16(1)), (Float16(0):Float16(0.05):Float16(12))),
         ((1f0, 1f0), (Float16(0):Float16(0.05):Float16(12))),
         ((2, 3), (0//1:12//1)),
+        ((3.0, 1e-310), (0.0:0.05:12.0))
     ])
 
     rmathcomp_tests("hyper", [