Skip to content

Commit 9dd5bf8

Browse files
andreasnoackandreasnoack
authored
Adjust accuracy tolerance in beta_inc_inv to accommodate Float64 (#376)
1 parent b16f79f commit 9dd5bf8

File tree

2 files changed

+72
-48
lines changed

2 files changed

+72
-48
lines changed

‎src/beta_inc.jl

Lines changed: 26 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -930,7 +930,8 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
930930
x = p
931931
r = sqrt(-2*log(p))
932932
p_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00)
933-
fpu = 1e-30
933+
fpu = floatmin(Float64)
934+
sae = log10(fpu)
934935
lb = logbeta(a, b)
935936

936937
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)
971972
x = .9999
972973
end
973974

974-
iex = max(-5.0/a^2 - 1.0/p^0.2 - 13.0, -30.0)
975-
acu = 10.0^iex
975+
# This first argument was proposed in
976+
#
977+
# K. J. Berry, P. W. Mielke, Jr and G. W. Cran (1990).
978+
# Algorithm AS R83: A Remark on Algorithm AS 109: Inverse of the
979+
# Incomplete Beta Function Ratio.
980+
# Journal of the Royal Statistical Society.
981+
# Series C (Applied Statistics), 39(2), 309–310. doi:10.2307/2347779
982+
#
983+
# but the last term has been changed from 13 to 34 since the
984+
# the original article
985+
#
986+
# Majumder, K. L., & Bhattacharjee, G. P. (1973).
987+
# Algorithm as 64: Inverse of the incomplete beta function ratio.
988+
# Journal of the Royal Statistical Society.
989+
# Series C (Applied Statistics), 22(3), 411-414.
990+
#
991+
# argues that the iex value should be set to -2r - 2 where r is the
992+
# required number of accurate digits.
993+
#
994+
# The idea with the -5.0/a^2 - 1.0/p^0.2 - 34.0 correction is to
995+
# use -2r - 2 (for 16 digits) for large values of a and p but use
996+
# a much smaller tolerance for small values of a and p.
997+
iex = max(-5.0/a^2 - 1.0/p^0.2 - 34.0, sae)
998+
acu = exp10(iex)
976999

9771000
#iterate
9781001
while true

‎test/beta_inc.jl

Lines changed: 46 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -234,52 +234,53 @@ end
234234

235235
@testset "inverse of incomplete beta" begin
236236
f(a,b,p) = beta_inc_inv(a,b,p)[1]
237-
@test f(.5, .5, 0.6376856085851985E-01) ≈ 0.01
238-
@test f(.5, .5, 0.20483276469913355) ≈ 0.1
239-
@test f(.5, .5, 1.0000) ≈ 1.0000
240-
@test f(1.0, .5, 0.0) ≈ 0.00
241-
@test f(1.0, .5, 0.5012562893380045E-02) ≈ 0.01
242-
@test f(1.0, .5, 0.5131670194948620E-01) ≈ 0.1
243-
@test f(1.0, .5, 0.2928932188134525) ≈ 0.5
244-
@test f(1.0, 1.0, .5) ≈ 0.5
245-
@test f(2.0, 2.0, .028) ≈ 0.1
246-
@test f(2.0, 2.0, 0.104) ≈ 0.2
247-
@test f(2.0, 2.0, .216) ≈ 0.3
248-
@test f(2.0, 2.0, .352) ≈ 0.4
249-
@test f(2.0, 2.0, .5) ≈ 0.5
250-
@test f(2.0, 2.0, 0.648) ≈ 0.6
251-
@test f(2.0, 2.0, 0.784) ≈ 0.7
252-
@test f(2.0, 2.0, 0.896) ≈ 0.8
253-
@test f(2.0, 2.0, .972) ≈ 0.9
254-
@test f(5.5, 5.0, 0.4361908850559777) ≈ .5
255-
@test f(10.0, .5, 0.1516409096347099) ≈ 0.9
256-
@test f(10.0, 5.0, 0.8978271484375000E-01) ≈ 0.5
257-
@test f(10.0, 5.0, 1.00) ≈ 1.00
258-
@test f(10.0, 10.0, .5) ≈ .5
259-
@test f(20.0, 5.0, 0.4598773297575791) ≈ 0.8
260-
@test f(20.0, 10.0, 0.2146816102371739) ≈ 0.6
261-
@test f(20.0, 10.0, 0.9507364826957875) ≈ 0.8
262-
@test f(20.0, 20.0, .5) ≈ .5
263-
@test f(20.0, 20.0, 0.8979413687105918) ≈ 0.6
264-
@test f(30.0, 10.0, 0.2241297491808366) ≈ 0.7
265-
@test f(30.0, 10.0, 0.7586405487192086) ≈ 0.8
266-
@test f(40.0, 20.0, 0.7001783247477069) ≈ 0.7
267-
@test f(1.0, 0.5, 0.5131670194948620E-01) ≈ 0.1
268-
@test f(1.0, 0.5, 0.1055728090000841) ≈ 0.2
269-
@test f(1.0, 0.5, 0.1633399734659245) ≈ 0.3
270-
@test f(1.0, 0.5, 0.2254033307585166) ≈ 0.4
271-
@test f(1.0, 2.0, .36) ≈ 0.2
272-
@test f(1.0, 3.0, .488) ≈ 0.2
273-
@test f(1.0, 4.0, .5904) ≈ 0.2
274-
@test f(1.0, 5.0, .67232) ≈ 0.2
275-
@test f(2.0, 2.0, .216) ≈ 0.3
276-
@test f(3.0, 2.0, 0.837e-1) ≈ 0.3
277-
@test f(4.0, 2.0, 0.3078e-1) ≈ 0.3
278-
@test f(5.0, 2.0, 0.10935e-1) ≈ 0.3
279-
@test f(1.30625000, 11.75620000, 0.9188846846205182) ≈ 0.225609
280-
@test f(1.30625000, 11.75620000, 0.21053116418502513) ≈ 0.033557
281-
@test f(1.30625000, 11.75620000, 0.18241165418408148) ≈ 0.029522
237+
@test f(.5, .5, 0.6376856085851985E-01) ≈ 0.01 rtol=8eps()
238+
@test f(.5, .5, 0.20483276469913355) ≈ 0.1 rtol=8eps()
239+
@test f(.5, .5, 1.0000) ≈ 1.0000 rtol=8eps()
240+
@test f(1.0, .5, 0.0) == 0.0
241+
@test f(1.0, .5, 0.5012562893380045E-02) ≈ 0.01 rtol=8eps()
242+
@test f(1.0, .5, 0.5131670194948620E-01) ≈ 0.1 rtol=8eps()
243+
@test f(1.0, .5, 0.2928932188134525) ≈ 0.5 rtol=8eps()
244+
@test f(1.0, 1.0, .5) ≈ 0.5 rtol=8eps()
245+
@test f(2.0, 2.0, .028) ≈ 0.1 rtol=8eps()
246+
@test f(2.0, 2.0, 0.104) ≈ 0.2 rtol=8eps()
247+
@test f(2.0, 2.0, .216) ≈ 0.3 rtol=8eps()
248+
@test f(2.0, 2.0, .352) ≈ 0.4 rtol=8eps()
249+
@test f(2.0, 2.0, .5) ≈ 0.5 rtol=8eps()
250+
@test f(2.0, 2.0, 0.648) ≈ 0.6 rtol=8eps()
251+
@test f(2.0, 2.0, 0.784) ≈ 0.7 rtol=8eps()
252+
@test f(2.0, 2.0, 0.896) ≈ 0.8 rtol=8eps()
253+
@test f(2.0, 2.0, .972) ≈ 0.9 rtol=8eps()
254+
@test f(5.5, 5.0, 0.4361908850559777) ≈ 0.5 rtol=8eps()
255+
@test f(10.0, .5, 0.1516409096347099) ≈ 0.9 rtol=8eps()
256+
@test f(10.0, 5.0, 0.8978271484375000E-01) ≈ 0.5 rtol=8eps()
257+
@test f(10.0, 5.0, 1.00) ≈ 1.0 rtol=8eps()
258+
@test f(10.0, 10.0, .5) ≈ 0.5 rtol=8eps()
259+
@test f(20.0, 5.0, 0.4598773297575791) ≈ 0.8 rtol=8eps()
260+
@test f(20.0, 10.0, 0.2146816102371739) ≈ 0.6 rtol=8eps()
261+
@test f(20.0, 10.0, 0.9507364826957875) ≈ 0.8 rtol=8eps()
262+
@test f(20.0, 20.0, .5) ≈ 0.5 rtol=8eps()
263+
@test f(20.0, 20.0, 0.8979413687105918) ≈ 0.6 rtol=8eps()
264+
@test f(30.0, 10.0, 0.2241297491808366) ≈ 0.7 rtol=8eps()
265+
@test f(30.0, 10.0, 0.7586405487192086) ≈ 0.8 rtol=8eps()
266+
@test f(40.0, 20.0, 0.7001783247477069) ≈ 0.7 rtol=8eps()
267+
@test f(1.0, 0.5, 0.5131670194948620E-01) ≈ 0.1 rtol=8eps()
268+
@test f(1.0, 0.5, 0.1055728090000841) ≈ 0.2 rtol=8eps()
269+
@test f(1.0, 0.5, 0.1633399734659245) ≈ 0.3 rtol=8eps()
270+
@test f(1.0, 0.5, 0.2254033307585166) ≈ 0.4 rtol=8eps()
271+
@test f(1.0, 2.0, .36) ≈ 0.2 rtol=8eps()
272+
@test f(1.0, 3.0, .488) ≈ 0.2 rtol=8eps()
273+
@test f(1.0, 4.0, .5904) ≈ 0.2 rtol=8eps()
274+
@test f(1.0, 5.0, .67232) ≈ 0.2 rtol=8eps()
275+
@test f(2.0, 2.0, .216) ≈ 0.3 rtol=8eps()
276+
@test f(3.0, 2.0, 0.837e-1) ≈ 0.3 rtol=8eps()
277+
@test f(4.0, 2.0, 0.3078e-1) ≈ 0.3 rtol=8eps()
278+
@test f(5.0, 2.0, 0.10935e-1) ≈ 0.3 rtol=8eps()
279+
@test f(1.30625000, 11.75620000, 0.9188846846205182) ≈ 0.225609 rtol=8eps()
280+
@test f(1.30625000, 11.75620000, 0.21053116418502513) ≈ 0.033557 rtol=8eps()
281+
@test f(1.30625000, 11.75620000, 0.18241165418408148) ≈ 0.029522 rtol=8eps()
282282
@test f(1000.0, 2.0, 9.0797754e-317) ≈ 0.48 # This one is a bit slow (but also hard)
283+
@test f(1.5, 5.0, beta_inc(1.5, 5.0, 0.2142857142857142)[1])[1] ≈ 0.2142857142857142 rtol=8eps()
283284

284285
for T in (Float16, Float32, Float64)
285286
p = rand(T)

0 commit comments

Comments
 (0)