qtU: 'uniroot()'-based Computing of t-Distribution Quantiles

Description

Currently, R's own qt() (aka qnt() in the non-central case) uses simple inversion of pt to compute quantiles in the case where ncp is specified.
That simple inversion (of pnt()) has seen to be deficient, even in cases where pnt(), i.e., R's pt(.., ncp=*) does not loose accuracy.

This uniroot()-based inversion does not suffer from these deficits in some cases.
qtU() is simply defined as

qtU <- Vectorize(qtU1, c("p","df","ncp"))

where in qtU1(p, df, ncp, *) each of (p, df, ncp) must be of length one.

Usage

qtU1(p, df, ncp, lower.tail = TRUE, log.p = FALSE, interval = c(-10, 10),
     tol = 1e-05, verbose = FALSE, ...)
qtU (p, df, ncp, lower.tail = TRUE, log.p = FALSE, interval = c(-10, 10),
     tol = 1e-05, verbose = FALSE, ...)

Value

numeric vector of t quantiles, properly recycled in (p, df, ncp).

Arguments

p, df, ncp: vectors of probabilities, degrees of freedom, and non-centrality parameter; see qt. As there, ncp may be missing which amounts to being zero.
lower.tail, log.p: logical; see qt.
interval: the interval in which quantiles should be searched; passed to uniroot(); the current default is arbitrary and suboptimal; when pt(q,*) is accurate enough and hence montone (increasing iff lower.tail), this interval is automatically correctly extended by uniroot().
tol: non-negative convergence tolerance passed to uniroot().
verbose: logical indicating if every call of the objective function should produce a line of console output.
...: optional further arguments passed to uniroot().

Author

Martin Maechler

Examples

Run this code

qtU1 # simple definition {with extras only for  'verbose = TRUE'}

## An example, seen to be deficient
## Stephen Berman to R-help, 13 June 2022,
## "Why does qt() return Inf with certain negative ncp values?"
q2 <- seq(-3/4, -1/4, by=1/128)
pq2 <- pt(q2, 35, ncp=-7, lower.tail=FALSE)
### ==> via qtU(), a simple uniroot() - based inversion of pt()
qpqU  <- qtU(pq2, 35, ncp=-7, lower.tail=FALSE, tol=1e-10)
stopifnot(all.equal(q2, qpqU, tol=1e-9)) # perfect!

## These two currently (2022-06-14) give Inf  whereas qtU() works fine
qt  (9e-12, df=35, ncp=-7, lower.tail=FALSE) # warnings; --> Inf
qntR(9e-12, df=35, ncp=-7, lower.tail=FALSE) #  (ditto)
## verbose = TRUE  shows all calls to pt():
qtU1(9e-12, df=35, ncp=-7, lower.tail=FALSE, verbose=TRUE)