Implement tdistinvcdf natively instead of inverting the incomplete beta - #229
Implement tdistinvcdf natively instead of inverting the incomplete beta#229andreasnoack wants to merge 3 commits into
Conversation
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).
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #229 +/- ##
==========================================
+ Coverage 75.50% 77.57% +2.06%
==========================================
Files 23 23
Lines 1135 1204 +69
==========================================
+ Hits 857 934 +77
+ Misses 278 270 -8 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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.
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.
|
The benchmark looks great, in particular with the latest commit 👍 |
|
@devmotion do you have any additional comments? |
| while true | ||
| tprob, f = _tdistcdf_pdf(ν, B, tp) | ||
| # the density underflows only where the start already carries full precision | ||
| f < floatmin(Float64) && break |
There was a problem hiding this comment.
For ν < 1 the exponent (ν + 1) / 2 is smaller than 1, so f is much larger than k2 and stays well above floatmin long after k2 itself has become subnormal. In that window beta_inc only sees a handful of significant bits of its argument, tprob comes back wrong in the third digit, and the step throws away a start that was already good to ~3e-14. I get a relative error of 3.4e-1 for (0.9, 1e-146) and 1.7e-1 for (0.35, 1e-57), and for eg (0.5, 1e-81) and (0.9, 1e-145) the loop runs into the iteration cap - 100 beta_inc evaluations, ~11 µs instead of ~200 ns. Per-df maxima over p ∈ [1e-300, 1/2]: 3.1e-8 at ν = 0.2, 1.7e-1 at 0.35, 1.5e-2 at 0.5, 8.9e-4 at 0.7, 3.4e-1 at 0.9. qt is fine (~5e-14) at all of these.
Could we guard on the incomplete beta argument rather than on the density, ie return k2 from _tdistcdf_pdf and break on k2 < floatmin(Float64)? With that the maximum for ν < 1 is back to 1.6e-13. The comment above assumes the two underflow together, which only holds for ν ≥ 1.
| 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)) |
There was a problem hiding this comment.
tol is essentially absolute for |tp| ≪ 1, and (tprob - pr) / pr cancels once tprob ≈ pr ≈ 1/2, so we lose relative accuracy close to the median. For ν = 1 I get 2.7e-13 at p = 0.5 - 1e-4, 2.6e-11 at 0.5 - 1e-6, 4.6e-8 at 0.5 - 1e-10 and 4.2e-4 at 0.5 - 1e-14, where the current implementation stays at ~1e-16 throughout (Rmath is as bad as this PR here). Same shape for every ν up to ~100. The absolute error is fine everywhere, so perhaps this is acceptable, but it is a regression and the description claims accuracy equivalence with qt without mentioning it.
It also looks cheap to avoid: beta_inc already returns 1/2 - cdf as its second component, and 0.5 - pr is exact for pr ≥ 1/4. Iterating on those two for pr > 1/4 gives me ~1e-16 over the whole range with the same number of beta_inc evaluations.
| small = 1.0e-14 | ||
| smalllpr = -small * log(pr) * pr | ||
|
|
||
| logB = logbeta(ν / 2, 0.5) |
There was a problem hiding this comment.
Neither logB nor B is used before the early return for large ν below, and together they cost ~65 ns. Could we move them into the branches that need them? That takes tdistinvcdf(1000, 0.975) and tdistinvcdf(1e5, 0.999) from ~105 ns to ~50 ns here.
| -exp(lν / 2 - (log(pr) + lν + logB) / ν) | ||
| else | ||
| return sqrt(fdistinvccdf(one(ν), ν, 2 * (1 - p))) | ||
| # invert the leading power-law term of the tail, pr ≈ f(0) * √ν * |t|^(-ν) |
There was a problem hiding this comment.
I don't quite follow this one: f(0) * √ν * |t|^(-ν) is |t|^(-ν) / B, and inverting that gives (B * pr)^(-1/ν) rather than √ν * (ν * B * pr)^(-1/ν). AFAICT the relation actually inverted here is the same as in the branch above, pr ≈ (ν / t^2)^(ν/2) / (ν * B), and the only difference is that sqrt(u - 1) * sqrt(u + 1) keeps x = ν / (ν + t^2) exact instead of dropping to the leading power of t. Could we use the same wording as above? The description has the same formula.
| logB = logbeta(ν / 2, 0.5) | ||
| B = exp(logB) | ||
|
|
||
| tp = if pr >= 0.5 || (ν >= 1.0 && _tdist_cf_start_is_better(pr, ν)) |
There was a problem hiding this comment.
Is pr >= 0.5 reachable? p == 0.5 returns earlier and 1 - p is exact for p > 1/2, so pr < 0.5 always as far as I can tell.
| tdistinvcdf(ν::Real, p::Real) = _tdistinvcdf(map(float, promote(ν, p))...) | ||
|
|
||
| tdistinvccdf(ν::Real, p::Real) = -tdistinvcdf(ν, p) | ||
| function tdistinvlogcdf(ν::T, logp::T) where {T <: Real} |
There was a problem hiding this comment.
Should tdistinvlogcdf go through the new kernel as well? Right now tdistinvcdf(0.5, 1e-8) is fixed but tdistinvlogcdf(0.5, log(1e-8)) still returns -Inf, and for (2.5, log(1e-12)) the two disagree by 7e-9. The tail start is computed in log space anyway, so it looks like log(pr) could just become logp there. Happy for it to be a follow-up.
| (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), |
There was a problem hiding this comment.
At both of these k2 underflows to exactly zero, so they exit through f == 0. The cases that go wrong are one decade earlier, where k2 is subnormal but nonzero - eg (0.9, 1e-144), (0.35, 1e-56) or (0.7, 1e-110), see my comment in tdist.jl. Could we add one of those? Nothing bounds the iteration count either, so the points that need 100 iterations pass silently.
| (1.5, 1.0e-200, -1.1245005997832135e133), | ||
| ) | ||
| @test tdistinvcdf(ν, p) ≈ t rtol = 1.0e-13 | ||
| @test tdistinvccdf(ν, p) == -tdistinvcdf(ν, p) |
There was a problem hiding this comment.
tdistinvccdf is defined as -tdistinvcdf, so this holds by construction - is it testing anything?
| end | ||
| end | ||
|
|
||
| @testset "round trips" begin |
There was a problem hiding this comment.
The round trips and the Rmath tests both stay away from p ≈ 1/2 - the Rmath grid only reaches p = 0.5 ± 0.03, and rmathcomp compares quantiles with atol = 2e-14 - so nothing here notices the loss of relative accuracy near the median. Would a couple of reference values at p = 0.5 - 1e-6 and similar be worth adding?
devmotion
left a comment
devmotion
left a comment
There was a problem hiding this comment.
I used Claude to check this against a 512-bit MPFR reference (own continued fraction for the incomplete beta, cross-checked against the closed forms for ν = 1 and ν = 2), and the numbers in the description reproduce for ν ≥ 1: the 15 reference values in the tests are all correct to ~1e-16, and over p ∈ [1e-300, 1/2] the maximum is 5.7e-14. Nice improvement.
Two things I ran into that I think are worth fixing, both commented inline. For ν < 1 there is a window of p where beta_inc is called with a subnormal argument, and the polish then makes the result much worse than the starting approximation it was given - up to 34% relative error, and in some cases 100 iterations, ie ~11 µs. And relative accuracy near the median is worse than what we have today.
On the remaining gap to qt for ν = 5: as far as I can tell it is not the number of iterations. Both implementations evaluate the cdf/pdf pair twice, and ours is the cheaper of the two, 76 ns for _tdistcdf_pdf against 84 ns for pt + dt here. The difference is almost exactly logbeta(ν / 2, 1/2) + exp, which Hill's algorithm does not need at all: 334 ns for the kernel, 270 ns with B passed in, and qt is 270 ns. So there are two levers if we care. B only scales the Newton step, so a few digits would do there (full precision is only needed for the tail start, where it enters as (ν B p)^(-1/ν)), and the second beta_inc at ν = 5 is only there because of the 64x margin - one step from the CF start already gives 1.0e-14 and the quadratic estimate predicts that well, 6.1e-15. A Halley step does not help, it gives 2.0e-14 from the same start.
2 participants
Summary
Addresses the quantile part of #228 (see the task list in my comment there): replace the
tdistinvcdfroute throughfdistinvccdf/beta_inc_inv— a Newton iteration where every step evaluates the full incomplete beta — with a pure Julia implementation based on VBAinvtdistby Ian Smith:pr ≈ f(0)·√ν·|t|^(−ν)with a closed-form next-order correction (crossover between the two starts follows Smith's empirical curves).log(cdf)— an ordinary Newton step near the center, nearly exact in the tails where the cdf is close to exponential in t. Each iteration costs onebeta_incevaluation; the pdf comes from the same argument reduction for one extrapow. Iterates are clamped to the half-line containing the solution.ν ≥ 250andxn² ≤ ν/100the Cornish–Fisher expansion has already converged to full precision and the polish is skipped (validated: relative error ≤ 2.2e-15 in this region, at the accuracy floor ofnorminvcdfitself).Only
Float16andFloat32are routed explicitly through the Float64 kernel (per review); wider types such asBigFloatthrow aMethodErrorrather than being silently computed at Float64 precision. The Newton polish stops without a verification-only cdf evaluation: after each step the remaining error is bounded by the quadratic-convergence estimate, with a 64x safety margin (per-df accuracy maxima re-validated unchanged).Accuracy
Against a 512-bit MPFR Newton reference on the incomplete beta representation:
qt(which shows up to 7.6e-9 at p = 1e-300 for small ν, and up to 5.5e-4 for ν < 1).tdistinvcdf(0.5, 1e-8)returned-Inf(true value ≈ −1.03e15), and (ν = 2.5, p = 1e-12) was only good to 7e-9.While testing round trips I also found that
tdistcdfitself loses precision in deep tails (e.g.tdistcdf(1, -3.18e7)returns 9.49e-9 instead of 1.0e-8) becausefdistccdfcollapses the tail intoinv(1 + ν/x²) ≈ 1. That is pre-existing and out of scope here, but the same paired argument reduction used by_tdistcdf_pdfin this PR would fix it — I can follow up.Performance
Apple Silicon, SpecialFunctions master (includes the merged but unreleased JuliaMath/SpecialFunctions.jl#542/#543):
qtWith the released SpecialFunctions 2.8.3 the
beta_inc-dependent rows are roughly 60-90 ns slower each.Test plan
qtis within the 2e-14 tolerance on the test grids) and thetdistinvcdf(0, 0.975) → NaNregression testt → cdf → invcdf(through the well-conditioned tail representation), edge cases (p ∈ {0, 0.5, 1},ν = Inf,NaNs, out-of-rangep), type stability forFloat16/32/64Disclosure: this PR was prepared by Claude Code at my direction; I have reviewed the implementation, the accuracy methodology, and the benchmarks.
🤖 Generated with Claude Code