Implement tdistinvcdf natively instead of inverting the incomplete beta

src/distrs/tdist.jl

@@ -36,16 +36,145 @@ tdistlogcdf(ν::Real, x::Real) = tdistlogcdf(map(float, promote(ν, x))...)
 
 tdistlogccdf(ν::Real, x::Real) = tdistlogcdf(ν, -x)
 
-function tdistinvcdf(ν::T, p::T) where {T <: Real}
-    if isinf(ν)
+# 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.5
-        return -sqrt(fdistinvccdf(one(ν), ν, 2 * 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
-        return sqrt(fdistinvccdf(one(ν), ν, 2 * (1 - p)))
+        # 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
+
+# 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(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}

test/misc.jl

@@ -98,3 +98,74 @@
         @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
+        # 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)
+        @test tdistinvcdf(3.5, 0.25) == -tdistinvcdf(3.5, 0.75)
+    end
+end