pnbeta {DPQ} | R Documentation |
Noncentral Beta Probabilities
Description
pnbetaAppr2()
and its inital version pnbetaAppr2v1()
provide the “approximation 2” of Chattamvelli and Shanmugam(1997)
to the noncentral Beta probability distribution.
pnbetaAS310()
is an R level interface to a C translation (and
“Rification”) of the AS 310
Fortran implementation.
Usage
pnbetaAppr2(x, a, b, ncp = 0, lower.tail = TRUE, log.p = FALSE)
pnbetaAS310(x, a, b, ncp = 0, lower.tail = TRUE, log.p = FALSE,
useAS226 = (ncp < 54.),
errmax = 1e-6, itrmax = 100)
Arguments
x |
numeric vector (of quantiles), typically from inside |
a , b |
the shape parameters of Beta, aka as |
ncp |
non-centrality parameter. |
log.p |
logical; if TRUE, probabilities p are given as log(p). |
lower.tail |
logical; if TRUE (default), probabilities are
|
useAS226 |
|
errmax |
non-negative number determining convergence for AS 310. |
itrmax |
positive integer number, only |
Value
a numeric vector of (log) probabilities of the same length as x
.
Note
The authors in the reference compare AS 310 with Lam(1995), Frick(1990) and Lenth(1987)
and state to be better than them. R's current (2019) noncentral beta
implementation builds on these, too, with some amendments though; still,
pnbetaAS310()
may potentially be better, at least in certain
corners of the 4-dimensional input space.
Author(s)
Martin Maechler; pnbetaAppr2()
in Oct 2007.
References
Chattamvelli, R., and Shanmugam, R. (1997) Algorithm AS 310: Computing the Non-Central Beta Distribution Function. Journal of the Royal Statistical Society. Series C (Applied Statistics) 46(1), 146–156, for “approximation 2” notably p.154; doi:10.1111/1467-9876.00055.
Lenth, R. V. (1987) Algorithm AS 226, ...,
Frick, H. (1990)'s AS R84, ..., and
Lam, M.L. (1995)'s AS R95 : See ‘References’ in R's pbeta
page.
See Also
R's own pbeta
.
Examples
## Same arguments as for Table 1 (p.151) of the reference
a <- 5*rep(1:3, each=3)
aargs <- cbind(a = a, b = a,
ncp = rep(c(54, 140, 170), 3),
x = 1e-4*c(8640, 9000, 9560, 8686, 9000, 9000, 8787, 9000, 9220))
aargs
pnbA2 <- apply(aargs, 1, function(aa) do.call(pnbetaAppr2, as.list(aa)))
pnA310<- apply(aargs, 1, function(aa) do.call(pnbetaAS310, as.list(aa)))
aar2 <- aargs; dimnames(aar2)[[2]] <- c(paste0("shape", 1:2), "ncp", "q")
pnbR <- apply(aar2, 1, function(aa) do.call(pbeta, as.list(aa)))
range(relD2 <- 1 - pnbA2 /pnbR)
range(relD310 <- 1 - pnA310/pnbR)
cbind(aargs, pnbA2, pnA310, pnbR,
relD2 = signif(relD2, 3), relD310 = signif(relD310, 3)) # <------> Table 1
stopifnot(abs(relD2) < 0.009) # max is 0.006286
stopifnot(abs(relD310) < 1e-5 ) # max is 6.3732e-6
## Arguments as for Table 2 (p.152) of the reference :
aarg2 <- cbind(a = c( 10, 10, 15, 20, 20, 20, 30, 30),
b = c( 20, 10, 5, 10, 30, 50, 20, 40),
ncp=c(150,120, 80,110, 65,130, 80,130),
x = c(868,900,880,850,660,720,720,800)/1000)
pnbA2 <- apply(aarg2, 1, function(aa) do.call(pnbetaAppr2, as.list(aa)))
pnA310<- apply(aarg2, 1, function(aa) do.call(pnbetaAS310, as.list(aa)))
aar2 <- aarg2; dimnames(aar2)[[2]] <- c(paste0("shape", 1:2), "ncp", "q")
pnbR <- apply(aar2, 1, function(aa) do.call(pbeta, as.list(aa)))
range(relD2 <- 1 - pnbA2 /pnbR)
range(relD310 <- 1 - pnA310/pnbR)
cbind(aarg2, pnbA2, pnA310, pnbR,
relD2 = signif(relD2, 3), relD310 = signif(relD310, 3)) # <------> Table 2
stopifnot(abs(relD2 ) < 0.006) # max is 0.00412
stopifnot(abs(relD310) < 1e-5 ) # max is 5.5953e-6
## Arguments as for Table 3 (p.152) of the reference :
aarg3 <- cbind(a = c( 10, 10, 10, 15, 10, 12, 30, 35),
b = c( 5, 10, 30, 20, 5, 17, 30, 30),
ncp=c( 20, 54, 80,120, 55, 64,140, 20),
x = c(644,700,780,760,795,560,800,670)/1000)
pnbA3 <- apply(aarg3, 1, function(aa) do.call(pnbetaAppr2, as.list(aa)))
pnA310<- apply(aarg3, 1, function(aa) do.call(pnbetaAS310, as.list(aa)))
aar3 <- aarg3; dimnames(aar3)[[2]] <- c(paste0("shape", 1:2), "ncp", "q")
pnbR <- apply(aar3, 1, function(aa) do.call(pbeta, as.list(aa)))
range(relD2 <- 1 - pnbA3 /pnbR)
range(relD310 <- 1 - pnA310/pnbR)
cbind(aarg3, pnbA3, pnA310, pnbR,
relD2 = signif(relD2, 3), relD310 = signif(relD310, 3)) # <------> Table 3
stopifnot(abs(relD2 ) < 0.09) # max is 0.06337
stopifnot(abs(relD310) < 1e-4) # max is 3.898e-5