Rmpfr (version 0.9-5)

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.

TODO (not yet): It's sufficient for one of \(a,b\) to be integer such that the result is a finite sum (but the coefficients will no longer be rational, see Abramowitz and Stegun, 26.5.6 and *.7, p.944).

Usage

pbetaI(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE,
       precBits = NULL,
       useRational = !log.p && !is.mpfr(q) && is.null(precBits) && int2,
       rnd.mode = c("N","D","U","Z","A"))

Value

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

Arguments

q

called \(x\), above; vector of quantiles, in \([0,1]\); can be numeric, or of class "mpfr" or also "bigq" (“big rational” from package gmp); in the latter case, if log.p = FALSE as by default, all computations are exact, using big rational arithmetic.

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

useRational

optional logical, specifying if we should try to do everything in exact rational arithmetic, i.e, using package gmp functionality only, and return bigq numbers instead of mpfr numbers.

rnd.mode

a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Author

Martin Maechler

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.

See Also

pbeta, sumBinomMpfr chooseZ.

Examples

x <- (0:12)/16 # not all the way up ..
a <- 7; b <- 788

p.  <- pbetaI(x, a, b) ## a bit slower:
system.time(
pp  <- pbetaI(x, a, b, precBits = 2048)
) # 0.23 -- 0.50 sec
## 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))
 )

## When 'q' is a  bigrational (i.e., class "bigq", package 'gmp'), everything
## is computed *exactly* with bigrational arithmetic:
(q4 <- as.bigq(1, 2^(0:4)))
pb4 <- pbetaI(q4, 10, 288, lower.tail=FALSE)
stopifnot( is.bigq(pb4) )
mpb4 <- as(pb4, "mpfr")
mpb4[1:2]
getPrec(mpb4) # 128 349 1100 1746 2362
(pb. <- pbeta(asNumeric(q4), 10, 288, lower.tail=FALSE))
stopifnot(mpb4[1] == 0,
          all.equal(mpb4, pb., tol=4e-15))

qbetaI. <- function(p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE,
    precBits = NULL, rnd.mode = c("N", "D", "U", "Z", "A"),
    tolerance = 1e-20, ...)
{
    if(is.na(a <- as.integer(shape1))) stop("a = shape1 is not coercable to finite integer")
    if(is.na(b <- as.integer(shape2))) stop("b = shape2 is not coercable to finite integer")
    unirootR(function(q) pbetaI(q, a, b, lower.tail=lower.tail, log.p=log.p,
                                precBits=precBits, rnd.mode=rnd.mode) - p,
             interval = if(log.p) c(-double.xmax, 0) else 0:1,
             tol = tolerance, ...)
} # end{qbetaI}

(p <- 1 - mpfr(1,128)/20) # 'p' must be high precision
q95.1.3 <- qbetaI.(p, 1,3, tolerance = 1e-29) # -> ~29 digits accuracy
str(q95.1.3) ; roundMpfr(q95.1.3$root, precBits = 29 * log2(10))
## relative error is really small:
(relE <- asNumeric(1 - pbetaI(q95.1.3$root, 1,3) / p))
stopifnot(abs(relE) < 1e-28)