Added inverse of incomplete beta Ix(a,b) and Iy(a,b) (#166)

Author: Sumegh-git

src/SpecialFunctions.jl

@@ -49,6 +49,7 @@ export
     trigamma,
     gamma_inc,
     beta_inc,
+    beta_inc_inv,
     gamma_inc_inv,
     hankelh1,
     hankelh1x,

src/beta_inc.jl

@@ -825,3 +825,120 @@ function beta_inc(a::T, b::T, x::T, y::T) where {T<:Union{Float16, Float32}}
 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,"")))
+
+#GW Cran, KJ Martin, GE Thomas, Remark AS R19 and Algorithm AS 109: A Remark on Algorithms AS 63: The Incomplete Beta Integral and AS 64: Inverse of the Incomplete Beta Integeral,
+#Applied Statistics,
+#Volume 26, Number 1, 1977, pages 111-114. 
+
+"""
+    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).
+"""
+function beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64; lb = logbeta(a,b)[1])
+    fpu = 1e-30
+    x = p
+    if p == 0.0
+        return (0.0, 1.0)
+    elseif p == 1.0
+        return (1.0, 0.0)
+    end
+
+    #change tail if necessary
+
+    if p > 0.5
+        pp = q 
+        aa = b
+        bb = a
+        indx = true
+    else
+        pp = p
+        aa = a
+        bb = b
+        indx = false
+    end
+
+    #Initial approx
+
+    r = sqrt(-log(pp^2))
+    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))
+    else
+        r = 2.0*bb
+        t = 1.0/(9.0*bb)
+        t = r*(1.0-t+pp_approx*sqrt(t))^3
+        if t <= 0.0
+            x = -expm1((log((1.0-pp)*bb)+lb)/bb)
+        else
+            t = (4.0*aa+r-2.0)/t
+            if t <= 1.0
+                x = exp((log(pp*aa)+lb)/aa)
+            else
+                x = 1.0 - 2.0/(t+1.0)
+            end
+        end
+    end
+
+    #solve x using modified newton-raphson iteration
+
+    r = 1.0 - aa
+    t = 1.0 - bb
+    pp_approx_prev = 0.0
+    sq = 1.0
+    prev = 1.0
+
+    if x < 0.0001
+        x = 0.0001
+    end
+    if x > .9999
+        x = .9999
+    end
+
+    iex = max(-5.0/aa^2 - 1.0/pp^0.2 - 13.0, -30.0)
+    acu = 10.0^iex
+
+    #iterate
+    while true
+        pp_approx = beta_inc(aa,bb,x)[1]
+        xin = x
+        pp_approx = (pp_approx-pp)*exp(lb+r*log(xin)+t*log1p(-xin))
+        if pp_approx * pp_approx_prev <= 0.0
+            prev = max(sq, fpu)
+        end
+        g = 1.0
+
+        tx = x - g*pp_approx
+        while true
+            adj = g*pp_approx
+            sq = adj^2
+            tx = x - adj
+            if (prev > sq && tx >= 0.0 && tx <= 1.0)
+                break
+            end
+            g /= 3.0
+        end
+
+        #check if current estimate is acceptable
+    
+        if prev <= acu || pp_approx^2 <= acu
+            x = tx
+            return indx ? (1.0 - x, x) : (x, 1.0-x)
+        end
+
+        if tx == x
+            return indx ? (1.0 - x, x) : (x, 1.0-x)
+        end
+
+        x = tx
+        pp_approx_prev = pp_approx
+    end
+end
+
+beta_inc_inv(a::Float64, b::Float64, p::Float64) = beta_inc_inv(a, b, p, 1.0-p)

test/runtests.jl

@@ -166,6 +166,54 @@ end
     @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 "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
+end
 @testset "elliptic integrals" begin
 #Computed using Wolframalpha EllipticK and EllipticE functions.
 	@test ellipk(0) ≈ 1.570796326794896619231322 rtol=2*eps()