Merge pull request #383 from JuliaMath/release-1.8-backports

Release 1.8 backports

Author: andreasnoack

src/beta_inc.jl

@@ -930,7 +930,8 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
     x = p
     r = sqrt(-2*log(p))
     p_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00)
-    fpu = 1e-30
+    fpu = floatmin(Float64)
+    sae = log10(fpu)
     lb = logbeta(a, b)
 
     if a > 1.0 && b > 1.0
@@ -971,8 +972,30 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
         x = .9999
     end
 
-    iex = max(-5.0/a^2 - 1.0/p^0.2 - 13.0, -30.0)
-    acu = 10.0^iex
+    # This first argument was proposed in
+    #
+    # K. J. Berry, P. W. Mielke, Jr and G. W. Cran (1990).
+    # Algorithm AS R83: A Remark on Algorithm AS 109: Inverse of the
+    #   Incomplete Beta Function Ratio.
+    # Journal of the Royal Statistical Society.
+    # Series C (Applied Statistics), 39(2), 309–310. doi:10.2307/2347779
+    #
+    # but the last term has been changed from 13 to 34 since the
+    # the original article
+    #
+    # Majumder, K. L., & Bhattacharjee, G. P. (1973).
+    # Algorithm as 64: Inverse of the incomplete beta function ratio.
+    # Journal of the Royal Statistical Society.
+    # Series C (Applied Statistics), 22(3), 411-414.
+    #
+    # argues that the iex value should be set to -2r - 2 where r is the
+    # required number of accurate digits.
+    #
+    # The idea with the -5.0/a^2 - 1.0/p^0.2 - 34.0 correction is to
+    # use -2r - 2 (for 16 digits) for large values of a and p but use
+    # a much smaller tolerance for small values of a and p.
+    iex = max(-5.0/a^2 - 1.0/p^0.2 - 34.0, sae)
+    acu = exp10(iex)
 
     #iterate
     while true

src/gamma_inc.jl

@@ -803,6 +803,10 @@ function _gamma_inc(a::Float64, x::Float64, ind::Integer)
         else
             return (1.0, 0.0)
         end
+    elseif isnan(a) || isnan(x)
+        return (a*x, a*x)
+    elseif isinf(x)
+        return (1.0, 0.0)
     end
 
     if a >= 1.0
@@ -1014,8 +1018,17 @@ function _gamma_inc_inv(a::Float64, p::Float64, q::Float64)
     return x
 end
 
-_gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
-_gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
+function _gamma_inc_inv(a::T, p::T, q::T) where {T <: Union{Float16, Float32}}
+    if p + q != one(T)
+        throw(ArgumentError("p + q must equal one but was $(p + q)"))
+    end
+    p64, q64 = if p < q
+        (Float64(p), 1 - Float64(p))
+    else
+        (1 - Float64(q), Float64(q))
+    end
+    return T(_gamma_inc_inv(Float64(a), p64, q64))
+end
 
 # like promote(x,y), but don't complexify real values
 promotereal(x::Real, y::Real) = promote(x, y)

test/beta_inc.jl

@@ -234,52 +234,53 @@ end
 
 @testset "inverse of incomplete beta" begin
     f(a,b,p) = beta_inc_inv(a,b,p)[1]
-    @test f(.5, .5, 0.6376856085851985E-01) ≈ 0.01
-    @test f(.5, .5, 0.20483276469913355) ≈ 0.1
-    @test f(.5, .5, 1.0000) ≈ 1.0000
-    @test f(1.0, .5, 0.0) ≈ 0.00
-    @test f(1.0, .5, 0.5012562893380045E-02) ≈ 0.01
-    @test f(1.0, .5, 0.5131670194948620E-01) ≈ 0.1
-    @test f(1.0, .5, 0.2928932188134525) ≈ 0.5
-    @test f(1.0, 1.0, .5) ≈ 0.5
-    @test f(2.0, 2.0, .028) ≈ 0.1
-    @test f(2.0, 2.0, 0.104) ≈ 0.2
-    @test f(2.0, 2.0, .216) ≈ 0.3
-    @test f(2.0, 2.0, .352) ≈ 0.4
-    @test f(2.0, 2.0, .5) ≈ 0.5
-    @test f(2.0, 2.0, 0.648) ≈ 0.6
-    @test f(2.0, 2.0, 0.784) ≈ 0.7
-    @test f(2.0, 2.0, 0.896) ≈ 0.8
-    @test f(2.0, 2.0, .972) ≈ 0.9
-    @test f(5.5, 5.0, 0.4361908850559777) ≈ .5
-    @test f(10.0, .5, 0.1516409096347099) ≈ 0.9
-    @test f(10.0, 5.0, 0.8978271484375000E-01) ≈ 0.5
-    @test f(10.0, 5.0, 1.00) ≈ 1.00
-    @test f(10.0, 10.0, .5) ≈ .5
-    @test f(20.0, 5.0, 0.4598773297575791) ≈ 0.8
-    @test f(20.0, 10.0, 0.2146816102371739) ≈ 0.6
-    @test f(20.0, 10.0, 0.9507364826957875) ≈ 0.8
-    @test f(20.0, 20.0, .5) ≈ .5
-    @test f(20.0, 20.0, 0.8979413687105918) ≈ 0.6
-    @test f(30.0, 10.0, 0.2241297491808366) ≈ 0.7
-    @test f(30.0, 10.0, 0.7586405487192086) ≈ 0.8
-    @test f(40.0, 20.0, 0.7001783247477069) ≈ 0.7
-    @test f(1.0, 0.5, 0.5131670194948620E-01) ≈ 0.1
-    @test f(1.0, 0.5, 0.1055728090000841) ≈ 0.2
-    @test f(1.0, 0.5, 0.1633399734659245) ≈ 0.3
-    @test f(1.0, 0.5, 0.2254033307585166) ≈ 0.4
-    @test f(1.0, 2.0, .36) ≈ 0.2
-    @test f(1.0, 3.0, .488) ≈ 0.2
-    @test f(1.0, 4.0, .5904) ≈ 0.2
-    @test f(1.0, 5.0, .67232) ≈ 0.2
-    @test f(2.0, 2.0, .216) ≈ 0.3
-    @test f(3.0, 2.0, 0.837e-1) ≈ 0.3
-    @test f(4.0, 2.0, 0.3078e-1) ≈ 0.3
-    @test f(5.0, 2.0, 0.10935e-1) ≈ 0.3
-    @test f(1.30625000, 11.75620000, 0.9188846846205182) ≈ 0.225609
-    @test f(1.30625000, 11.75620000, 0.21053116418502513) ≈ 0.033557
-    @test f(1.30625000, 11.75620000, 0.18241165418408148) ≈ 0.029522
+    @test f(.5, .5, 0.6376856085851985E-01)    ≈ 0.01   rtol=8eps()
+    @test f(.5, .5, 0.20483276469913355)       ≈ 0.1    rtol=8eps()
+    @test f(.5, .5, 1.0000)                    ≈ 1.0000 rtol=8eps()
+    @test f(1.0, .5, 0.0)                     == 0.0
+    @test f(1.0, .5, 0.5012562893380045E-02)   ≈ 0.01   rtol=8eps()
+    @test f(1.0, .5, 0.5131670194948620E-01)   ≈ 0.1    rtol=8eps()
+    @test f(1.0, .5, 0.2928932188134525)       ≈ 0.5    rtol=8eps()
+    @test f(1.0, 1.0, .5)                      ≈ 0.5    rtol=8eps()
+    @test f(2.0, 2.0, .028)                    ≈ 0.1    rtol=8eps()
+    @test f(2.0, 2.0, 0.104)                   ≈ 0.2    rtol=8eps()
+    @test f(2.0, 2.0, .216)                    ≈ 0.3    rtol=8eps()
+    @test f(2.0, 2.0, .352)                    ≈ 0.4    rtol=8eps()
+    @test f(2.0, 2.0, .5)                      ≈ 0.5    rtol=8eps()
+    @test f(2.0, 2.0, 0.648)                   ≈ 0.6    rtol=8eps()
+    @test f(2.0, 2.0, 0.784)                   ≈ 0.7    rtol=8eps()
+    @test f(2.0, 2.0, 0.896)                   ≈ 0.8    rtol=8eps()
+    @test f(2.0, 2.0, .972)                    ≈ 0.9    rtol=8eps()
+    @test f(5.5, 5.0, 0.4361908850559777)      ≈ 0.5    rtol=8eps()
+    @test f(10.0, .5, 0.1516409096347099)      ≈ 0.9    rtol=8eps()
+    @test f(10.0, 5.0, 0.8978271484375000E-01) ≈ 0.5    rtol=8eps()
+    @test f(10.0, 5.0, 1.00)                   ≈ 1.0    rtol=8eps()
+    @test f(10.0, 10.0, .5)                    ≈ 0.5    rtol=8eps()
+    @test f(20.0, 5.0, 0.4598773297575791)     ≈ 0.8    rtol=8eps()
+    @test f(20.0, 10.0, 0.2146816102371739)    ≈ 0.6    rtol=8eps()
+    @test f(20.0, 10.0, 0.9507364826957875)    ≈ 0.8    rtol=8eps()
+    @test f(20.0, 20.0, .5)                    ≈ 0.5    rtol=8eps()
+    @test f(20.0, 20.0, 0.8979413687105918)    ≈ 0.6    rtol=8eps()
+    @test f(30.0, 10.0, 0.2241297491808366)    ≈ 0.7    rtol=8eps()
+    @test f(30.0, 10.0, 0.7586405487192086)    ≈ 0.8    rtol=8eps()
+    @test f(40.0, 20.0, 0.7001783247477069)    ≈ 0.7    rtol=8eps()
+    @test f(1.0, 0.5, 0.5131670194948620E-01)  ≈ 0.1    rtol=8eps()
+    @test f(1.0, 0.5, 0.1055728090000841)      ≈ 0.2    rtol=8eps()
+    @test f(1.0, 0.5, 0.1633399734659245)      ≈ 0.3    rtol=8eps()
+    @test f(1.0, 0.5, 0.2254033307585166)      ≈ 0.4    rtol=8eps()
+    @test f(1.0, 2.0, .36)                     ≈ 0.2    rtol=8eps()
+    @test f(1.0, 3.0, .488)                    ≈ 0.2    rtol=8eps()
+    @test f(1.0, 4.0, .5904)                   ≈ 0.2    rtol=8eps()
+    @test f(1.0, 5.0, .67232)                  ≈ 0.2    rtol=8eps()
+    @test f(2.0, 2.0, .216)                    ≈ 0.3    rtol=8eps()
+    @test f(3.0, 2.0, 0.837e-1)                ≈ 0.3    rtol=8eps()
+    @test f(4.0, 2.0, 0.3078e-1)               ≈ 0.3    rtol=8eps()
+    @test f(5.0, 2.0, 0.10935e-1)              ≈ 0.3    rtol=8eps()
+    @test f(1.30625000, 11.75620000, 0.9188846846205182)  ≈ 0.225609 rtol=8eps()
+    @test f(1.30625000, 11.75620000, 0.21053116418502513) ≈ 0.033557 rtol=8eps()
+    @test f(1.30625000, 11.75620000, 0.18241165418408148) ≈ 0.029522 rtol=8eps()
     @test f(1000.0, 2.0, 9.0797754e-317) ≈ 0.48 # This one is a bit slow (but also hard)
+    @test f(1.5, 5.0, beta_inc(1.5, 5.0, 0.2142857142857142)[1])[1] ≈ 0.2142857142857142 rtol=8eps()
 
     for T in (Float16, Float32, Float64)
         p = rand(T)

test/gamma_inc.jl

@@ -44,10 +44,17 @@
     @test gamma_inc(1.1, 1e3)[2] == 0.0
     @test gamma_inc(24.0, 1e-323)[1] == 0.0
     @test gamma_inc(6311.0, 6311.0*0.59999)[1] < 1e-300
+    @test gamma_inc(0.5, Inf) == (1.0, 0.0)
+    @test all(isnan, gamma_inc(NaN, 1.0))
+    @test all(isnan, gamma_inc(1.0, NaN))
+    @test all(isnan, gamma_inc(NaN, Inf))
     @test_throws DomainError gamma_inc(-1, 2, 2)
     @test_throws DomainError gamma_inc(0, 0, 1)
     @test_throws DomainError gamma_inc(7.098843361278083e33, 7.09884336127808e33)
     @test_throws DomainError gamma_inc(6.693195169205051e28, 6.693195169205049e28)
+    @test_throws DomainError gamma_inc(NaN, -1.0)
+    @test_throws DomainError gamma_inc(-1.0, NaN)
+    @test_throws DomainError gamma_inc(1.0, -Inf)
 end
 
 @testset "inverse of incomplete gamma ratios" begin
@@ -164,7 +171,15 @@ end
     for x = .5:5.0:100.0
         @test SpecialFunctions.stirling_error(x) ≈ log(gamma(x)) - (x-.5)*log(x)+x- log(2*pi)/2.0
     end
+
     @test_throws ArgumentError("p + q must equal one but is 0.5") gamma_inc_inv(0.4, 0.2, 0.3)
+
+    @testset "Low precision with Float64(p) + Float64(q) != 1" for T in (Float16, Float32)
+        @test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.1), T(0.9)))[1]::T ≈ T(0.1)
+        @test gamma_inc(T(1.0), gamma_inc_inv(T(1.0), T(0.9), T(0.1)))[2]::T ≈ T(0.1)
+        @test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.1), T(0.92))
+        @test_throws ArgumentError("p + q must equal one but was 1.02") gamma_inc_inv(T(1.0), T(0.92), T(0.1))
+    end
 end
 
 double(x::Real) = Float64(x)