From cd255e33e34c5c9fd4893f3c0614e451f827dcfc Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 16:31:21 +0200 Subject: [PATCH 01/11] Some fixes and simplications to the beta function(s). The handling of large magnitude differences was only applied to logabsbeta but not logbeta. They now share code so the special handling applies to both functions. The cutoff for special handling has also been lowered to avoid loss of precision for intermediate argument size differences. Finally, helper functions for the beta function now have assertions to ensure that they are only applied for the ranges they are written for. --- src/beta_inc.jl | 115 ++++++++++++++++++++++++++--------------------- src/gamma.jl | 82 ++++++++++++++------------------- src/gamma_inc.jl | 7 ++- test/gamma.jl | 45 ++++++++++++------- 4 files changed, 132 insertions(+), 117 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 89dc4bde..6a9183f4 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -13,6 +13,7 @@ loggammadiv(a::Number, b::Number) = _loggammadiv(float(a), float(b)) _loggammadiv(a::T, b::T) where {T<:Base.IEEEFloat} = T(_loggammadiv(Float64(a), Float64(b))) # handle Float16, Float32 function _loggammadiv(a::Float64, b::Float64) + @assert b >= 8 if a > b h = b/a c = 1.0/(1.0 + h) @@ -46,11 +47,13 @@ end stirling_corr(a0,b0) Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) -for a0,b0 >= 8 +for a0, b0 >= 8 """ function stirling_corr(a0::Float64, b0::Float64) - a = min(a0,b0) - b = max(a0,b0) + @assert a0 >= 0 + @assert b0 >= 0 + a = min(a0, b0) + b = max(a0, b0) h = a/b c = h/(1.0 + h) @@ -90,11 +93,11 @@ function esum(mu::Float64, x::Float64) end """ - beta_integrand(a,b,x,y,mu=0.0) + beta_integrand(a, b, x, y, mu=0.0) Compute ``e^{mu} * x^{a}y^{b}/B(a,b)`` """ -function beta_integrand(a::Float64,b::Float64,x::Float64,y::Float64,mu::Float64=0.0) +function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Float64=0.0) a0, b0 = minmax(a,b) if a0 >= 8.0 if a > b @@ -175,7 +178,7 @@ end """ beta_inc_cont_fraction(a,b,x,y,lambda,epps) -Compute ``I_{x}(a,b)`` using continued fraction expansion when `a,b > 1`. +Compute ``I_{x}(a,b)`` using continued fraction expansion when `a, b > 1`. It is assumed that ``\\lambda = (a+b)*y - b`` External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) @@ -186,6 +189,8 @@ See also: [`beta_inc`](@ref) `BFRAC(A,B,X,Y,LAMBDA,EPS)` from Didonato and Morris (1982) """ function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, lambda::Float64, epps::Float64) + @assert a > 1 + @assert b > 1 ans = beta_integrand(a,b,x,y) if ans == 0.0 return 0.0 @@ -241,7 +246,7 @@ end """ beta_inc_asymptotic_symmetric(a,b,lambda,epps) -Compute ``I_{x}(a,b)`` using asymptotic expansion for `a,b >= 15`. +Compute ``I_{x}(a,b)`` using asymptotic expansion for `a, b >= 15`. It is assumed that ``\\lambda = (a+b)*y - b``. External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) @@ -252,6 +257,8 @@ See also: [`beta_inc`](@ref) `BASYM(A,B,LAMBDA,EPS)` from Didonato and Morris (1982) """ function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64) + @assert a >= 15 + @assert b >= 15 a0 =zeros(22) b0 = zeros(22) c = zeros(22) @@ -338,9 +345,9 @@ function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, end """ - beta_inc_asymptotic_asymmetric(a,b,x,y,w,epps) + beta_inc_asymptotic_asymmetric(a, b, x, y, w, epps) -Evaluation of ``I_{x}(a,b)`` when `b < min(epps,epps*a)` and `x <= 0.5` using asymptotic expansion. +Evaluation of ``I_{x}(a,b)`` using asymptotic expansion. It is assumed `a >= 15` and `b <= 1`, and epps is tolerance used. External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) @@ -348,6 +355,8 @@ External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wi See also: [`beta_inc`](@ref) """ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64) + @assert a >= 15 + @assert b <= 1 c = zeros(31) d = zeros(31) bm1 = b - 1.0 @@ -424,6 +433,8 @@ See also: [`beta_inc`](@ref) `FPSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64) + @assert b < eps(Float64)*min(1.0, a) + @assert x <= 0.5 ans = 1.0 if a > 1.0e-3 * epps ans = 0.0 @@ -466,6 +477,9 @@ See also: [`beta_inc`](@ref) `APSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64) + @assert a <= eps(Float64)*min(1.0, b) + @assert b*x <= 1 + @assert x <= 0.5 g = Base.MathConstants.eulergamma bx = b*x t = x - bx @@ -492,7 +506,7 @@ end #B .LE. 1 OR B*X .LE. 0.7 """ - beta_inc_power_series(a,b,x,epps) + beta_inc_power_series(a, b, x, epps) Computes ``I_x(a,b)`` using power series : ```math @@ -506,6 +520,7 @@ See also: [`beta_inc`](@ref) `BPSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64) + @assert b <= 1 || b*x <= 0.7 ans = 0.0 if x == 0.0 return 0.0 @@ -592,7 +607,7 @@ function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64 end """ - beta_inc_diff(a,b,x,y,n,epps) + beta_inc_diff(a, b, x, y, n, epps) Compute ``I_{x}(a,b) - I_{x}(a+n,b)`` where `n` is positive integer and epps is tolerance. A generalised version of [DLMF](https://dlmf.nist.gov/8.17.20). @@ -694,7 +709,7 @@ end #Wikipedia : https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function """ - beta_inc(a,b,x) + beta_inc(a, b, x[, y]) Returns a tuple ``(I_{x}(a,b),1.0-I_{x}(a,b))`` where the Regularized Incomplete Beta Function is given by: ```math @@ -704,40 +719,40 @@ where ``B(a,b) = \\Gamma(a)\\Gamma(b)/\\Gamma(a+b)``. External links: [DLMF](https://dlmf.nist.gov/8.17.1), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) -See also: [`beta_inc_inv(a,b,p,q)`](@ref SpecialFunctions.beta_inc_inv) +See also: [`beta_inc_inv(a, b, p, q)`](@ref SpecialFunctions.beta_inc_inv) """ function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) p = 0.0 q = 0.0 # lambda = a - (a+b)*x if a < 0.0 || b < 0.0 - return error("a or b is negative") + return throw(DomainError((a, b), "a or b is negative")) elseif a == 0.0 && b == 0.0 - return error("a and b are 0.0") + return throw(DomainError((a, b), "a and b are 0.0")) elseif x < 0.0 || x > 1.0 - return error("x < 0 or x > 1") + return throw(DomainError(x, "x < 0 or x > 1")) elseif y < 0.0 || y > 1.0 - return error("y < 0 or y > 1") + return throw(DomainError(y, "y < 0 or y > 1")) else z = x + y - 1.0 if abs(z) > 3.0*eps() - return error("x + y != 1.0") # ERROR HANDLING + return throw(DomainError((x, y), "x + y != 1.0")) # ERROR HANDLING end end if x == 0.0 - return (0.0,1.0) + return (0.0, 1.0) elseif y == 0.0 - return (1.0,0.0) + return (1.0, 0.0) elseif a == 0.0 - return (1.0,0.0) + return (1.0, 0.0) elseif b == 0.0 - return (0.0,1.0) + return (0.0, 1.0) end #EVALUATION OF ALGOS FOR PROPER SUB-DOMAINS ABOVE epps = max(eps(), 1.0e-15) - if max(a,b) < 1.0E-3 * epps - return (b/(a+b), a/(a+b)) + if max(a, b) < 1.0E-3 * epps + return (b/(a + b), a/(a + b)) end ind = false a0 = a @@ -745,9 +760,9 @@ function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) x0 = x y0 = y - if min(a0,b0) > 1.0 + if min(a0, b0) > 1.0 #PROCEDURE FOR A0>1 AND B0>1 - lambda = a > b ? (a+b)*y - b : a - (a+b)*x + lambda = a > b ? (a + b)*y - b : a - (a + b)*x if lambda < 0.0 ind = true a0 = b @@ -757,41 +772,41 @@ function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) lambda = abs(lambda) end if b0 < 40.0 && b0*x0 <= 0.7 - p = beta_inc_power_series(a0,b0,x0,epps) + p = beta_inc_power_series(a0, b0, x0, epps) q = 1.0 - p elseif b0 < 40.0 - n = trunc(Int, b0) - b0 -= float(n) - if b0 == 0.0 - n-=1 - b0=1.0 - end - p = beta_inc_diff(b0,a0,y0,x0,n,epps) - if x0 <= 0.7 - p += beta_inc_power_series(a0,b0,x0,epps) - q = 1.0 - p - else - if a0 <= 15.0 - n = 20 - p += beta_inc_diff(a0, b0, x0, y0, n, epps) - a0 += n - end - p = beta_inc_asymptotic_asymmetric(a0,b0,x0,y0,p,15.0*eps()) - q = 1.0 - p + n = trunc(Int, b0) + b0 -= n + if b0 == 0.0 + n -= 1 + b0 = 1.0 + end + p = beta_inc_diff(b0, a0, y0, x0, n, epps) + if x0 <= 0.7 + p += beta_inc_power_series(a0, b0, x0, epps) + q = 1.0 - p + else + if a0 <= 15.0 + n = 20 + p += beta_inc_diff(a0, b0, x0, y0, n, epps) + a0 += n end + p = beta_inc_asymptotic_asymmetric(a0, b0, x0, y0, p, 15.0*eps()) + q = 1.0 - p + end elseif a0 > b0 if b0 <= 100.0 || lambda > 0.03*b0 - p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps()) + p = beta_inc_cont_fraction(a0, b0, x0, y0, lambda, 15.0*eps()) q = 1.0 - p else - p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps()) + p = beta_inc_asymptotic_symmetric(a0, b0, lambda, 100.0*eps()) q = 1.0 - p end elseif a0 <= 100.0 || lambda > 0.03*a0 - p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps()) + p = beta_inc_cont_fraction(a0, b0, x0, y0, lambda, 15.0*eps()) q = 1.0 - p else - p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps()) + p = beta_inc_asymptotic_symmetric(a0, b0, lambda, 100.0*eps()) q = 1.0 - p end return ind ? (q, p) : (p, q) @@ -872,7 +887,7 @@ beta_inc(a::T, b::T, x::T, y::T) where {T<:AbstractFloat} = throw(MethodError(be #Volume 26, Number 1, 1977, pages 111-114. """ - beta_inc_inv(a,b,p,q,lb=logbeta(a,b)[1]) + beta_inc_inv(a, b, p, q, lb=logbeta(a,b)[1]) Computes inverse of incomplete beta function. Given `a`,`b` and ``I_{x}(a,b) = p`` find `x` and return tuple `(x,y)`. diff --git a/src/gamma.jl b/src/gamma.jl index a86ea640..db467f00 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -772,38 +772,9 @@ Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y) """ function beta end -function beta(a::Real, b::Real) - #Special case for negative integer argument - if a <= 0.0 - if isinteger(a) && 1-a-b > 0 - sgn = isinteger(b/2) ? 1 : -1 - return sgn* beta(1-a-b,b) - end - end - if b <= 0.0 - if isinteger(b) && 1-a-b > 0 - sgn = isinteger(a/2) ? 1 : -1 - return sgn* beta(1-a-b,a) - end - end - if a < b - a,b = b,a - end - #asymptotic expansion for log(B(a,b)) for |a| >> |b| - if abs(a) > 1e5*abs(b) && abs(a) > 1e5 - return exp(loggammadiv(b,a) + loggamma(b)) - end - ya, sa = logabsgamma(a) - yb, sb = logabsgamma(b) - yab, sab = logabsgamma(a+b) - return exp(ya + yb - yab) * (sa*sb*sab) -end - function beta(a::Number, b::Number) - ya = loggamma(a) - yb = loggamma(b) - yab = loggamma(a+b) - return exp(ya + yb - yab) + lab, sign = logabsbeta(a, b) + return sign*exp(lab) end """ @@ -813,7 +784,13 @@ Natural logarithm of the [`beta`](@ref) function ``\\log(|\\operatorname{B}(x,y) See also [`logabsbeta`](@ref). """ -logbeta(a::Number, b::Number) = loggamma(a)+loggamma(b)-loggamma(a+b) +function logbeta(a::Number, b::Number) + lab, sign = logabsbeta(a, b) + if sign < 0 + throw(DomainError((a, b), "`beta(a, b)` must be non-negative")) + end + return lab +end """ logabsbeta(x, y) @@ -822,31 +799,38 @@ Compute the natural logarithm of the absolute value of the [`beta`](@ref) functi See also [`logbeta`](@ref). """ -function logabsbeta(a::Real, b::Real) - if a <= 0.0 - if isinteger(a) && 1-a-b > 0 - sgn = isinteger(b/2) ? 1 : -1 - return logabsbeta(1-a-b,b) - end +function logabsbeta(a::T, b::T) where T<:Real + # ensure that a <= b + if a > b + return logabsbeta(b, a) end - if b <= 0.0 - if isinteger(b) && 1-a-b > 0 - sgn = isinteger(a/2) ? 1 : -1 - return logabsbeta(1-a-b,a) + + if a <= 0.0 && isinteger(a) + if a + b <= 0 + r = logbeta(1 - a - b, b) + sgn = isinteger(b/2) ? 1 : -1 + return r, sgn + else + return -log(zero(first(promote(a, b)))), 1 end end - if a < b - a,b = b,a - end - #asymptotic expansion for log(B(a,b)) for |a| >> |b| - if abs(a) > 1e5*abs(b) && abs(a) > 1e5 - return (loggammadiv(b,a) + loggamma(b)), 1 + + # asymptotic expansion for large `b` and positive arguments + # FIXME! We probably lose precision for negative arguments. However, `loggammadiv` + # is based on the NSWC code which doesn't bother with negative arguments + if a > 0 && b > 8 + return (loggammadiv(a, b) + loggamma(a)), 1 end + ya, sa = logabsgamma(a) yb, sb = logabsgamma(b) - yab, sab = logabsgamma(a+b) + yab, sab = logabsgamma(a + b) (ya + yb - yab), (sa*sb*sab) end +logabsbeta(a::Real, b::Real) = logabsbeta(promote(a, b)...) + +logabsbeta(a::Number, b::Number) = loggamma(a) + loggamma(b) - loggamma(a + b), 1 + ## from base/mpfr.jl # log of absolute value of gamma function diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 797c56db..d416dbd5 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -37,6 +37,7 @@ const d80 = -.652623918595309E-03 Uses the relation `gamma(a+1) = a*gamma(a)`. """ function rgamma1pm1(a::Float64) + @assert -0.5 <= a <= 1.5 t = a rangereduce = a > 0.5 t = rangereduce ? a-1 : a #-0.5<= t <= 0.5 @@ -86,6 +87,7 @@ end Compute function `g` in ``1/\\Gamma(x+1) = 1+x*(x-1)*g(x)``, `-1 <= x <= 1`. """ function auxgam(x::Float64) + @assert -1 <= x <= 1 if x < 0 return -(1.0 + (1.0 + x)*(1.0 + x)*auxgam(1.0 + x))/(1.0 - x) else @@ -95,11 +97,12 @@ function auxgam(x::Float64) end """ - loggamma1p(a) + loggamma1p(x) -Compute ``log(\\Gamma(1+a))`` for `-1 < a <= 1`. +Compute ``log(\\Gamma(1+x))`` for `-1 < x <= 1`. """ function loggamma1p(x::Float64) + @assert -1 < x <= 1 return -log1p(x*(x - 1.0)*auxgam(x)) end diff --git a/test/gamma.jl b/test/gamma.jl index d3c4f681..fff103d3 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -248,22 +248,35 @@ end @test first.(logabsbinomial.(200, 0:200)) ≈ log.(binomial.(BigInt(200), (0:200))) end -@testset "beta, lbeta" begin - @test beta(3/2,7/2) ≈ 5π/128 - @test beta(3,5) ≈ 1/105 - @test logbeta(5,4) ≈ log(beta(5,4)) - @test beta(5,4) ≈ beta(4,5) - @test beta(-2.0,2.0) ≈ 0.5 - @test logabsbeta(-2.0,2.0)[1] ≈ -0.69314718055994529 - @test beta(1e8,0.5) ≈ 0.00017724538531210809 - @test logabsbeta(1e8,0.5)[1] ≈ -8.637975427801484 - @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3 - @test logabsbeta(-1/2, 3)[1] ≈ log(16/3) - @test beta(Float32(5),Float32(4)) == beta(Float32(4),Float32(5)) - @test beta(3,5) ≈ beta(3+0im,5+0im) - @test(beta(3.2+0.1im,5.3+0.3im) ≈ exp(logbeta(3.2+0.1im,5.3+0.3im)) ≈ +@testset "beta, logbeta, and logabsbeta" begin + @test beta(3/2, 7/2) ≈ 5π/128 + @test beta(3, 5) ≈ 1/105 + @test logbeta(5, 4) ≈ log(beta(5,4)) + @test beta(5, 4) ≈ beta(4,5) + + @testset "negative integer argument" begin + @test beta(-2.0, 2.0) ≈ 1/2 rtol=1e-14 + @test beta(-5.0, 2.0) ≈ 1/20 rtol=1e-14 + @test logabsbeta(-2.0, 2.0)[1] ≈ -log(2) rtol=1e-14 + @test beta(-2.0, -2.0) == Inf + @test logbeta(-2.0, -2.0) == Inf + @test beta(-2.0, 2.1) == Inf + @test logbeta(-2.0, 2.1) == Inf + end + + @testset "large difference in magnitude" begin + @test beta(1e4, 1.5) ≈ 8.86193693673874630607029e-7 rtol=1e-14 + @test logabsbeta(1e4, 1.5)[1] ≈ log(8.86193693673874630607029e-7) rtol=1e-14 + @test beta(1e8, 0.5) ≈ 0.00017724538531210809 rtol=1e-14 + @test logabsbeta(1e8, 0.5)[1] ≈ log(0.00017724538531210809) rtol=1e-14 + end + + @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3 + @test logabsbeta(-1/2, 3)[1] ≈ log(16/3) + @test beta(Float32(5), Float32(4)) == beta(Float32(4), Float32(5)) + @test beta(3, 5) ≈ beta(3 + 0im, 5 + 0im) + @test(beta(3.2 + 0.1im, 5.3 + 0.3im) ≈ exp(logbeta(3.2 + 0.1im, 5.3 + 0.3im)) ≈ 0.00634645247782269506319336871208405439180447035257028310080 - 0.00169495384841964531409376316336552555952269360134349446910im) - - @test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps() + @test beta(big(1.0), big(1.2)) ≈ beta(1.0, 1.2) rtol=4*eps() end From 226b28a4140056e93a09434a52f1a9995b5f0ed9 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 17:35:54 +0200 Subject: [PATCH 02/11] Update src/beta_inc.jl Co-authored-by: David Widmann --- src/beta_inc.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 6a9183f4..28e0785f 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -50,10 +50,9 @@ Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) for a0, b0 >= 8 """ function stirling_corr(a0::Float64, b0::Float64) - @assert a0 >= 0 - @assert b0 >= 0 a = min(a0, b0) b = max(a0, b0) + @assert a >= 8 h = a/b c = h/(1.0 + h) From 01f55e5882311e23c6ba63ecd222d1e32e78c201 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 17:41:56 +0200 Subject: [PATCH 03/11] Update src/beta_inc.jl Co-authored-by: David Widmann --- src/beta_inc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 28e0785f..8554b4a5 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -432,7 +432,7 @@ See also: [`beta_inc`](@ref) `FPSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64) - @assert b < eps(Float64)*min(1.0, a) + @assert b < epps*min(1.0, a) @assert x <= 0.5 ans = 1.0 if a > 1.0e-3 * epps From 7bcfd01a15fe5e02e77ee6de57d43d32dd561751 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 17:42:09 +0200 Subject: [PATCH 04/11] Update src/beta_inc.jl Co-authored-by: David Widmann --- src/beta_inc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 8554b4a5..fec18135 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -476,7 +476,7 @@ See also: [`beta_inc`](@ref) `APSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64) - @assert a <= eps(Float64)*min(1.0, b) + @assert a <= epps*min(1.0, b) @assert b*x <= 1 @assert x <= 0.5 g = Base.MathConstants.eulergamma From ec68aaa62c4f173b039de4033a18dd75677d5ad3 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 17:43:27 +0200 Subject: [PATCH 05/11] Update src/gamma.jl Co-authored-by: David Widmann --- src/gamma.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gamma.jl b/src/gamma.jl index db467f00..b2685361 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -805,7 +805,7 @@ function logabsbeta(a::T, b::T) where T<:Real return logabsbeta(b, a) end - if a <= 0.0 && isinteger(a) + if a <= 0 && isinteger(a) if a + b <= 0 r = logbeta(1 - a - b, b) sgn = isinteger(b/2) ? 1 : -1 From 443376d9dfebe70e4f616e63aa952ca2c8caca71 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 18:01:30 +0200 Subject: [PATCH 06/11] Update src/gamma.jl Co-authored-by: David Widmann --- src/beta_inc.jl | 32 ++++++++++++++++---------------- src/gamma.jl | 6 +++--- test/gamma.jl | 7 +++++-- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index fec18135..fbfe690c 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -630,7 +630,7 @@ function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Intege d = exp(-t) end - ans = beta_integrand(a,b,x,y,mu)/a + ans = beta_integrand(a, b, x, y, mu)/a if n == 1 || ans == 0.0 return ans end @@ -820,53 +820,53 @@ function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64) end if b0 < min(epps, epps*a0) - p = beta_inc_power_series2(a0,b0,x0,epps) + p = beta_inc_power_series2(a0, b0, x0, epps) q = 1.0 - p elseif a0 < min(epps, epps*b0) && b0*x0 <= 1.0 - q = beta_inc_power_series1(a0,b0,x0,epps) + q = beta_inc_power_series1(a0, b0, x0, epps) p = 1.0 - q elseif max(a0,b0) > 1.0 if b0 <= 1.0 - p = beta_inc_power_series(a0,b0,x0,epps) + p = beta_inc_power_series(a0, b0, x0, epps) q = 1.0 - p elseif x0 >= 0.3 - q = beta_inc_power_series(b0,a0,y0,epps) + q = beta_inc_power_series(b0, a0, y0, epps) p = 1.0 - q elseif x0 >= 0.1 if b0 > 15.0 - q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps()) + q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps()) p = 1.0 - q else n = 20 - q = beta_inc_diff(b0,a0,y0,x0,n,epps) + q = beta_inc_diff(b0, a0, y0, x0, n, epps) b0 += n - q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps()) + q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps()) p = 1.0 - q end elseif (x0*b0)^(a0) <= 0.7 - p = beta_inc_power_series(a0,b0,x0,epps) + p = beta_inc_power_series(a0, b0, x0, epps) q = 1.0 - p else n = 20 - q = beta_inc_diff(b0,a0,y0,x0,n,epps) + q = beta_inc_diff(b0, a0, y0, x0, n, epps) b0 += n - q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps()) + q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps()) p = 1.0 - q end elseif a0 >= min(0.2, b0) - p = beta_inc_power_series(a0,b0,x0,epps) + p = beta_inc_power_series(a0, b0, x0, epps) q = 1.0 - p elseif x0^a0 <= 0.9 - p = beta_inc_power_series(a0,b0,x0,epps) + p = beta_inc_power_series(a0, b0, x0, epps) q = 1.0 - p elseif x0 >= 0.3 - q = beta_inc_power_series(b0,a0,y0,epps) + q = beta_inc_power_series(b0, a0, y0, epps) p = 1.0 - q else n = 20 - q = beta_inc_diff(b0,a0,y0,x0,n,epps) + q = beta_inc_diff(b0, a0, y0, x0, n, epps) b0 += n - q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps()) + q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps()) p = 1.0 - q end diff --git a/src/gamma.jl b/src/gamma.jl index b2685361..8aaa7e70 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -806,12 +806,12 @@ function logabsbeta(a::T, b::T) where T<:Real end if a <= 0 && isinteger(a) - if a + b <= 0 + if a + b <= 0 && isinteger(b) r = logbeta(1 - a - b, b) - sgn = isinteger(b/2) ? 1 : -1 + sgn = iseven(Int(b)) ? 1 : -1 return r, sgn else - return -log(zero(first(promote(a, b)))), 1 + return -log(zero(a)), 1 end end diff --git a/test/gamma.jl b/test/gamma.jl index fff103d3..790b6718 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -255,11 +255,14 @@ end @test beta(5, 4) ≈ beta(4,5) @testset "negative integer argument" begin - @test beta(-2.0, 2.0) ≈ 1/2 rtol=1e-14 - @test beta(-5.0, 2.0) ≈ 1/20 rtol=1e-14 + @test beta(-2.0, 1.0) ≈ -1/2 rtol=1e-14 + @test beta(-2.0, 2.0) ≈ 1/2 rtol=1e-14 + @test beta(-5.0, 2.0) ≈ 1/20 rtol=1e-14 @test logabsbeta(-2.0, 2.0)[1] ≈ -log(2) rtol=1e-14 @test beta(-2.0, -2.0) == Inf @test logbeta(-2.0, -2.0) == Inf + @test beta(-2.0, 1.9) == Inf + @test logbeta(-2.0, 1.9) == Inf @test beta(-2.0, 2.1) == Inf @test logbeta(-2.0, 2.1) == Inf end From 5dd89c86a01888cf6d602f908726eced2f51c90e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 21:48:52 +0200 Subject: [PATCH 07/11] Improve test coverage for beta_inc --- test/beta_inc.jl | 223 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 203 insertions(+), 20 deletions(-) diff --git a/test/beta_inc.jl b/test/beta_inc.jl index e3d70d97..ab51c6dd 100644 --- a/test/beta_inc.jl +++ b/test/beta_inc.jl @@ -1,30 +1,213 @@ @testset "incomplete beta" begin -#Compared with scipy.special.betainc - ans1 = [0.063768560858519868, 0.090334470601733122, 0.11082468660445946, 0.12818843369794991, 0.14356629312870628, 0.15754242542404437, 0.17046342838319906, 0.18255489099382846, 0.19397336804135662, 0.2048327646991335, 0.21521902554364258, 0.22519890061842537, 0.23482546905318047, 0.24414177534616024, 0.253183311106635, 0.26197976086890934, 0.27055626226830082, 0.27893433611357099, 0.28713258625741256, 0.29516723530086664, 0.30305254089347, 0.31080112365240059, 0.3184242286319004, 0.32593193612447635, 0.33333333333333337, 0.34063665547711719, 0.34784940276303394, 0.35497843812588953, 0.36203006950165795, 0.36901011956554541, 0.375923985234415, 0.3827766887550389, 0.38957292183289544, 0.39631708397254195, 0.40301331597932183, 0.409665529398267, 0.41627743252839544, 0.42285255354095019, 0.42939426114224771, 0.43590578315102518, 0.44239022330290334, 0.44885057654807836, 0.45528974307049164, 0.46171054122585353, 0.46811571957074011, 0.47450796813452323, 0.48088992906936007, 0.48726420680024518, 0.49363337778673011, 0.50000000000000011, 0.50636662221327011, 0.51273579319975504, 0.5191100709306401, 0.52549203186547688, 0.53188428042926006, 0.53828945877414658, 0.54471025692950836, 0.55114942345192175, 0.55760977669709688, 0.56409421684897509, 0.57060573885775256, 0.57714744645905003, 0.58372256747160478, 0.59033447060173316, 0.59698668402067834, 0.60368291602745827, 0.61042707816710473, 0.61722331124496133, 0.62407601476558505, 0.63098988043445459, 0.63796993049834205, 0.64502156187411064, 0.65215059723696633, 0.65936334452288281, 0.66666666666666663, 0.67406806387552365, 0.68157577136809966, 0.68919887634759947, 0.69694745910653, 0.70483276469913336, 0.71286741374258733, 0.72106566388642879, 0.72944373773169924, 0.73802023913109061, 0.74681668889336494, 0.75585822465383967, 0.76517453094681931, 0.7748010993815746, 0.78478097445635719, 0.79516723530086619, 0.8060266319586431, 0.81744510900617107, 0.82953657161680039, 0.84245757457595483, 0.8564337068712925, 0.87181156630205003, 0.88917531339554046, 0.90966552939826684, 0.9362314391414801] - ans2 = [0.087041654488676928, 0.12317839372423395, 0.15096439050993934, 0.1744376469592186, 0.19516124171507579, 0.21393635267239447, 0.23123866443564164, 0.24737796264883508, 0.26256933221933759, 0.27696931578459572, 0.29069605231718387, 0.30384130015177169, 0.31647801578481272, 0.32866534497198974, 0.34045203033834726, 0.35187880931248239, 0.36298014555988484, 0.37378550722502435, 0.38432032902829749, 0.39460674882922503, 0.40466418009101346, 0.41450976283569702, 0.42415872320430908, 0.43362466329673527, 0.44291979714289581, 0.4520551445675528, 0.46104069179265272, 0.46988552550762908, 0.47859794558683333, 0.48718556048052553, 0.49565536843945707, 0.50401382707500741, 0.51226691325208229, 0.52042017492130455, 0.52847877619206929, 0.53644753670811252, 0.54433096619706123, 0.55213329491367624, 0.55985850057457809, 0.56751033228373404, 0.57509233186792286, 0.58260785297598572, 0.59006007824198947, 0.59745203476815989, 0.60478660814681062, 0.61206655521005915, 0.61929451567078364, 0.62647302279711903, 0.63360451324511835, 0.64069133615944129, 0.64773576163960578, 0.65473998865910954, 0.66170615251627152, 0.66863633188876026, 0.67553255555825598, 0.68239680886742871, 0.68923103996827928, 0.69603716591885445, 0.70281707868434762, 0.70957265109866863, 0.71630574284372639, 0.72301820650601145, 0.72971189377369172, 0.73638866184254259, 0.74305038010582547, 0.74969893721203962, 0.75633624858569781, 0.76296426452045429, 0.76958497897174627, 0.77620043919849924, 0.78281275643163062, 0.78942411778262311, 0.79603679965050977, 0.80265318294305021, 0.80927577050162136, 0.81590720721472276, 0.82255030342945812, 0.82920806243432144, 0.83588371300491959, 0.84258074829830698, 0.84930297278272149, 0.85605455944435893, 0.86284012029219526, 0.86966479429526578, 0.876534358507564, 0.8834553705447894, 0.89043535424294118, 0.89748304605629159, 0.90460872897483335, 0.91182469611459804, 0.91914591281205937, 0.92659099459483496, 0.93418371187592941, 0.94195542513481079, 0.94994928861995065, 0.95822815489579449, 0.96689132900746566, 0.97611719389699281, 0.98631372146716156] - ans3 = [0.012873483750773953, 0.024045778495719504, 0.034668868563965893, 0.044957919792848362, 0.05501104929272229, 0.064884481926849022, 0.074614610511183635, 0.084226884888052606, 0.093740088419712889, 0.10316864948642543, 0.11252399935922629, 0.12181542160322313, 0.13105061024399095, 0.14023605142518913, 0.14937729307408915, 0.15847914076041122, 0.16754580333263841, 0.1765810034326174, 0.18558806286120144, 0.19456996956034567, 0.20352943091060927, 0.21246891667645043, 0.22139069400685801, 0.23029685625993337, 0.23918934697038266, 0.24806997995710448, 0.2569404563342797, 0.26580237901716475, 0.27465726518533057, 0.28350655706915756, 0.29235163135146713, 0.30119380741925028, 0.31003435465621887, 0.31887449893226855, 0.32771542841860729, 0.33655829883559579, 0.34540423822300514, 0.3542543513084796, 0.36310972353877363, 0.37197142482927481, 0.38084051307999506, 0.38971803750029138, 0.39860504177981332, 0.40750256713936639, 0.41641165529238072, 0.42533335134536215, 0.43426870666399331, 0.44321878173037438, 0.45218464901620103, 0.46116739589644162, 0.47016812762827714, 0.47918797042070466, 0.48822807462128687, 0.49728961804808525, 0.50637380949686728, 0.51548189245630038, 0.5246151490670925, 0.53377490436500419, 0.5429625308524717, 0.55217945344936725, 0.56142715488038319, 0.57070718156489009, 0.58002115008515442, 0.58937075432090513, 0.59875777335283087, 0.60818408025528292, 0.61765165191996785, 0.62716258007870107, 0.63671908372557429, 0.6463235231787503, 0.65597841607162255, 0.66568645562500284, 0.67545053162996704, 0.68527375466991214, 0.6951594842368839, 0.70511136156043297, 0.71513334817976937, 0.72522977156950186, 0.73540537950094409, 0.74566540532108994, 0.75601564701303459, 0.76646256384385103, 0.77701339572859029, 0.78767631232750368, 0.79846060164096611, 0.80937691194811789, 0.82043756714268778, 0.83165698521328602, 0.84305224521672972, 0.85464387408595721, 0.86645696969790587, 0.87852285861624935, 0.89088164472968334, 0.90358633052105619, 0.91670992507941651, 0.93035879743560423, 0.94470095224765227, 0.96003789419656471, 0.97705809868034199] - ans4 = [5.6598501882086638e-163, 1.2793886161240335e-138, 2.2576938967074611e-124, 2.8979946053659182e-114, 2.0091646864167303e-106, 5.1249028123107484e-100, 1.3384365465929538e-94, 6.5928463236275696e-90, 9.0834048016126672e-86, 4.5811711646109305e-82, 1.024712535445947e-78, 1.1712904358478081e-75, 7.61899830618238e-73, 3.0663999960291821e-70, 8.1593401517362853e-68, 1.5142302711017038e-65, 2.0475200955891378e-63, 2.0916713536958276e-61, 1.6640409038082517e-59, 1.057759838815703e-57, 5.4913266744369942e-56, 2.3725786887367975e-54, 8.6718335278405774e-53, 2.7197933177013249e-51, 7.4118251807120739e-50, 1.7744725579765609e-48, 3.7689144668459763e-47, 7.1639177579230802e-46, 1.2281691523716354e-44, 1.9124082898174377e-43, 2.7218251674439459e-42, 3.5610505888183843e-41, 4.3051319713293141e-40, 4.8320967998003304e-39, 5.0570559392737623e-38, 4.954351648811439e-37, 4.5601168072923055e-36, 3.9565121052116655e-35, 3.2458692695921472e-34, 2.5250201231776632e-33, 1.8674800424552049e-32, 1.31632743396051e-31, 8.862853181709709e-31, 5.71217922566864e-30, 3.5310511291576422e-29, 2.0973893941833731e-28, 1.1991528375926732e-27, 6.6098385509764045e-27, 3.5179140837714274e-26, 1.8104027085942885e-25, 9.0207305335112628e-25, 4.3574464645305935e-24, 2.0429719208160492e-23, 9.307222587504958e-23, 4.124448596486994e-22, 1.7796548980836483e-21, 7.4841565165421277e-21, 3.0702859704148574e-20, 1.2297464356858357e-19, 4.8128959875339663e-19, 1.8419917152708652e-18, 6.8988896144258444e-18, 2.5303817306725186e-17, 9.094925675449239e-17, 3.2055039961949189e-16, 1.1085187441597278e-15, 3.7635118688411301e-15, 1.255134520239119e-14, 4.1140267831674448e-14, 1.3260121388841605e-13, 4.2048092809091946e-13, 1.3124119186938038e-12, 4.0338344386888606e-12, 1.2214673719698131e-11, 3.6454138300247904e-11, 1.0727343115700782e-10, 3.1138059245591031e-10, 8.9189474914175862e-10, 2.5218693952601578e-09, 7.0417217701699225e-09, 1.9424040709629306e-08, 5.2949202359208791e-08, 1.4268930041189607e-07, 3.8026633846818863e-07, 1.0025399986934807e-06, 2.6156990305099061e-06, 6.756237677059056e-06, 1.7282977043885323e-05, 4.3803033456977776e-05, 0.00011004048858461744, 0.0002741427949192326, 0.00067767979451775295, 0.0016633950407366552, 0.0040576338171432955, 0.009848744173397107, 0.023828559323224537, 0.057638830150965534, 0.14020121761365209, 0.34813332055564178] - ans5 = [0.013962049747018378, 0.042738447399963485, 0.080235508245346435, 0.12332127072222762, 0.16988388360150664, 0.21838766458819614, 0.26767957944384452, 0.31688247132531383, 0.36532771975839284, 0.41250897644193646, 0.45804848436293449, 0.50167153949877097, 0.54318655272778527, 0.58246914958502727, 0.61944929138313776, 0.65410072456419888, 0.68643226670884316, 0.71648056882768607, 0.74430408220004141, 0.76997801988794912, 0.79359014748064294, 0.8152372703488121, 0.83502230934244948, 0.85305187581929465, 0.86943427172671017, 0.88427785226197508, 0.89768969816152788, 0.90977455245768279, 0.92063398298239707, 0.93036573727946092, 0.93906326112402183, 0.94681535570471753, 0.95370595182410001, 0.95981398231488424, 0.96521333632950712, 0.96997288129947024, 0.97415654022800602, 0.97782341361450698, 0.9810279367445075, 0.98382006434176061, 0.98624547569157095, 0.98834579432593017, 0.99015881722713328, 0.99171874927107939, 0.993056439306095, 0.99419961485797559, 0.99517311297582334, 0.99599910519385948, 0.99669731498846048, 0.99728522646317797, 0.99777828330274299, 0.9981900773047403, 0.99853252602897213, 0.99881603930330098, 0.99904967449435056, 0.99924128059493644, 0.99939763130024395, 0.99952454734409613, 0.99962700844743568, 0.99970925529544441, 0.99977488200944087, 0.99982691961654058, 0.99986791104559958, 0.99989997819362175, 0.99992488161389415, 0.9999440733768028, 0.99995874364766302, 0.99996986151394762, 0.9999782105779107, 0.99998441981059161, 0.99998899014027587, 0.99999231722334625, 0.99999471081867697, 0.99999641115882809, 0.99999760268277504, 0.99999842546616124, 0.99999898465647652, 0.99999935819244778, 0.99999960305956725, 0.99999976030730964, 0.99999985902839927, 0.99999991947664768, 0.99999995547750509, 0.9999999762646713, 0.99999998785693356, 0.99999999407190587, 0.99999999725752275, 0.99999999880800261, 0.99999999951849838, 0.99999999982174348, 0.99999999994062116, 0.99999999998264033, 0.99999999999569866, 0.9999999999991418, 0.99999999999987266, 0.99999999999998768, 0.99999999999999944, 0.99999999999999989, 0.99999999999999989] - ans6 = [3.3560143825993687e-143, 2.213518371644144e-113, 4.0145526832291446e-96, 5.2109131251766219e-84, 1.0203891990294805e-74, 3.2677379646480312e-67, 6.1015737258099584e-61, 1.4185802923209925e-55, 6.6985847646944376e-51, 8.9719726417923444e-47, 4.3251435416063278e-43, 8.9466478285748974e-40, 9.0742566193142145e-37, 5.0057011672920701e-34, 1.6303369052296386e-31, 3.3493886577780779e-29, 4.5813387723702873e-27, 4.3628327684456486e-25, 3.0028168459240925e-23, 1.5416424964057036e-21, 6.0647292131758379e-20, 1.8708373622428274e-18, 4.6165844488973474e-17, 9.2725535948439169e-16, 1.5390751973611367e-14, 2.1394034328449769e-13, 2.5200592719451178e-12, 2.5418366766245469e-11, 2.215842537778695e-10, 1.6834088907362684e-09, 1.122874450436064e-08, 6.6201731659826321e-08, 3.4707638415115294e-07, 1.6269341185131871e-06, 6.8526269592546404e-06, 2.6052419808499182e-05, 8.97713031375554e-05, 0.0002814332703651891, 0.00080554169029503023, 0.0021120197810150015, 0.0050879707571202049, 0.01129547597461162, 0.023174688420508269, 0.044064890651698461, 0.077870360082158169, 0.12827033082474237, 0.19756362321699061, 0.28548625344477141, 0.38850225052293164, 0.50000000000002465, 0.61149774947706836, 0.71451374655522859, 0.80243637678300939, 0.87172966917525962, 0.92212963991784314, 0.95593510934830217, 0.97682531157949204, 0.98870452402538833, 0.99491202924287991, 0.99788798021898495, 0.99919445830970499, 0.99971856672963477, 0.99991022869686241, 0.99997394758019154, 0.99999314737304079, 0.99999837306588146, 0.99999965292361581, 0.99999993379826835, 0.99999998877125551, 0.99999999831659114, 0.99999999977841569, 0.99999999997458167, 0.9999999999974799, 0.99999999999978606, 0.99999999999998457, 0.99999999999999911, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989] - @testset "a=0.5, b=0.5, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans1) - @test beta_inc(0.5, 0.5, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(0.5,0.5,x) + # Compared with Wolfram Engine's BetaRegularized + @testset "a=0.5, b=0.5, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-25,0.06376856085851985,0.0903344706017331,0.1108246866044594,0.1281884336979499, + 0.1435662931287063,0.1575424254240443,0.170463428383199,0.1825548909938284,0.1939733680413566, + 0.2048327646991335,0.2152190255436425,0.2251989006184253,0.2348254690531804,0.2441417753461602, + 0.253183311106635,0.2619797608689092,0.2705562622683007,0.2789343361135709,0.2871325862574125, + 0.2951672353008665,0.3030525408934699,0.3108011236524005,0.3184242286319003,0.3259319361244763, + 0.3333333333333333,0.3406366554771171,0.3478494027630338,0.3549784381258894,0.3620300695016579, + 0.3690101195655454,0.375923985234415,0.3827766887550388,0.3895729218328954,0.3963170839725418, + 0.4030133159793217,0.4096655293982669,0.4162774325283954,0.42285255354095,0.4293942611422476, + 0.4359057831510251,0.4423902233029032,0.4488505765480783,0.4552897430704916,0.4617105412258534, + 0.4681157195707401,0.4745079681345231,0.48088992906936,0.4872642068002451,0.49363337778673,0.5, + 0.50636662221327,0.5127357931997549,0.51911007093064,0.5254920318654769,0.5318842804292599, + 0.5382894587741466,0.5447102569295084,0.5511494234519217,0.5576097766970968,0.5640942168489749, + 0.5706057388577524,0.57714744645905,0.5837225674716046,0.5903344706017331,0.5969866840206783, + 0.6036829160274582,0.6104270781671046,0.6172233112449612,0.624076014765585,0.6309898804344546, + 0.6379699304983421,0.6450215618741106,0.6521505972369662,0.6593633445228829,0.6666666666666667, + 0.6740680638755237,0.6815757713680997,0.6891988763475995,0.6969474591065301,0.7048327646991335, + 0.7128674137425875,0.7210656638864291,0.7294437377316993,0.7380202391310908,0.746816688893365, + 0.7558582246538398,0.7651745309468196,0.7748010993815747,0.7847809744563575,0.7951672353008665, + 0.8060266319586434,0.8174451090061716,0.829536571616801,0.8424575745759557,0.8564337068712937, + 0.8718115663020501,0.8891753133955406,0.9096655293982669,0.9362314391414802,1.]) + @test beta_inc(0.5, 0.5, x, 1.0 - x)[1] ≈ val rtol=1e-14 end - @testset "a=0.5, b=0.8, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans2) - @test beta_inc(0.5, 0.8, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(0.5,0.8,x) + @testset "a=0.5, b=0.8, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-24,0.08704165448867688,0.1231783937242339,0.1509643905099393,0.1744376469592185, + 0.1951612417150757,0.2139363526723944,0.2312386644356416,0.247377962648835,0.2625693322193375, + 0.2769693157845956,0.2906960523171838,0.3038413001517716,0.3164780157848126,0.3286653449719896, + 0.3404520303383471,0.3518788093124823,0.3629801455598847,0.3737855072250243,0.3843203290282974, + 0.394606748829225,0.4046641800910133,0.4145097628356969,0.424158723204309,0.4336246632967352, + 0.4429197971428956,0.4520551445675527,0.4610406917926526,0.4698855255076289,0.4785979455868333, + 0.4871855604805254,0.4956553684394569,0.5040138270750072,0.5122669132520822,0.5204201749213043, + 0.5284787761920692,0.5364475367081124,0.5443309661970611,0.5521332949136762,0.5598585005745779, + 0.5675103322837338,0.5750923318679227,0.5826078529759856,0.5900600782419893,0.5974520347681598, + 0.6047866081468104,0.612066555210059,0.6192945156707836,0.6264730227971188,0.6336045132451182, + 0.6406913361594411,0.6477357616396057,0.6547399886591094,0.6617061525162714,0.6686363318887601, + 0.6755325555582559,0.6823968088674286,0.6892310399682792,0.6960371659188544,0.7028170786843475, + 0.7095726510986684,0.7163057428437263,0.7230182065060113,0.7297118937736915,0.7363886618425424, + 0.7430503801058253,0.7496989372120396,0.7563362485856976,0.7629642645204542,0.7695849789717461, + 0.7762004391984992,0.7828127564316306,0.7894241177826232,0.7960367996505097,0.8026531829430502, + 0.8092757705016214,0.8159072072147228,0.8225503034294581,0.8292080624343213,0.8358837130049195, + 0.8425807482983071,0.8493029727827216,0.8560545594443592,0.8628401202921954,0.869664794295266, + 0.8765343585075641,0.8834553705447895,0.8904353542429415,0.8974830460562918,0.9046087289748337, + 0.9118246961145984,0.9191459128120601,0.9265909945948355,0.9341837118759302,0.941955425134812, + 0.9499492886199519,0.9582281548957946,0.9668913290074657,0.9761171938969928,0.9863137214671616,1.]) + @test beta_inc(0.5, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-14 end - @testset "a=0.9, b=0.8, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans3) - @test beta_inc(0.9, 0.8, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(0.9,0.8,x) + @testset "a=0.9, b=0.8, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-45,0.01287348375077395,0.0240457784957195,0.0346688685639659,0.04495791979284836, + 0.0550110492927223,0.06488448192684902,0.07461461051118363,0.08422688488805261,0.0937400884197129, + 0.1031686494864254,0.1125239993592263,0.1218154216032231,0.1310506102439909,0.1402360514251891, + 0.1493772930740891,0.1584791407604112,0.1675458033326384,0.1765810034326174,0.1855880628612014, + 0.1945699695603456,0.2035294309106093,0.2124689166764504,0.2213906940068579,0.2302968562599333, + 0.2391893469703827,0.2480699799571045,0.2569404563342797,0.2658023790171647,0.2746572651853306, + 0.2835065570691575,0.2923516313514671,0.3011938074192502,0.3100343546562188,0.3188744989322686, + 0.3277154284186073,0.3365582988355958,0.3454042382230052,0.3542543513084796,0.3631097235387736, + 0.3719714248292748,0.3808405130799951,0.3897180375002914,0.3986050417798133,0.4075025671393664, + 0.4164116552923807,0.4253333513453621,0.4342687066639933,0.4432187817303744,0.452184649016201, + 0.4611673958964416,0.4701681276282771,0.4791879704207046,0.4882280746212869,0.4972896180480853, + 0.5063738094968672,0.5154818924563003,0.5246151490670925,0.5337749043650042,0.5429625308524718, + 0.5521794534493672,0.5614271548803832,0.5707071815648901,0.5800211500851545,0.5893707543209051, + 0.5987577733528309,0.608184080255283,0.6176516519199679,0.6271625800787011,0.6367190837255744, + 0.6463235231787504,0.6559784160716227,0.6656864556250031,0.6754505316299672,0.6852737546699122, + 0.695159484236884,0.7051113615604331,0.7151333481797695,0.7252297715695021,0.7354053795009442, + 0.74566540532109,0.7560156470130347,0.7664625638438512,0.7770133957285906,0.7876763123275041, + 0.7984606016409664,0.8093769119481182,0.8204375671426882,0.8316569852132864,0.8430522452167303, + 0.8546438740859578,0.8664569696979065,0.8785228586162501,0.8908816447296843,0.9035863305210574, + 0.9167099250794181,0.9303587974356043,0.9447009522476523,0.9600378941965647,0.977058098680342,1.]) + @test beta_inc(0.9, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-14 end - @testset "a=80.9, b=0.8, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans4) - @test beta_inc(80.9, 0.8, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(80.9,0.8,x) + @testset "a=80.9, b=0.8, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-4037,5.659850188208874e-163,1.279388616124076e-138,2.257693896707542e-124, + 2.897994605366003e-114,2.009164686416781e-106,5.124902812310911e-100,1.338436546592981e-94, + 6.592846323627736e-90,9.083404801612933e-86,4.581171164611029e-82,1.024712535445972e-78, + 1.171290435847841e-75,7.618998306182543e-73,3.066399996029232e-70,8.159340151736502e-68, + 1.514230271101736e-65,2.047520095589173e-63,2.091671353695881e-61,1.664040903808287e-59, + 1.057759838815721e-57,5.49132667443713e-56,2.372578688736847e-54,8.671833527840728e-53, + 2.719793317701389e-51,7.411825180712226e-50,1.774472557976592e-48,3.768914466846032e-47, + 7.163917757923168e-46,1.228169152371666e-44,1.912408289817481e-43,2.721825167444e-42, + 3.561050588818447e-41,4.305131971329379e-40,4.832096799800394e-39,5.057055939273883e-38, + 4.954351648811546e-37,4.560116807292394e-36,3.956512105211734e-35,3.245869269592197e-34, + 2.525020123177697e-33,1.867480042455247e-32,1.316327433960537e-31,8.862853181709876e-31, + 5.712179225668737e-30,3.531051129157696e-29,2.097389394183402e-28,1.199152837592699e-27, + 6.609838550976534e-27,3.517914083771492e-26,1.810402708594319e-25,9.020730533511399e-25, + 4.357446464530653e-24,2.042971920816074e-23,9.307222587505059e-23,4.124448596487034e-22, + 1.779654898083663e-21,7.4841565165423e-21,3.070285970414923e-20,1.22974643568586e-19, + 4.812895987534055e-19,1.841991715270897e-18,6.898889614425955e-18,2.530381730672556e-17, + 9.094925675449364e-17,3.205503996194959e-16,1.108518744159741e-15,3.763511868841169e-15, + 1.255134520239131e-14,4.114026783167532e-14,1.326012138884187e-13,4.204809280909274e-13, + 1.312411918693827e-12,4.033834438688928e-12,1.221467371969832e-11,3.645413830024844e-11, + 1.072734311570093e-10,3.113805924559142e-10,8.918947491417688e-10,2.521869395260185e-9, + 7.041721770169991e-9,1.942404070962948e-8,5.294920235920981e-8,1.426893004118986e-7, + 3.802663384681951e-7,1.002539998693497e-6,2.615699030509946e-6,6.756237677059152e-6, + 0.00001728297704388556,0.00004380303345697832,0.0001100404885846187,0.0002741427949192357, + 0.0006776797945177599,0.001663395040736671,0.004057633817143372,0.009848744173397283, + 0.02382855932322502,0.05763883015096727,0.1402012176136536,0.348133320555652,1.]) + @test beta_inc(80.9, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-13 end - @testset "a=1.7, b=10.5, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans5) - @test beta_inc(1.7, 10.5, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(1.7,10.5,x) + @testset "a=1.7, b=10.5, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-83,0.01396204974701841,0.04273844739996358,0.0802355082453466,0.1233212707222279, + 0.169883883601507,0.2183876645881966,0.2676795794438451,0.3168824713253145,0.3653277197583937, + 0.4125089764419372,0.4580484843629353,0.5016715394987723,0.5431865527277862,0.5824691495850261, + 0.6194492913831362,0.654100724564198,0.6864322667088425,0.7164805688276856,0.7443040822000412, + 0.7699780198879489,0.7935901474806426,0.8152372703488115,0.8350223093424488,0.8530518758192943, + 0.8694342717267099,0.8842778522619746,0.8976896981615276,0.9097745524576825,0.9206339829823968, + 0.9303657372794607,0.9390632611240216,0.9468153557047174,0.9537059518240998,0.9598139823148841, + 0.9652133363295071,0.9699728812994702,0.974156540228006,0.9778234136145069,0.9810279367445075, + 0.9838200643417606,0.9862454756915709,0.9883457943259302,0.9901588172271333,0.9917187492710794, + 0.993056439306095,0.9941996148579756,0.9951731129758233,0.9959991051938595,0.9966973149884605, + 0.997285226463178,0.997778283302743,0.9981900773047402,0.9985325260289722,0.998816039303301, + 0.9990496744943506,0.9992412805949365,0.9993976313002439,0.9995245473440961,0.9996270084474357, + 0.9997092552954444,0.9997748820094408,0.9998269196165406,0.9998679110455995,0.9998999781936217, + 0.9999248816138942,0.9999440733768028,0.9999587436476631,0.9999698615139476,0.9999782105779107, + 0.9999844198105916,0.9999889901402758,0.9999923172233462,0.9999947108186769,0.9999964111588281, + 0.9999976026827751,0.9999984254661612,0.9999989846564765,0.9999993581924477,0.9999996030595673, + 0.9999997603073096,0.9999998590283993,0.9999999194766476,0.9999999554775051,0.9999999762646713, + 0.9999999878569335,0.9999999940719059,0.9999999972575227,0.9999999988080026,0.999999999518498, + 0.999999999821743,0.999999999940621,0.99999999998264,0.999999999995699,0.999999999999142, + 0.999999999999873,0.999999999999988,0.999999999999999,1.,1.,1.]) + @test beta_inc(1.7, 10.5, x, 1.0 - x)[1] ≈ val rtol=1e-14 end - @testset "a=100.5, b=100.5, x=$x" for (x, val) in zip(0.01:0.01:0.99, ans6) - @test beta_inc(100.5, 100.5, x, 1.0 - x)[1] ≈ val # scipy.special.betainc(100.5,100.5,x) + @testset "a=100.5, b=100.5, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-4966,3.356014382599029e-143,2.21351837164392e-113,4.014552683228967e-96, + 5.21091312517625e-84,1.020389199029434e-74,3.267737964647878e-67,6.101573725809464e-61, + 1.418580292320946e-55,6.698584764694289e-51,8.971972641791594e-47,4.325143541605987e-43, + 8.94664782857427e-40,9.074256619313614e-37,5.005701167291761e-34,1.630336905229558e-31, + 3.349388657777907e-29,4.581338772370011e-27,4.362832768445309e-25,3.002816845923933e-23, + 1.541642496405598e-21,6.064729213175483e-20,1.870837362242701e-18,4.616584448897111e-17, + 9.27255359484341e-16,1.539075197361062e-14,2.139403432844863e-13,2.520059271944931e-12, + 2.541836676624415e-11,2.215842537778625e-10,1.683408890736215e-9,1.122874450435998e-8, + 6.620173165982372e-8,3.470763841511403e-7,1.626934118513085e-6,6.852626959254337e-6, + 0.00002605241980849797,0.00008977130313754974,0.0002814332703651759,0.0008055416902949821, + 0.002112019781014849,0.00508797075711991,0.01129547597461094,0.0231746884205069, + 0.04406489065169599,0.07787036008215303,0.1282703308247337,0.1975636232169813,0.2854862534447546, + 0.3885022505229126,0.5,0.6114977494770874,0.7145137465552454,0.8024363767830187,0.8717296691752663, + 0.922129639917847,0.955935109348304,0.9768253115794931,0.9887045240253891,0.9949120292428801, + 0.9978879802189852,0.999194458309705,0.9997185667296348,0.9999102286968625,0.9999739475801915, + 0.9999931473730407,0.9999983730658815,0.9999996529236158,0.9999999337982683,0.9999999887712555, + 0.9999999983165911,0.999999999778416,0.999999999974582,0.99999999999748,0.999999999999786, + 0.999999999999985,0.999999999999999,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1., + 1.,1.,1.]) + @test beta_inc(100.5, 100.5, x, 1.0 - x)[1] ≈ val rtol=1e-13 end - @test beta_inc(1.5,200.5,0.07,0.93)[1] ≈ 0.99999790408564 + @testset "a=1.5, b=1e-20, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [0.e-95,6.706954621511613e-24,1.908573860989784e-23,3.527823559800854e-23, + 5.465108108164382e-23,7.685476511624814e-23,1.016654428420256e-22,1.289283267454361e-22, + 1.585315312526521e-22,1.903920840622343e-22,2.244476844084102e-22,2.606514996691386e-22, + 2.989687056036453e-22,3.393740920821702e-22,3.818503532521727e-22,4.26386835711828e-22, + 4.729786038720361e-22,5.21625731553983e-22,5.72332759176175e-22,6.251082749838933e-22, + 6.799645911929102e-22,7.369174942192775e-22,7.95986053851056e-22,8.571924801968895e-22, + 9.205620200883313e-22,9.861228866810969e-22,1.053906217532698e-21,1.12394605758909e-21, + 1.196279364399263e-21,1.270946033569195e-21,1.347988943018207e-21,1.427454015049574e-21, + 1.50939029562088e-21,1.593850050518841e-21,1.680888878423608e-21,1.77056584110116e-21, + 1.862943611198906e-21,1.958088638346952e-21,2.056071334492481e-21,2.156966279623312e-21, + 2.260852449274566e-21,2.367813465464866e-21,2.477937872981194e-21,2.591319443230187e-21, + 2.708057508204485e-21,2.828257327482878e-21,2.952030491600198e-21,3.079495365596476e-21, + 3.210777577095268e-21,3.346010553881064e-21,3.48533611665991e-21,3.628905133513816e-21, + 3.776878243519303e-21,3.929426658119426e-21,4.086733050147924e-21,4.248992541941733e-21, + 4.416413805790527e-21,4.589220292116253e-21,4.767651603322342e-21,4.951965034288911e-21, + 5.142437304125938e-21,5.339366508167657e-21,5.543074324470629e-21,5.753908515483088e-21, + 5.972245773362194e-21,6.19849496698491e-21,6.433100860486935e-21,6.67654838776915e-21, + 6.929367585613005e-21,7.192139310873532e-21,7.465501896037407e-21,7.75015893407293e-21, + 8.046888430431106e-21,8.356553620653897e-21,8.680115830952776e-21,9.018649862807561e-21, + 9.373362520215628e-21,9.745615082692113e-21,1.013695077726479e-20,1.05491286460553e-20, + 1.098416568357789e-20,1.14443897916644e-20,1.193250706559674e-20,1.245168833304345e-20, + 1.300568195962486e-20,1.359896310974523e-20,1.423693458276194e-20,1.492620220780857e-20, + 1.567496069508049e-20,1.649354777878116e-20,1.739526322363106e-20,1.839762082220938e-20, + 1.952434090962006e-20,2.080868037972059e-20,2.229934708861127e-20,2.407185551639959e-20, + 2.625271544895813e-20,2.907909078985462e-20,3.308342534794925e-20,3.996470818039522e-20,1.]) + @test beta_inc(1.5, 1e-20, x, 1.0 - x)[1] ≈ val rtol=1e-14 + end + @testset "a=1e-20, b=0.99, x=$x" for (x, val) in zip( + 0.0:0.01:1.0, + [1.0,4.621640571507914e-20,3.928392625456373e-20,3.522826232996826e-20,3.235042350168417e-20, + 3.011796455118398e-20,2.82937201372135e-20,2.675117900735502e-20,2.541482518522368e-20, + 2.423594928780016e-20,2.318129286271936e-20,2.222713398379447e-20,2.135595723375987e-20, + 2.05544611885328e-20,1.981230641875294e-20,1.912129648210598e-20,1.847482377915803e-20, + 1.786748370079641e-20,1.729479923207514e-20,1.675302011465343e-20,1.623897358450671e-20, + 1.574995156483639e-20,1.528362412553572e-20,1.483797219538297e-20,1.441123460536074e-20, + 1.400186594951488e-20,1.360850271543896e-20,1.322993581024732e-20,1.286508808543862e-20, + 1.251299580739944e-20,1.217279327043997e-20,1.184369993372369e-20,1.152500960103253e-20, + 1.121608126599496e-20,1.091633132431086e-20,1.062522691510309e-20,1.034228020045673e-20, + 1.006704342884576e-20,9.799104656964201e-21,9.538084027304864e-21,9.283630517029258e-21, + 9.035419088275798e-21,8.793148181840232e-21,8.556537505729393e-21,8.325326077896178e-21, + 8.099270488866332e-21,7.878143355244009e-21,7.661731939451283e-21,7.449836914688898e-21, + 7.242271257138218e-21,7.03885924996661e-21,6.839435585837694e-21,6.643844556434574e-21, + 6.451939319035191e-21,6.263581231480807e-21,6.078639247988984e-21,5.896989369212464e-21, + 5.718514140760424e-21,5.543102195099824e-21,5.370647832359387e-21,5.201050636081688e-21, + 5.034215120421851e-21,4.870050405684712e-21,4.708469919434276e-21,4.549391120707029e-21, + 4.392735245120113e-21,4.238427068891405e-21,4.086394689985318e-21,3.936569324769163e-21, + 3.788885118712737e-21,3.643278969790933e-21,3.499690363356917e-21,3.358061217343224e-21, + 3.21833573672026e-21,3.080460276196264e-21,2.944383210178798e-21,2.810054809033594e-21, + 2.677427120668968e-21,2.546453856438284e-21,2.41709028028191e-21,2.289293099913057e-21, + 2.163020358672793e-21,2.038231326414587e-21,1.914886387391528e-21,1.79294692255393e-21, + 1.67237518283336e-21,1.553134148749398e-21,1.435187369793683e-21,1.318498774123754e-21, + 1.20303243443583e-21,1.088752268197901e-21,9.756216372555088e-22,8.63602788212576e-22, + 7.526560302762232e-22,6.427384567004279e-22,5.338018164784628e-22,4.25788652801385e-22, + 3.186244273953916e-22,2.121983888512568e-22,1.063002842170579e-22,0]) + @test beta_inc(1e-20, 0.99, x, 1.0 - x)[2] ≈ val rtol=1e-14 + end + @test beta_inc(1.5, 200.5, 0.07,0.93)[1] ≈ 0.99999790408564 + @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89))) @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1) end From 8863a5ffd2fd97eac6990bf6ef16d49aab26efb2 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 May 2021 23:05:38 +0200 Subject: [PATCH 08/11] Fix a bug in beta_inc --- src/beta_inc.jl | 36 +++++++++++++++++++----------------- test/beta_inc.jl | 2 ++ 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index fbfe690c..ecbacda8 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -105,10 +105,10 @@ function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Floa y0 = h/(1.0 + h) lambda = (a+b)*y - b else - h = a/b - x0 = h/(1.0 + h) - y0 = 1.0/(1.0 + h) - lambda = a - (a+b)*x + h = a/b + x0 = h/(1.0 + h) + y0 = 1.0/(1.0 + h) + lambda = a - (a+b)*x end e = -lambda/a u = abs(e) > 0.6 ? e - log(x/x0) : - LogExpFunctions.log1pmx(e) @@ -155,23 +155,25 @@ function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Floa t = 1.0 + rgamma1pm1(apb) end return a0*(esum(mu,z))*(1.0 + rgamma1pm1(b0))/t + else + ans = esum(mu, z) + if ans == 0.0 + return 0.0 + end + apb = a + b + if apb > 1.0 + z = (1.0 + rgamma1pm1(apb - 1.0))/apb + else + z = 1.0 + rgamma1pm1(apb) + end + c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z + return ans*(a0*c)/(1.0 + a0/b0) end else - z -= logbeta(a,b) - ans = esum(mu,z) + z -= logbeta(a, b) + ans = esum(mu, z) return ans end - if ans == 0.0 - return 0.0 - end - apb = a + b - if apb > 1.0 - z = (1.0 + rgamma1pm1(apb - 1.0))/apb - else - z = 1.0 + rgamma1pm1(apb) - end - c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z - return ans*(a0*c)/(1.0 + a0/b0) end """ diff --git a/test/beta_inc.jl b/test/beta_inc.jl index ab51c6dd..402212ce 100644 --- a/test/beta_inc.jl +++ b/test/beta_inc.jl @@ -206,7 +206,9 @@ 3.186244273953916e-22,2.121983888512568e-22,1.063002842170579e-22,0]) @test beta_inc(1e-20, 0.99, x, 1.0 - x)[2] ≈ val rtol=1e-14 end + @test beta_inc(1.5, 200.5, 0.07,0.93)[1] ≈ 0.99999790408564 + @test beta_inc(1e-20, 0.000001, 0.2)[2] ≈ 1.0000013862929421e-14 @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89))) @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1) From e41adcfbfaa9b86dceb2acd8b9e5cda716c5f2a2 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 18 May 2021 17:56:09 +0200 Subject: [PATCH 09/11] Avoid an infinite recursion in beta_inc_inv by ensuring that search direction is finite. --- src/beta_inc.jl | 28 ++++++++------- test/beta_inc.jl | 91 ++++++++++++++++++++++++------------------------ 2 files changed, 61 insertions(+), 58 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index ecbacda8..7d8c336d 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -919,28 +919,28 @@ function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbe #Initial approx - r = sqrt(-log(pp^2)) + r = sqrt(-2*log(pp)) pp_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00) if a > 1.0 && b > 1.0 r = (pp_approx^2 - 3.0)/6.0 s = 1.0/(2*aa - 1.0) t = 1.0/(2*bb - 1.0) - h = 2.0/(s+t) - w = pp_approx*sqrt(h+r)/h - (t-s)*(r + 5.0/6.0 - 2.0/(3.0*h)) - x = aa/ (aa+bb*exp(w^2)) + h = 2.0/(s + t) + w = pp_approx*sqrt(h + r)/h - (t - s)*(r + 5.0/6.0 - 2.0/(3.0*h)) + x = aa/(aa + bb*exp(w + w)) else r = 2.0*bb t = 1.0/(9.0*bb) - t = r*(1.0-t+pp_approx*sqrt(t))^3 + t = r*(1.0 - t + pp_approx*sqrt(t))^3 if t <= 0.0 - x = -expm1((log((1.0-pp)*bb)+lb)/bb) + x = -expm1((log((1.0 - pp)*bb) + lb)/bb) else - t = (4.0*aa+r-2.0)/t + t = (4.0*aa + r - 2.0)/t if t <= 1.0 - x = exp((log(pp*aa)+lb)/aa) + x = exp((log(pp*aa) + lb)/aa) else - x = 1.0 - 2.0/(t+1.0) + x = 1.0 - 2.0/(t + 1.0) end end end @@ -965,9 +965,10 @@ function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbe #iterate while true - pp_approx = beta_inc(aa,bb,x)[1] + pp_approx = beta_inc(aa, bb, x)[1] xin = x - pp_approx = (pp_approx-pp)*exp(lb+r*log(xin)+t*log1p(-xin)) + pp_approx = (pp_approx - pp)*min(floatmax(), exp(lb + r*log(xin) + t*log1p(-xin))) + if pp_approx * pp_approx_prev <= 0.0 prev = max(sq, fpu) end @@ -978,6 +979,7 @@ function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbe adj = g*pp_approx sq = adj^2 tx = x - adj + prev if (prev > sq && tx >= 0.0 && tx <= 1.0) break end @@ -988,11 +990,11 @@ function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbe if prev <= acu || pp_approx^2 <= acu x = tx - return indx ? (1.0 - x, x) : (x, 1.0-x) + return indx ? (1.0 - x, x) : (x, 1.0 - x) end if tx == x - return indx ? (1.0 - x, x) : (x, 1.0-x) + return indx ? (1.0 - x, x) : (x, 1.0 - x) end x = tx diff --git a/test/beta_inc.jl b/test/beta_inc.jl index 402212ce..3badd737 100644 --- a/test/beta_inc.jl +++ b/test/beta_inc.jl @@ -216,49 +216,50 @@ 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 + @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(1000.0, 2.0, 9.0797754e-317) ≈ 0.48 # This one is a bit slow (but also hard) end From cc8c06d7e53224e7f1aec0c1f34a649c4a7cb0d0 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 19 May 2021 13:03:16 +0200 Subject: [PATCH 10/11] Use float literals in comparisons in some methods restricted to Float64 arguments. --- src/beta_inc.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 7d8c336d..a16c48a3 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -52,7 +52,7 @@ for a0, b0 >= 8 function stirling_corr(a0::Float64, b0::Float64) a = min(a0, b0) b = max(a0, b0) - @assert a >= 8 + @assert a >= 8.0 h = a/b c = h/(1.0 + h) @@ -190,8 +190,8 @@ See also: [`beta_inc`](@ref) `BFRAC(A,B,X,Y,LAMBDA,EPS)` from Didonato and Morris (1982) """ function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, lambda::Float64, epps::Float64) - @assert a > 1 - @assert b > 1 + @assert a > 1.0 + @assert b > 1.0 ans = beta_integrand(a,b,x,y) if ans == 0.0 return 0.0 @@ -258,8 +258,8 @@ See also: [`beta_inc`](@ref) `BASYM(A,B,LAMBDA,EPS)` from Didonato and Morris (1982) """ function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64) - @assert a >= 15 - @assert b >= 15 + @assert a >= 15.0 + @assert b >= 15.0 a0 =zeros(22) b0 = zeros(22) c = zeros(22) @@ -356,8 +356,8 @@ External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wi See also: [`beta_inc`](@ref) """ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64) - @assert a >= 15 - @assert b <= 1 + @assert a >= 15.0 + @assert b <= 1.0 c = zeros(31) d = zeros(31) bm1 = b - 1.0 @@ -521,7 +521,7 @@ See also: [`beta_inc`](@ref) `BPSER(A,B,X,EPS)` from Didonato and Morris (1982) """ function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64) - @assert b <= 1 || b*x <= 0.7 + @assert b <= 1.0 || b*x <= 0.7 ans = 0.0 if x == 0.0 return 0.0 @@ -691,7 +691,7 @@ function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Intege end kp1 = k + 1 - for i = kp1:nm1 + for i in kp1:nm1 l = i - 1 d *= ((apb + l)/(ap1 + l))*x w += d From d87a95ed477c61dca54501335acb75882e084f06 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 19 May 2021 20:30:40 +0200 Subject: [PATCH 11/11] Update src/beta_inc.jl Co-authored-by: David Widmann --- src/beta_inc.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index a16c48a3..b6b75cb3 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -979,7 +979,6 @@ function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbe adj = g*pp_approx sq = adj^2 tx = x - adj - prev if (prev > sq && tx >= 0.0 && tx <= 1.0) break end