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