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.
1 2 
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 1letter string specifying how rounding
should happen at Clevel 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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40  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 = 1e230),
all.equal( Ip, pbetaI(x, a, b, precBits = 2048, lower.tail=FALSE))
)
}
rErr < function(approx, true, eps = 1e200) {
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=1e15) 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))
)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.