Fix Julia 1.0 compatibility (#138)

* Revert "Use julia implementations for pdfs and some cdf-like functions (#113)"

This reverts commit 28d732a.

* Fix test errors on Julia 1.0

* Run CI on Julia 1.0

* Bump version

* Update Project.toml

.github/workflows/ci.yml

@@ -26,6 +26,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
+          - '1.0'
           - '1'
         os:
           - ubuntu-latest

Project.toml

@@ -1,10 +1,9 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.17"
+version = "0.9.18"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
@@ -14,7 +13,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 ChainRulesCore = "1"
-HypergeometricFunctions = "0.3"
 InverseFunctions = "0.1"
 IrrationalConstants = "0.1"
 LogExpFunctions = "0.3.2"

src/distrs/beta.jl

@@ -1,17 +1,14 @@
 # functions related to beta distributions
 
-using HypergeometricFunctions: _₂F₁
-
 # R implementations
+# For pdf and logpdf we use the Julia implementation
 using .RFunctions:
-    # betapdf,
-    # betalogpdf,
-    # betacdf,
-    # betaccdf,
-    # betalogcdf,
-    # betalogccdf,
-    # betainvcdf,
-    # betainvccdf,
+    betacdf,
+    betaccdf,
+    betalogcdf,
+    betalogccdf,
+    betainvcdf,
+    betainvccdf,
     betainvlogcdf,
     betainvlogccdf
 
@@ -25,60 +22,3 @@ function betalogpdf(α::T, β::T, x::T) where {T<:Real}
     val = xlogy(α - 1, y) + xlog1py(β - 1, -y) - logbeta(α, β)
     return x < 0 || x > 1 ? oftype(val, -Inf) : val
 end
-
-function betacdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
-    if iszero(α) && β > 0
-        return float(last(promote(α, β, x, x >= 0)))
-    end
-
-    return first(beta_inc(α, β, clamp(x, 0, 1)))
-end
-
-function betaccdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
-    if iszero(α) && β > 0
-        return float(last(promote(α, β, x, x < 0)))
-    end
-
-    last(beta_inc(α, β, clamp(x, 0, 1)))
-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
-    if iszero(α) && β > 0
-        return log(last(promote(x, x >= 0)))
-    end
-
-    _x = clamp(x, 0, 1)
-    p, q = beta_inc(α, β, _x)
-    if p < floatmin(p)
-        # see https://dlmf.nist.gov/8.17#E7
-        return -log(α) + xlogy(α, _x) + log(_₂F₁(promote(α, 1 - β, α + 1, _x)...)) - logbeta(α, β)
-    elseif p <= 0.7
-        return log(p)
-    else
-        return log1p(-q)
-    end
-end
-betalogcdf(α::Real, β::Real, x::Real) = betalogcdf(promote(α, β, x)...)
-
-function betalogccdf(α::Real, β::Real, x::Real)
-    # Handle a degenerate case
-    if iszero(α) && β > 0
-        return log(last(promote(α, β, x, x < 0)))
-    end
-
-    p, q = beta_inc(α, β, clamp(x, 0, 1))
-    if q < 0.7
-        return log(q)
-    else
-        return log1p(-p)
-    end
-end
-
-betainvcdf(α::Real, β::Real, p::Real) = first(beta_inc_inv(α, β, p))
-
-betainvccdf(α::Real, β::Real, p::Real) = last(beta_inc_inv(β, α, p))

src/distrs/binom.jl

@@ -1,18 +1,18 @@
 # functions related to binomial distribution
 
 # R implementations
+# For pdf and logpdf we use the Julia implementation
 using .RFunctions:
-    # binompdf,
-    # binomlogpdf,
-    # binomcdf,
-    # binomccdf,
-    # binomlogcdf,
-    # binomlogccdf,
+    binomcdf,
+    binomccdf,
+    binomlogcdf,
+    binomlogccdf,
     binominvcdf,
     binominvccdf,
     binominvlogcdf,
     binominvlogccdf
 
+
 # Julia implementations
 binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k))
 
@@ -22,23 +22,3 @@ function binomlogpdf(n::T, p::T, k::T) where {T<:Real}
     val = min(0, betalogpdf(m + 1, n - m + 1, p) - log(n + 1))
     return 0 <= k <= n && isinteger(k) ? val : oftype(val, -Inf)
 end
-
-for l in ("", "log"), compl in (false, true)
-    fbinom = Symbol(string("binom", l, ifelse(compl, "c", "" ), "cdf"))
-    fbeta  = Symbol(string("beta" , l, ifelse(compl,  "", "c"), "cdf"))
-    @eval function ($fbinom)(n::Real, p::Real, k::Real)
-        if isnan(k)
-            return last(promote(n, p, k))
-        end
-        res = ($fbeta)(max(0, floor(k) + 1), max(0, n - floor(k)), p)
-
-        # When p == 1 the betaccdf doesn't return the correct result
-        # so these cases have to be special cased
-        if isone(p)
-            newres = oftype(res, $compl ? k < n : k >= n)
-            return $(l === "" ? :newres : :(log(newres)))
-        else
-            return res
-        end
-    end
-end

src/distrs/chisq.jl

@@ -1,11 +1,22 @@
 # functions related to chi-square distribution
 
-# Just use the Gamma definitions
-for f in ("pdf", "logpdf", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccdf", "invlogcdf", "invlogccdf")
-    _chisqf = Symbol("chisq"*f)
-    _gammaf = Symbol("gamma"*f)
-    @eval begin
-        $(_chisqf)(k::Real, x::Real) = $(_chisqf)(promote(k, x)...)
-        $(_chisqf)(k::T, x::T) where {T<:Real} = $(_gammaf)(k/2, 2, x)
-    end
-end
+# R implementations
+# For pdf and logpdf we use the Julia implementation
+using .RFunctions:
+    chisqcdf,
+    chisqccdf,
+    chisqlogcdf,
+    chisqlogccdf,
+    chisqinvcdf,
+    chisqinvccdf,
+    chisqinvlogcdf,
+    chisqinvlogccdf
+
+# Julia implementations
+# promotion ensures that we do forward e.g. `chisqpdf(::Int, ::Float32)` to
+# `gammapdf(::Float32, ::Int, ::Float32)` but not `gammapdf(::Float64, ::Int, ::Float32)`
+chisqpdf(k::Real, x::Real) = chisqpdf(promote(k, x)...)
+chisqpdf(k::T, x::T) where {T<:Real} = gammapdf(k / 2, 2, x)
+
+chisqlogpdf(k::Real, x::Real) = chisqlogpdf(promote(k, x)...)
+chisqlogpdf(k::T, x::T) where {T<:Real} = gammalogpdf(k / 2, 2, x)

src/distrs/fdist.jl

@@ -1,5 +1,17 @@
 # functions related to F distribution
 
+# R implementations
+# For pdf and logpdf we use the Julia implementation
+using .RFunctions:
+    fdistcdf,
+    fdistccdf,
+    fdistlogcdf,
+    fdistlogccdf,
+    fdistinvcdf,
+    fdistinvccdf,
+    fdistinvlogcdf,
+    fdistinvlogccdf
+
 # Julia implementations
 fdistpdf(ν1::Real, ν2::Real, x::Real) = exp(fdistlogpdf(ν1, ν2, x))
 
@@ -11,19 +23,3 @@ function fdistlogpdf(ν1::T, ν2::T, x::T) where {T<:Real}
     val = (xlogy(ν1, ν1ν2) + xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, 1 + ν1ν2 * y)) / 2 - logbeta(ν1 / 2, ν2 / 2)
     return x < 0 ? oftype(val, -Inf) : val
 end
-
-for f in ("cdf", "ccdf", "logcdf", "logccdf")
-    ff = Symbol("fdist"*f)
-    bf = Symbol("beta"*f)
-    @eval $ff(ν1::T, ν2::T, x::T) where {T<:Real} = $bf(ν1/2, ν2/2, inv(1 + ν2/(ν1*max(0, x))))
-    @eval $ff(ν1::Real, ν2::Real, x::Real) = $ff(promote(ν1, ν2, x)...)
-end
-for f in ("invcdf", "invccdf", "invlogcdf", "invlogccdf")
-    ff = Symbol("fdist"*f)
-    bf = Symbol("beta"*f)
-    @eval function $ff(ν1::T, ν2::T, y::T) where {T<:Real}
-        x = $bf(ν1/2, ν2/2, y)
-        return x/(1 - x)*ν2/ν1
-    end
-    @eval $ff(ν1::Real, ν2::Real, y::Real) = $ff(promote(ν1, ν2, y)...)
-end

src/distrs/gamma.jl

@@ -1,17 +1,14 @@
 # functions related to gamma distribution
 
-using HypergeometricFunctions: drummond1F1
-
 # R implementations
+# For pdf and logpdf we use the Julia implementation
 using .RFunctions:
-    # gammapdf,
-    # gammalogpdf,
-    # gammacdf,
-    # gammaccdf,
-    # gammalogcdf,
-    # gammalogccdf,
-    # gammainvcdf,
-    # gammainvccdf,
+    gammacdf,
+    gammaccdf,
+    gammalogcdf,
+    gammalogccdf,
+    gammainvcdf,
+    gammainvccdf,
     gammainvlogcdf,
     gammainvlogccdf
 
@@ -25,80 +22,3 @@ function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
     val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
     return x < 0 ? oftype(val, -Inf) : val
 end
-
-function gammacdf(k::T, θ::T, x::T) where {T<:Real}
-    # Handle the degenerate case
-    if iszero(k)
-        return float(last(promote(x, x >= 0)))
-    end
-    return first(gamma_inc(k, max(0, x)/θ))
-end
-gammacdf(k::Real, θ::Real, x::Real) = gammacdf(map(float, promote(k, θ, x))...)
-
-function gammaccdf(k::T, θ::T, x::T) where {T<:Real}
-    # Handle the degenerate case
-    if iszero(k)
-        return float(last(promote(x, x < 0)))
-    end
-    return last(gamma_inc(k, max(0, x)/θ))
-end
-gammaccdf(k::Real, θ::Real, x::Real) = gammaccdf(map(float, promote(k, θ, x))...)
-
-gammalogcdf(k::Real, θ::Real, x::Real) = _gammalogcdf(map(float, promote(k, θ, x))...)
-
-# Implemented via the non-log version. For tiny values, we recompute the result with
-# loggamma. In that situation, there is little risk of significant cancellation.
-function _gammalogcdf(k::Float64, θ::Float64, x::Float64)
-    # Handle the degenerate case
-    if iszero(k)
-        return log(x >= 0)
-    end
-
-    xdθ = max(0, x)/θ
-    l, u = gamma_inc(k, xdθ)
-    if l < floatmin(Float64) && isfinite(k) && isfinite(xdθ)
-        return -log(k) + k*log(xdθ) - xdθ + log(drummond1F1(1.0, 1 + k, xdθ)) - loggamma(k)
-    elseif l < 0.7
-        return log(l)
-    else
-        return log1p(-u)
-    end
-end
-function _gammalogcdf(k::T, θ::T, x::T) where {T<:Union{Float16,Float32}}
-    return T(_gammalogcdf(Float64(k), Float64(θ), Float64(x)))
-end
-
-gammalogccdf(k::Real, θ::Real, x::Real) = _gammalogccdf(map(float, promote(k, θ, x))...)
-
-# Implemented via the non-log version. For tiny values, we recompute the result with
-# loggamma. In that situation, there is little risk of significant cancellation.
-function _gammalogccdf(k::Float64, θ::Float64, x::Float64)
-    # Handle the degenerate case
-    if iszero(k)
-        return log(x < 0)
-    end
-
-    xdθ = max(0, x)/θ
-    l, u = gamma_inc(k, xdθ)
-    if u < floatmin(Float64)
-        return loggamma(k, xdθ) - loggamma(k)
-    elseif u < 0.7
-        return log(u)
-    else
-        return log1p(-l)
-    end
-end
-
-function _gammalogccdf(k::T, θ::T, x::T) where {T<:Union{Float16,Float32}}
-    return T(_gammalogccdf(Float64(k), Float64(θ), Float64(x)))
-end
-
-function gammainvcdf(k::Real, θ::Real, p::Real)
-    _k, _θ, _p = promote(k, θ, p)
-    return _θ*gamma_inc_inv(_k, _p, 1 - _p)
-end
-
-function gammainvccdf(k::Real, θ::Real, p::Real)
-    _k, _θ, _p = promote(k, θ, p)
-    return _θ*gamma_inc_inv(_k, 1 - _p, _p)
-end

src/distrs/pois.jl

@@ -1,13 +1,12 @@
 # functions related to Poisson distribution
 
 # R implementations
+# For pdf and logpdf we use the Julia implementation
 using .RFunctions:
-    # poispdf,
-    # poislogpdf,
-    # poiscdf,
-    # poisccdf,
-    # poislogcdf,
-    # poislogccdf,
+    poiscdf,
+    poisccdf,
+    poislogcdf,
+    poislogccdf,
     poisinvcdf,
     poisinvccdf,
     poisinvlogcdf,
@@ -21,12 +20,3 @@ function poislogpdf(λ::T, x::T) where {T <: Real}
     val = xlogy(x, λ) - λ - loggamma(x + 1)
     return x >= 0 && isinteger(x) ? val : oftype(val, -Inf)
 end
-
-# Just use the Gamma definitions
-poiscdf(λ::Real, x::Real) = gammaccdf(max(0, floor(x + 1)), 1, λ)
-
-poisccdf(λ::Real, x::Real) = gammacdf(max(0, floor(x + 1)), 1, λ)
-
-poislogcdf(λ::Real, x::Real) = gammalogccdf(max(0, floor(x + 1)), 1, λ)
-
-poislogccdf(λ::Real, x::Real) = gammalogcdf(max(0, floor(x + 1)), 1, λ)

src/distrs/tdist.jl

@@ -1,9 +1,8 @@
 # functions related to student's T distribution
 
 # R implementations
+# For pdf and logpdf we use the Julia implementation
 using .RFunctions:
-    # tdistpdf,
-    # tdistlogpdf,
     tdistcdf,
     tdistccdf,
     tdistlogcdf,

test/misc.jl

@@ -42,7 +42,12 @@ end
             #  B₁(a, b) = B(a, b)
             a = rand(eltya)
             b = rand(eltyb)
-            @test logmvbeta(1, a, b) ≈ logbeta(a, b)
+            # Older SpecialFunctions + Julia versions require slightly larger tolerance
+            if VERSION < v"1.3" && promote_type(eltya, eltyb) === Float64
+                @test logmvbeta(1, a, b) ≈ logbeta(a, b) rtol = 10 * sqrt(eps(Float64))
+            else
+                @test logmvbeta(1, a, b) ≈ logbeta(a, b)
+            end
         end
     end