Use _₁F₁ instead of drummond1F1 in gammalogcdf (#151)

After the native Julia implementation for the Gamma distribution
was implemented, HypergeometricFunctions added the _₁F₁ function
which will decide on the algorithm based on the arguments. This
ensures good behavior in a broader part of the parameter space.

src/distrs/gamma.jl

@@ -1,6 +1,6 @@
 # functions related to gamma distribution
 
-using HypergeometricFunctions: drummond1F1
+using HypergeometricFunctions: _₁F₁
 
 # R implementations
 using .RFunctions:
@@ -57,7 +57,7 @@ function _gammalogcdf(k::Float64, θ::Float64, x::Float64)
     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)
+        return -log(k) + k*log(xdθ) - xdθ + log(_₁F₁(1.0, 1 + k, xdθ)) - loggamma(k)
     elseif l < 0.7
         return log(l)
     else

test/misc.jl

@@ -85,8 +85,13 @@ end
 
 # https://github.com/JuliaStats/StatsFuns.jl/issues/143
 # https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/47
-@testset "logbetacdf: numerical issue" begin
+@testset "betalogcdf: 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
+
+# https://github.com/JuliaStats/StatsFuns.jl/issues/150
+@testset "gammalogcdf: numerical issue" begin
+    @test gammalogcdf(42648.50647826826, 2.2498007956420723e-5, 0.6991377135675367) ≈ -1933.2698959040617410
+end