pbetaI: Accurate Incomplete Beta / Beta Probabilities For Integer...

In Rmpfr: Interface R to MPFR - Multiple Precision Floating-Point Reliable

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

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 (from gmp) numbers instead of mpfr numbers.
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

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
   system.time(
   stopifnot(exprs = {
     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),
               tolerance = 1e-230)
     all.equal( Ip, pbetaI(x, a, b, precBits = 2048, lower.tail=FALSE))
   })
   ) # 0.375 sec -- "slow" ???
}

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

cbind(x
    , pb    = rErr(pbeta(x, a, b), pp)
    , pbUp  = rErr(pbeta(x, a, b, lower.tail=FALSE),  Ip)
    , ln.p  = rErr(pbeta(x, a, b, log.p = TRUE    ),  lp)
    , ln.pUp= rErr(pbeta(x, a, b, lower.tail=FALSE, log.p=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., tolerance = 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)) # -5.877e-39
stopifnot(abs(relE) < 1e-28)

Rmpfr documentation built on Nov. 18, 2024, 3:01 p.m.