dgamma-utils {DPQ} | R Documentation |
dgamma()
– Pure R VersionsMostly, 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.
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)
x, n |
|
lambda, np, M |
each |
log |
logical indicating if the log-density should be returned,
otherwise the density at |
verbose |
logical indicating if some information about the computations are to be printed. |
delta, bd0.delta |
a positive number, a cutoff for |
tol_logcf, eps2, minL1, trace.lcf, logCF |
optional tuning arguments passed to |
maxit |
the number of |
s0 |
the very small s_0 determining that |
version |
a |
scheme |
a |
cutoffs |
an increasing numeric vector, required to start with
with |
use.halves |
|
nalgm |
number of terms to use for Chebyshev polynomial approxmation
in |
xbig |
a large positive number; if |
bd0()
:Loader's “Binomial Deviance” function; for
x, M > 0 (where the limit x -> 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 * D_0(x/M),
where D_0(u) := u \log(u) + 1-u = u(\log(u) - 1) + 1. Hence
bd_0(x,M) = M *((x/M)*(log(x/M) - 1) +1) = x log(x/M) - x + M.
A different way to rewrite this from Martyn Plummer, notably for important situation when |x-M| << M, is using t := (x-M)/M (and |t| << 1 for that situation), equivalently, x/M = 1+t. Using t,
bd_0(x,M) = log(1+t) - t * M = M * ((t+1)(log(1+t) - 1) + 1) = M * ((t+1) log(1+t) - t) = M * p1l1(t),
and
p1l1 (t) := (t+1)*log(1+t) - t = t^2/2 - t^3/6 ...
where the Taylor series expansion is useful for small |t|.
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 + ε_B, 1].
Martin Maechler
C. Loader (2000), see dbinom
's documentation.
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)