From 0424b0420a42383198927e5cf8ae8def9719197f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Felix=20N=C3=B6=C3=9Fler?= <felixnoessler@aol.com>
Date: Mon, 22 Apr 2024 18:43:45 +0200
Subject: [PATCH] try to fix numerical instability with high values of the
 degrees of freedom of the t-distribution

---
 src/distrs/fdist.jl | 12 ++++++++++--
 src/distrs/tdist.jl |  4 +++-
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl
index 17dfe17..4e6cd74 100644
--- a/src/distrs/fdist.jl
+++ b/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
diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl
index e2673af..f58edb2 100644
--- a/src/distrs/tdist.jl
+++ b/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