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)) )