dnt {DPQ} | R Documentation |
Non-central t-Distribution Density - Algorithms and Approximations
Description
dntJKBf1
implements the summation formulas
of Johnson, Kotz and Balakrishnan (1995),
(31.15) on page 516 and (31.15') on p.519, the latter being typo-corrected
for a missing factor 1 / j!
.
dntJKBf()
is Vectorize(dntJKBf1,
c("x","df","ncp"))
, i.e., works vectorized in all three main
arguments x
, df
and ncp
.
The functions .dntJKBch1()
and .dntJKBch()
are only there
for didactical reasons allowing to check that indeed formula (31.15')
in the reference is missing a j!
factor in the denominator.
The dntJKBf*()
functions are written to also work with
arbitrary precise numbers of class
"mpfr"
(from package Rmpfr)
as arguments.
Usage
dntJKBf1(x, df, ncp, log = FALSE, M = 1000)
dntJKBf (x, df, ncp, log = FALSE, M = 1000)
## The "checking" versions, only for proving correctness of formula:
.dntJKBch1(x, df, ncp, log = FALSE, M = 1000, check=FALSE, tol.check = 1e-7)
.dntJKBch (x, df, ncp, log = FALSE, M = 1000, check=FALSE, tol.check = 1e-7)
Arguments
x , df , ncp |
|
log |
as in |
M |
the number of terms to be used, a positive integer. |
check |
logical indicating if checks of the formula equalities should be done. |
tol.check |
tolerance to be used for |
Details
How to choose M
optimally has not been investigated yet.
Note that relatedly,
R's source code ‘R/src/nmath/dnt.c’ has claimed from 2003 till 2014
but wrongly that the noncentral t density f(x, *)
is
f(x, df, ncp) = df^(df/2) * exp(-.5*ncp^2) / (sqrt(pi)*gamma(df/2)*(df+x^2)^((df+1)/2)) * sum_{k=0}^Inf gamma((df + k + df)/2)*ncp^k / prod(1:k)*(2*x^2/(df+x^2))^(k/2) .
These functions (and this help page) prove that it was wrong.
Value
a number for dntJKBf1()
and .dntJKBch1()
.
a numeric vector of the same length as the maximum of the lengths of
x, df, ncp
for dntJKBf()
and .dntJKBch()
.
Author(s)
Martin Maechler
References
Johnson, N.L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions Vol~2, 2nd ed.; Wiley; chapter 31, Section 5 Distribution Function, p.514 ff
See Also
R's dt
;
(an improved version of) Viechtbauer's proposal: dtWV
.
Examples
tt <- seq(0, 10, len = 21)
ncp <- seq(0, 6, len = 31)
dt3R <- outer(tt, ncp, dt, df = 3)
dt3JKB <- outer(tt, ncp, dntJKBf, df = 3)
all.equal(dt3R, dt3JKB) # Lnx(64-b): 51 NA's in dt3R
x <- seq(-1,12, by=1/16)
fx <- dt(x, df=3, ncp=5)
re1 <- 1 - .dntJKBch(x, df=3, ncp=5) / fx ; summary(warnings()) # slow, with warnings
op <- options(warn = 2) # (=> warning == error, for now)
re2 <- 1 - dntJKBf (x, df=3, ncp=5) / fx # faster, no warnings
stopifnot(all.equal(re1[!is.na(re1)], re2[!is.na(re1)], tol=1e-6))
head( cbind(x, fx, re1, re2) , 20)
matplot(x, log10(abs(cbind(re1, re2))), type = "o", cex = 1/4)
## One of the numerical problems in "base R"'s non-central t-density:
options(warn = 0) # (factory def.)
x <- 2^seq(-12, 32, by=1/8) ; df <- 1/10
dtm <- cbind(dt(x, df=df, log=TRUE),
dt(x, df=df, ncp=df/2, log=TRUE),
dt(x, df=df, ncp=df, log=TRUE),
dt(x, df=df, ncp=df*2, log=TRUE)) #.. quite a few warnings:
summary(warnings())
matplot(x, dtm, type="l", log = "x", xaxt="n",
main = "dt(x, df=1/10, log=TRUE) central and noncentral")
sfsmisc::eaxis(1)
legend("right", legend=c("", paste0("ncp = df",c("/2","","*2"))),
lty=1:4, col=1:4, bty="n")
# ---- using MPFR high accuracy arithmetic (too slow for routine testing) ---
## no such kink here:
x. <- if(requireNamespace("Rmpfr")) Rmpfr::mpfr(x, 256) else x
system.time(dtJKB <- dntJKBf(x., df=df, ncp=df, log=TRUE)) # 21s (!) was only 7s ???
lines(x, dtJKB, col=adjustcolor(3, 1/2), lwd=3)
options(op) # reset to prev.
## Relative Difference / Approximation errors :
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x")
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x", xaxt="n", ylim=c(-1,1)*1e-3); sfsmisc::eaxis(1)
plot(x, 1 - dtJKB / dtm[,3], type="l", log="x", xaxt="n", ylim=c(-1,1)*1e-7); sfsmisc::eaxis(1)
plot(x, abs(1 - dtJKB / dtm[,3]), type="l", log="xy", axes=FALSE, main =
"dt(*, 1/10, 1/10, log=TRUE) relative approx. error",
sub= paste("Copyright © 2019 Martin Mächler --- ", R.version.string))
for(j in 1:2) sfsmisc::eaxis(j)