stirlerrM: Stirling Formula Approximation Error
In DPQmpfr: DPQ (Density, Probability, Quantile) Distribution Computations using MPFR
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
Compute the log() of the error of Stirling's formula for n!.
Used in certain accurate approximations of (negative) binomial and Poisson probabilities.
stirlerrM() currently simply uses the direct mathematical formula,
based on lgamma(), adapted for use with mpfr-numbers.
Usage
1
2
stirlerrM(n, minPrec = 128L)
stirlerrSer(n, k)
Arguments
n
numeric or “numeric-alike” vector, typically
“large” positive integer or half integer valued, here typically an
"mpfr"-number vector.
k
integer scalar, now in 1:11.
minPrec
minimal precision (in bits) to be used when coercing
number-alikes, say, biginteger (bigz) to "mpfr".
Details
Stirling's approximation to n! has been
n! ~= (n/e)^n * sqrt(2*pi*n)
,
where by definition the error is the difference of the left and right
hand side of this formula, in \log-scale,
delta(n) = logΓ(n + 1) - n* log(n) + n - log(2*pi*n)/2.
See the vignette log1pmx, bd0, stirlerr, ... from package
DPQ, where the series expansion of δ(n) is used with
11 terms, starting with
delta(n) = 1/(12 n) - 1/(360 n^3) + 1/(1260 n^5) +/- O(n^{-7}).
Value
a numeric or other “numeric-alike” class vector, e.g.,
mpfr, of the same length as n.
Note
In principle, the direct formula should be replaced by a few terms of the
series in powers of 1/n for large n, but we assume using
high enough precision for n should be sufficient and “easier”.
Author(s)
Martin Maechler
References
Catherine Loader, see dbinom;
Martin Maechler (2021)
log1pmx(), bd0(), stirlerr() – Computing Poisson, Binomial, Gamma Probabilities in R.
https://CRAN.R-project.org/package=DPQ/vignettes/log1pmx-etc.pdf
See Also
dbinom, stirlerr() in package
DPQ which is a pure R version R's mathlib-internal C function.
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
### ---------------- Regular R double precision -------------------------------
n <- n. <- c(1:10, 15, 20, 30, 50*(1:6), 100*(4:9), 10^(3:12))
(stE <- stirlerrM(n)) # direct formula is *not* good when n is large:
require(graphics)
plot(stirlerrM(n) ~ n, log = "x", type = "b", xaxt="n")
sfsmisc::eaxis(1, sub10=3)
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## for larger n, current values are even *negative* ==> dbl prec *not* sufficient
## y in log-scale [same conclusion]
plot (stirlerrM(n) ~ n, log = "xy", type = "b", ylim = c(1e-13, 0.08))
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("topright", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, ncol=2, bty="n")
## the numbers:
options(digits=4, width=111)
stEmat. <- cbind(sM = stirlerrM(n),
sapply(setNames(1:8, paste0("k=",1:8)),
function(k) stirlerrSer(n=n, k=k)))
stEmat.
## for printing n=<nice>:
N <- Rmpfr::asNumeric
dfm <- function(n, mm) data.frame(n=formatC(N(n)), N(mm), check.names=FALSE)
## relative differences:
dfm(n, stEmat.[,-1]/stEmat.[,1] - 1)
# => stirlerrM() {with dbl prec} deteriorates after ~ n = 200--500
dfm(n, stEmat.[,-(1+8)]/stEmat.[,1+8] - 1)
### ---------------- MPFR High Accuracy -------------------------------
stopifnot(require(gmp),
require(Rmpfr))
n <- as.bigz(n.)
## now repeat everything .. from above ... FIXME shows bugs !
## fully accurate using big rational arithmetic
class(stEserQ <- sapply(setNames(1:8, paste0("k=",1:8)),
function(k) stirlerrSer(n=n, k=k))) # list ..
stopifnot(sapply(stEserQ, class) == "bigq") # of exact big rationals
str(stEsQM <- lapply(stEserQ, as, Class="mpfr"))# list of 8; each prec. 128..702
stEsQM. <- lapply(stEserQ, .bigq2mpfr, precB = 512) # constant higher precision
stEsQMm <- sapply(stEserQ, asNumeric) # a matrix
stEM <- stirlerrM(mpfr(n, 128)) # now ok (loss of precision, but still ~ 10 digits correct)
stEM4k <- stirlerrM(mpfr(n, 4096))# assume "perfect"
## ==> what's the accuracy of the 128-bit 'stEM'?
N <- asNumeric # short
dfm(n, stEM/stEM4k - 1)
## 29 1e+06 4.470e-25
## 30 1e+07 -7.405e-23
## 31 1e+08 -4.661e-21
## 32 1e+09 -7.693e-20
## 33 1e+10 3.452e-17 (still ok)
## 34 1e+11 -3.472e-15 << now start losing
## 35 1e+12 -3.138e-13 <<<<
## same conclusion via number of correct (decimal) digits:
dfm(n, log10(abs(stEM/stEM4k - 1)))
plot(N(-log10(abs(stEM/stEM4k - 1))) ~ N(n), type="o", log="x",
xlab = quote(n), main = "#{correct digits} of 128-bit stirlerrM(n)")
ubits <- c(128, 52) # above 128-bit and double precision
abline(h = ubits* log10(2), lty=2)
text(1, ubits* log10(2), paste0(ubits,"-bit"), adj=c(0,0))
stopifnot(identical(stirlerrM(n), stEM)) # for bigz & bigq, we default to precBits = 128
all.equal(roundMpfr(stEM4k, 64),
stirlerrSer (n., 8)) # 0.00212 .. because of 1st few n. ==> drop these
all.equal(roundMpfr(stEM4k,64)[n. >= 3], stirlerrSer (n.[n. >= 3], 8)) # 6.238e-8
plot(asNumeric(abs(stirlerrSer(n., 8) - stEM4k)) ~ n.,
log="xy", type="b", main="absolute error of stirlerrSer(n, 8) & (n, 5)")
abline(h = 2^-52, lty=2); text(1, 2^-52, "52-bits", adj=c(1,-1)/8)
lines(asNumeric(abs(stirlerrSer(n., 5) - stEM4k)) ~ n., col=2)
plot(asNumeric(stirlerrM(n)) ~ n., log = "x", type = "b")
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## y in log-scale
plot(asNumeric(stirlerrM(n)) ~ n., log = "xy", type = "b", ylim = c(1e-13, 0.08))
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## all "looks" perfect (so we could skip this)
## the numbers ...
## %% FIXME a list instead of mpfrMatrix ... FIXME _____________
## FIXME ... asNumeric() needed or as(*, "mpfr") or ...
ks <- 1:8 ## k <= 5 === FIXME --- use DPQ's version !!
stirlS.l <- lapply(setNames(ks, paste0("k=",ks)),
function(k) stirlerrSer(n=n, k=k))
## ==> an mpfrMatrix of dim 35 x 5 :
mss <- do.call(cbind, lapply(stirlS.l, mpfr, precBits=256))
stEmat <- cbind(sM = stEM4k, mss)
signif(asNumeric(stEmat), 6) # so it prints nicely
## print *relative errors* nicely :
## simple double precision version of direct formula (cancellation for n >> 1 !):
stE <- stirlerrM(n.)
dfm(n , cbind(stEmat[,-1], dbl=stE)/stEM4k - 1)
## relative differences:
dfm(n, stEmat[,-1] / stEmat[,1] - 1)
dfm(n., stEmat[,-(1+8)]/ stEmat[,1+8] - 1)
DPQmpfr documentation built on May 18, 2021, 1:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
Description
Compute the log() of the error of Stirling's formula for n!.
Used in certain accurate approximations of (negative) binomial and Poisson probabilities.
stirlerrM() currently simply uses the direct mathematical formula,
based on lgamma(), adapted for use with mpfr-numbers.
Usage
1 2 | stirlerrM(n, minPrec = 128L)
stirlerrSer(n, k)
|
Arguments
n |
numeric or “numeric-alike” vector, typically
“large” positive integer or half integer valued, here typically an
|
k |
integer scalar, now in |
minPrec |
minimal precision (in bits) to be used when coercing
number-alikes, say, biginteger ( |
Details
Stirling's approximation to n! has been
n! ~= (n/e)^n * sqrt(2*pi*n)
, where by definition the error is the difference of the left and right hand side of this formula, in \log-scale,
delta(n) = logΓ(n + 1) - n* log(n) + n - log(2*pi*n)/2.
See the vignette log1pmx, bd0, stirlerr, ... from package DPQ, where the series expansion of δ(n) is used with 11 terms, starting with
delta(n) = 1/(12 n) - 1/(360 n^3) + 1/(1260 n^5) +/- O(n^{-7}).
Value
a numeric or other “numeric-alike” class vector, e.g.,
mpfr, of the same length as n.
Note
In principle, the direct formula should be replaced by a few terms of the
series in powers of 1/n for large n, but we assume using
high enough precision for n should be sufficient and “easier”.
Author(s)
Martin Maechler
References
Catherine Loader, see dbinom;
Martin Maechler (2021) log1pmx(), bd0(), stirlerr() – Computing Poisson, Binomial, Gamma Probabilities in R. https://CRAN.R-project.org/package=DPQ/vignettes/log1pmx-etc.pdf
See Also
dbinom, stirlerr() in package
DPQ which is a pure R version R's mathlib-internal C function.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | ### ---------------- Regular R double precision -------------------------------
n <- n. <- c(1:10, 15, 20, 30, 50*(1:6), 100*(4:9), 10^(3:12))
(stE <- stirlerrM(n)) # direct formula is *not* good when n is large:
require(graphics)
plot(stirlerrM(n) ~ n, log = "x", type = "b", xaxt="n")
sfsmisc::eaxis(1, sub10=3)
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## for larger n, current values are even *negative* ==> dbl prec *not* sufficient
## y in log-scale [same conclusion]
plot (stirlerrM(n) ~ n, log = "xy", type = "b", ylim = c(1e-13, 0.08))
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("topright", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, ncol=2, bty="n")
## the numbers:
options(digits=4, width=111)
stEmat. <- cbind(sM = stirlerrM(n),
sapply(setNames(1:8, paste0("k=",1:8)),
function(k) stirlerrSer(n=n, k=k)))
stEmat.
## for printing n=<nice>:
N <- Rmpfr::asNumeric
dfm <- function(n, mm) data.frame(n=formatC(N(n)), N(mm), check.names=FALSE)
## relative differences:
dfm(n, stEmat.[,-1]/stEmat.[,1] - 1)
# => stirlerrM() {with dbl prec} deteriorates after ~ n = 200--500
dfm(n, stEmat.[,-(1+8)]/stEmat.[,1+8] - 1)
### ---------------- MPFR High Accuracy -------------------------------
stopifnot(require(gmp),
require(Rmpfr))
n <- as.bigz(n.)
## now repeat everything .. from above ... FIXME shows bugs !
## fully accurate using big rational arithmetic
class(stEserQ <- sapply(setNames(1:8, paste0("k=",1:8)),
function(k) stirlerrSer(n=n, k=k))) # list ..
stopifnot(sapply(stEserQ, class) == "bigq") # of exact big rationals
str(stEsQM <- lapply(stEserQ, as, Class="mpfr"))# list of 8; each prec. 128..702
stEsQM. <- lapply(stEserQ, .bigq2mpfr, precB = 512) # constant higher precision
stEsQMm <- sapply(stEserQ, asNumeric) # a matrix
stEM <- stirlerrM(mpfr(n, 128)) # now ok (loss of precision, but still ~ 10 digits correct)
stEM4k <- stirlerrM(mpfr(n, 4096))# assume "perfect"
## ==> what's the accuracy of the 128-bit 'stEM'?
N <- asNumeric # short
dfm(n, stEM/stEM4k - 1)
## 29 1e+06 4.470e-25
## 30 1e+07 -7.405e-23
## 31 1e+08 -4.661e-21
## 32 1e+09 -7.693e-20
## 33 1e+10 3.452e-17 (still ok)
## 34 1e+11 -3.472e-15 << now start losing
## 35 1e+12 -3.138e-13 <<<<
## same conclusion via number of correct (decimal) digits:
dfm(n, log10(abs(stEM/stEM4k - 1)))
plot(N(-log10(abs(stEM/stEM4k - 1))) ~ N(n), type="o", log="x",
xlab = quote(n), main = "#{correct digits} of 128-bit stirlerrM(n)")
ubits <- c(128, 52) # above 128-bit and double precision
abline(h = ubits* log10(2), lty=2)
text(1, ubits* log10(2), paste0(ubits,"-bit"), adj=c(0,0))
stopifnot(identical(stirlerrM(n), stEM)) # for bigz & bigq, we default to precBits = 128
all.equal(roundMpfr(stEM4k, 64),
stirlerrSer (n., 8)) # 0.00212 .. because of 1st few n. ==> drop these
all.equal(roundMpfr(stEM4k,64)[n. >= 3], stirlerrSer (n.[n. >= 3], 8)) # 6.238e-8
plot(asNumeric(abs(stirlerrSer(n., 8) - stEM4k)) ~ n.,
log="xy", type="b", main="absolute error of stirlerrSer(n, 8) & (n, 5)")
abline(h = 2^-52, lty=2); text(1, 2^-52, "52-bits", adj=c(1,-1)/8)
lines(asNumeric(abs(stirlerrSer(n., 5) - stEM4k)) ~ n., col=2)
plot(asNumeric(stirlerrM(n)) ~ n., log = "x", type = "b")
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## y in log-scale
plot(asNumeric(stirlerrM(n)) ~ n., log = "xy", type = "b", ylim = c(1e-13, 0.08))
for(k in 1:8) lines(n, stirlerrSer(n, k), col = k+1)
legend("top", c("stirlerrM(n)", paste0("stirlerrSer(n, k=", 1:8, ")")),
pch=c(1,rep(NA,8)), col=1:(8+1), lty=1, bty="n")
## all "looks" perfect (so we could skip this)
## the numbers ...
## %% FIXME a list instead of mpfrMatrix ... FIXME _____________
## FIXME ... asNumeric() needed or as(*, "mpfr") or ...
ks <- 1:8 ## k <= 5 === FIXME --- use DPQ's version !!
stirlS.l <- lapply(setNames(ks, paste0("k=",ks)),
function(k) stirlerrSer(n=n, k=k))
## ==> an mpfrMatrix of dim 35 x 5 :
mss <- do.call(cbind, lapply(stirlS.l, mpfr, precBits=256))
stEmat <- cbind(sM = stEM4k, mss)
signif(asNumeric(stEmat), 6) # so it prints nicely
## print *relative errors* nicely :
## simple double precision version of direct formula (cancellation for n >> 1 !):
stE <- stirlerrM(n.)
dfm(n , cbind(stEmat[,-1], dbl=stE)/stEM4k - 1)
## relative differences:
dfm(n, stEmat[,-1] / stEmat[,1] - 1)
dfm(n., stEmat[,-(1+8)]/ stEmat[,1+8] - 1)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.