Avoid an infinite loop in beta_inc for NaN input (#445)

This follows suit with `gamma_inc`, which returns `NaN` if either of its
two arguments are `NaN`, by returning `NaN` if any of the four arguments
are `NaN`.

src/beta_inc.jl

@@ -754,7 +754,9 @@ function _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64=1-x)
         end
     end
 
-    if x == 0.0
+    if isnan(x) || isnan(y) || isnan(a) || isnan(b)
+        return (NaN, NaN)
+    elseif x == 0.0
         return (0.0, 1.0)
     elseif y == 0.0
         return (1.0, 0.0)

test/beta_inc.jl

@@ -304,6 +304,8 @@ end
         # See https://github.com/JuliaStats/StatsFuns.jl/issues/133#issuecomment-1069602721
         y = 2.0e-280
         @test beta_inc(2.0, 1.0, beta_inc_inv(2.0, 1.0, y, 1.0)[1])[1] ≈ y
+        # See https://github.com/JuliaStats/GLM.jl/issues/538#issuecomment-1603447448
+        @test isequal(beta_inc(6.0, 112.5, NaN), (NaN, NaN))
     end
 
     @testset "StatsFuns#145" begin