lgamma1p {DPQ} | R Documentation |
log(gamma(a+1))
Compute
lG[1](a) := log(Gamma(a+1)) = log(a * Gamma(a)) = log(a) + log(Gamma(a)),
which is “in principle” the same as
log(gamma(a+1))
or lgamma(a+1)
,
accurately also for (very) small a (0 < a < 0.5).
lgamma1p (a, tol_logcf = 1e-14) lgamma1p.(a, cutoff.a = 1e-6, k = 3) lgamma1p_series(x, k)
a, x |
a numeric vector. |
tol_logcf |
for |
cutoff.a |
for |
k |
an integer, the number of terms in the series expansion used internally. |
lgamma1p()
is an R translation of the function (in Fortran) in
Didonato and Morris (1992) which uses a 40-degree polynomial approximation.
lgamma1p_series(x, k)
is Taylor series approximation of order k
,
(derived via Maple), which is gamma*x + pi^2 * x^2/ 12 + O(x^3), where gamma
is Euler's constant 0.5772156649....
a numeric vector with the same attributes as a
.
Morten Welinder (C code of Jan 2005, see R's bug issue
PR#7307) for lgamma1p()
.
Martin Maechler, notably for lgamma1p_series()
which works
with package Rmpfr but otherwise may be much less
accurate than Morten's 40 term series!
Didonato, A. and Morris, A., Jr, (1992)
Algorithm 708: Significant digit computation of the incomplete beta function ratios.
ACM Transactions on Mathematical Software, 18, 360–373;
see also pbeta
.
curve(-log(x*gamma(x)), 1e-30, .8, log="xy", col="gray50", lwd = 3, axes = FALSE, ylim = c(1e-30,1)) sfsmisc::eaxis(1); sfsmisc::eaxis(2) at <- 10^(1-4*(0:8)) abline(h = at, v = at, col = "lightgray", lty = "dotted") curve(-lgamma( 1+x), add=TRUE, col="red2", lwd=1/2)# underflows even earlier curve(-lgamma1p (x), add=TRUE, col="blue") curve(-lgamma1p.(x), add=TRUE, col=adjustcolor("forest green",1/4), lwd = 5, lty = 2) for(k in 1:7) curve(-lgamma1p_series(x, k=k), add=TRUE, col=paste0("gray",30+k*8), lty = 3)