try to fix numerical instability with high values of the degrees of f…

…reedom of the t-distribution

src/distrs/fdist.jl

@@ -21,9 +21,17 @@ end
 for f in ("invcdf", "invccdf", "invlogcdf", "invlogccdf")
     ff = Symbol("fdist"*f)
     bf = Symbol("beta"*f)
+    cf = Symbol("chisq"*f)
     @eval function $ff(ν1::T, ν2::T, y::T) where {T<:Real}
-        x = $bf(ν1/2, ν2/2, y)
-        return x/(1 - x)*ν2/ν1
+        # the evaluation of the CDF with the t-distribution combined with the
+        # beta distribution is unstable for high ν2,
+        # therefore, we switch from the beta to the chi-squared distribution at ν > 1.0e7
+        if ν2 > 1.0e7
+            return $cf(ν1, y)/ν1
+        else
+            x = $bf(ν1/2, ν2/2, y)
+            return x/(1 - x)*ν2/ν1
+        end
     end
     @eval $ff(ν1::Real, ν2::Real, y::Real) = $ff(promote(ν1, ν2, y)...)
 end

src/distrs/tdist.jl

@@ -4,7 +4,9 @@ tdistpdf(ν::Real, x::Real) = exp(tdistlogpdf(ν, x))
 
 tdistlogpdf(ν::Real, x::Real) = tdistlogpdf(promote(ν, x)...)
 function tdistlogpdf(ν::T, x::T) where {T<:Real}
-    isinf(ν) && return normlogpdf(x)
+    # the logpdf equation of the t-distribution is numerically unstable for high ν,
+    # therefore we switch at ν > 1.0e7 to the normal distribution
+    ν > 1.0e7 && return normlogpdf(x)
     νp12 = (ν + 1) / 2
     return loggamma(νp12) - (logπ + log(ν)) / 2 - loggamma(ν / 2) - νp12 * log1p(x^2 / ν)
 end