dtWV: Asymptotic Noncentral t Distribution Density by Viechtbauer

In DPQ: Density, Probability, Quantile ('DPQ') Computations

Asymptotic Noncentral t Distribution Density by Viechtbauer

Description

Compute the density function f(x) of the t distribution with df degrees of freedom and non-centrality parameter ncp, according to Wolfgang Viechtbauer's proposal in 2002. This is an asymptotic formula for “large” df = \nu, or mathematically \nu \to \infty.

Usage

dtWV(x, df, ncp = 0, log = FALSE)

Arguments

x	numeric vector.
df	degrees of freedom (> 0, maybe non-integer). df = Inf is allowed.
ncp	non-centrality parameter \delta; If omitted, use the central t distribution.
log	logical; if TRUE, log(f(x)) is returned instead of f(x).

Details

The formula used is “asymptotic”: Resnikoff and Lieberman (1957), p.1 and p.25ff, proposed to use recursive polynomials for (integer !) degrees of freedom f = 1,2,\dots, 20, and then, for df = f > 20, use the asymptotic approximation which Wolfgang Viechtbauer proposed as a first version of a non-central t density for R (when dt() did not yet have an ncp argument).

Value

numeric vector of density values, properly recycled in (x, df, ncp).

Author(s)

Wolfgang Viechtbauer (2002) post to R-help (https://stat.ethz.ch/pipermail/r-help/2002-October/026044.html), and Martin Maechler (log argument; tweaks, notably recycling).

References

Resnikoff, George J. and Lieberman, Gerald J. (1957) Tables of the non-central t-distribution; Technical report no. 32 (LIE ONR 32), April 1, 1957; Applied Math. and Stat. Lab., Stanford University. https://statistics.stanford.edu/technical-reports/tables-non-central-t-distribution-density-function-cumulative-distribution

See Also

dt, R's (C level) implementation of the (non-central) t density; dntJKBf, for Johnson et al.'s summation formula approximation.

Examples

tt <- seq(0, 10, len = 21)
ncp <- seq(0, 6, len = 31)
dt3R  <- outer(tt, ncp, dt  , df = 3)
dt3WV <- outer(tt, ncp, dtWV, df = 3)
all.equal(dt3R, dt3WV) # rel.err 0.00063
dt25R  <- outer(tt, ncp, dt  , df = 25)
dt25WV <- outer(tt, ncp, dtWV, df = 25)
all.equal(dt25R, dt25WV) # rel.err 1.1e-5

x <- -10:700
fx  <- dt  (x, df = 22, ncp =100)
lfx <- dt  (x, df = 22, ncp =100, log=TRUE)
lfV <- dtWV(x, df = 22, ncp =100, log=TRUE)

head(lfx, 20) # shows that R's dt(*, log=TRUE) implementation is "quite suboptimal"

## graphics
opa <- par(no.readonly=TRUE)
par(mar=.1+c(5,4,4,3), mgp = c(2, .8,0))
plot(fx ~ x, type="l")
par(new=TRUE) ; cc <- c("red", adjustcolor("orange", 0.4))
plot(lfx ~ x, type = "o", pch=".", col=cc[1], cex=2, ann=FALSE, yaxt="n")
sfsmisc::eaxis(4, col=cc[1], col.axis=cc[1], small.args = list(col=cc[1]))
lines(x, lfV, col=cc[2], lwd=3)
dtt1 <- "      dt"; dtt2 <- "(x, df=22, ncp=100"; dttL <- paste0(dtt2,", log=TRUE)")
legend("right", c(paste0(dtt1,dtt2,")"), paste0(c(dtt1,"dtWV"), dttL)),
       lty=1, lwd=c(1,1,3), col=c("black", cc), bty = "n")
par(opa) # reset

DPQ documentation built on Nov. 3, 2024, 3 a.m.