Rmpfr (version 0.8-1)

pbetaI: Accurate Incomplete Beta / Beta Probabilities For Integer Shapes

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

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 \le 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.

See Also

pbeta, sumBinomMpfr chooseZ.

Examples

# NOT RUN {
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))
 )
# }