diff --git a/docs/src/index.md b/docs/src/index.md index 734c7b45..b5490b3f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -20,6 +20,7 @@ libraries. | [`trigamma(x)`](@ref SpecialFunctions.trigamma) | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) | | [`polygamma(m,x)`](@ref SpecialFunctions.polygamma) | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `lgamma` function at `x`) | | [`gamma_inc(a,x,IND)`](@ref SpecialFunctions.gamma_inc) | [incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates P(a,x) and Q(a,x)for accuracy specified by IND and returns tuple (p,q)) | +| [`beta_inc(a,b,x,y)`](@ref SpecialFunctions.beta_inc) | [incomplete beta function ratio Ix(a,b) and Iy(a,b)](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q)) | | [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv) | [inverse of incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates x given P(a,x)=p and Q(a,x)=q | | [`ellipk(x)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 1st kind at `x`) | | [`ellipe(x)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 2nd kind at `x`) | diff --git a/docs/src/special.md b/docs/src/special.md index 3a236591..cac0587c 100644 --- a/docs/src/special.md +++ b/docs/src/special.md @@ -51,6 +51,7 @@ SpecialFunctions.zeta SpecialFunctions.gamma SpecialFunctions.gamma_inc SpecialFunctions.gamma_inc_inv +SpecialFunctions.beta_inc SpecialFunctions.loggamma SpecialFunctions.logabsgamma SpecialFunctions.logfactorial diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 2f8aff05..db2deb12 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -48,6 +48,7 @@ export polygamma, trigamma, gamma_inc, + beta_inc, gamma_inc_inv, hankelh1, hankelh1x, @@ -64,6 +65,7 @@ include("ellip.jl") include("sincosint.jl") include("gamma.jl") include("gamma_inc.jl") +include("beta_inc.jl") include("deprecated.jl") for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, diff --git a/src/beta_inc.jl b/src/beta_inc.jl new file mode 100644 index 00000000..0a2bceb2 --- /dev/null +++ b/src/beta_inc.jl @@ -0,0 +1,827 @@ +using Base.Math: @horner +using Base.MPFR: ROUNDING_MODE +const exparg_n = log(nextfloat(floatmin(Float64))) +const exparg_p = log(prevfloat(floatmax(Float64))) + +#COMPUTE log(gamma(b)/gamma(a+b)) when b >= 8 +""" + loggammadiv(a,b) + +Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8 +""" +function loggammadiv(a::Float64, b::Float64) + if a > b + h = b/a + c = 1.0/(1.0 + h) + x = h/(1.0 + h) + d = a + (b - 0.5) + else + h = a/b + c = h/(1.0 + h) + x = 1.0/(1.0 + h) + d = b + a - 0.5 + end + x² = x*x + s₃ = 1.0 + (x + x²) + s₅ = 1.0 + (x + x²*s₃) + s₇ = 1.0 + (x + x²*s₅) + s₉ = 1.0 + (x + x²*s₇) + s₁₁ = 1.0 + (x + x²*s₉) + + # SET W = stirling(b) - stirling(a+b) + t = inv(b)^2 + w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁) + w *= c/b + + #COMBINING + u = d*log1p(a/b) + v = a*(log(b) - 1.0) + return u <= v ? w - u - v : w - v - u +end + +""" + stirling_corr(a0,b0) + +Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) +for a0,b0 >= 8 +""" +function stirling_corr(a0::Float64, b0::Float64) + a = min(a0,b0) + b = max(a0,b0) + + h = a/b + c = h/(1.0 + h) + x = 1.0/(1.0 + h) + x² = x*x + #SET SN = (1-X^N)/(1-X) + s₃ = 1.0 + (x + x²) + s₅ = 1.0 + (x + x²*s₃) + s₇ = 1.0 + (x + x²*s₅) + s₉ = 1.0 + (x + x²*s₇) + s₁₁ = 1.0 + (x + x²*s₉) + t = inv(b)^2 + w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁) + w *= c/b + # COMPUTE stirling(a) + w + t = inv(a)^2 + return @horner(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w +end + +""" + esum(mu,x) + +Compute ``e^{mu+x}`` +""" +function esum(mu::Float64, x::Float64) + if x > 0.0 + if mu > 0.0 || mu + x < 0.0 + return exp(mu)*exp(x) + else + return exp(mu + x) + end + elseif mu < 0.0 || mu + x > 0.0 + return exp(mu)*exp(x) + else + return exp(mu + x) + end +end + +""" + 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) + a0, b0 = minmax(a,b) + if a0 >= 8.0 + if a > b + h = b/a + x0 = 1.0/(1.0 + h) + 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 + end + e = -lambda/a + u = abs(e) > 0.6 ? u = e - log(x/x0) : - log1pmx(e) + e = lambda/b + v = abs(e) > 0.6 ? e - log(y/y0) : - log1pmx(e) + z = esum(mu, -(a*u + b*v)) + return (1.0/sqrt(2*pi))*sqrt(b*x0)*z*exp(-stirling_corr(a,b)) + elseif x > 0.375 + if y > 0.375 + lnx = log(x) + lny = log(y) + else + lnx = log1p(-y) + lny = log(y) + end + else + lnx = log(x) + lny = log1p(-x) + end + z = a*lnx + b*lny + if a0 < 1.0 + b0 = max(a,b) + if b0 >= 8.0 + u = loggamma1p(a0) + loggammadiv(a0,b0) + return a0*(esum(mu, z-u)) + elseif b0 > 1.0 + u = loggamma1p(a0) + n = trunc(Int,b0 - 1.0) + if n >= 1 + c = 1.0 + for i = 1:n + b0 -= 1.0 + c *= (b0/(a0+b0)) + end + u += log(c) + end + z -= u + b0 -= 1.0 + apb = a0 + b0 + if apb > 1.0 + u = a0 + b0 - 1.0 + t = (1.0 + rgamma1pm1(u))/apb + else + t = 1.0 + rgamma1pm1(apb) + end + return a0*(esum(mu,z))*(1.0 + rgamma1pm1(b0))/t + end + else + 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 + +""" + beta_inc_cont_fraction(a,b,x,y,lambda,epps) + +Compute ``I_{x}(a,b)`` using continued fraction expansion when a,b > 1. +It is assumed that ``\\lambda = (a+b)*y - b`` +DLMF : https://dlmf.nist.gov/8.17#E22 +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) + ans = beta_integrand(a,b,x,y) + if ans == 0.0 + return 0.0 + end + c = 1.0 + lambda + c0 = b/a + c1 = 1.0 + 1.0/a + yp1 = y + 1.0 + + n = 0.0 + p = 1.0 + s = a + 1.0 + an = 0.0 + bn = 1.0 + anp1 = 1.0 + bnp1 = c/c1 + r = c1/c + #CONT FRACTION + + while true + n += 1.0 + t = n/a + w = n*(b - n)*x + e = a/s + alpha = (p*(p+c0)*e*e)*(w*x) + e = (1.0 + t)/(c1 + 2*t) + beta = n + w/s +e*(c + n*yp1) + p = 1.0 + t + s += 2.0 + + #update an, bn, anp1, bnp1 + t = alpha*an + beta*anp1 + an = anp1 + anp1 = t + t = alpha*bn + beta*bnp1 + bn = bnp1 + bnp1 = t + + r0 = r + r = anp1/bnp1 + if abs(r - r0) <= epps*r + break + end + #rescale + an /= bnp1 + bn /= bnp1 + anp1 = r + bnp1 = 1.0 + end + return ans*r +end + +""" + beta_inc_asymptotic_symmetric(a,b,lambda,epps) + +Compute ``I_{x}(a,b)`` using asymptotic expansion for a,b >= 15. +It is assumed that ``\\lambda = (a+b)*y - b`` +BASYM(A,B,LAMBDA,EPS) from Didonato and Morris (1982) +""" +function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64) + a0 =zeros(22) + b0 = zeros(22) + c = zeros(22) + d = zeros(22) + e0 = 2/sqrt(pi) + e1 = 2^(-1.5) + sm = 0.0 + ans = 0.0 + if a > b + h = b/a + r0 = 1.0/(1.0 + h) + r1 = (b-a)/a + w0 = 1.0/sqrt(b*(1.0+h)) + else + h = a/b + r0 = 1.0/(1.0 + h) + r1 = (b-a)/b + w0 = 1.0/sqrt(a*(1.0+h)) + end + f = -a*log1pmx(-(lambda/a)) - b*log1pmx((lambda/b)) + t = exp(-f) + if t == 0.0 + return ans + end + z0 = sqrt(f) + z = 0.5*(z0/e1) + z² = 2.0*f + + a0[1] = (2.0/3.0)*r1 + c[1] = -0.5*a0[1] + d[1] = - c[1] + j0 = (0.5/e0)*erfcx(z0) + j1 = e1 + sm = j0 + d[1]*w0*j1 + + s = 1.0 + h² = h*h + hn = 1.0 + w = w0 + znm1 = z + zn = z² + + for n = 2: 2: 20 + hn *= h² + a0[n] = 2.0*r0*(1.0 + h*hn)/(n + 2.0) + s += hn + a0[n+1] = 2.0*r1*s/(n+3.0) + + for i = n: n+1 + r = -0.5*(i + 1.0) + b0[1] = r*a0[1] + for m = 2:i + bsum = 0.0 + for j =1: m-1 + bsum += (j*r - (m-j))*a0[j]*b0[m-j] + end + b0[m] = r*a0[m] + bsum/m + end + c[i] = b0[i]/(i+1.0) + dsum = 0.0 + for j = 1: i-1 + imj = i - j + dsum += d[imj]*c[j] + end + d[i] = -(dsum + c[i]) + end + + j0 = e1*znm1 + (n - 1)*j0 + j1 = e1*zn + n*j1 + znm1 *= z² + zn *= z² + w *= w0 + t0 = d[n]*w*j0 + w *= w0 + t1 = d[n+1]*w*j1 + sm += (t0 + t1) + if (abs(t0) + abs(t1)) <= epps*sm + break + end + end + + u = exp(-stirling_corr(a,b)) + return e0*t*u*sm +end + +""" + 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. +It is assumed a >= 15 and b <= 1, and epps is tolerance used. +""" +function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64) + c = zeros(31) + d = zeros(31) + bm1 = b - 1.0 + nu = a + 0.5*bm1 + if y > 0.375 + lnx = log(x) + else + lnx = log1p(-y) + end + z = -nu*lnx + if b*z == 0.0 + return error("expansion can't be computed") + end + + # COMPUTATION OF THE EXPANSION + #SET R = EXP(-Z)*Z**B/GAMMA(B) + r = b*(1.0 + rgamma1pm1(b))*exp(b*log(z)) + r *= exp(a*lnx)*exp(0.5*bm1*lnx) + u = loggammadiv(b,a) + b*log(nu) + u = r*exp(-u) + if u == 0.0 + return error("expansion can't be computed") + end + (p, q) = gamma_inc(b,z,0) + v = inv(nu)^2/4 + t2 = lnx^2/4 + l = w/u + j = q/r + sm = j + t = 1.0 + cn = 1.0 + n2 = 0.0 + for n = 1:30 + bp2n = b + n2 + j = (bp2n*(bp2n + 1.0)*j + (z + bp2n + 1.0)*t)*v + n2 += 2.0 + t *= t2 + cn /= n2*(n2 + 1.0) + c[n] = cn + s = 0.0 + if n != 1 + nm1 = n -1 + coef = b - n + for i = 1:nm1 + s += coef*c[i]*d[n-i] + coef += b + end + end + d[n] = bm1*cn + s/n + dj = d[n] * j + sm += dj + if sm <= 0.0 + return error("expansion can't be computed") + end + if abs(dj) <= epps*(sm+l) + break + end + end + return w + u*sm +end + + +#For b < min(eps, eps*a) and x <= 0.5 +""" + beta_inc_power_series2(a,b,x,epps) + +Variant of BPSER(A,B,X,EPS). +FPSER(A,B,X,EPS) from Didonato and Morris (1982) +""" +function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64) + ans = 1.0 + if a > 1.0e-3 * epps + ans = 0.0 + t = a*log(x) + if t < exparg_n + return ans + end + ans = exp(t) + end + ans *= b/a + tol = epps/a + an = a + 1.0 + t = x + s = t/an + an += 1.0 + t *= x + c = t/an + s += c + while abs(c) > tol + an += 1.0 + t *= x + c = t/an + s += c + end + ans *= (1.0 + a*s) + return ans +end + +#A <= MIN(EPS,EPS*B), B*X <= 1, AND X <= 0.5., A is small +""" + beta_inc_power_series1(a,b,x,epps) + +Another variant of BPSER(A,B,X,EPS) +APSER(A,B,X,EPS) from Didonato and Morris (1982) +""" +function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64) + g = Base.MathConstants.eulergamma + bx = b*x + t = x - bx + if b*epps > 2e-2 + c = log(bx) + g + t + else + c = log(x) + digamma(b) + g + t + end + tol = 5.0*epps*abs(c) + j = 1.0 + s = 0.0 + j += 1.0 + t *= (x - bx/j) + aj = t/j + s += aj + while abs(aj) > tol + j += 1.0 + t *= (x - bx/j) + aj = t/j + s += aj + end + return -a*(c + s) +end + +#B .LE. 1 OR B*X .LE. 0.7 +""" + beta_inc_power_series(a,b,x,epps) + +Computes Ix(a,b) using power series : +```math +I_{x}(a,b) = G(a,b)x^{a}/a (1 + a\\sum_{1}^{\\infty}((1-b)(2-b)...(j-b)/j!(a+j)) x^{j}) +``` +BPSER(A,B,X,EPS) from Didonato and Morris (1982) +""" +function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64) + ans = 0.0 + if x == 0.0 + return 0.0 + end + a0 = min(a,b) + b0 = max(a,b) + if a0 >= 1.0 + z = a*log(x) - logbeta(a,b) + ans = exp(z)/a + else + + if b0 >= 8.0 + u = loggamma1p(a0) + loggammadiv(a0,b0) + z = a*log(x) - u + ans = (a0/a)*exp(z) + if ans == 0.0 || a <= 0.1*epps + return ans + end + elseif b0 > 1.0 + u = loggamma1p(a0) + m = b0 - 1.0 + if m >= 1.0 + c = 1.0 + for i = 1:m + b0 -= 1.0 + c *= (b0/(a0+b0)) + end + u += log(c) + end + z = a*log(x) - u + b0 -= 1.0 + apb = a0 + b0 + if apb > 1.0 + u = a0 + b0 - 1.0 + t = (1.0 + rgamma1pm1(u))/apb + else + t = 1.0 + rgamma1pm1(apb) + end + ans = exp(z)*(a0/a)*(1.0 + rgamma1pm1(b0))/t + if ans == 0.0 || a <= 0.1*epps + return ans + end + else + #PROCEDURE FOR A0 < 1 && B0 < 1 + ans = x^a + if ans == 0.0 + return ans + end + apb = a + b + if apb > 1.0 + u = a + b - 1.0 + z = (1.0 + rgamma1pm1(u))/apb + else + z = 1.0 + rgamma1pm1(apb) + end + c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z + ans *= c*(b/apb) + #label l70 start + if ans == 0.0 || a <= 0.1*epps + return ans + end + end + end + if ans == 0.0 || a <= 0.1*epps + return ans + end + # COMPUTE THE SERIES + + sm = 0.0 + n = 0.0 + c = 1.0 + tol = epps/a + n += 1.0 + c *= x*(1.0 - b/n) + w = c/(a + n) + sm += w + while abs(w) > tol + n += 1.0 + c *= x*(1.0 - b/n) + w = c/(a+n) + sm += w + end + return ans*(1.0 + a*sm) +end + +""" + 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 more generalised version of https://dlmf.nist.gov/8.17#E20 +""" +function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Integer, epps::Float64) + apb = a + b + ap1 = a + 1.0 + mu = 0.0 + d = 1.0 + if n != 1 && a >= 1.0 && apb >= 1.1*ap1 + mu = abs(exparg_n) + k = exparg_p + if k < mu + mu = k + end + t = mu + d = exp(-t) + end + + ans = beta_integrand(a,b,x,y,mu)/a + if n == 1 || ans == 0.0 + return ans + end + nm1 = n -1 + w = d + + k = 0 + if b <= 1.0 + kp1 = k + 1 + for i = kp1:nm1 + l = i - 1 + d *= ((apb + l)/(ap1 + l))*x + w += d + if d <= epps*w + break + end + end + return ans*w + elseif y > 1.0e-4 + r = trunc(Int,(b - 1.0)*x/y - a) + if r < 1.0 + kp1 = k + 1 + for i = kp1:nm1 + l = i - 1 + d *= ((apb + l)/(ap1 + l))*x + w += d + if d <= epps*w + break + end + end + return ans*w + end + k = t = nm1 + if r < t + k = r + end + # ADD INC TERMS OF SERIES + for i = 1:k + l = i -1 + d *= ((apb + l)/(ap1 + l))*x + w += d + end + if k == nm1 + return ans*w + end + else + k = nm1 + for i = 1:k + l = i -1 + d *= ((apb + l)/(ap1 + l))*x + w += d + end + if k == nm1 + return ans*w + end + end + + kp1 = k + 1 + for i = kp1:nm1 + l = i - 1 + d *= ((apb + l)/(ap1 + l))*x + w += d + if d <= epps*w + break + end + end + return ans*w +end + + +#SIGNIFICANT DIGIT COMPUTATION OF INCOMPLETE BETA FUNCTION RATIOS +#by ARMIDO R. DIDONATO AND ALFRED H. MORRIS, JR. +#ACM Transactions on Mathematical Software. Vol18, No3, September1992, Pages360-373 +#DLMF : https://dlmf.nist.gov/8.17#E1 +#Wikipedia : https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function + +""" + beta_inc(a,b,x,y) + +Computes Incomplete Beta Function Ratios given by: +```math +I_{x}(a,b) = G(a,b) \\int_{0}^{x} t^{a-1}(1-t)^{b-1} dt, +``` +and ``I_{y}(a,b) = 1.0 - I_{x}(a,b)``. +given +``B(a,b) = 1/G(a,b) = \\Gamma(a)\\Gamma(b)/\\Gamma(a+b)`` and ``x+y = 1``. +""" +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") + elseif a == 0.0 && b == 0.0 + return error("a and b are 0.0") + elseif x < 0.0 || x > 1.0 + return error("x < 0 or x > 1") + elseif y < 0.0 || y > 1.0 + return error("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 + end + end + + if x == 0.0 + return (0.0,1.0) + elseif y == 0.0 + return (1.0,0.0) + elseif a == 0.0 + return (1.0,0.0) + elseif b == 0.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)) + end + ind = false + a0 = a + b0 = b + x0 = x + y0 = y + + if min(a0,b0) > 1.0 + #PROCEDURE FOR A0>1 AND B0>1 + lambda = a > b ? (a+b)*y - b : a - (a+b)*x + if lambda < 0.0 + ind = true + a0 = b + b0 = a + x0 = y + y0 = x + lambda = abs(lambda) + end + if b0 < 40.0 && b0*x0 <= 0.7 + 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 + 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()) + q = 1.0 - p + else + 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()) + q = 1.0 - p + else + p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps()) + q = 1.0 - p + end + return ind ? (q, p) : (p, q) + end +#PROCEDURE FOR A0<=1 OR B0<=1 + if x > 0.5 + ind = true + a0 = b + b0 = a + y0 = x + x0 = y + end + + if b0 < min(epps, epps*a0) + 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) + p = 1.0 - q + elseif max(a0,b0) > 1.0 + if b0 <= 1.0 + 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) + 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()) + p = 1.0 - q + else + n = 20 + 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()) + p = 1.0 - q + end + elseif (x0*b0)^(a0) <= 0.7 + 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) + b0 += n + 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) + q = 1.0 - p + elseif x0^a0 <= 0.9 + 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) + p = 1.0 - q + else + n = 20 + 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()) + p = 1.0 - q + end + +#TERMINATION + return ind ? (q, p) : (p, q) +end + +beta_inc(a::Float64, b::Float64, x::Float64) = beta_inc(a, b, x, 1.0 - x) +function beta_inc(a::T, b::T, x::T, y::T) where {T<:Union{Float16, Float32}} + T.(beta_inc(Float64(a), Float64(b), Float64(x), Float64(y))) +end +beta_inc(a::Real, b::Real, x::Real, y::Real) = beta_inc(promote(float(a), float(b), float(x), float(y))...) +beta_inc(a::T, b::T, x::T, y::T) where {T<:AbstractFloat} = throw(MethodError(beta_inc,(a, b, x, y,""))) diff --git a/test/runtests.jl b/test/runtests.jl index a8fcefce..b2ad1ac1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -124,6 +124,48 @@ end @test SpecialFunctions.stirling(x) ≈ log(gamma(x)) - (x-.5)*log(x)+x- log(2*pi)/2.0 end end +@testset "incomplete beta" begin +#Compared with scipy.special.betainc + ctr = 1 + 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] + for x = 0.01:0.01:0.99 + @test beta_inc(0.5,0.5,x,1.0-x)[1] ≈ ans1[ctr]#scipy.special.betainc(0.5,0.5,x) + ctr += 1 + end + ctr = 1 + for x = 0.01:0.01:0.99 + @test beta_inc(0.5,0.8,x,1.0-x)[1] ≈ ans2[ctr]#scipy.special.betainc(0.5,0.8,x) + ctr += 1 + end + ctr = 1 + for x = 0.01:0.01:0.99 + @test beta_inc(0.9,0.8,x,1.0-x)[1] ≈ ans3[ctr]#scipy.special.betainc(0.9,0.8,x) + ctr += 1 + end + ctr = 1 + for x = 0.01:0.01:0.99 + @test beta_inc(80.9,0.8,x,1.0-x)[1] ≈ ans4[ctr]#scipy.special.betainc(80.9,0.8,x) + ctr += 1 + end + ctr = 1 + for x = 0.01:0.01:0.99 + @test beta_inc(1.7,10.5,x,1.0-x)[1] ≈ ans5[ctr]#scipy.special.betainc(1.7,10.5,x) + ctr += 1 + end + ctr = 1 + for x = 0.01:0.01:0.99 + @test beta_inc(100.5,100.5,x,1.0-x)[1] ≈ ans6[ctr]#scipy.special.betainc(100.5,100.5,x) + ctr += 1 + 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(11.99) + SpecialFunctions.stirling(100.1) - SpecialFunctions.stirling(11.99 + 100.1) +end @testset "elliptic integrals" begin #Computed using Wolframalpha EllipticK and EllipticE functions. @test ellipk(0) ≈ 1.570796326794896619231322 rtol=2*eps()