Fix a bug in beta_inc

Author: andreasnoack

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
 
 """

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)