dbinom_raw: R's C Mathlib (Rmath) dbinom_raw() Binomial Probability pure...

In DPQ: Density, Probability, Quantile ('DPQ') Computations

R's C Mathlib (Rmath) dbinom_raw() Binomial Probability pure R Function

Description

A pure R implementation of R's C API (‘Mathlib’ specifically) dbinom_raw() function which computes binomial probabilities and is continuous in x, i.e., also “works” for non-integer x.

Usage

dbinom_raw (x, n, p, q = 1-p, log = FALSE,
            version = c("2008", "R4.4"),
            verbose = getOption("verbose"))

Arguments

x	vector with values typically in 0:n, but here allowed to non-integer values.
n	called size in R's dbinom().
p	called prob in R's dbinom(), the success probability, hence in [0, 1].
q	mathemtically the same as 1 - p, but may be (much) more accurate, notably when small.
log	logical indicating if the log() of the resulting probability should be returned; useful notably in case the probability itself would underflow to zero.
version	a character string; originally, "2008" was the only option. Still the default currently, this may change in the future.
verbose	integer indicating the amount of verbosity of diagnostic output, 0 means no output, 1 more, etc.

Value

numeric vector of the same length as x which may have to be thought of recycled along n, p and/or q.

Author(s)

R Core and Martin Maechler

See Also

Note that our CRAN package Rmpfr provides dbinom, an mpfr-accurate function to be used used instead of R's or this pure R version relying bd0() and stirlerr() where the latter currently only provides accurate double precision accuracy.

Examples

for(n in c(3, 10, 27, 100, 500, 2000, 5000, 1e4, 1e7, 1e10)) {
 x <- if(n <= 2000) 0:n else round(seq(0, n, length.out=2000))
 p <- 3/4
 stopifnot(all.equal(dbinom_raw(x, n, p, q=1-p) -> dbin,
                     dbinom    (x, n, p), tolerance = 1e-13))# 1.636e-14 (Apple clang 14.0.3)
 stopifnot(all.equal(dbin, dbinom_raw(x, n, p, q=1-p, version = "R4.4") -> dbin44,
                     tolerance = 1e-13))
 cat("n = ", n, ": ", (aeq <- all.equal(dbin44, dbin, tolerance = 0)), "\n")
 if(n < 3000) stopifnot(is.character(aeq)) # check that dbin44 is "better" ?!
}

n <- 1024 ; x <- 0:n
plot(x, dbinom_raw(x, n, p, q=1-p) - dbinom(x, n, p), type="l", main = "|db_r(x) - db(x)|")
plot(x, dbinom_raw(x, n, p, q=1-p) / dbinom(x, n, p) - 1, type="b", log="y",
     main = "rel.err.  |db_r(x / db(x) - 1)|")

DPQ documentation built on Sept. 11, 2024, 8:37 p.m.