DPQ (version 0.4-3)

dgamma-utils: Utility Functions for dgamma() -- Pure R Versions

Description

Mostly, pure R transcriptions of the C code utility functions for dgamma() and similar “base” density functions by Catherine Loader.

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

These have extra arguments with defaults that correspond to the C Mathlib code hardwired cutoffs and tolerances.

Usage

dpois_raw(x, lambda, log=FALSE,
          version, 
          ## the defaults for version will probably change in the future
          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)
bd0(x, np,
    delta = 0.1, maxit = 1000L,
    s0 = .Machine$double.xmin,
    verbose = getOption("verbose"))
bd0C(x, np, delta = 0.1, maxit = 1000L, version = "R4.0", verbose = getOption("verbose"))
ebd0(x, M, verbose = getOption("verbose"))
stirlerr(n, scheme = c("R3", "R4.1"),
         cutoffs = switch(scheme
                        , R3   = c(15, 35, 80, 500)
                        , R4.1 = c(7.5, 8.5, 10.625, 12.125, 20, 26, 55, 200, 3300)
                          ),
         use.halves = missing(cutoffs),
         verbose = FALSE)
lgammacor(x, nalgm = 5, xbig = 2^26.5)

Arguments

x, n

numeric (or number-alike such as "mpfr").

lambda, np, M

each numeric; distrubution 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.

delta, bd0.delta

a positive number, a cutoff for bd0() where the logcf() series expansion is used when \(|x - M| < delta*(x+M)\).

tol_logcf, eps2, minL1, trace.lcf, logCF

optional tuning arguments passed to log1pmx().

maxit

the number of logcf() terms to be used in bd0() when \(|x-M|\) is small.

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

scheme

a character string specifying the cutoffs scheme.

cutoffs

an increasing numeric vector, required to start with with cutoffs[1] <= 15 specifying the cutoffs to switch from 2 to 3 to ..., up to 10 term approximations for non-small n, where the direct formula loses precision. When missing (as by default), scheme is used, where scheme = "R3" chooses (15, 35, 80, 500), the cutoffs in use in R versions up to (and including) 4.0.z.

use.halves

logical indicating if the full-accuracy prestored values should be use when \(2n \in \{0,1,\dots,30\}\), i.e., \(n \le 15\) and n is integer or integer + \(\frac{1}{2}\). Turn this off to judge the underlying approximation accuracy by comparison with MPFR. However, keep the default TRUE for back-compatibility.

nalgm

number of terms to use for Chebyshev polynomial approxmation in lgammacor(). The default, 5, is the value hard wired in R's C Mathlib.

xbig

a large positive number; if x >= xbig, the simple asymptotic approximation lgammacor(x) := 1/(12*x) is used. The default, \(2^{26.5} = 94906265.6\), is the value hard wired in R's C Mathlib.

Value

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

lgammacor(x) originally returned NaN for all \(|x| < 10\), as its Chebyshev polynomial approximation has been constructed for \(x \in [10, xbig]\), specifically for \(u \in [-1,1]\) where \(t := 10/x \in [1/x_B, 1]\) and \(u := 2t^2 -1 \in [-1 + \epsilon_B, 1]\).

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

% item

References

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

See Also

dgamma, dpois.

Examples

# NOT RUN {
n <- seq(1, 50, by=1/4)
st.n <- stirlerr(n) # now vectorized
stopifnot(identical(st.n, sapply(n, stirlerr)))
plot(n, st.n, type = "b", log="xy", ylab = "stirlerr(n)")

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)
stopifnot(all.equal(bd0x1kC, bd0x1k, tol=1e-15)) # even tol=0 currently ..
if(FALSE) ## FIXME !
ebd0x1k <- ebd0(x, 1000) # FIXME: subscript out of bout
if(FALSE) # the bug is  only here:
    ebd0(1001, 1000)
## so this works:
ex. <- ebd0(x[x != 1001], 1000)
## but there is more wrong currently:
lines(x[x!=1001], colSums(ex.), col=3)

x <- 1:80
dp   <- dpois    (x, 48, log=TRUE)# R's 'stats' pkg function
dp.r <- dpois_raw(x, 48, log=TRUE)
all.equal(dp, dp.r, tol = 0) # on Linux 64b, see TRUE
stopifnot(all.equal(dp, dp.r, tol = 1e-14))
# matplot(x, cbind(dp, dp.r), type = "b") # looks "ok", .. but not if you look closely:
# plot(x, dp.r - dp, type = "b",
#      main = "dpois_raw(x, 48, log=T)  -  dpois(..)  -- something's wrong!")
# abline(h=0, lty=3)
# }