Fix betalogcdf with large parameters (#144)

* Fix `betalogcdf` with large parameters

* Handle special case of infinite parameter values in `logbetacdf`

* Revert some unrelated changes

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.0.0"
+version = "1.0.1"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -14,7 +14,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 ChainRulesCore = "1"
-HypergeometricFunctions = "0.3"
+HypergeometricFunctions = "0.3.10"
 InverseFunctions = "0.1"
 IrrationalConstants = "0.1"
 LogExpFunctions = "0.3.2"

src/distrs/beta.jl

@@ -55,8 +55,9 @@ function betalogcdf(α::T, β::T, x::T) where {T<:Real}
     _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(α, β)
+        # see https://dlmf.nist.gov/8.17#E8
+        # we use E8 instead of E7 due to https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/47
+        return -log(α) + xlogy(α, _x) + xlog1py(β, -_x) + log(_₂F₁(promote(α + β, 1, α + 1, _x)...; method=:positive)) - logbeta(α, β)
     elseif p <= 0.7
         return log(p)
     else

test/misc.jl

@@ -82,3 +82,11 @@ end
     # for 64.0f0 the expansion is used but for 64.0 the BigFloat value is rounded
     @test Float32(lstirling_asym(64.0)) ≈ @inferred lstirling_asym(64.0f0)
 end
+
+# https://github.com/JuliaStats/StatsFuns.jl/issues/143
+# https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/47
+@testset "logbetacdf: numerical issue" begin
+    # Mathematica: N[Log[CDF[BetaDistribution[6041, 2496], 1/10]], 10]
+    @test betalogcdf(6041, 2496, 0.1) ≈ -9020.029401
+    @test betainvlogcdf(6041, 2496, betalogcdf(6041, 2496, 0.1)) ≈ 0.1
+end