lgamma1p {DPQ} | R Documentation |
Accurate log(gamma(a+1))
Description
Compute
l\Gamma_1(a) := \log\Gamma(a+1) = \log(a\cdot \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)
.
Usage
lgamma1p (a, tol_logcf = 1e-14, f.tol = 1, ...)
lgamma1p.(a, cutoff.a = 1e-6, k = 3)
lgamma1p_series(x, k)
lgamma1pC(x)
Arguments
a , x |
a numeric vector. |
tol_logcf |
for |
f.tol |
numeric (factor) used in
|
... |
further optional arguments passed on to |
cutoff.a |
for |
k |
an integer, the number of terms in the series expansion used internally. |
Details
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. ...
lgamma1pC()
is an interface to R C API (‘Mathlib’ / ‘Rmath.h’) function.
Value
a numeric vector with the same attributes as a
.
Author(s)
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!
References
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
.
See Also
Examples
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") -> lgxy
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)
stopifnot(with(lgxy, all.equal(y, -lgamma1pC(x))))