stirlerr {DPQ} | R Documentation |
Stirling's Error Function - Auxiliary for Gamma, Beta, etc
Description
Stirling's approximation (to the factorial or \Gamma
function)
error in \log
scale is the difference of the left and right hand
side of Stirling's approximation to n!
,
n! \approx \bigl(\frac{n}{e}\bigr)^n \sqrt{2\pi n},
i.e., stirlerr(n) :=
\delta(n)
,
where
\delta(n) = \log\Gamma(n + 1) - n\log(n) + n - \log(2 \pi n)/2.
Partly, pure R transcriptions of the C code utility functions for
dgamma()
, dbinom()
, dpois()
, dt()
,
and similar “base” density functions by Catherine Loader.
These DPQ versions typically have extra arguments with defaults that correspond to R's Mathlib C code hardwired cutoffs and tolerances.
lgammacor(x)
is “the same” as stirlerr(x)
, both
computing delta(x)
accurately, however is only defined for x
\ge 10
, and has been crucially used for R's own lgamma()
and lbeta()
computations.
Usage
stirlerr(n, scheme = c("R3", "R4.4_0"),
cutoffs = switch(scheme
, R3 = c(15, 35, 80, 500)
, R4.4_0 = c(5.25, rep(6.5, 4), 7.1, 7.6, 8.25, 8.8, 9.5, 11,
14, 19, 25, 36, 81, 200, 3700, 17.4e6)
),
use.halves = missing(cutoffs),
direct.ver = c("R3", "lgamma1p", "MM2", "n0"),
order = NA,
verbose = FALSE)
stirlerrC(n, version = c("R3", "R4..1", "R4.4_0"))
stirlerr_simpl(n, version = c("R3", "lgamma1p", "MM2", "n0"), minPrec = 128L)
lgammacor(x, nalgm = 5, xbig = 2^26.5)
Arguments
x , n |
|
verbose |
logical indicating if some information about the computations are to be printed. |
version |
a |
scheme |
a |
cutoffs |
an increasing numeric vector, required to start with
with |
use.halves |
|
direct.ver |
a |
order |
approximation order, |
minPrec |
a positive integer; for |
nalgm |
number of terms to use for Chebyshev polynomial approximation
in |
xbig |
a large positive number; if |
Details
stirlerr()
:Stirling's error,
stirlerr(n):=
\delta(n)
has asymptotic (n \to\infty
) expansion\delta(n) = \frac 1{12 n} - \frac 1{360 n^3} + \frac 1{1260 n^5} \pm O(n^{-7}),
and this expansion is used up to remainder
O(n^{-35})
in current (package DPQ)stirlerr(n)
; different numbers of terms between different cutoffs forn
, and using the direct formula forn <= c_1
, wherec_1
is the first cutoff,cutoff[1]
.Note that (new in 2024-01)
stirlerr(n, order = k)
will not usecutoffs
nor the direct formula (with itsdirect.ver
), nor halves (use.halves=TRUE
), and allowsk \le 20
. Tests seem to indicate that for current double precision arithmetic, onlyk \le 17
seem to make sense.lgammacor(x)
:The “same” Stirling's error, but only defined for
x \ge 10
, returningNaN
otherwise. The example below suggests that R's hardwired default ofnalgm = 5
loses more than one digit accuracy, andnalgm = 6
seems much better. OTOH, the use oflgammacor()
in R's (Mathlib/‘libRmath’) C code is always in conjunction with considerably larger terms such that small inaccuracies inlgammacor()
will not become visible in the values of the functions usinglgammacor()
internally, notablylbeta()
andlgamma()
.
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]
.
Author(s)
Martin Maechler
References
C. Loader (2000), see dbinom
's documentation.
Our package vignette log1pmx, bd0, stirlerr - Probability Computations in R.
See Also
dgamma
,
dpois
.
High precision versions stirlerrM(n)
and
stirlerrSer(n,k)
in package DPQmpfr (via the
Rmpfr and gmp packages).
Examples
n <- seq(1, 50, by=1/4)
st.n <- stirlerr(n) # now vectorized
stopifnot(identical(st.n, sapply(n, stirlerr)))
st3. <- stirlerr(n, "R3", direct.ver = "R3") # previous default
st3 <- stirlerr(n, "R3", direct.ver = "lgamma1p") # new? default
## for these n, there is *NO* difference:
stopifnot(st3 == st3.)
plot(n, st.n, type = "b", log="xy", ylab = "stirlerr(n)")
st4 <- stirlerr(n, "R4.4_0", verbose = TRUE) # verbose: give info on cases
## order = k = 1:20 terms in series approx:
k <- 1:20
stirlOrd <- sapply(k, function(k) stirlerr(n, order = k))
matlines(n, stirlOrd)
matplot(n, stirlOrd - st.n, type = "b", cex=1/2, ylim = c(-1,1)/10, log = "x",
main = substitute(list(stirlerr(n, order=k) ~~"error", k == 1:mK), list(mK = max(k))))
matplot(n, abs(stirlOrd - st.n), type = "b", cex=1/2, log = "xy",
main = "| stirlerr(n, order=k) error |")
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
colnames(stirlOrd) <- paste0("k=", k)
stCn <- stirlerrC(n)
all.equal(st.n, stCn, tolerance = 0) # see 6.7447e-14
stopifnot(all.equal(st.n, stCn, tolerance = 1e-12))
stC2 <- stirlerrC(n, version = "R4..1")
stC4 <- stirlerrC(n, version = "R4.4_0")
## lgammacor(n) : only defined for n >= 10
lgcor <- lgammacor(n)
lgcor6 <- lgammacor(n, nalgm = 6) # more accurate?
all.equal(lgcor[n >= 10], st.n[n >= 10], tolerance=0)# .. rel.diff.: 4.687e-14
stopifnot(identical(is.na(lgcor), n < 10),
all.equal(lgcor[n >= 10],
st.n [n >= 10], tolerance = 1e-12))
## look at *relative* errors -- need "Rmpfr" for "truth" % Rmpfr / DPQmpfr in 'Suggests'
if(requireNamespace("Rmpfr") && requireNamespace("DPQmpfr")) {
## stirlerr(n) uses DPQmpfr::stirlerrM() automagically when n is <mpfr>
relErrV <- sfsmisc::relErrV; eaxis <- sfsmisc::eaxis
mpfr <- Rmpfr::mpfr; asNumeric <- Rmpfr::asNumeric
stM <- stirlerr(mpfr(n, 512))
relE <- asNumeric(relErrV(stM, cbind(st3, st4, stCn, stC4,
lgcor, lgcor6, stirlOrd)))
matplot(n, pmax(abs(relE),1e-20), type="o", cex=1/2, log="xy", ylim =c(8e-17, 0.1),
xaxt="n", yaxt="n", main = quote(abs(relErr(stirlerr(n)))))
## mark "lgcor*" -- lgammacor() particularly !
col.lgc <- adjustcolor(c(2,4), 2/3)
matlines(n, abs(relE[,c("lgcor","lgcor6")]), col=col.lgc, lwd=3)
lines(n, abs(relE[,"lgcor6"]), col=adjustcolor(4, 2/3), lwd=3)
eaxis(1, sub10=2); eaxis(2); abline(h = 2^-(53:51), lty=3, col=adjustcolor(1, 1/2))
axis(1, at=15, col=NA, line=-1); abline(v=c(10,15), lty=2, col=adjustcolor(1, 1/4))
legend("topright", legend=colnames(relE), cex = 3/4,
col=1:6, lty=1:5, pch= c(1L:9L, 0L, letters)[seq_len(ncol(relE))])
legend("topright", legend=colnames(relE)[1:6], cex = 3/4, lty=1:5, lwd=3,
col=c(rep(NA,4), col.lgc), bty="n")
## Note that lgammacor(x) {default, n=5} is clearly inferior,
## but lgammacor(x, 6) is really good in [10, 50]
## ===> Larger n's:
nL <- c(seq(50, 99, by = 1/2), 100*2^seq(0,8, by = 1/4))
stMl <- stirlerr(mpfr(nL, 512))
lgc5 <- lgammacor(nL, nalgm = 5)
lgc6 <- lgammacor(nL, nalgm = 6)
stir7 <- stirlerr(nL, order = 7)
relEl <- asNumeric(relErrV(stMl,
cbind(lgammacor.5 = lgc5, lgammacor.6 = lgc6, 'stirlerr_k=7' = stir7)))
matplot(nL, pmax(abs(relEl),2^-55), type="o", cex=2/3, log="xy",
ylim = c(2^-54.5, max(abs(relEl))), ylab = quote(abs(rE)),
xaxt="n", yaxt="n", main = quote(abs(relErr(stirlerr(n)))))
eaxis(1, sub10=2); eaxis(2, cex.axis=.8)
abline(h = 2^-(54:51), lty=3, col=adjustcolor(1, 1/2))
legend("center", legend=colnames(relEl), lwd=2, pch = paste(1:3), col=1:3, lty=1:3)
}# end if( <Rmpfr> )