| pbetaI {Rmpfr} | R Documentation |
For integers a, b, I(x; a,b) aka
pbeta(x, a,b) is a polynomial in x with rational coefficients,
and hence arbitarily accurately computable.
pbetaI(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE,
precBits = NULL, rnd.mode = c("N","D","U","Z","A"))
q |
called x, above; vector of quantiles, in [0,1]. |
shape1, shape2 |
the positive Beta “shape” parameters, called a, b, above. Must be integer valued for this function. |
ncp |
unused, only for compatibility with |
lower.tail |
logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x]. |
log.p |
logical; if TRUE, probabilities p are given as log(p). |
precBits |
the precision (in number of bits) to be used in
|
rnd.mode |
a 1-letter string specifying how rounding
should happen at C-level conversion to MPFR, see |
an "mpfr" vector of the same length as q.
For upper tail probabilities, i.e., when lower.tail=FALSE,
we may need large precBits, because the implicit or explicit
1 - P computation suffers from severe cancellation.
Martin Maechler
x <- (0:12)/16 # not all the way up ..
a <- 7; b <- 788
p. <- pbetaI(x, a, b) ## still slow: %% TOO slow -- FIXME
pp <- pbetaI(x, a, b, precBits = 2048)
## Currently, the lower.tail=FALSE are computed "badly":
lp <- log(pp) ## = pbetaI(x, a, b, log.p=TRUE)
lIp <- log1p(-pp) ## = pbetaI(x, a, b, lower.tail=FALSE, log.p=TRUE)
Ip <- 1 - pp ## = pbetaI(x, a, b, lower.tail=FALSE)
if(Rmpfr:::doExtras()) { ## somewhat slow
stopifnot(
all.equal(lp, pbetaI(x, a, b, precBits = 2048, log.p=TRUE)),
all.equal(lIp, pbetaI(x, a, b, precBits = 2048, lower.tail=FALSE, log.p=TRUE),
tol = 1e-230),
all.equal( Ip, pbetaI(x, a, b, precBits = 2048, lower.tail=FALSE))
)
}
rErr <- function(approx, true, eps = 1e-200) {
true <- as.numeric(true) # for "mpfr"
ifelse(Mod(true) >= eps,
## relative error, catching '-Inf' etc :
ifelse(true == approx, 0, 1 - approx / true),
## else: absolute error (e.g. when true=0)
true - approx)
}
rErr(pbeta(x, a, b), pp)
rErr(pbeta(x, a, b, lower=FALSE), Ip)
rErr(pbeta(x, a, b, log = TRUE), lp)
rErr(pbeta(x, a, b, lower=FALSE, log = TRUE), lIp)
a.EQ <- function(..., tol=1e-15) all.equal(..., tolerance=tol)
stopifnot(
a.EQ(pp, pbeta(x, a, b)),
a.EQ(lp, pbeta(x, a, b, log.p=TRUE)),
a.EQ(lIp, pbeta(x, a, b, lower.tail=FALSE, log.p=TRUE)),
a.EQ( Ip, pbeta(x, a, b, lower.tail=FALSE))
)