Skip to content

Commit

Permalink
Merge pull request #383 from JuliaMath/release-1.8-backports
Browse files Browse the repository at this point in the history 
Release 1.8 backports
  • Loading branch information
@andreasnoack
andreasnoack authored Feb 18, 2022
2 parents 0597dd0 + 677c261 commit 9767e2b
Show file tree
Hide file tree
Showing 4 changed files with 102 additions and 50 deletions.
29 changes: 26 additions & 3 deletions src/beta_inc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
17 changes: 15 additions & 2 deletions src/gamma_inc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
91 changes: 46 additions & 45 deletions test/beta_inc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
15 changes: 15 additions & 0 deletions test/gamma_inc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 9767e2b

Please sign in to comment.