DPQ (version 0.6-1)

pbetaRv1: Pure R Implementation of Old pbeta()

Description

pbetaRv1() is an implementation of the original (“version 1” pbeta() function in R (versions <= 2.2.x), before we started using TOMS 708 bratio() instead, see the current pbeta help page also for references.

pbetaRv1() is basically a manual translation from C to R of the underlying pbeta_raw() C function, see in R's source tree at https://svn.r-project.org/R/branches/R-2-2-patches/src/nmath/pbeta.c

For consistency within R, we are using R's argument names (q, shape1, shape2) instead of C code's (x, pin, qin ).

It is only for the central beta distribution.

Usage

pbetaRv1(q, shape1, shape2, lower.tail = TRUE,
         eps = 0.5 * .Machine$double.eps,
         sml = .Machine$double.xmin,
         verbose = 0)

Value

a number.

Arguments

q, shape1, shape2

non-negative numbers, q in \([0,1]\), see pbeta.

lower.tail

indicating if \(F(q; *)\) should be returned or the upper tail probability \(1 - F(q)\).

eps

the tolerance used to determine congerence. eps has been hard coded in C code to 0.5 * .Machine$double.eps which is equal to \(2^{-53}\) or 1.110223e-16.

sml

the smallest positive number on the typical platform. The default .Machine$double.xmin is hard coded in the C code (as DBL_MIN), and this is equal to \(2^{-1022}\) or 2.225074e-308 on all current platforms.

verbose

integer indicating the amount of verbosity of diagnostic output, 0 means no output, 1 more, etc.

Author

Martin Maechler

References

(From the C code:)

Nancy E. Bosten and E.L. Battiste (1974). Remark on Algorithm 179 (S14): Incomplete Beta Ratio. Communications of the ACM, 17(3), 156--7.

See Also

pbeta; our approximations to the non-central beta pnbetaAppr2.

Examples

all.equal(pbetaRv1(1/4, 2, 3),
          pbeta   (1/4, 2, 3))
set.seed(101)

N <- 1000
x <- sample.int(7, N, replace=TRUE) / 8
a <-   rlnorm(N)
b <- 5*rlnorm(N)
pbt <- pbeta(x, a, b)
for(i in 1:N) {
   stopifnot(all.equal(pbetaRv1(x[i], a[i], b[i]), pbt[i]))
   cat(".", if(i %% 20 == 0) paste0(i, "\n")) #%
}

if(requireNamespace("Rmpfr"))
    pbetaRv1(.1, 1e-3, Rmpfr::mpfr(1e7, 128), sml = 2^Rmpfr::mpfr(-1e20, 128),
             lower.tail = FALSE, verbose=2)