From 192904016f5addc267ce20f71a492838f32f73e9 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 17 Aug 2026 22:46:59 +0200 Subject: [PATCH 1/3] Implement tdistinvcdf natively instead of inverting the incomplete beta MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Based on VBA invtdist by Ian Smith. The quantile is started from a Cornish-Fisher expansion around the normal quantile in the central region, or from inverting the leading power-law term of the tail, and polished with Newton's method applied to log(cdf), where each iteration costs a single evaluation of the incomplete beta function; the pdf is recovered from the same argument reduction at negligible cost. For ν >= 250 with xn^2 <= ν / 100 the expansion is already fully converged and the polish is skipped. Compared to the previous implementation, which inverted the incomplete beta function numerically through fdistinvccdf, this is 2-30x faster and fixes deep-tail failures for small ν, e.g. tdistinvcdf(0.5, 1e-8) previously returned -Inf instead of about -1.03e15. Verified against a 512-bit reference: max relative error 3.8e-14 for ν in [1, 10000] over p in [1e-300, 1 - 1e-12] (mean about 1e-15). --- src/distrs/tdist.jl | 136 ++++++++++++++++++++++++++++++++++++++++++++ test/misc.jl | 68 ++++++++++++++++++++++ 2 files changed, 204 insertions(+) diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl index 8b1fef9..5e7b852 100644 --- a/src/distrs/tdist.jl +++ b/src/distrs/tdist.jl @@ -36,6 +36,142 @@ tdistlogcdf(ν::Real, x::Real) = tdistlogcdf(map(float, promote(ν, x))...) tdistlogccdf(ν::Real, x::Real) = tdistlogcdf(ν, -x) +# Pure Julia implementation of the inverse CDF, based on VBA invtdist by Ian Smith: +# start from a Cornish-Fisher expansion of the quantile around the normal quantile +# (central region) or from inverting the leading power-law term of the tail (tail +# region), then polish with Newton's method applied to log(cdf), where each +# iteration costs a single evaluation of the incomplete beta function. + +# Crossover curves between the Cornish-Fisher and tail starting approximations. +# Based on VBA BetterThanTailApprox by Ian Smith. +function _tdist_cf_start_is_better(pr::Float64, ν::Float64) + if ν <= 2.0 + return pr > 0.25 * exp((1.0 - ν) * 1.78514841051368) + elseif ν <= 5.0 + return pr > 0.045 * exp((2.0 - ν) * 1.30400766847605) + elseif ν <= 20.0 + return pr > 0.0009 * exp((5.0 - ν) * 0.921034037197618) + else + return pr > 9.0e-10 * exp((20.0 - ν) * 0.690775527898214) + end +end + +# cdf and pdf of the t distribution at x <= 0, with B = beta(ν / 2, 1 / 2) precomputed +function _tdistcdf_pdf(ν::Float64, B::Float64, x::Float64) + if abs(x) >= min(1.0, ν) + # this form of k2 = ν / (ν + x^2) and x2 = x^2 / (ν + x^2) avoids + # premature overflow of x^2 + k2 = ν / x + t = x + k2 + k2 = k2 / t + x2 = x / t + else + x² = x * x + t = ν + x² + x2 = x² / t + k2 = ν / t + end + p = first(beta_inc(ν / 2, 0.5, k2, x2)) / 2 + f = k2^((ν + 1) / 2) / (sqrt(ν) * B) + return p, f +end + +function _tdistinvcdf(ν::Float64, p::Float64) + if isnan(ν) || isnan(p) || !(0.0 <= p <= 1.0) || !(ν > 0.0) + return NaN + elseif isinf(ν) + return norminvcdf(p) + elseif p == 0.0 + return -Inf + elseif p == 1.0 + return Inf + elseif p == 0.5 + return 0.0 + end + + # work with the smaller tail; the result is negated for p > 1/2 on return + pr = p > 0.5 ? 1.0 - p : p + small = 1.0e-14 + smalllpr = -small * log(pr) * pr + + logB = logbeta(ν / 2, 0.5) + B = exp(logB) + + local tp::Float64 + tprob = 0.0 + tpdif = 0.0 + if pr >= 0.5 || (ν >= 1.0 && _tdist_cf_start_is_better(pr, ν)) + # Cornish-Fisher expansion of the t quantile in powers of 1/ν around the + # normal quantile (Fisher & Cornish, Technometrics 2 (1960), 209-225) + xn = norminvcdf(pr) + x = xn * xn + tp = (((((27.0 * x + 339.0) * x + 930.0) * x - 1782.0) * x - 765.0) * x + 17955.0) / (368640.0 * ν) + tp = (tp + ((((79.0 * x + 776.0) * x + 1482.0) * x - 1920.0) * x - 945.0) / 92160.0) / ν + tp = (tp + (((3.0 * x + 19.0) * x + 17.0) * x - 15.0) / 384.0) / ν + tp = (tp + ((5.0 * x + 16.0) * x + 3.0) / 96.0) / ν + tp = (tp + (x + 1.0) / 4.0) / ν + tp = xn * (1.0 + tp) + # for large ν the expansion has already converged to full precision and the + # Newton polish can be skipped (validated relative error <= 2.2e-15 in this + # region against a high-precision reference) + if ν >= 250.0 && x <= ν / 100.0 + return p > 0.5 ? -tp : tp + end + tprob = 0.0 + tpdif = 1.0 + abs(tp) + elseif ν < 1.0 + # leading power-law term of the tail solved in log space: + # pr ≈ (ν / t^2)^(ν / 2) / (ν * B) + lν = log(ν) + tp = -exp(lν / 2 - (log(pr) + lν + logB) / ν) + isfinite(tp) || return p > 0.5 ? -tp : tp + tprob, f = _tdistcdf_pdf(ν, B, tp) + if f < floatmin(Float64) + tpdif = 0.0 + else + tpdif = tprob / f * log1p((tprob - pr) / pr) + tpnew = tp - tpdif + tp = tpnew < 0.0 ? tpnew : tp / 2 + end + else + # invert the leading power-law term of the tail: pr ≈ f(0) * √ν * |t|^(-ν) + # where f(0) = 1 / (√ν * B) is the density at zero, map back to the t scale, + # apply one closed-form next-order correction, and one Newton step + u = exp(-log(ν * B * pr) / ν) + tp = -sqrt(ν) * sqrt(u - 1.0) * sqrt(u + 1.0) + isfinite(tp) || return p > 0.5 ? -tp : tp + tpdif = tp / ν + tpdif = -log1p((0.5 - 1.0 / (ν + 2.0)) / (1.0 + tpdif * tp)) * (tpdif + 1.0 / tp) + tp -= tpdif + tprob, f = _tdistcdf_pdf(ν, B, tp) + if f < floatmin(Float64) + tpdif = 0.0 + else + tpdif = tprob / f * log1p((tprob - pr) / pr) + tpnew = tp - tpdif + tp = tpnew < 0.0 ? tpnew : tp / 2 + end + end + + # Newton iteration applied to log(cdf): step = cdf / pdf * log(cdf / pr). + # Near the center this is an ordinary Newton step; in the tails, where the cdf + # is nearly exponential in t, the log transform makes the problem nearly linear. + iter = 0 + while abs(tprob - pr) > smalllpr && abs(tpdif) > small * (1.0 + abs(tp)) + (iter += 1) > 100 && break + tprob, f = _tdistcdf_pdf(ν, B, tp) + f < floatmin(Float64) && break + tpdif = tprob / f * log1p((tprob - pr) / pr) + tpnew = tp - tpdif + # keep the iterate in the negative half-line, where the solution lies + tp = tpnew < 0.0 ? tpnew : tp / 2 + end + return p > 0.5 ? -tp : tp +end + +function tdistinvcdf(ν::T, p::T) where {T <: Union{Float16, Float32, Float64}} + return T(_tdistinvcdf(Float64(ν), Float64(p))) +end function tdistinvcdf(ν::T, p::T) where {T <: Real} if isinf(ν) return norminvcdf(p) diff --git a/test/misc.jl b/test/misc.jl index 378b454..83d5c41 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -98,3 +98,71 @@ @test isnan(@inferred(tdistinvcdf(0, 0.975))) end end + +@testitem "tdistinvcdf" begin + using StatsFuns + using Test + + @testset "reference values" begin + # reference values computed with a 512-bit MPFR Newton iteration on the + # regularized incomplete beta representation of the cdf + for (ν, p, t) in ( + (0.35, 0.499, -0.0041388393417127554), + (0.5, 0.3, -1.0095258786071661), + (0.5, 1.0e-8, -1.02849115631634e15), + (1.0, 1.0e-100, -3.1830988618379064e99), + (2.5, 1.0e-12, -55306.174076515817), + (5.0, 1.0e-8, -62.40450611096729), + (5.0, 0.95, 2.0150483733330233), + (12.0, 0.3, -0.53861766820191637), + (50.0, 0.999, 3.261409055798318), + (300.0, 0.025, -1.9679030112610869), + (1000.0, 0.975, 1.9623390808264081), + (1.0e4, 1.0e-20, -9.282474153254304), + (1.0e6, 1 - 1.0e-12, 7.0345756932732169), + ) + @test tdistinvcdf(ν, p) ≈ t rtol = 1.0e-13 + @test tdistinvccdf(ν, p) == -tdistinvcdf(ν, p) + end + end + + @testset "round trips" begin + # t -> cdf -> invcdf on a grid where the cdf itself is accurate; deep + # tails are covered by the reference values above since tdistcdf + # currently loses precision there + # the quantile is always represented through its own tail: going through + # the complementary probability is ill-conditioned for any implementation + for ν in (0.5, 1.0, 2.5, 5.0, 20.0, 100.0, 1.0e3, 1.0e6), + t in (-8.0, -3.0, -0.5, 0.0, 1.0, 6.0) + + if t <= 0 + p = tdistcdf(ν, t) + 0 < p < 1 || continue + @test tdistinvcdf(ν, p) ≈ t atol = 1.0e-14 rtol = 1.0e-11 + else + q = tdistccdf(ν, t) + 0 < q < 1 || continue + @test tdistinvccdf(ν, q) ≈ t atol = 1.0e-14 rtol = 1.0e-11 + end + end + end + + @testset "edge cases and types" begin + @test tdistinvcdf(5, 0.0) == -Inf + @test tdistinvcdf(5, 1.0) == Inf + @test tdistinvcdf(5, 0.5) === 0.0 + @test isnan(tdistinvcdf(5.0, NaN)) + @test isnan(tdistinvcdf(NaN, 0.5)) + @test isnan(tdistinvcdf(5.0, -0.1)) + @test isnan(tdistinvcdf(5.0, 1.1)) + @test isnan(tdistinvcdf(-1.0, 0.5)) + @test tdistinvcdf(Inf, 0.975) == norminvcdf(0.975) + @test @inferred(tdistinvcdf(1, 0.75f0)) isa Float32 + @test @inferred(tdistinvcdf(Float16(1), Float16(0.75))) isa Float16 + @test @inferred(tdistinvcdf(1, 0.75)) isa Float64 + # StatsFuns#228: deep tails for small ν used to underflow to -Inf + @test isfinite(tdistinvcdf(0.5, 1.0e-8)) + # symmetry (1 - 0.75 is exact in binary) + @test tdistinvcdf(3.5, 0.25) == -tdistinvcdf(3.5, 0.75) + end +end From d31e6e8733a91cda37adb108ee45890b70a0952a Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 18 Aug 2026 06:22:21 +0200 Subject: [PATCH 2/3] Route all Real argument types of tdistinvcdf through the Float64 kernel Previously only mixed-type calls were promoted to floats; a direct same-type call such as tdistinvcdf(1//2, 1//2) dispatched to the old fdistinvccdf-based method and inherited its deep-tail failures. There is now a single generic method that converts to Float64, evaluates the kernel, and converts back to float(T), like the Rmath-based functions elsewhere in the package. As a side effect, BigFloat arguments now return a BigFloat computed at Float64 precision where they previously threw a MethodError. --- src/distrs/tdist.jl | 15 ++++----------- test/misc.jl | 3 +++ 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl index 5e7b852..6699b81 100644 --- a/src/distrs/tdist.jl +++ b/src/distrs/tdist.jl @@ -169,19 +169,12 @@ function _tdistinvcdf(ν::Float64, p::Float64) return p > 0.5 ? -tp : tp end -function tdistinvcdf(ν::T, p::T) where {T <: Union{Float16, Float32, Float64}} - return T(_tdistinvcdf(Float64(ν), Float64(p))) -end +# The kernel operates in Float64, like the Rmath-based functions elsewhere in the +# package; argument types with more precision than Float64 are truncated. function tdistinvcdf(ν::T, p::T) where {T <: Real} - if isinf(ν) - return norminvcdf(p) - elseif p < 0.5 - return -sqrt(fdistinvccdf(one(ν), ν, 2 * p)) - else - return sqrt(fdistinvccdf(one(ν), ν, 2 * (1 - p))) - end + return convert(float(T), _tdistinvcdf(Float64(ν), Float64(p))) end -tdistinvcdf(ν::Real, p::Real) = tdistinvcdf(map(float, promote(ν, p))...) +tdistinvcdf(ν::Real, p::Real) = tdistinvcdf(promote(ν, p)...) tdistinvccdf(ν::Real, p::Real) = -tdistinvcdf(ν, p) function tdistinvlogcdf(ν::T, logp::T) where {T <: Real} diff --git a/test/misc.jl b/test/misc.jl index 83d5c41..b29956a 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -160,6 +160,9 @@ end @test @inferred(tdistinvcdf(1, 0.75f0)) isa Float32 @test @inferred(tdistinvcdf(Float16(1), Float16(0.75))) isa Float16 @test @inferred(tdistinvcdf(1, 0.75)) isa Float64 + # same-type non-IEEE arguments use the same kernel + @test tdistinvcdf(1 // 2, 1 // 100) == tdistinvcdf(0.5, 0.01) + @test @inferred(tdistinvcdf(5, 1)) == Inf # StatsFuns#228: deep tails for small ν used to underflow to -Inf @test isfinite(tdistinvcdf(0.5, 1.0e-8)) # symmetry (1 - 0.75 is exact in binary) From e0b0cee05606f7dcd9a201191ef99db095b920ff Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 18 Aug 2026 08:58:51 +0200 Subject: [PATCH 3/3] Address review comments and avoid the verification-only cdf evaluation Dispatch: only Float16 and Float32 are routed explicitly through the Float64 kernel, as suggested in review; wider types such as BigFloat now throw a MethodError instead of being silently computed at Float64 precision. Stopping rule: the Newton polish previously certified convergence by evaluating the cdf once more after the converging step. The loop is restructured so that after each step the remaining error is bounded by the quadratic-convergence estimate |g''/g'| * step^2 / 2 for g = log(cdf), with a 64x safety margin; the extra evaluation is skipped when the bound is below the tolerance. Re-validated against the 512-bit reference: per-df maxima are unchanged (3.8e-14 max for df in [1, 10000], 1.3e-13 down to df = 0.05). The restructuring also merges the duplicated density-underflow guards into a single loop guard and removes the uncovered local declaration flagged by codecov; two reference tests where the density underflows cover the remaining early exit. --- src/distrs/tdist.jl | 85 +++++++++++++++++++++------------------------ test/misc.jl | 4 +++ 2 files changed, 43 insertions(+), 46 deletions(-) diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl index 6699b81..3dbed97 100644 --- a/src/distrs/tdist.jl +++ b/src/distrs/tdist.jl @@ -97,84 +97,77 @@ function _tdistinvcdf(ν::Float64, p::Float64) logB = logbeta(ν / 2, 0.5) B = exp(logB) - local tp::Float64 - tprob = 0.0 - tpdif = 0.0 - if pr >= 0.5 || (ν >= 1.0 && _tdist_cf_start_is_better(pr, ν)) + tp = if pr >= 0.5 || (ν >= 1.0 && _tdist_cf_start_is_better(pr, ν)) # Cornish-Fisher expansion of the t quantile in powers of 1/ν around the # normal quantile (Fisher & Cornish, Technometrics 2 (1960), 209-225) xn = norminvcdf(pr) x = xn * xn - tp = (((((27.0 * x + 339.0) * x + 930.0) * x - 1782.0) * x - 765.0) * x + 17955.0) / (368640.0 * ν) - tp = (tp + ((((79.0 * x + 776.0) * x + 1482.0) * x - 1920.0) * x - 945.0) / 92160.0) / ν - tp = (tp + (((3.0 * x + 19.0) * x + 17.0) * x - 15.0) / 384.0) / ν - tp = (tp + ((5.0 * x + 16.0) * x + 3.0) / 96.0) / ν - tp = (tp + (x + 1.0) / 4.0) / ν - tp = xn * (1.0 + tp) + t = (((((27.0 * x + 339.0) * x + 930.0) * x - 1782.0) * x - 765.0) * x + 17955.0) / (368640.0 * ν) + t = (t + ((((79.0 * x + 776.0) * x + 1482.0) * x - 1920.0) * x - 945.0) / 92160.0) / ν + t = (t + (((3.0 * x + 19.0) * x + 17.0) * x - 15.0) / 384.0) / ν + t = (t + ((5.0 * x + 16.0) * x + 3.0) / 96.0) / ν + t = (t + (x + 1.0) / 4.0) / ν + t = xn * (1.0 + t) # for large ν the expansion has already converged to full precision and the # Newton polish can be skipped (validated relative error <= 2.2e-15 in this # region against a high-precision reference) if ν >= 250.0 && x <= ν / 100.0 - return p > 0.5 ? -tp : tp + return p > 0.5 ? -t : t end - tprob = 0.0 - tpdif = 1.0 + abs(tp) + t elseif ν < 1.0 # leading power-law term of the tail solved in log space: # pr ≈ (ν / t^2)^(ν / 2) / (ν * B) lν = log(ν) - tp = -exp(lν / 2 - (log(pr) + lν + logB) / ν) - isfinite(tp) || return p > 0.5 ? -tp : tp - tprob, f = _tdistcdf_pdf(ν, B, tp) - if f < floatmin(Float64) - tpdif = 0.0 - else - tpdif = tprob / f * log1p((tprob - pr) / pr) - tpnew = tp - tpdif - tp = tpnew < 0.0 ? tpnew : tp / 2 - end + -exp(lν / 2 - (log(pr) + lν + logB) / ν) else - # invert the leading power-law term of the tail: pr ≈ f(0) * √ν * |t|^(-ν) - # where f(0) = 1 / (√ν * B) is the density at zero, map back to the t scale, - # apply one closed-form next-order correction, and one Newton step + # invert the leading power-law term of the tail, pr ≈ f(0) * √ν * |t|^(-ν) + # where f(0) = 1 / (√ν * B) is the density at zero, map back to the t + # scale, and apply one closed-form next-order correction u = exp(-log(ν * B * pr) / ν) - tp = -sqrt(ν) * sqrt(u - 1.0) * sqrt(u + 1.0) - isfinite(tp) || return p > 0.5 ? -tp : tp - tpdif = tp / ν - tpdif = -log1p((0.5 - 1.0 / (ν + 2.0)) / (1.0 + tpdif * tp)) * (tpdif + 1.0 / tp) - tp -= tpdif - tprob, f = _tdistcdf_pdf(ν, B, tp) - if f < floatmin(Float64) - tpdif = 0.0 - else - tpdif = tprob / f * log1p((tprob - pr) / pr) - tpnew = tp - tpdif - tp = tpnew < 0.0 ? tpnew : tp / 2 + t = -sqrt(ν) * sqrt(u - 1.0) * sqrt(u + 1.0) + if isfinite(t) + d = t / ν + t -= -log1p((0.5 - 1.0 / (ν + 2.0)) / (1.0 + d * t)) * (d + 1.0 / t) end + t end + # the true quantile may overflow in the extreme tails + isfinite(tp) || return p > 0.5 ? -tp : tp # Newton iteration applied to log(cdf): step = cdf / pdf * log(cdf / pr). # Near the center this is an ordinary Newton step; in the tails, where the cdf # is nearly exponential in t, the log transform makes the problem nearly linear. iter = 0 - while abs(tprob - pr) > smalllpr && abs(tpdif) > small * (1.0 + abs(tp)) - (iter += 1) > 100 && break + while true tprob, f = _tdistcdf_pdf(ν, B, tp) + # the density underflows only where the start already carries full precision f < floatmin(Float64) && break tpdif = tprob / f * log1p((tprob - pr) / pr) + # second-order coefficient of the iteration, |d²log(cdf)/dt² / dlog(cdf)/dt|, + # for the quadratic-convergence estimate of the error remaining after the step + curv = abs(-(ν + 1.0) * tp / (ν + tp * tp) - f / tprob) tpnew = tp - tpdif # keep the iterate in the negative half-line, where the solution lies tp = tpnew < 0.0 ? tpnew : tp / 2 + tol = small * (1.0 + abs(tp)) + # converged if the cdf already matched, the step was negligible, or the + # estimated error remaining after the step is negligible (the exact + # second-order coefficient is 1/2; the factor 32 is a 64x safety margin) + abs(tprob - pr) <= smalllpr && break + abs(tpdif) <= tol && break + 32.0 * curv * tpdif * tpdif <= tol && break + (iter += 1) >= 100 && break end return p > 0.5 ? -tp : tp end -# The kernel operates in Float64, like the Rmath-based functions elsewhere in the -# package; argument types with more precision than Float64 are truncated. -function tdistinvcdf(ν::T, p::T) where {T <: Real} - return convert(float(T), _tdistinvcdf(Float64(ν), Float64(p))) -end -tdistinvcdf(ν::Real, p::Real) = tdistinvcdf(promote(ν, p)...) +# Only Float16 and Float32 are routed explicitly through the Float64 kernel; +# wider types such as BigFloat are unsupported rather than silently computed +# at Float64 precision. +_tdistinvcdf(ν::Float16, p::Float16) = convert(Float16, _tdistinvcdf(Float64(ν), Float64(p))) +_tdistinvcdf(ν::Float32, p::Float32) = convert(Float32, _tdistinvcdf(Float64(ν), Float64(p))) +tdistinvcdf(ν::Real, p::Real) = _tdistinvcdf(map(float, promote(ν, p))...) tdistinvccdf(ν::Real, p::Real) = -tdistinvcdf(ν, p) function tdistinvlogcdf(ν::T, logp::T) where {T <: Real} diff --git a/test/misc.jl b/test/misc.jl index b29956a..007224f 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -120,6 +120,10 @@ end (1000.0, 0.975, 1.9623390808264081), (1.0e4, 1.0e-20, -9.282474153254304), (1.0e6, 1 - 1.0e-12, 7.0345756932732169), + # the density underflows at these two points, exercising the + # early exit from the Newton polish + (0.9, 1.0e-150, -1.2783541486158767e166), + (1.5, 1.0e-200, -1.1245005997832135e133), ) @test tdistinvcdf(ν, p) ≈ t rtol = 1.0e-13 @test tdistinvccdf(ν, p) == -tdistinvcdf(ν, p)