dgamma-utils: 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"))

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

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!

Author

Martin Maechler

Details

% ../vignettes/log1pmx-etc.Rnw <<<<<<<<<

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().

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.

Examples

Run this code

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

Run the code above in your browser using DataLab