From 4d64e3ce44835482bdf1bfcc6629cb6d41657a14 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 30 Jan 2022 10:05:05 +0100 Subject: [PATCH 1/4] Handle the case when Float64(p) + Float64(q) != 1 in gamma_inc_inv for lower precision p and q. --- src/gamma_inc.jl | 18 ++++++++++++++++-- test/gamma_inc.jl | 8 ++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 36b72ac9..30dc5c0b 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -1014,8 +1014,22 @@ 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}} + p64, q64 = if p < q + _q = 1 - Float64(p) + if q != T(_q) + throw(ArgumentError("p + q must equal one but was $(p + q)")) + end + Float64(p), _q + else + _p = 1 - Float64(q) + if p != T(_p) + throw(ArgumentError("p + q must equal one but was $(p + q)")) + end + _p, 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) diff --git a/test/gamma_inc.jl b/test/gamma_inc.jl index b4e7c6c1..1ae42e9c 100644 --- a/test/gamma_inc.jl +++ b/test/gamma_inc.jl @@ -164,7 +164,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) From 6b84662da26d53d4abd1f7b821f2a67e48c0cdb7 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 31 Jan 2022 13:35:06 +0100 Subject: [PATCH 2/4] Use Speve's version --- src/gamma_inc.jl | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 30dc5c0b..eb278ed1 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -1015,18 +1015,13 @@ function _gamma_inc_inv(a::Float64, p::Float64, q::Float64) end 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 - _q = 1 - Float64(p) - if q != T(_q) - throw(ArgumentError("p + q must equal one but was $(p + q)")) - end - Float64(p), _q + (Float64(p), 1 - Float64(p)) else - _p = 1 - Float64(q) - if p != T(_p) - throw(ArgumentError("p + q must equal one but was $(p + q)")) - end - _p, Float64(q) + (1 - Float64(q), Float64(q)) end return T(_gamma_inc_inv(Float64(a), p64, q64)) end From 9866ac95c29468d4a38a264eb3f77dd76ee9e081 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 31 Jan 2022 23:30:37 +0100 Subject: [PATCH 3/4] Adjust accuracy tolerance in beta_inc_inv to accommodate Float64 --- src/beta_inc.jl | 29 +++++++++++++-- test/beta_inc.jl | 91 ++++++++++++++++++++++++------------------------ 2 files changed, 72 insertions(+), 48 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 67c0d350..8ce4adc9 100644 --- a/src/beta_inc.jl +++ b/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 diff --git a/test/beta_inc.jl b/test/beta_inc.jl index 8f9fcbb6..2b0d04d5 100644 --- a/test/beta_inc.jl +++ b/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) From 677c2619f43e1ad0f224d19a382a13aec4f04bfe Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 6 Feb 2022 22:09:49 +0100 Subject: [PATCH 4/4] Handle NaNs and Infs in gamma_inc --- src/gamma_inc.jl | 4 ++++ test/gamma_inc.jl | 7 +++++++ 2 files changed, 11 insertions(+) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index eb278ed1..3f5cdc20 100644 --- a/src/gamma_inc.jl +++ b/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 diff --git a/test/gamma_inc.jl b/test/gamma_inc.jl index 1ae42e9c..07d3089c 100644 --- a/test/gamma_inc.jl +++ b/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