betaD94 {DPQmpfr}

Ding(1994) (non-central) Beta Distribution Functions

Description

The three functions "p" (cumulative distribution, CDF), "d" (density (PDF)), and "q" (quantile) use Ding(1994)'s algorithm A, B, and C, respectively, each of which implements a recursion formula using only simple arithmetic and log and exp.

These are particularly useful also for using with high precision "mpfr" numbers from the Rmpfr CRAN package.

Usage

dbetaD94(x, shape1, shape2, ncp = 0, log = FALSE,
         eps = 1e-10, itrmax = 100000L, verbose = FALSE)
pbetaD94(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE,
	 log_scale = (a * b > 0) && (a + b > 100 || c >= 500),
         eps = 1e-10, itrmax = 100000L, verbose = FALSE)

qbetaD94(p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE,
	 log_scale = (a * b > 0) && (a + b > 100 || c >= 500),
         delta = 1e-6,
         eps = delta^2,
         itrmax = 100000L,
         iterN = 1000L,
         verbose = FALSE)

Arguments

x, q

numeric vector of values in [0,1] as beta variates.

shape1, shape2

the two shape parameters of the beta distribution, must be positive.

ncp

the noncentrality parameter; by default zero for the (central) beta distribution; if positive, we have a noncentral beta distribution.

p

numeric vector of probabilities, log()ged in case log.p is true.

log, log.p

logical indicating if the density or probability values should be log()ged.

lower.tail

logical indicating if the lower or upper tail probability should be computed, or for qbeta*() are provided.

eps

a non-negative number specifying the desired accuracy for computing F() and f().

itrmax

the maximal number of steps for computing F() and f().

delta

[For qbeta*():] non-negative number indicating the desired accuracy for computing x_p (the root of pbeta*() == p), i.e., the convergence tolerance for the Newton iterations. This sets default eps = delta^2 which is sensible but may be too small, such that eps should be specified in addition to delta.

iterN

[For qbeta*():] The maximal number of Newton iterations.

log_scale

logical indicating if most of the computations should happen in log scale, which protects from “early” overflow and underflow but takes more computations. The current default is somewhat arbitrary, still derived from the facts that gamma(172) overflows to Inf already and exp(-750) underflows to 0 already.

verbose

logical (or integer) indicating the amount of diagnostic output during computation; by default none.

Value

In all three cases, a numeric vector with the same attributes as x (or q respectively), containing (an approximation) to the correponding beta distribution function.

Author(s)

Martin Maechler, notably log_scale was not part of Ding's proposals.

References

Cherng G. Ding (1994) On the computation of the noncentral beta distribution. Computational Statistics & Data Analysis 18, 449–455.

See Also

pbeta. Package Rmpfr's pbetaI() needs both shape1 and shape2 to be integer but is typically more efficient than the current pbetaD94() implementation.

Examples

          
## Low precision (eps, delta) values as "e.g." in Ding(94): ------------------

## Compare with  Table 3  of  Baharev_et_al 2017 %% ===> ./qbBaha2017.Rd <<<<<<<<<<<
aa <- c(0.5, 1, 1.5, 2, 2.5, 3, 5, 10, 25)
bb <- c(1:15, 10*c(2:5, 10, 25, 50))
utime <-
 qbet <- matrix(NA_real_, length(aa), length(bb),
                dimnames = list(a = formatC(aa), b = formatC(bb)))
(doExtras <- DPQmpfr:::doExtras())
if(doExtras) qbetL <- utimeL <- utime

p <- 0.95
delta <- 1e-4
eps   <- 1e-6
system.t.usr <- function(expr)
  system.time(gcFirst = FALSE, expr)[["user.self"]]

system.time(
for(ia in seq_along(aa)) {
    a <- aa[ia]; cat("\n--==--\na=",a,":\n")
    for(ib in seq_along(bb)) {
        b <- bb[ib]; cat("\n>> b=",b,"\n")
        utime [ia, ib] <- system.t.usr(
          qbet[ia, ib] <-   qbetaD94(p, a, b, ncp = 0, delta=delta, eps=eps, verbose = 2))
        if(doExtras)
          utimeL[ia, ib] <- system.t.usr(
           qbetL[ia, ib] <-   qbetaD94(p, a, b, ncp = 0, delta=delta, eps=eps,
                                       verbose = 2, log_scale=TRUE))
    }
    cat("\n")
}
)# system.time(.): ~ 1 sec (lynne i7-7700T, Fedora 32, 2020)
sum(print(table(round(1000*utime)))) # lynne .. :
##  0  1  2  3  4  5  6  7  8  9 10 11 14 15 16 29
## 53 94 15  3  3 12  2  2  2  2  1  2  3  1  2  1
## [1] 198
if(doExtras) print(sum(print(table(round(1000*utimeL))))) # lynne .. :