Student-t CDF and quantile are significantly slower than Rmath

[Santymax98]

While optimizing Student pair-copula conditionals in VineCopulas.jl, I found that the native Student-t CDF and especially the quantile can be substantially slower than the corresponding Rmath routines.

Using:

• StatsFuns v2.2.1
• Rmath v0.9.0

a small benchmark gives:

using BenchmarkTools, StatsFuns, Rmath

ν = 5.0
x = 1.5

@btime StatsFuns.tdistcdf($ν, $x)
# 355.423 ns (0 allocations)

@btime Rmath.pt($x, $ν)
# 230.705 ns (0 allocations)

For the quantile:

p = 0.95

@btime StatsFuns.tdistinvcdf($ν, $p)
# 1.396 μs (0 allocations)

@btime Rmath.qt($p, $ν)
# 561.773 ns (0 allocations)

and in the tail:

p = 1e-8

@btime StatsFuns.tdistinvcdf($ν, $p)
# 1.962 μs (0 allocations)

@btime Rmath.qt($p, $ν)
# 302.579 ns (0 allocations)

The numerical results agree essentially to machine precision.

For these examples, Rmath is about 1.5× faster for the CDF, 2.5× faster for a central quantile, and 6.5× faster for the tail quantile.

This matters downstream in Copulas.jl and VineCopulas.jl, where Student conditional CDFs and quantiles are repeatedly evaluated during Rosenblatt transforms and simulation. A related discussion is tracked in lrnv/Copulas.jl#381.

Would there be room to optimize the native tdistcdf / tdistinvcdf implementations, particularly the inverse CDF, while keeping the native Julia implementation?

[andreasnoack]

Thanks for the report. I looked into where the time goes. StatsFuns has no dedicated Student-t code: tdistcdf/tdistinvcdf route through the F distribution into SpecialFunctions' beta_inc/beta_inc_inv, so essentially all of the time is spent in the incomplete beta function, and the quantile additionally pays for inverting it with a Newton iteration where every step evaluates the full beta_inc.

For the CDF, profiling turned up two fixable issues in SpecialFunctions' translation of the NSWC incomplete beta routines, now opened as PRs:

• Restore integer truncation semantics from BPSER and BUP in the beta_inc translation JuliaMath/SpecialFunctions.jl#542 — a missed Fortran integer truncation that makes the power-series branch construct a TwicePrecision float range on every call
• Use the NSWC GAMLN1 rational approximation in the incomplete beta routines JuliaMath/SpecialFunctions.jl#543 — use the original's cheap GAMLN1 rational approximation instead of the ~8× more expensive loggamma1p where it is safe

Benchmarks (Apple Silicon, ν = 5; results are bit-identical throughout):

case	StatsFuns today	with both PRs	Rmath
tdistcdf(5, 1.5)	223 ns	129 ns	141 ns (pt)
tdistinvcdf(5, 0.95)	1004 ns	802 ns	378 ns (qt)
tdistinvcdf(5, 1e-8)	1642 ns	810 ns	229 ns (qt)

So the PRs fully close the CDF gap (slightly faster than Rmath afterwards), and they help the quantile since each Newton step gets cheaper — but the quantile remains 2–3.5× slower. That part is algorithmic: qt-style implementations use a dedicated closed-form Student-t quantile expansion (Hill's Algorithm 396, CACM 1970) with a one-step correction instead of numerically inverting the incomplete beta, so matching it requires a dedicated tdistinvcdf implementation rather than further tuning of beta_inc_inv.

Splitting the issue accordingly:

• tdistcdf performance: addressed by Restore integer truncation semantics from BPSER and BUP in the beta_inc translation JuliaMath/SpecialFunctions.jl#542 and Use the NSWC GAMLN1 rational approximation in the incomplete beta routines JuliaMath/SpecialFunctions.jl#543; done once they are merged, released, and the compat bound allows them here
• tdistinvcdf performance: implement a dedicated Student-t quantile (e.g. Hill's Algorithm 396 with a Taylor/Newton correction) instead of going through fdistinvccdf/beta_inc_inv — addressed by Implement tdistinvcdf natively instead of inverting the incomplete beta #229