Avoid an infinite recursion in beta_inc_inv by ensuring that search

direction is finite.

Author: andreasnoack

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

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