dgamma-utils: Binomial Deviance - Auxiliary Functions for 'dgamma()' Etc

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

Binomial Deviance – Auxiliary Functions for dgamma() Etc

Description

The “binomial deviance” function bd0(x, M) := D_0(x,M) := M \cdot d_0(x/M), where d_0(r) := r\log(r) + 1-r \stackrel{!}{=} p_1l_1(r-1). Mostly, pure R transcriptions of the C code utility functions for dgamma(), dbinom(), dpois(), dt(), and similar “base” density functions by Catherine Loader.
These have extra arguments with defaults that correspond to R's Mathlib C code hardwired cutoffs and tolerances.

Usage

dpois_raw(x, lambda, log=FALSE,
          version, 
          ## the defaults for `version` will probably change in the future
          small.x__lambda = .Machine$double.eps,
          bd0.delta = 0.1,
          ## optional arguments of log1pmx() :
          tol_logcf = 1e-14, eps2 = 0.01, minL1 = -0.79149064, trace.lcf = verbose,
          logCF = if (is.numeric(x)) logcf else logcfR,
          verbose = FALSE)

dpois_simpl (x, lambda, log=FALSE)
dpois_simpl0(x, lambda, log=FALSE)

bd0(x, np, delta = 0.1,
    maxit = max(1L, as.integer(-1100 / log2(delta))),
    s0 = .Machine$double.xmin,
    verbose = getOption("verbose"))
bd0C(x, np, delta = 0.1, maxit = 1000L, version = c("R4.0", "R_2025_0510"),
     verbose = getOption("verbose"))
# "simple" log1pmx() based versions :
bd0_p1l1d1(x, M, tol_logcf = 1e-14, ...)
bd0_p1l1d (x, M, tol_logcf = 1e-14, ...)
# p1l1() based version {possibly using log1pmx(), too}:
bd0_p1l1  (x, M, ...)
# using faster formula for large |M-x|/x
bd0_l1pm  (x, M, tol_logcf = 1e-14, ...)

ebd0 (x, M, verbose = getOption("verbose"), ...) # experimental, may disappear !!
ebd0C(x, M, verbose = getOption("verbose"))

Arguments

x	numeric (or number-alike such as "mpfr").
lambda, np, M	each numeric (or number-alike ..); distribution parameters.
log	logical indicating if the log-density should be returned, otherwise the density at x.
verbose	logical indicating if some information about the computations are to be printed.
small.x__lambda	small positive number; for dpois_raw(x, lambda), when x/lambda is not larger than small.x__lambda, the direct log poisson formula is used instead of ebd0() or bd0() and stirlerr().
delta, bd0.delta	a non-negative number \delta < 1, a cutoff for bd0() where a continued fraction series expansion is used when |x - M| \le \delta\cdot(x+M).
tol_logcf, eps2, minL1, trace.lcf, logCF, ...	optional tuning arguments passed to log1pmx(), and to its options passed to logcf().
maxit	the number of series expansion terms to be used in bd0() when |x-M| is small. The default is k such that \delta^{2k} \le 2^{-1022-52}, i.e., will underflow to zero.
s0	the very small s_0 determining that bd0() = s already before the locf series expansion.
version	a character string specifying the version of bd0() to use. The versions available and the current default are partly experimental!

Details

bd0():

Loader's “Binomial Deviance” function; for x, M > 0 (where the limit x \to 0 is allowed). In the case of dbinom, x are integers (and M = n p), but in general x is real.

bd_0(x,M) := M \cdot D_0\bigl(\frac{x}{M}\bigr),

where D_0(u) := u \log(u) + 1-u = u(\log(u) - 1) + 1. Hence

bd_0(x,M) = M \cdot \bigl(\frac{x}{M}(\log(\frac{x}{M}) -1) +1 \bigr) = x \log(\frac{x}{M}) - x + M.

A different way to rewrite this from Martyn Plummer, notably for important situation when \left|x-M \right| \ll M, is using t := (x-M)/M (and \left|t \right| \ll 1 for that situation), equivalently, \frac{x}{M} = 1+t. Using t,

bd_0(x,M) = \log(1+t) - t \cdot M = M \cdot [(t+1)(\log(1+t) - 1) + 1] = M \cdot [(t+1) \log(1+t) - t] = M \cdot p_1l_1(t),

and

p_1l_1(t) := (t+1)\log(1+t) - t = \frac{t^2}{2} - \frac{t^3}{6} ...

where the Taylor series expansion is useful for small |t|, see p1l1.

Note that bd0(x, M) now also works when x and/or M are arbitrary-accurate mpfr-numbers (package Rmpfr).

bd0C() interfaces to C code which corresponds to R's C Mathlib (Rmath) bd0().

Value

a numeric vector “like” x; in some cases may also be an (high accuracy) "mpfr"-number vector, using CRAN package Rmpfr.

ebd0() (R code) and ebd0C() (interface to C code) are experimental, meant to be precision-extended version of bd0(), returning (yh, yl) (high- and low-part of y, the numeric result). In order to work for long vectors x, yh, yl need to be list components; hence we return a two-column data.frame with column names "yh" and "yl".

Author(s)

Martin Maechler

References

C. Loader (2000), see dbinom's documentation.

Martin Mächler (2021 ff) log1pmx, ... Computing ... Probabilities in R. DPQ package vignette https://CRAN.R-project.org/package=DPQ/vignettes/log1pmx-etc.pdf.

See Also

p1l1() and its variants, e.g., p1l1ser(); log1pmx() which is called from bd0_p1l1d*() and bd0_l1pm(). Then, stirlerr for Stirling's error function, complementing bd0() for computation of Gamma, Beta, Binomial and Poisson probabilities. R's own dgamma, dpois.

Examples

x <- 800:1200
bd0x1k <- bd0(x, np = 1000)
plot(x, bd0x1k, type="l", ylab = "bd0(x, np=1000)")
bd0x1kC <- bd0C(x, np = 1000)
lines(x, bd0x1kC, col=2)
bd0.1d1 <- bd0_p1l1d1(x, 1000)
bd0.1d  <- bd0_p1l1d (x, 1000)
bd0.1p1 <- bd0_p1l1  (x, 1000)
bd0.1pm <- bd0_l1pm  (x, 1000)
stopifnot(exprs = {
    all.equal(bd0x1kC, bd0x1k,  tol=1e-14) # even tol=0 currently ..
    all.equal(bd0x1kC, bd0.1d1, tol=1e-14)
    all.equal(bd0x1kC, bd0.1d , tol=1e-14)
    all.equal(bd0x1kC, bd0.1p1, tol=1e-14)
    all.equal(bd0x1kC, bd0.1pm, tol=1e-14)
})

str(log1pmx) #--> play with { tol_logcf, eps2, minL1, trace.lcf, logCF }  --> vignette

ebd0x1k <- ebd0 (x, 1000)
exC     <- ebd0C(x, 1000)
stopifnot(all.equal(exC, ebd0x1k,          tol = 4e-16),
          all.equal(bd0x1kC, rowSums(exC), tol = 2e-15))
lines(x, rowSums(ebd0x1k), col=adjustcolor(4, 1/2), lwd=4)

## very large args --> need new version "R_2025..."
np <- 117e306;  np / .Machine$double.xmax # 0.65
x <- (1:116)*1e306
tail(bd0L <-  bd0C(x, np, version = "R4")) # underflow to 0
tail(bd0Ln <- bd0C(x, np, version = "R_2025_0510"))
bd0L.R <- bd0(x, np)
all.equal(bd0Ln, bd0L.R, tolerance = 0) # see TRUE
stopifnot(exprs = {
    is.finite(bd0Ln)
    bd0Ln > 4e303
    tail(bd0L, 54) == 0
    all.equal(bd0Ln, np*p1l1((x-np)/np), tolerance = 5e-16)
    all.equal(bd0Ln, bd0L.R, tolerance = 5e-16)
})

DPQ documentation built on July 18, 2025, 3:01 p.m.