# pbetaI: Accurate Incomplete Beta / Beta Probabilities For Integer... In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

## Description

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.

## Usage

 ```1 2``` ```pbetaI(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE, precBits = NULL, rnd.mode = c("N","D","U","Z","A")) ```

## Arguments

 `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 `pbeta`, must be kept at its default, 0. `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 `sumBinomMpfr()`. `rnd.mode` a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see `mpfr`.

## Value

an `"mpfr"` vector of the same length as `q`.

## Note

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.

## Author(s)

Martin Maechler

`pbeta`, `sumBinomMpfr` `chooseZ`.
 ``` 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 = 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)) ) ```