Nothing
tests/stirlerr-tst.R
In DPQ: Density, Probability, Quantile ('DPQ') Computations
#### Testing stirlerr()
#### =================== {previous 2nd part of this, now -->>> ./bd0-tst.R <<<---
require(DPQ)
for(pkg in c("Rmpfr", "DPQmpfr"))
if(!requireNamespace(pkg)) {
cat("no CRAN package", sQuote(pkg), " ---> no tests here.\n")
q("no")
}
require("Rmpfr")
source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE))
## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS()
##_ options(conflicts.policy = list(depends.ok=TRUE, error=FALSE, warn=FALSE))
require(sfsmisc) # masking 'list_' *and* gmp's factorize(), is.whole()
##_ options(conflicts.policy = NULL)o
## plot1cuts() , etc: ---> ../inst/extraR/relErr-plots.R <<<<<<<
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
source(system.file(package="DPQ", "extraR", "relErr-plots.R", mustWork=TRUE))
do.pdf <- TRUE # (manually)
do.pdf <- !dev.interactive(orNone = TRUE)
do.pdf
if(do.pdf) {
pdf.options(width = 9, height = 6.5) # for all pdf plots {9/6.5 = 1.38 ; 4/3 < 1.38 < sqrt(2) [A4]
pdf("stirlerr-tst.pdf")
}
showProc.time()
(doExtras <- DPQ:::doExtras())
(noLdbl <- (.Machine$sizeof.longdouble <= 8)) ## TRUE when --disable-long-double or 'M1mac' ..
(M.mac <- grepl("aarch64-apple", R.version$platform)) # Mac with M1, M2, ... proc
abs19 <- function(r) pmax(abs(r), 1e-19) # cut |err| to positive {for log-plots}
options(width = 100, nwarnings = 1e5, warnPartialMatchArgs = FALSE)
##=== Really, dpois_raw() and dbinom_raw() *both* use stirlerr(x) for "all" 'x > 0'
## ~~~~~~~~~~~ ~~~~~~~~~~~~ =========== ===== !
## below, 6 "it's okay, but *far* from perfect:" ===> need more terms in stirlerr() [
## April 20: MM added more terms up to S10; 2024-01: up to S12 ..helps a little only
x <- lseq(1/16, 6, length=2048)
system.time(stM <- DPQmpfr::stirlerrM(Rmpfr::mpfr(x,2048))) # 1.7 sec elapsed
plot(x, stirlerr(x, use.halves=FALSE) - stM, type="l", log="x", main="absolute Error")
plot(x, stirlerr(x, use.halves=FALSE) / stM - 1, type="l", log="x", main="relative Error")
plot(x, abs(stirlerr(x, use.halves=FALSE) / stM - 1), type="l", log="xy",main="|relative Error|")
drawEps.h(-(52:50))
## lgammacor() does *NOT* help, as it is *designed* for x >= 10 ... (but is interesting there)
##
## ==> Need another chebyshev() or rational-approx. for x in [.1, 7] or so !!
##=============> see also ../Misc/stirlerr-trms.R <===============
## ~~~~~~~~~~~~~~~~~~~~~~~
cutoffs <- c(15,35,80,500) # cut points, n=*, in the stirlerr() "algorithm"
##
n <- c(seq(1,15, by=1/4),seq(16, 25, by=1/2), 26:30, seq(32,50, by=2), seq(55,1000, by=5),
20*c(51:99), 50*(40:80), 150*(27:48), 500*(15:20))
st.n <- stirlerr(n, "R3")# rather use.halves=TRUE; but here , use.halves=FALSE
plot(st.n ~ n, log="xy", type="b") ## looks good now (straight line descending {in log-log !}
nM <- mpfr(n, 2048)
st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose
all.equal(asNumeric(st.nM), st.n)# TRUE
all.equal(st.nM, as(st.n,"mpfr"))# .. difference: 3.381400e-14 was 1.05884.........e-15
all.equal(roundMpfr(st.nM, 64), as(st.n,"mpfr"), tolerance=1e-16)# (ditto)
## --- Look at the direct formula -- why is it not good for n ~= 5 ?
##
## Preliminary Conclusions :
## 1. there is *some* cancellation even for small n (how much?)
## 2. lgamma1p(n) does really not help much compared to lgamma(n+1) --- but a tiny bit in some cases
### 1. Investigating lgamma1p(n) vs lgamma(n+1) for n < 1 =============================================
##' @title Relative Error of lgamma(n+1) vs lgamma1p() vs MM's stirlerrD2():
##' @param n numeric, typically n << 1
##' @param precBits
##' @return relative error WRT mpfr(n, precBits)
##' @author Martin Maechler
relE.lgam1 <- function(n, precBits = if(doExtras) 1024 else 320) {
M_LN2PI <- 1.837877066409345483 # ~ log(2*Const("pi",60)); very slightly more accurate than log(2*pi)
st <- lgamma(n +1) - (n +0.5)*log(n) + n - M_LN2PI/2
st. <- lgamma1p(n) - (n +0.5)*log(n) + n - M_LN2PI/2 # "lgamma1p"
st2 <- lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - M_LN2PI)/2 # "MM2"
st0 <- -(l.n + M_LN2PI)/2 # "n0"
nM <- mpfr(n, precBits)
stM <- lgamma(nM+1) - (nM+0.5)*log(nM) + nM - log(2*Const("pi", precBits))/2
## stM <- roundMpfr(stM, 128)
cbind("R3" = asNumeric(relErrV(stM, st))
, "lgamma1p"= asNumeric(relErrV(stM, st.))
, "MM2" = asNumeric(relErrV(stM, st2))
, "n0" = asNumeric(relErrV(stM, st0))
)
}
n <- 2^-seq.int(1022, 1, by = -1/4)
relEx <- relE.lgam1(n)
showProc.time()
## Is *equivalent* to 'new' stirlerr_simpl(n version = *) [not for <mpfr> though, see 'relEmat']:
(simpVer <- eval(formals(stirlerr_simpl)$version))
if(is.null(simpVer)) { warning("got wrong old version of package 'DPQ':")
print(packageDescription("DPQ"))
stop("invalid outdated version package 'DPQ'")
}
stir.allS <- function(n) sapply(simpVer, function(v) stirlerr_simpl(n, version=v))
stirS <- stir.allS(n) # matrix
nM <- mpfr(n, 256) # "high" precision = 256 should suffice!
stirM <- stirlerr(nM)
releS <- asNumeric(relErrV(stirM, stirS))
all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux
stopifnot(all.equal(relEx, releS, tolerance = if(noLdbl) 2e-15 else 1e-15))
simpVer3 <- simpVer[simpVer != "lgamma1p"] # have no mpfr-ified lgamma1p()!
## stirlerr_simpl(<mpfr>, *) :
stirM2 <- sapplyMpfr(simpVer3, function(v) stirlerr_simpl(nM, version=v))
## TODO ?:
## apply(stirM2, 2, function(v) all.equal(v, stirS, check.class=FALSE))
## releS2 <- asNumeric(relErrV(stirM2, stirS))
## all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux
## stopifnot(all.equal(relEx, releS, tolerance = 1e-15))
relEmat <- matrix(NA, ncol(stirM2), ncol(stirS),
dimnames = list(simpVer3, colnames(stirS)))
for(j in seq_len(ncol(stirM2)))
for(k in seq_len(ncol(stirS)))
relEmat[j,k] <- asNumeric(relErr(stirM2[,j], stirS[,k]))
relEmat
round(-log10(relEmat), 2) # well .. {why? / expected ?}
cols <- c("gray30", adjustcolor(c(2,3,4), 1/2)); lwd <- c(1, 3,3,3)
stopifnot((k <- length(cols)) == ncol(relEx), k == length(lwd))
matplot(n, relEx, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
mtext("really small errors are dominated by small (< 2^-53) errors of log(n)")
## very interesting: there are different intervals <---> log(n) Qpattern !!
## -- but very small difference, only for n >~= 1/1000 but not before
drawEps.h(negative=TRUE) # abline(h= c(-4,-2:2, 4)*2^-53, lty=c(2,2,2, 1, 2,2,2), col="gray")
legend("topleft", legend = colnames(relEx), col=cols, lwd=3)
## zoomed in a bit:
n. <- 2^-seq.int(400,2, by = -1/4)
relEx. <- relE.lgam1(n.)
matplot(n., relEx., type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
drawEps.h(negative=TRUE)
legend("topleft", legend = colnames(relEx.), col=cols, lwd=3)
##====> Absolute errors (and look at "n0") --------------------------------------
matplot(n., abs19(relEx.), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
drawEps.h(); legend("top", legend = colnames(relEx.), col=cols, lwd=3)
lines(n., abs19(relEx.[,"n0"]), type = "o", cex=1/4, col=cols[4], lwd=2)
## more zooom-in
n.2 <- 2^-seq.int(85, 50, by= -1/100)
stirS.2 <- sapply(c("R3", "lgamma1p", "n0"), function(v) stirlerr_simpl(n.2, version=v))
releS.2 <- asNumeric(relErrV(stirlerr(mpfr(n.2, 320)), stirS.2))
matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3)
abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3)
axis(1, at=5e-17, col.axis=cb, line=-1/4, cex = 3/4)
matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
xaxt="n", xlim = c(8e-18, 1e-15), ## <<<<<<<<<<<<<<<<<<< Zoom-in
xlab = quote(n), main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
eaxis(1); drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3)
abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3)
mtext('stirlerr_simpl(*, "n0") is as good as others for n <= 5e-17', col=adjustcolor(cols[3], 2))
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ===> "all *but* "n0" approximations for "larger" small n:
n2 <- 2^-seq.int(20,0, length.out=1000)
relEx2 <- relE.lgam1(n2)[,c("R3", "lgamma1p", "MM2")] # "n0" is "bad": |relE| >= 2.2e-6 !
cols <- c("gray30", adjustcolor(c(2,4), 1/2)); lwd <- c(1, 3,3)
stopifnot((k <- length(cols)) == ncol(relEx2), k == length(lwd))
matplot(n2, relEx2, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-3,3)*1e-15, xaxt="n",
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
eaxis(1, sub10=c(-3,0)); drawEps.h(negative=TRUE)
legend("topleft", legend = colnames(relEx2), col=cols, lwd=3)
##==> "MM" is *NOT* good for n < 1 *but*
## "the same" -- even larger small n:
n3 <- seq(.01, 5, length=1000)
relEx3 <- relE.lgam1(n3)[,c("R3", "lgamma1p", "MM2")] # "no" is "bad" ..
stopifnot((k <- length(cols)) == ncol(relEx3), k == length(lwd))
matplot(n3, relEx3, type = "l", col=cols, lwd=lwd,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
legend("topleft", legend = colnames(relEx3), col=cols, lwd=3)
drawEps.h(negative=TRUE)
matplot(n3, abs19(relEx3), type = "l", col=cols, lwd=lwd,
log="y", ylim = 2^-c(54, 44), yaxt = "n", ylab = quote(abs(relE)), xlab=quote(n),
main = "|relative errors| of direct (approx.) formula for stirlerr(n), small n")
eaxis(2, cex.axis=0.9); legend("topleft", legend = colnames(relEx3), col=cols, lwd=3)
drawEps.h()
lines(n3, smooth.spline(abs(relEx3)[,1], df=12)$y, lwd=3, col=cols[1])
lines(n3, smooth.spline(abs(relEx3)[,2], df=12)$y, lwd=3, col=adjustcolor(cols[2], 1/2))
lines(n3, smooth.spline(abs(relEx3)[,3], df=12)$y, lwd=4, col=adjustcolor(cols[3], offset = rep(.2,4)))
## ===> from n >~= 1, "MM2" is definitely better up to n = 5 !!
## Check log() only :
plot(n, asNumeric(relErrV(log(mpfr(n, 256)), log(n))), ylim = c(-1,1)*2^-53,
log="x", type="l", xaxt="n") ## ===> indeed --- log(n) approximation pattern !!
eaxis(1) ; drawEps.h(negative=TRUE)
showProc.time()
## =========== "R3" vs "lgamma1p" -------------------------- which is better?
## really for the very small n, all is dominated by -(n+0.5)*log(n); and lgamma1p() is unnecessary!
i <- 1:20; ni <- n[i]
lgamma1p(ni)
- (ni +0.5)*log(ni) + ni
## much less extreme:
n2 <- lseq(2^-12, 1/2, length=1000)
relE2 <- relE.lgam1(n2)[,-4]
cols <- c("gray30", adjustcolor(2:3, 1/2)); lwd <- c(1,3,3)
matplot(n2, relE2, type = "l", log="x", col=cols, lwd=lwd)
legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3)
drawEps.h(negative=TRUE)
matplot(n2, abs19(relE2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(6e-17, 1e-15),
xaxt = "n"); eaxis(1, sub10=c(-2,0))
legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3)
drawEps.h()
## "MM2" is *worse* here, n < 1/2
for(j in 1:3) lines(n2, smooth.spline(abs(relE2[,j]), df=10)$y, lwd=3,
col=adjustcolor(cols[j], 1.5, offset = rep(-1/4, 4)))
## "lgammap very slightly better in [0.002, 0.05] ...
## "TODO": draw 0.90-quantile curves {--> cobs::cobs() ?} instead of mean-curves?
## which is better? ... "random difference"
d.absrelE <- abs(relE2[,"R3"]) - abs(relE2[,"lgamma1p"])
plot (n2, d.absrelE, type = "l", log="x", # no clear picture ...
main = "|relE_R3| - |relE_lgamma1p|", axes=FALSE, frame.plot=TRUE)
eaxis(1, sub10=c(-2,1)); eaxis(2); axis(3, at=max(n2)); abline(v = max(n2), lty=3, col="gray")
## 'lgamma1p' very slightly better:
lines(n2, smooth.spline(d.absrelE, df=12)$y, lwd=3, col=2)
## not really small n at all == here see, how "bad" the direct formula gets for 1 < n < 10 or so
n3 <- lseq(2^-14, 2^2, length=800)
relE3 <- relE.lgam1(n3)[, -4]
matplot(n3, relE3, type = "l", log="x", col=cols, lty=1, lwd = c(1,3),
main = quote(rel.lgam1(n)), xlab=quote(n))
matplot(n3, abs19(relE3), type = "l", log="xy", col=cols, lwd = c(1,3), xaxt="n",
main = quote(abs(rel.lgam1(n))), xlab=quote(n), ylim = c(2e-17, 4e-14))
drawEps.h(); eaxis(1, sub10=c(-2,3))
legend("topleft", legend=colnames(relE3), col=cols, lwd=2)
## very small difference --- draw the 3 smoothers :
for(j in 1:3) {
ll <- lowess(log(n3), abs19(relE3[, j]), f= 1/12)
with(ll, lines(exp(x), y, col=adjustcolor(cols[j], 1.5), lwd=3))
}
## ==> lgamma1p(.) very slightly in n ~ 10^-4 -- 10^-2 --- but not where it matters: n ~ 0.1 -- 1 !!
## "MM2" gets best from n >~ 1 !
abline(v=1, lty=3, col = adjustcolor(1, 3/4))
showProc.time()
### 2. relErr( stirlerr(.) ) ============================================================
##' Very revealing plot showing the *relative* approximation error of stirlerr(<dblprec>)
##'
p.stirlerrDev <- function(n, precBits = if(doExtras) 2048L else 512L,
stnM = stirlerr(mpfr(n, precBits), use.halves=use.halves, verbose=verbose),
abs = FALSE,
## cut points, n=*, in the stirlerr() algorithm; "FIXME": sync with ../R/dgamma.R <<<<
scheme = c("R3", "R4.4_0"),
cutoffs = switch(match.arg(scheme)
, R3 = c(15, 35, 80, 500)
, R4.4_0 = c(4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.7,
6.1, 6.5, 7, 7.9, 8.75, 10.5, 13,
20, 26, 60, 200, 3300, 17.4e6)
## {FIXME: need to sync} <==> ../man/stirlerr.Rd <==> ../R/dgamma.R
),
use.halves = missing(cutoffs),
direct.ver = c("R3", "lgamma1p", "MM2", "n0"),
verbose = getOption("verbose"),
type = "b", cex = 1,
col = adjustcolor(1, 3/4), colnB = adjustcolor("orange4", 1/3),
log = if(abs) "xy" else "x",
xlim=NULL, ylim = if(abs) c(8e-18, max(abs(N(relE)))))
{
op <- par(las = 1, mgp=c(2, 0.6, 0))
on.exit(par(op))
require("Rmpfr"); require("sfsmisc")
st <- stirlerr(n, scheme=scheme, cutoffs=cutoffs, use.halves=use.halves, direct.ver=direct.ver,
verbose=verbose)
relE <- relErrV(stnM, st) # eps0 = .Machine$double.xmin
N <- asNumeric
form <- if(abs) abs(N(relE)) ~ n else N(relE) ~ n
plot(form, log=log, type=type, cex=cex, col=col, xlim=xlim, ylim=ylim,
ylab = quote(relErrV(stM, st)), axes=FALSE, frame.plot=TRUE,
main = sprintf("stirlerr(n, cutoffs) rel.error [wrt stirlerr(Rmpfr::mpfr(n, %d))]",
precBits))
eaxis(1, sub10=3)
eaxis(2)
mtext(paste("cutoffs =", deparse1(cutoffs)))
ylog <- par("ylog")
## FIXME: improve this ---> drawEps.h() above
if(ylog) {
epsC <- c(1,2,4,8)*2^-52
epsCxp <- expression(epsilon[C],2*epsilon[C], 4*epsilon[C], 8*epsilon[C])
} else {
epsC <- (-2:2)*2^-52
epsCxp <- expression(-2*epsilon[C],-epsilon[C], 0, +epsilon[C], +2*epsilon[C])
}
dy <- diff(par("usr")[3:4])
if(diff(range(if(ylog) log10(epsC) else epsC)) > dy/50) {
lw <- rep(1/2, 5); lw[if(ylog) 1 else 3] <- 2
abline( h=epsC, lty=3, lwd=lw)
axis(4, at=epsC, epsCxp, las=2, cex.axis = 3/4, mgp=c(3/4, 1/4, 0), tck=0)
} else ## only x-axis
abline(h=if(ylog) epsC else 0, lty=3, lwd=2)
abline(v = cutoffs, col=colnB)
axis(3, at=cutoffs, col=colnB, col.axis=colnB,
labels = formatC(cutoffs, digits=3, width=1))
invisible(relE)
} ## p.stirlerrDev()
showProc.time()
n <- lseq(2^-10, 1e10, length=4096)
n <- lseq(2^-10, 5000, length=4096)
## store "expensive" stirlerr() result, and re-use many times below:
nM <- mpfr(n, if(doExtras) 2048 else 512)
st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose
p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE) # default cutoffs= c(15, 40, 85, 600)
p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE, ylim = c(-1,1)*1e-12) # default cutoffs= c(15, 40, 85, 600)
## show the zoom-in region in next plot
yl2 <- 3e-14*c(-1,1)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
if(do.pdf) { dev.off() ; pdf("stirlerr-relErr_1.pdf") }
## drop n < 7:
p.stirlerrDev(n=n, stnM=st.nM, xlim = c(7, max(n)), use.halves=FALSE) # default cutoffs= c(15, 40, 85, 600)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
## The first plot clearly shows we should do better:
## Current code is switching to less terms too early, loosing up to 2 decimals precision
if(FALSE) # no visible difference {use.halves = T / F }:
p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = FALSE)
p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = TRUE)# exact at n/2 (n <= ..)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
showProc.time()
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-1.pdf") }
### ~19.April 2021: "This is close to *the* solution" (but see 'cuts' below)
cuts <- c(7, 12, 20, 26, 60, 200, 3300)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
st. <- stirlerr(n=n , cutoffs = cuts, verbose=TRUE)
st.nM <- stirlerr(n=nM, cutoffs = cuts, use.halves=FALSE) ## << on purpose
relE <- asNumeric(relErrV(st.nM, st.))
head(cbind(n, relE), 20)
## nice printout :
print(cbind(n = format(n, drop0trailing = TRUE),
stirlerr= format(st.,scientific=FALSE, digits=4),
relErr = signif(relE, 4))
, quote=FALSE)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts)
## and zoom in:
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2/20)
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-2.pdf") }
## zoom in ==> {good for n >= 10}
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", ylim = 2e-15*c(-1,1),
cutoffs = cuts)## old default cutoffs = c(15,35, 80, 500)
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-3.pdf") }
showProc.time()
##-- April 20: have more terms up to S10 in stirlerr() --> can use more cutoffs
n <- n5m <- lseq(1/64, 5000, length=4096)
nM <- mpfr(n, if(doExtras) 2048L # a *lot* accuracy for stirlerr(nM,*)
else 512L)
ct10.1 <- c( 5.4, 7.5, 8.5, 10.625, 12.125, 20, 26, 60, 200, 3300)# till 2024-01-19
ct10.2 <- c( 5.4, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300)
cuts <-
ct12.1 <- c(5.22, 6.5, 7.0, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## 5.25 is "too small" but the direct formula is already really bad there, ...
st.nM <- roundMpfr(stirlerr(nM, use.halves=FALSE, ## << on purpose;
verbose=TRUE), precBits = 128)
## NB: for x=xM <mpfr>; `cutoffs` are *not* used.
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0")
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*3e-15)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*1e-15)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", abs=TRUE)
axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4)
abline(v = 2:6, lty=3, col=cola)
if(FALSE)## using exact values sferr_halves[] *instead* of MPFR ones: ==> confirmation they lay on top
lines((0:30)/2, abs(stirlerr((0:30)/2, cutoffs=cuts, verbose=TRUE)/DPQ:::sferr_halves - 1), type="o", col=2,lwd=2)
if(FALSE) ## nice (but unneeded) printout :
print(cbind(n = format(n, drop0trailing = TRUE),
stirlerr= format(st.,scientific=FALSE, digits=4),
relErr = signif(relE, 4))
, quote=FALSE)
showProc.time()
## ========== where should the cutoffs be ? ===================================================
.stirl.cutoffs <- function(scheme)
eval(do.call(substitute, list(formals(stirlerr)$cutoffs, list(scheme = scheme))))
drawCuts <- function(scheme, axis=NA, lty = 3, col = "skyblue", ...) {
abline(v = (ct <- .stirl.cutoffs(scheme)), lty=lty, col=col, ...)
if(is.finite(axis)) axisCuts(side = axis, at = ct, col=col, ...)
}
axisCuts <- function(scheme, side = 3, at = .stirl.cutoffs(scheme), col = "skyblue", line = -3/4, ...)
axis(side, at=at, labels=formatC(at), col.axis = col, col=NA, col.ticks=NA, line=line, ...)
mtextCuts <- function(cutoffs, scheme, ...) {
if(!missing(scheme)) cutoffs <- .stirl.cutoffs(scheme)
mtext(paste("cutoffs =", deparse1(cutoffs)), ...)
}
if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k.pdf") }
mK <- 20L # := max(k)
## order = k = 1:mK terms in series approx:
k <- 1:mK
n <- 2^seq(1, 28, by=1/16)
nM <- mpfr(n, 1024)
stnM <- stirlerr(nM) # the "true" values
stirlOrd <- sapply(k, function(k.) stirlerr(n, order = k.))
stirlO_lgcor <- cbind(stirlOrd, sapply(5:6, function(nal) lgammacor(n, nalgm = nal)))
relE <- asNumeric(stirlO_lgcor/stnM -1) # "true" relative error
## use a "smooth" but well visible polette :
palROBG <- colorRampPalette(c("red", "darkorange2", "blue", "seagreen"), space = "Lab")
palette(adjustcolor(palROBG(mK+2), 3/4))
## -- 2 lgamcor()'s
(tit.k <- substitute(list( stirlerr(n, order=k) ~~"error", k == 1:mK), list(mK = mK)))
(tit.kA <- substitute(list(abs(stirlerr(n, order=k) ~~"error"), k == 1:mK), list(mK = mK)))
lgammacorTit <- function(...) mtext("+ lgammacor(x, 5) [bad] + lgammacor(x, 6) [good]", col=1, ...)
matplotB(n, relE, cex=2/3, ylim = c(-1,1)*1e-13, col=k,
log = "x", xaxt="n", main = tit.k)
lgammacorTit()
eaxis(1, nintLog = 20)
drawCuts("R4.4_0")
## zoom in (ylim)
matplotB(n, relE, cex=2/3, ylim = c(-1,1)*5e-15, col=k,
log = "x", xaxt="n", main = tit.k)
lgammacorTit()
eaxis(1, nintLog = 20); abline(h = (-2:2)*2^-53, lty=3, lwd=1/2)
drawCuts("R4.4_0", axis = 3)
## log-log |rel.Err| -- "linear"
matplotB(n, abs19(relE), cex=2/3, col=k, ylim = c(8e-17, 1e-7), log = "xy", main=tit.kA)
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
lgammacorTit(line=-1)
drawCuts("R4.4_0", axis = 3)
## zoom in -- still "large" {no longer carry the lgammacor() .. }:
n2c <- 2^seq(2, 8, by=1/256)
nMc <- mpfr(n2c, 1024)
stnMc <- stirlerr(nMc) # the "true" values
stirlOrc <- sapply(k, function(k.) stirlerr(n2c, order = k.))
relEc <- asNumeric(stirlOrc/stnMc -1) # "true" relative error
matplotB(n2c, relEc, cex=2/3, ylim = c(-1,1)*1e-13, col=k,
log = "x", xaxt="n", main = tit.k)
eaxis(1, sub10 = 2)
drawCuts("R4.4_0", axis=3)
## log-log |rel.Err| -- "linear"
matplotB(n2c, abs19(relEc), cex=2/3, col=k, ylim = c(8e-17, 1e-3), log = "xy", main=tit.kA)
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
drawCuts("R4.4_0", axis = 3)
## zoom into the critical n region
nc <- seq(3.5, 11, by=1/128)
ncM <- mpfr(nc, 256)
stncM <- stirlerr(ncM) # the "true" values
stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k))
relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error
## log-log |rel.Err| -- "linear"
matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(2e-17, 1e-8),
log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2)
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6))
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4))
## lines(nc, abs19(asNumeric(stirlerr_simpl(nc,"MM2")/stncM - 1)), lwd=3, col=adjustcolor("purple", 2/3))
legend(10^par("usr")[1], 1e-9, legend=paste0("k=", k), bty="n", lwd=2,
col=k, lty=1:5, pch= c(1L:9L, 0L, letters)[seq_along(k)])
drawCuts("R4.4_0", axis=3)
## Zoom-in [only]
matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(4e-17, 1e-11), xlim = c(4.8, 6.5),
log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2)
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6))
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4))
k. <- k[-(1:6)]
legend("bottomleft", legend=paste0("k=", k.), bty="n", lwd=2,
col=k., lty=1:5, pch= c(1L:9L, 0L, letters)[k.])
drawCuts("R4.4_0", axis=3)
showProc.time()
##--- Accuracy of "R4.4_0" -------------------------------------------------------
for(nc in list(seq(4.75, 28, by=1/512), # for a bigger pix
seq(4.75, 9, by=1/1024)))
{
ncM <- mpfr(nc, 1024)
stncM <- stirlerr(ncM) # the "true" values
stirl.440 <- stirlerr(nc, scheme = "R4.4_0")
stirl.3 <- stirlerr(nc, scheme = "R3")
relE440 <- asNumeric(relErrV(stncM, stirl.440))
relE3 <- asNumeric(relErrV(stncM, stirl.3 ))
##
plot(nc, abs19(relE440), xlab=quote(n), main = quote(abs(relErr(stirlerr(n, '"R4.4_0"')))),
type = "l", log = "xy", ylim = c(4e-17, 1e-13))
mtextCuts(scheme="R4.4_0", cex=4/5)
drawCuts("R4.4_0", lty=2, lwd=2, axis=4)
drawEps.h()
if(max(nc) <= 10) abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(4, 1/2))
if(TRUE) { # but just so ...
c3 <- adjustcolor("royalblue", 1/2)
lines(nc, pmax(abs(relE3), 1e-18), col=c3)
title(quote(abs(relErr(stirlerr(n, '"R3"')))), adj=1, col.main = c3)
drawCuts("R3", lty=4, col=c3); mtextCuts(scheme="R3", adj=1, col=c3)
}
addOrd <- TRUE
addOrd <- dev.interactive(orNone=TRUE)
if(addOrd) {
if(max(nc) >= 100) {
i <- (15 <= nc & nc <= 85) # (no-op in [4.75, 7] !)
ni <- nc[i]
} else { i <- TRUE; ni <- nc }
for(k in 8:17) lines(ni, abs19(asNumeric(relErrV(stncM[i], stirlerr(ni, order=k)))), col=adjustcolor(k, 1/3))
title(sub = "stirlerr(*, order k = 8:17)")
}
} ## for(nc ..)
if(FALSE)
lines(nc, abs19(relE440))# *re* draw!
showProc.time()
palette("Tableau")
## Focus more:
stirlerrPlot <- function(nc, k, res=NULL, legend.xy = "left", full=TRUE, precB = 1024,
ylim = c(3e-17, 2e-13), cex = 5/4) {
stopifnot(require("Rmpfr"), require("graphics"))
if(is.list(res) && all(c("nc", "k", "relEc","splarelE") %in% names(res))) { ## do *not* recompute
list2env(res, envir = environment())
} else { ## compute
stopifnot(is.finite(nc), nc > 0, length(nc) >= 100, k == as.integer(k), 0 <= k, k <= 20)
ncM <- mpfr(nc, precB)
stncM <- stirlerr(ncM) # the "true" values
stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k))
relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error
## log |rel.Err| -- "linear"
## smooth() on log-scale {and transform back}:
splarelE <- apply(log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y))
## the direct formulas (default "R3", "MM2"):
arelEs0 <- abs(asNumeric(stirlerr_simpl(nc )/stncM - 1))
arelEs2 <- abs(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1))
}
pch. <- c(1L:9L, 0L, letters)[k]
if(full)
matplotB(nc, abs19(relEc), col=k, pch = pch., cex=cex, ylim=ylim, log = "y",
xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
else ## smooth only
matplotB(nc, splarelE, col=adjustcolor(k,2/3), pch=pch., lwd=2, cex=cex, ylim=ylim, log = "y",
xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
legend(legend.xy, legend=paste0("k=", k), bty="n", lwd=2, col=k, lty=1:5, pch = pch.)
abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(10, 1/2))
drawCuts("R4.4_0", axis=3)
if(full) {
matlines(nc, splarelE, col=adjustcolor(k,2/3), lwd=4)
lines(nc, pmax(arelEs0, 1e-19), lwd=1.5, col=adjustcolor( 2, 0.2))
lines(nc, pmax(arelEs2, 1e-19), lwd=1, col=adjustcolor(10, 0.2))
}
lines(nc, smooth.spline(arelEs0, df=12)$y, lwd=3, col= adjustcolor( 2, 1/2))
lines(nc, smooth.spline(arelEs2, df=12)$y, lwd=3, col= adjustcolor(10, 1/2))
invisible(list(nc=nc, k=k, relEc = relEc, splarelE = splarelE, arelEs0=arelEs0, arelEs2=arelEs2))
}
rr1 <- stirlerrPlot(nc = seq(4.75, 9.0, by=1/1024),
k = 7:20)
stirlerrPlot(res = rr1, full=FALSE, ylim = c(8e-17, 1e-13))
if(interactive())
stirlerrPlot(res = rr1)
rr <- stirlerrPlot(nc = seq(5, 6.25, by=1/2048), k = 9:18)
stirlerrPlot(res = rr, full=FALSE, ylim = c(8e-17, 1e-13))
showProc.time()
palette("default")
if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k-vs-k1.pdf") }
##' Find 'cuts', i.e., a region c(k) +/- s(k) i.e. intervals [c(k) - s(k), c(k) + s(k)]
##' where c(k) is such that relE(n=c(k), k) ~= eps)
##' 1. Find the c1(k) such that |relE(n, k)| ~= c1(k) * n^{-2k}
findC1 <- function(n, ks, e1 = 1e-15, e2 = 1e-5, res=NULL, precBits = 1024, do.plot = TRUE, ...)
{
if(is.list(res) && all(c("n", "ks", "arelE") %in% names(res))) { ## do *not* recompute, take from 'res':
list2env(res, envir = environment())
} else { ## compute
stopifnot(require("Rmpfr"),
is.numeric(ks), ks == (k. <- as.integer(ks)), length(ks <- k.) >= 1,
length(e1) == 1L, length(e2) == 1L, is.finite(c(e1,e2)), e1 >= 0, e2 >= e1,
0 <= ks, ks <= 20, is.numeric(n), n > 0, is.finite(n), length(n) >= 100)
nM <- mpfr(n, precBits)
stirM <- stirlerr(nM) # the "true" values
stirlOrd <- sapply(ks, function(k) stirlerr(n, order = k))
arelE <- abs(asNumeric(stirlOrd/stirM -1)) # "true" relative error
}
arelE19 <- pmax(arelE, 1e-19)
## log |rel.Err| -- "linear"
## on log-scale {and transform back; for linear fit, only use values inside [e1, e2]
if(do.plot) # experi
matplotB(n, arelE19, log="xy", ...)
## matplot(n, arelE19, type="l", log="xy", xlim = c(min(n), 20))
## re-compute these, as they *also* depend on (e1, e2)
c1 <- vapply(seq_along(ks), function(i) {
k <- ks[i]
y <- arelE19[,i]
iUse <- e1 <= y & y <= e2
if(sum(iUse) < 10) stop("only", sum(iUse), "values in [e1,e2]")
## .lm.fit(cbind(1, log(n[iUse])), log(y[iUse]))$coefficients
## rather, we *know* the error is c* n^{-2k} , i.e.,
## log |relE| = log(c) - 2k * log(n)
## <==> c = exp( log|relE| + 2k * log(n))
exp(mean(log(y[iUse]) + 2*k * log(n[iUse])))
}, numeric(1))
if(do.plot) {
drawEps.h()
for(i in seq_along(ks))
lines(n, c1[i] * n^(-2*ks[i]), col=adjustcolor(i, 1/3), lwd = 4, lty = 2)
}
invisible(list(n=n, ks=ks, arelE = arelE, c1 = c1))
} ## findC1()
c1.Res <- findC1(n = 2^seq(2, 26, by=1/128), ks = 1:18)
(s.c1.fil <- paste0("stirlerr-c1Res-", myPlatform(), ".rds"))
saveRDS(c1.Res, file = s.c1.fil)
if(!exists("c1.Res")) {
c1.Res <- readRDS(s.c1.fil)
## re-do the "default" plot of findC1():
findC1(res = c1.Res, xaxt="n"); eaxis(1, sub10=2)
}
## the same, zoomed in:
findC1(res = c1.Res, xlim = c(4, 40), ylim = c(2e-17, 1e-12))
ks <- c1.Res$ks; pch. <- c(1L:9L, 0L, letters)[ks]
legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.)
## smaller set : larger e1 :
c1.r2 <- findC1(res = c1.Res, xlim = c(4, 30), ylim = c(4e-17, 1e-13), e1 = 4e-15)
legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.)
print(digits = 4,
cbind(ks, c1. = c1.Res$c1, c1.2 = c1.r2$c1,
relD = round(relErrV(c1.r2$c1, c1.Res$c1), 4)))
c1.. <- c1.r2[c("ks", "c1")] # just the smallest part is needed here:
(s.c1.fil <- paste0("stirlerr-c1r2-", myPlatform(), ".rds"))
saveRDS(c1.., file = s.c1.fil) # was "stirlerr-c1.rds"
## 2. Now, find the n(k) +/- se.n(k) intervals
## Use these c1 from above
##' Given c1-results relErr(n,k); |relE(n,k) ~= c1 * n^{-2k} , find n such that
##' |relE(n,k)| line {in log-log scale} cuts y = eps, i.e., n* such that |relE(n,k)| <= eps for all n >= n*
n.ep <- function(eps, c1Res, ks = c1Res$ks, c1 = c1Res$c1, ...) {
stopifnot(is.finite(eps), length(eps) == 1L, eps > 0,
length(ks) == length(c1), is.numeric(c1), is.integer(ks), ks >= 1)
## n: given k, the location where the |relE(n,k)| line {in log-log} cuts y = eps
## |relE(n,k) ~= c1 * n^{-2k} <==>
## log|relE(n,k)| ~= log(c1) - 2k* log(n) <==>
## c := mean{ exp( log|relE(n,k)| + 2k* log(n) ) } ------- see findC1()
## now, solve for n :
## c1 * n^{-2k} == eps
## log(c1) - 2k* log(n) == log(eps)
## log(n) == (log(eps) - log(c1)) / (-2k) <==>
## n == exp((log(c1) - log(eps)) / 2k)
exp((log(c1) - log(eps))/(2*ks))
}
## get c1..
if(!exists("c1..")) c1.. <- readRDS("stirlerr-c1.rds")
ne2 <- n.ep(2^-51, c1Res = c1..) ## ok
ne1 <- n.ep(2^-52, c1Res = c1..)
ne. <- n.ep(2^-53, c1Res = c1..)
form <- function(n) format(signif(n, 3), scientific=FALSE)
data.frame(k = ks, ne2 = form(ne2), ne1 = form(ne1), ne. = form(ne.),
cutoffs = form(rev(.stirl.cutoffs("R4.4_0")[-1])))
## ------- Linux F 36/38 x86_64 (nb-mm5|v-lynne)
## k ne2 ne1 ne. cutoffs
## 1 8660000.00 12200000.00 17300000.00 17400000.00
## 2 2150.00 2560.00 3040.00 3700.00
## 3 159.00 178.00 200.00 200.00
## 4 46.60 50.80 55.40 81.00
## 5 23.40 25.10 26.90 36.00
## 6 15.20 16.10 17.10 25.00
## 7 11.50 12.10 12.70 19.00
## 8 9.43 9.85 10.30 14.00
## 9 8.20 8.53 8.86 11.00
## 10 7.42 7.68 7.95 9.50
## 11 6.89 7.11 7.34 8.80
## 12 6.53 6.72 6.92 8.25
## 13 6.27 6.44 6.62 7.60
## 14 6.10 6.25 6.41 7.10
## 15 5.98 6.12 6.26 6.50 * (not used)
## 16 5.90 6.03 6.16 6.50 * " "
## 17 5.85 5.97 6.10 6.50 * " "
## 18 5.83 5.95 6.06 6.50 << used all the way down to 5.25
## ok --- correct order of magnitude ! --- good!
## 2b. find *interval* around the 'n(eps)' values
## -- Try simply
d.k <- ne. - ne1
## interval
int.k <- cbind(ne1 - d.k,
ne1 + d.k)
## look at e.g.
data.frame(k=ks, `n(k)` = form(ne1), int = form(int.k))
## k n.k. int.1 int.2
## 1 12200000.00 7180000.00 17300000.00
## 2 2560.00 2070.00 3040.00
## 3 178.00 156.00 200.00
## 4 50.80 46.20 55.40
## 5 25.10 23.30 26.90
## 6 16.10 15.20 17.10
## 7 12.10 11.40 12.70
## 8 9.85 9.41 10.30
## 9 8.53 8.19 8.86
## 10 7.68 7.41 7.95
## 11 7.11 6.88 7.34
## 12 6.72 6.52 6.92
## 13 6.44 6.27 6.62
## 14 6.25 6.10 6.41
## 15 6.12 5.98 6.26
## 16 6.03 5.90 6.16
## 17 5.97 5.85 6.10
## 18 5.95 5.83 6.06
##' as function {well, *not* computing c1.k from scratch
nInt <- function(k, c1.k, ep12 = 2^-(52:53)) {
if(length(k) == 1L) { # special convention to call for *one* k, with c1.k vector
stopifnot(k == (k <- as.integer(k)), k >= 1, length(c1.k) >= k)
c1.k <- c1.k[k]
}
## see n.ep() above
n_ <- function(eps, k, c1) {
stopifnot(is.finite(eps), length(eps) == 1L, eps > 0,
length(k) == length(c1), is.numeric(c1), is.integer(k), k >= 1)
exp((log(c1) - log(eps))/(2*k))
}
ne1 <- n_(ep12[1], k, c1.k)
ne. <- n_(ep12[2], k, c1.k)
d.k <- ne. - ne1
stopifnot(d.k > 0)
## interval: {"fudge" 0.5 / 2.5} from results -- also 'noLdbl' gives *quite* different pic!:
odd <- k %% 2 == 1
cbind(ne1 - ifelse( odd & noLdbl, 1.5, 0.5) * d.k,
ne1 + ifelse(!odd , 8, 2.5) * d.k)
}
nInt(k= 1, c1..$c1)
nInt(k= 2, c1..$c1)
nInt(k=18, c1..$c1)
nints.k <- nInt(ks, c1..$c1)
## for printing
form(as.data.frame( nints.k ))
## (-.5 , +2.5) ## originally ( -1, +1)
## 1 9710000.00 24900000.00 # 7180000.00 17300000.00
## 2 2320.00 3770.00 # 2070.00 3040.00
## 3 167.00 232.00 # 156.00 200.00
## 4 48.50 62.30 # 46.20 55.40
## 5 24.20 29.60 # 23.30 26.90
## 6 15.70 18.50 # 15.20 17.10
## 7 11.70 13.60 # 11.40 12.70
## 8 9.63 10.90 # 9.41 10.30
## 9 8.36 9.36 # 8.19 8.86
## 10 7.54 8.36 # 7.41 7.95
## 11 7.00 7.68 # 6.88 7.34
## 12 6.62 7.21 # 6.52 6.92
## 13 6.36 6.88 # 6.27 6.62
## 14 6.17 6.64 # 6.10 6.41
## 15 6.05 6.48 # 5.98 6.26
## 16 5.96 6.36 # 5.90 6.16
## 17 5.91 6.28 # 5.85 6.10
## 18 5.89 6.24 # 5.83 6.06
## 3. Then for each of the intervals, compare order k vs k+1
## -- ==> optimal cutoff { how much platform dependency ?? }
### Here, compute only
find1cuts <- function(k, c1,
n = nInt(k, c1), # the *set* of n's or the 'range'
len.n = 1000,
precBits = 1024, nM = mpfr(n, precBits),
stnM = stirlerr(nM),
stirlOrd = sapply(k+(0:1), function(.k.) stirlerr(n, order = .k.)),
relE = asNumeric(stirlOrd/stnM -1), # "true" relative error for the {k, k+1}
do.spl=TRUE, df.spline = 9, # df = 5 gives 2 cutpoints for k==1
do.low=TRUE, f.lowess = 0.2,
do.cobs = require("cobs", quietly=TRUE), tau = 0.90)
{
## check relErrV( stirlerr(n, order=k ) vs
## stirlerr(n, order=k+1)
if(length(n) == 2L)
if(n[1] < n[2]) n <- seq(n[1], n[2], length.out = len.n)
else stop("'n' must be *increasing")
force(relE)
y <- abs19(relE)
## NB: all smoothing --- as in stirlerrPlot() above -- should happen in log-space
## (log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y))
ly <- log(y) # == log(abs19(relE)) == log(max(|r|, 1e-19))
if(do.spl) {##
s1 <- exp(smooth.spline(ly[,1], df=df.spline)$y)
s2 <- exp(smooth.spline(ly[,2], df=df.spline)$y)
}
if(do.low) { ## lowess
s1l <- exp(lowess(ly[,1], f=f.lowess)$y)
s2l <- exp(lowess(ly[,2], f=f.lowess)$y)
}
## also use cobs() splines for the 90% quantile !!
EE <- environment() # so can set do.cobs to FALSE in case of error
if(do.cobs) { ## <==> require("cobs") # yes, this is in tests/
cobsF <- function(Y) cobs(n, Y, tau=tau, nknots = 6, lambda = -1,
print.warn=FALSE, print.mesg=FALSE)
## sparseM::chol(<matrix.csr>) now gives error when it gave warning {about singularity}
cobsF <- function(Y) {
r <- tryCatch(cobs(n, Y, tau=tau, nknots = 6, lambda = -1,
print.warn=FALSE, print.mesg=FALSE),
error = identity)
if(inherits(r, "error")) {
assign("do.cobs", FALSE, envir = EE) # and return
list(fitted = FALSE)
}
else
r
}
cs1 <- exp(cobsF(ly[,1])$fitted)
cs2 <- exp(cobsF(ly[,2])$fitted)
}
smooths <- list(spl = if(do.spl ) cbind(s1, s2 ),
low = if(do.low ) cbind(s1l,s2l),
cobs= if(do.cobs) cbind(cs1,cs2))
## diffL <- list(spl = if(do.spl ) s2 -s1,
## low = if(do.low ) s2l-s1l,
## cobs= if(do.cobs) cs2-cs1)
### FIXME: simplification does not always *work* -- (i, n.) are not always ok
### ------ notably within R-devel-no-ldouble {hence probably macOS M1 ... ..}
sapply(smooths, function(s12) { d <- s12[,2] - s12[,1]
## typically a (almost or completely) montone increasing function, crossing zero *once*
## compute cutpoint:
## i := the first n[i] with d(n[i]) >= 0
i <- which(d >= 0)[1]
if(length(i) == 1L && !is.na(i) && i > 1L) {
i_ <- i - 1L # ==> d(n[i_]) < 0
## cutpoint must be in [n[i_], n[i]] --- do linear interpolation
n. <- n[i_] - (n[i] - n[i_])* d[i_] / (d[i] - d[i_])
} else {
if(length(i) != 1L) i <- -length(i)
n. <- NA_integer_
}
c(i=i, n.=n.)
}) -> n.L
list(k=k, n=n, relE = unname(relE), smooths=smooths, i.n = n.L)
} ## find1cuts()
k. <- 1:15
system.time(
## Failed on lynne [2024-06-04, R 4.4.1 beta] with
## Error in .local(x, ...) : insufficient space ---> SparseM :: chol(<..>)
## ==> now we catch this inside find1cuts():
resL <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1))
) ## -- warnings, notably from cobs() not converging
## needs 12 sec (!!) user system elapsed =
(s.find15.fil <- paste0("stirlerr-find1_1-15_", myPlatform(), ".rds"))
## now we catch cobs() errors {from SparseM::chol},
## okCuts <- !inherits(resL, "error")
## if(okCuts) {
saveRDS(resL, file = s.find15.fil) # was "stirlerr-find1_1-15.rds"
## } else traceback()
## 11: stop(mess)
## 10: .local(x, ...)
## 9: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax)
## 8: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax)
## 7: rq.fit.sfnc(Xeq, Yeq, Xieq, Yieq, tau = tau, rhs = rhs, control = rqCtrl)
## 6: drqssbc2(x, y, w, pw = pw, knots = knots, degree = degree, Tlambda = if (select.lambda) lambdaSet else lambda,
## constraint = constraint, ptConstr = ptConstr, maxiter = maxiter,
## trace = trace - 1, nrq, nl1, neqc, niqc, nvar, tau = tau,
## select.lambda = select.lambda, give.pseudo.x = keep.x.ps,
## rq.tol = rq.tol, tol.0res = tol.0res, print.warn = print.warn)
## 5: cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE,
## print.mesg = FALSE) at stirlerr-tst.R!udBzpT#32
## 4: cobsF(ly[, 1]) at stirlerr-tst.R!udBzpT#34
## 3: find1cuts(k = k, c1 = c1..$c1) at #1
## 2: FUN(X[[i]], ...)
## 1: lapply(setNames(, k.), function(k) find1cuts(k = k, c1 = c1..$c1))
ok1cutsLst <- function(res) {
stopifnot(is.list(res), sapply(res, is.list)) # must be list of lists
cobsL <- lapply(lapply(res, `[[`, "smooths"), `[[`, "cobs")
vapply(cobsL, is.array, NA)
}
(resLok <- ok1cutsLst(resL))
## 1 2 3 4 5 6 ...... 15
## FALSE TRUE TRUE TRUE TRUE TRUE ...... TRUE
if(FALSE) {
## e.g. in R-devel-no-ldouble:
(r1 <- find1cuts(k=1, c1=c1..$c1))$i.n # list --- no longer ok {SparseM::chol -> insufficient space}
(r2 <- find1cuts(k=2, c1=c1..$c1))$i.n # "good"
(r3 <- find1cuts(k=3, c1=c1..$c1))$i.n # 3-vector: only 3x 'i' == 1
}
## if(okCuts) {
mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1")$old.par -> opar
invisible(lapply(resL[resLok], plot1cuts))
## plus the "last" ones {also showing that k=15 is worse here anyway than k=17}
str(r17 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = 1500), c1=c1..$c1))
plot1cuts(r17) # no-ldouble is *very* different than normal: k *much better* than k+1
(s.find17.fil <- paste0("stirlerr-find1_17_", myPlatform(), ".rds"))
saveRDS(r17, file = s.find17.fil) # was "stirlerr-find1_17.rds"
## }
## ditto for tau = 0.8 {quantile for cobs()}:
system.time(
resL.8 <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1, tau = 0.80))
) # warnings from cobs {again 12.1 sec}
(resL8ok <- ok1cutsLst(resL.8))
mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1 -- tau = 0.8")
invisible(lapply(resL.8[resL8ok], plot1cuts))
## plus "last" one
str(r17.8 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = 1500), c1=c1..$c1, tau = 0.80))
plot1cuts(r17.8)
par(opar)
n.mn <- sapply(resL, `[[`, "i.n", simplify = "array")
n.mn8 <- sapply(resL.8, `[[`, "i.n", simplify = "array")
## of course, only the cobs part differs:
## ... when I later see errors here ... what's going on??
str(n.mn) # all NULL for "no-double" !!
dim(n.mn["n.",,])
str(n.mn8["n.", "cobs",])
sessionInfo()
form(data.frame(k = k., cutoffs = rev(.stirl.cutoffs("R4.4_0")[-1])[k.],
t(n.mn ["n.",,]), cobs.80 = n.mn8["n.", "cobs",]))
## k cutoffs spl low cobs cobs.80
## 1 1 17400000.00 15600000.00 15800000.00 15800000.00 15700000.00
## 2 2 3700.00 NA NA NA NA == from tau=0.90 I'd use ~ 5000
## 3 3 200.00 200.00 201.00 206.00 205.00 205
## 4 4 81.00 NA NA 86.60 86.20 86
## 5 5 36.00 27.20 27.10 27.00 27.10 27
## 6 6 25.00 23.50 NA NA NA 23.5
## 7 7 19.00 12.80 12.80 12.80 12.80 12.8
## 8 8 14.00 NA NA 12.30 12.80 12.3
## 9 9 11.00 8.91 8.95 8.89 8.86 8.9
## 10 10 9.50 9.69 NA 9.72 NA --skip--
## 11 11 8.80 7.26 7.31 7.32 7.29 7.3
## 12 12 8.25 8.16 NA 8.10 7.76 --skip--
## 13 13 7.60 6.51 6.53 6.51 6.52 6.52
## 14 14 7.10 7.43 NA 7.47 7.43 --skip--
## 15 15 6.50 6.10 6.10 6.08 6.09 6.10
## 16 --skip--
## 17 5.25 or something all the way dow
##---- In the no-ldouble case --- this is very different:
## k cutoffs spl low cobs cobs.80
## 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00
## 2 3700.00 NA NA 5890.00 6180.00
## 3 200.00 159.00 159.00 160.00 159.00
## 4 81.00 87.10 87.10 84.80 86.00
## 5 36.00 23.50 23.50 23.50 23.50
## 6 25.00 23.70 23.70 NA NA
## 7 19.00 11.50 11.50 11.50 11.50
## 8 14.00 NA NA 12.90 13.00
## 9 11.00 8.20 8.20 8.21 8.21
## 10 9.50 NA NA NA 9.69
## 11 8.80 6.84 6.84 6.85 6.83
## 12 8.25 NA NA 8.27 7.95
## 13 7.60 NA NA NA NA
## 14 7.10 NA NA 7.43 7.40
## 15 6.50 NA NA NA NA
## M1 macOS is similar to no-ldouble : "==" if equaln
## M1 mac | aarch64-apple-darwin20 | macOS Ventura 13.3.1 | R-devel (2024-01-29 r85841) ; LAPACK version 3.12.0
## k cutoffs spl low cobs cobs.80
## 1 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00 ==
## 2 2 3700.00 NA NA 5900.00 NA ~=
## 3 3 200.00 159.00 159.00 160.00 159.00 ==
## 4 4 81.00 87.10 87.10 84.80 86.00 ==
## 5 5 36.00 23.50 23.50 23.50 23.50
## 6 6 25.00 23.70 23.70 NA NA ==
## 7 7 19.00 11.50 11.50 11.50 11.50
## 8 8 14.00 NA NA 12.90 13.00 ==
## 9 9 11.00 8.20 8.20 8.21 8.21
## 10 10 9.50 NA NA NA 9.69 ==
## 11 11 8.80 6.84 6.84 6.85 6.83
## 12 12 8.25 NA NA 8.27 7.95 ==
## 13 13 7.60 NA NA NA NA ==
## 14 14 7.10 NA NA 7.43 7.40 ==
## 15 15 6.50 NA NA NA NA ==
## ========== Try a slightly better direct formula ======================================================
## after some trial error: gamma(n+1) = n*gamma(n) ==> lgamma(n+1) = lgamma(n) + log(n)
## stirlerrD2 <- function(n) lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - log(2*pi))/2
palette("default")
if(do.pdf) { dev.off(); pdf("stirlerr-tst_others.pdf") }
n <- n5m # (1/64 ... 5000) -- goes with st.nM
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, abs=TRUE)
axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4)
abline(v = 2:6, lty=3, col=cola)
i.n <- 1 <= n & n <= 15 ; cr2 <- adjustcolor(2, 0.6)
lines(n[i.n], abs(relErrV(st.nM[i.n], stirlerr_simpl(n[i.n], "MM2" ))), col=cr2, lwd=2)
legend(20, 1e-13, legend = quote( relErr(stirlerr_simpl(n, '"MM2"'))), col=cr2, lwd=3, bty="n")
n <- seq(1,6, by= 1/200)
stM <- stirlerr(mpfr(n, 512), use.halves = FALSE)
relE <- asNumeric(relErrV(stM, cbind(stD = stirlerr_simpl(n), stD2 = stirlerr_simpl(n, "MM2"))))
signif(apply(abs(relE), 2, summary), 4)
## stD stD2
## Min. 3.093e-17 7.081e-18
## 1st Qu. 2.777e-15 1.936e-15
## Median 8.735e-15 5.238e-15
## Mean 1.670e-14 1.017e-14 .. well "67% better"
## 3rd Qu. 2.380e-14 1.408e-14
## Max. 1.284e-13 9.382e-14
c2 <- adjustcolor(1:2, 0.6)
matplotB(n, abs19(relE), type="o", cex=3/4, log="y", ylim = c(8e-17, 1.3e-13), yaxt="n", col=c2)
eaxis(2)
abline(h = 2^-53, lty=1, col="gray")
smrelE <- apply(abs(relE), 2, \(y) lowess(n, y, f = 0.1)$y)
matlines(n, smrelE, lwd=3, lty=1)
legend("topleft", legend = expression(stirlerr_simpl(n), stirlerr_simple(n, "MM2")),
bty='n', col=c2, lwd=3, lty=1)
drawEps.h(-(53:48))
## should we e.g., use interpolation spline through sfserr_halves[] for n <= 7.5
## -- doing the interpolation on the log(1 - 12*x*stirlerr(x)) vs log2(x) scale -- maybe ?
stirL <- curve(1-12*x*stirlerr(x, cutoffs = cuts, verbose=TRUE), 1/64, 8, log="xy", n=2048)
## just need "true" values for x = 2^-(6,5,4,3,2) in addition to those we already have at x = 1/2, 1.5, 2, 2.5, ..., 7.5, 8
xM <- mpfr(stirL$x, 2048)
stilEM <- roundMpfr(1 - 12*xM*stirlerr(xM, verbose=TRUE), 128)
relEsml <- relErrV(stilEM, stirL$y)
plot(stirL$x, relEsml) # again: stirlerr() is "limited.." for ~ [2, 6]
## The function to interpolate:
plot(asNumeric( (stilEM)) ~ x, data = stirL, type = "l", log="x")
plot(asNumeric(sqrt(stilEM)) ~ x, data = stirL, type = "l", log="x")
plot(asNumeric( log(stilEM)) ~ x, data = stirL, type = "l", log="x")
y <- asNumeric( log(stilEM))
x <- stirL$x
spl <- splinefun(log(x), y)
plot(asNumeric(log(stilEM)) ~ x, data = stirL, type = "l", log="x")
lines(x, spl(log(x)), col=2)
summary(rE <- relErrV(target = y, current = spl(log(x)))) # all 0 {of course: interpolation at *these* x}
ssmpl <- smooth.spline(log(x), y)
str(ssmpl$fit)# only 174 knots, 170 coefficients
methods(class="smooth.spline.fit") #-> ..
str(yp <- predict(ssmpl$fit, x=log(x)))
lines(x, yp$y, col=adjustcolor(3, 1/3), lwd=3) # looks fine
## but
summary(reSpl <- relErrV(target = y, current = yp$y))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.35e-10 -7.55e-11 1.10e-12 0.00e+00 7.32e-11 5.53e-10
## which is *NOT* good enough, of course ....
showProc.time()
Try the DPQ package in your browser
Any scripts or data that you put into this service are public.
DPQ documentation built on Jan. 29, 2025, 3:01 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.
#### Testing stirlerr()
#### =================== {previous 2nd part of this, now -->>> ./bd0-tst.R <<<---
require(DPQ)
for(pkg in c("Rmpfr", "DPQmpfr"))
if(!requireNamespace(pkg)) {
cat("no CRAN package", sQuote(pkg), " ---> no tests here.\n")
q("no")
}
require("Rmpfr")
source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE))
## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS()
##_ options(conflicts.policy = list(depends.ok=TRUE, error=FALSE, warn=FALSE))
require(sfsmisc) # masking 'list_' *and* gmp's factorize(), is.whole()
##_ options(conflicts.policy = NULL)o
## plot1cuts() , etc: ---> ../inst/extraR/relErr-plots.R <<<<<<<
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
source(system.file(package="DPQ", "extraR", "relErr-plots.R", mustWork=TRUE))
do.pdf <- TRUE # (manually)
do.pdf <- !dev.interactive(orNone = TRUE)
do.pdf
if(do.pdf) {
pdf.options(width = 9, height = 6.5) # for all pdf plots {9/6.5 = 1.38 ; 4/3 < 1.38 < sqrt(2) [A4]
pdf("stirlerr-tst.pdf")
}
showProc.time()
(doExtras <- DPQ:::doExtras())
(noLdbl <- (.Machine$sizeof.longdouble <= 8)) ## TRUE when --disable-long-double or 'M1mac' ..
(M.mac <- grepl("aarch64-apple", R.version$platform)) # Mac with M1, M2, ... proc
abs19 <- function(r) pmax(abs(r), 1e-19) # cut |err| to positive {for log-plots}
options(width = 100, nwarnings = 1e5, warnPartialMatchArgs = FALSE)
##=== Really, dpois_raw() and dbinom_raw() *both* use stirlerr(x) for "all" 'x > 0'
## ~~~~~~~~~~~ ~~~~~~~~~~~~ =========== ===== !
## below, 6 "it's okay, but *far* from perfect:" ===> need more terms in stirlerr() [
## April 20: MM added more terms up to S10; 2024-01: up to S12 ..helps a little only
x <- lseq(1/16, 6, length=2048)
system.time(stM <- DPQmpfr::stirlerrM(Rmpfr::mpfr(x,2048))) # 1.7 sec elapsed
plot(x, stirlerr(x, use.halves=FALSE) - stM, type="l", log="x", main="absolute Error")
plot(x, stirlerr(x, use.halves=FALSE) / stM - 1, type="l", log="x", main="relative Error")
plot(x, abs(stirlerr(x, use.halves=FALSE) / stM - 1), type="l", log="xy",main="|relative Error|")
drawEps.h(-(52:50))
## lgammacor() does *NOT* help, as it is *designed* for x >= 10 ... (but is interesting there)
##
## ==> Need another chebyshev() or rational-approx. for x in [.1, 7] or so !!
##=============> see also ../Misc/stirlerr-trms.R <===============
## ~~~~~~~~~~~~~~~~~~~~~~~
cutoffs <- c(15,35,80,500) # cut points, n=*, in the stirlerr() "algorithm"
##
n <- c(seq(1,15, by=1/4),seq(16, 25, by=1/2), 26:30, seq(32,50, by=2), seq(55,1000, by=5),
20*c(51:99), 50*(40:80), 150*(27:48), 500*(15:20))
st.n <- stirlerr(n, "R3")# rather use.halves=TRUE; but here , use.halves=FALSE
plot(st.n ~ n, log="xy", type="b") ## looks good now (straight line descending {in log-log !}
nM <- mpfr(n, 2048)
st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose
all.equal(asNumeric(st.nM), st.n)# TRUE
all.equal(st.nM, as(st.n,"mpfr"))# .. difference: 3.381400e-14 was 1.05884.........e-15
all.equal(roundMpfr(st.nM, 64), as(st.n,"mpfr"), tolerance=1e-16)# (ditto)
## --- Look at the direct formula -- why is it not good for n ~= 5 ?
##
## Preliminary Conclusions :
## 1. there is *some* cancellation even for small n (how much?)
## 2. lgamma1p(n) does really not help much compared to lgamma(n+1) --- but a tiny bit in some cases
### 1. Investigating lgamma1p(n) vs lgamma(n+1) for n < 1 =============================================
##' @title Relative Error of lgamma(n+1) vs lgamma1p() vs MM's stirlerrD2():
##' @param n numeric, typically n << 1
##' @param precBits
##' @return relative error WRT mpfr(n, precBits)
##' @author Martin Maechler
relE.lgam1 <- function(n, precBits = if(doExtras) 1024 else 320) {
M_LN2PI <- 1.837877066409345483 # ~ log(2*Const("pi",60)); very slightly more accurate than log(2*pi)
st <- lgamma(n +1) - (n +0.5)*log(n) + n - M_LN2PI/2
st. <- lgamma1p(n) - (n +0.5)*log(n) + n - M_LN2PI/2 # "lgamma1p"
st2 <- lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - M_LN2PI)/2 # "MM2"
st0 <- -(l.n + M_LN2PI)/2 # "n0"
nM <- mpfr(n, precBits)
stM <- lgamma(nM+1) - (nM+0.5)*log(nM) + nM - log(2*Const("pi", precBits))/2
## stM <- roundMpfr(stM, 128)
cbind("R3" = asNumeric(relErrV(stM, st))
, "lgamma1p"= asNumeric(relErrV(stM, st.))
, "MM2" = asNumeric(relErrV(stM, st2))
, "n0" = asNumeric(relErrV(stM, st0))
)
}
n <- 2^-seq.int(1022, 1, by = -1/4)
relEx <- relE.lgam1(n)
showProc.time()
## Is *equivalent* to 'new' stirlerr_simpl(n version = *) [not for <mpfr> though, see 'relEmat']:
(simpVer <- eval(formals(stirlerr_simpl)$version))
if(is.null(simpVer)) { warning("got wrong old version of package 'DPQ':")
print(packageDescription("DPQ"))
stop("invalid outdated version package 'DPQ'")
}
stir.allS <- function(n) sapply(simpVer, function(v) stirlerr_simpl(n, version=v))
stirS <- stir.allS(n) # matrix
nM <- mpfr(n, 256) # "high" precision = 256 should suffice!
stirM <- stirlerr(nM)
releS <- asNumeric(relErrV(stirM, stirS))
all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux
stopifnot(all.equal(relEx, releS, tolerance = if(noLdbl) 2e-15 else 1e-15))
simpVer3 <- simpVer[simpVer != "lgamma1p"] # have no mpfr-ified lgamma1p()!
## stirlerr_simpl(<mpfr>, *) :
stirM2 <- sapplyMpfr(simpVer3, function(v) stirlerr_simpl(nM, version=v))
## TODO ?:
## apply(stirM2, 2, function(v) all.equal(v, stirS, check.class=FALSE))
## releS2 <- asNumeric(relErrV(stirM2, stirS))
## all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux
## stopifnot(all.equal(relEx, releS, tolerance = 1e-15))
relEmat <- matrix(NA, ncol(stirM2), ncol(stirS),
dimnames = list(simpVer3, colnames(stirS)))
for(j in seq_len(ncol(stirM2)))
for(k in seq_len(ncol(stirS)))
relEmat[j,k] <- asNumeric(relErr(stirM2[,j], stirS[,k]))
relEmat
round(-log10(relEmat), 2) # well .. {why? / expected ?}
cols <- c("gray30", adjustcolor(c(2,3,4), 1/2)); lwd <- c(1, 3,3,3)
stopifnot((k <- length(cols)) == ncol(relEx), k == length(lwd))
matplot(n, relEx, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
mtext("really small errors are dominated by small (< 2^-53) errors of log(n)")
## very interesting: there are different intervals <---> log(n) Qpattern !!
## -- but very small difference, only for n >~= 1/1000 but not before
drawEps.h(negative=TRUE) # abline(h= c(-4,-2:2, 4)*2^-53, lty=c(2,2,2, 1, 2,2,2), col="gray")
legend("topleft", legend = colnames(relEx), col=cols, lwd=3)
## zoomed in a bit:
n. <- 2^-seq.int(400,2, by = -1/4)
relEx. <- relE.lgam1(n.)
matplot(n., relEx., type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
drawEps.h(negative=TRUE)
legend("topleft", legend = colnames(relEx.), col=cols, lwd=3)
##====> Absolute errors (and look at "n0") --------------------------------------
matplot(n., abs19(relEx.), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
drawEps.h(); legend("top", legend = colnames(relEx.), col=cols, lwd=3)
lines(n., abs19(relEx.[,"n0"]), type = "o", cex=1/4, col=cols[4], lwd=2)
## more zooom-in
n.2 <- 2^-seq.int(85, 50, by= -1/100)
stirS.2 <- sapply(c("R3", "lgamma1p", "n0"), function(v) stirlerr_simpl(n.2, version=v))
releS.2 <- asNumeric(relErrV(stirlerr(mpfr(n.2, 320)), stirS.2))
matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3)
abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3)
axis(1, at=5e-17, col.axis=cb, line=-1/4, cex = 3/4)
matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16),
xaxt="n", xlim = c(8e-18, 1e-15), ## <<<<<<<<<<<<<<<<<<< Zoom-in
xlab = quote(n), main = quote(abs(relErr(stirlerr_simpl(n, '*')))))
eaxis(1); drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3)
abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3)
mtext('stirlerr_simpl(*, "n0") is as good as others for n <= 5e-17', col=adjustcolor(cols[3], 2))
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ===> "all *but* "n0" approximations for "larger" small n:
n2 <- 2^-seq.int(20,0, length.out=1000)
relEx2 <- relE.lgam1(n2)[,c("R3", "lgamma1p", "MM2")] # "n0" is "bad": |relE| >= 2.2e-6 !
cols <- c("gray30", adjustcolor(c(2,4), 1/2)); lwd <- c(1, 3,3)
stopifnot((k <- length(cols)) == ncol(relEx2), k == length(lwd))
matplot(n2, relEx2, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-3,3)*1e-15, xaxt="n",
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
eaxis(1, sub10=c(-3,0)); drawEps.h(negative=TRUE)
legend("topleft", legend = colnames(relEx2), col=cols, lwd=3)
##==> "MM" is *NOT* good for n < 1 *but*
## "the same" -- even larger small n:
n3 <- seq(.01, 5, length=1000)
relEx3 <- relE.lgam1(n3)[,c("R3", "lgamma1p", "MM2")] # "no" is "bad" ..
stopifnot((k <- length(cols)) == ncol(relEx3), k == length(lwd))
matplot(n3, relEx3, type = "l", col=cols, lwd=lwd,
main = "relative errors of direct (approx.) formula for stirlerr(n), small n")
legend("topleft", legend = colnames(relEx3), col=cols, lwd=3)
drawEps.h(negative=TRUE)
matplot(n3, abs19(relEx3), type = "l", col=cols, lwd=lwd,
log="y", ylim = 2^-c(54, 44), yaxt = "n", ylab = quote(abs(relE)), xlab=quote(n),
main = "|relative errors| of direct (approx.) formula for stirlerr(n), small n")
eaxis(2, cex.axis=0.9); legend("topleft", legend = colnames(relEx3), col=cols, lwd=3)
drawEps.h()
lines(n3, smooth.spline(abs(relEx3)[,1], df=12)$y, lwd=3, col=cols[1])
lines(n3, smooth.spline(abs(relEx3)[,2], df=12)$y, lwd=3, col=adjustcolor(cols[2], 1/2))
lines(n3, smooth.spline(abs(relEx3)[,3], df=12)$y, lwd=4, col=adjustcolor(cols[3], offset = rep(.2,4)))
## ===> from n >~= 1, "MM2" is definitely better up to n = 5 !!
## Check log() only :
plot(n, asNumeric(relErrV(log(mpfr(n, 256)), log(n))), ylim = c(-1,1)*2^-53,
log="x", type="l", xaxt="n") ## ===> indeed --- log(n) approximation pattern !!
eaxis(1) ; drawEps.h(negative=TRUE)
showProc.time()
## =========== "R3" vs "lgamma1p" -------------------------- which is better?
## really for the very small n, all is dominated by -(n+0.5)*log(n); and lgamma1p() is unnecessary!
i <- 1:20; ni <- n[i]
lgamma1p(ni)
- (ni +0.5)*log(ni) + ni
## much less extreme:
n2 <- lseq(2^-12, 1/2, length=1000)
relE2 <- relE.lgam1(n2)[,-4]
cols <- c("gray30", adjustcolor(2:3, 1/2)); lwd <- c(1,3,3)
matplot(n2, relE2, type = "l", log="x", col=cols, lwd=lwd)
legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3)
drawEps.h(negative=TRUE)
matplot(n2, abs19(relE2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(6e-17, 1e-15),
xaxt = "n"); eaxis(1, sub10=c(-2,0))
legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3)
drawEps.h()
## "MM2" is *worse* here, n < 1/2
for(j in 1:3) lines(n2, smooth.spline(abs(relE2[,j]), df=10)$y, lwd=3,
col=adjustcolor(cols[j], 1.5, offset = rep(-1/4, 4)))
## "lgammap very slightly better in [0.002, 0.05] ...
## "TODO": draw 0.90-quantile curves {--> cobs::cobs() ?} instead of mean-curves?
## which is better? ... "random difference"
d.absrelE <- abs(relE2[,"R3"]) - abs(relE2[,"lgamma1p"])
plot (n2, d.absrelE, type = "l", log="x", # no clear picture ...
main = "|relE_R3| - |relE_lgamma1p|", axes=FALSE, frame.plot=TRUE)
eaxis(1, sub10=c(-2,1)); eaxis(2); axis(3, at=max(n2)); abline(v = max(n2), lty=3, col="gray")
## 'lgamma1p' very slightly better:
lines(n2, smooth.spline(d.absrelE, df=12)$y, lwd=3, col=2)
## not really small n at all == here see, how "bad" the direct formula gets for 1 < n < 10 or so
n3 <- lseq(2^-14, 2^2, length=800)
relE3 <- relE.lgam1(n3)[, -4]
matplot(n3, relE3, type = "l", log="x", col=cols, lty=1, lwd = c(1,3),
main = quote(rel.lgam1(n)), xlab=quote(n))
matplot(n3, abs19(relE3), type = "l", log="xy", col=cols, lwd = c(1,3), xaxt="n",
main = quote(abs(rel.lgam1(n))), xlab=quote(n), ylim = c(2e-17, 4e-14))
drawEps.h(); eaxis(1, sub10=c(-2,3))
legend("topleft", legend=colnames(relE3), col=cols, lwd=2)
## very small difference --- draw the 3 smoothers :
for(j in 1:3) {
ll <- lowess(log(n3), abs19(relE3[, j]), f= 1/12)
with(ll, lines(exp(x), y, col=adjustcolor(cols[j], 1.5), lwd=3))
}
## ==> lgamma1p(.) very slightly in n ~ 10^-4 -- 10^-2 --- but not where it matters: n ~ 0.1 -- 1 !!
## "MM2" gets best from n >~ 1 !
abline(v=1, lty=3, col = adjustcolor(1, 3/4))
showProc.time()
### 2. relErr( stirlerr(.) ) ============================================================
##' Very revealing plot showing the *relative* approximation error of stirlerr(<dblprec>)
##'
p.stirlerrDev <- function(n, precBits = if(doExtras) 2048L else 512L,
stnM = stirlerr(mpfr(n, precBits), use.halves=use.halves, verbose=verbose),
abs = FALSE,
## cut points, n=*, in the stirlerr() algorithm; "FIXME": sync with ../R/dgamma.R <<<<
scheme = c("R3", "R4.4_0"),
cutoffs = switch(match.arg(scheme)
, R3 = c(15, 35, 80, 500)
, R4.4_0 = c(4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.7,
6.1, 6.5, 7, 7.9, 8.75, 10.5, 13,
20, 26, 60, 200, 3300, 17.4e6)
## {FIXME: need to sync} <==> ../man/stirlerr.Rd <==> ../R/dgamma.R
),
use.halves = missing(cutoffs),
direct.ver = c("R3", "lgamma1p", "MM2", "n0"),
verbose = getOption("verbose"),
type = "b", cex = 1,
col = adjustcolor(1, 3/4), colnB = adjustcolor("orange4", 1/3),
log = if(abs) "xy" else "x",
xlim=NULL, ylim = if(abs) c(8e-18, max(abs(N(relE)))))
{
op <- par(las = 1, mgp=c(2, 0.6, 0))
on.exit(par(op))
require("Rmpfr"); require("sfsmisc")
st <- stirlerr(n, scheme=scheme, cutoffs=cutoffs, use.halves=use.halves, direct.ver=direct.ver,
verbose=verbose)
relE <- relErrV(stnM, st) # eps0 = .Machine$double.xmin
N <- asNumeric
form <- if(abs) abs(N(relE)) ~ n else N(relE) ~ n
plot(form, log=log, type=type, cex=cex, col=col, xlim=xlim, ylim=ylim,
ylab = quote(relErrV(stM, st)), axes=FALSE, frame.plot=TRUE,
main = sprintf("stirlerr(n, cutoffs) rel.error [wrt stirlerr(Rmpfr::mpfr(n, %d))]",
precBits))
eaxis(1, sub10=3)
eaxis(2)
mtext(paste("cutoffs =", deparse1(cutoffs)))
ylog <- par("ylog")
## FIXME: improve this ---> drawEps.h() above
if(ylog) {
epsC <- c(1,2,4,8)*2^-52
epsCxp <- expression(epsilon[C],2*epsilon[C], 4*epsilon[C], 8*epsilon[C])
} else {
epsC <- (-2:2)*2^-52
epsCxp <- expression(-2*epsilon[C],-epsilon[C], 0, +epsilon[C], +2*epsilon[C])
}
dy <- diff(par("usr")[3:4])
if(diff(range(if(ylog) log10(epsC) else epsC)) > dy/50) {
lw <- rep(1/2, 5); lw[if(ylog) 1 else 3] <- 2
abline( h=epsC, lty=3, lwd=lw)
axis(4, at=epsC, epsCxp, las=2, cex.axis = 3/4, mgp=c(3/4, 1/4, 0), tck=0)
} else ## only x-axis
abline(h=if(ylog) epsC else 0, lty=3, lwd=2)
abline(v = cutoffs, col=colnB)
axis(3, at=cutoffs, col=colnB, col.axis=colnB,
labels = formatC(cutoffs, digits=3, width=1))
invisible(relE)
} ## p.stirlerrDev()
showProc.time()
n <- lseq(2^-10, 1e10, length=4096)
n <- lseq(2^-10, 5000, length=4096)
## store "expensive" stirlerr() result, and re-use many times below:
nM <- mpfr(n, if(doExtras) 2048 else 512)
st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose
p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE) # default cutoffs= c(15, 40, 85, 600)
p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE, ylim = c(-1,1)*1e-12) # default cutoffs= c(15, 40, 85, 600)
## show the zoom-in region in next plot
yl2 <- 3e-14*c(-1,1)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
if(do.pdf) { dev.off() ; pdf("stirlerr-relErr_1.pdf") }
## drop n < 7:
p.stirlerrDev(n=n, stnM=st.nM, xlim = c(7, max(n)), use.halves=FALSE) # default cutoffs= c(15, 40, 85, 600)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
## The first plot clearly shows we should do better:
## Current code is switching to less terms too early, loosing up to 2 decimals precision
if(FALSE) # no visible difference {use.halves = T / F }:
p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = FALSE)
p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = TRUE)# exact at n/2 (n <= ..)
abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2)
showProc.time()
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-1.pdf") }
### ~19.April 2021: "This is close to *the* solution" (but see 'cuts' below)
cuts <- c(7, 12, 20, 26, 60, 200, 3300)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
st. <- stirlerr(n=n , cutoffs = cuts, verbose=TRUE)
st.nM <- stirlerr(n=nM, cutoffs = cuts, use.halves=FALSE) ## << on purpose
relE <- asNumeric(relErrV(st.nM, st.))
head(cbind(n, relE), 20)
## nice printout :
print(cbind(n = format(n, drop0trailing = TRUE),
stirlerr= format(st.,scientific=FALSE, digits=4),
relErr = signif(relE, 4))
, quote=FALSE)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts)
## and zoom in:
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2/20)
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-2.pdf") }
## zoom in ==> {good for n >= 10}
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", ylim = 2e-15*c(-1,1),
cutoffs = cuts)## old default cutoffs = c(15,35, 80, 500)
if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-3.pdf") }
showProc.time()
##-- April 20: have more terms up to S10 in stirlerr() --> can use more cutoffs
n <- n5m <- lseq(1/64, 5000, length=4096)
nM <- mpfr(n, if(doExtras) 2048L # a *lot* accuracy for stirlerr(nM,*)
else 512L)
ct10.1 <- c( 5.4, 7.5, 8.5, 10.625, 12.125, 20, 26, 60, 200, 3300)# till 2024-01-19
ct10.2 <- c( 5.4, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300)
cuts <-
ct12.1 <- c(5.22, 6.5, 7.0, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## 5.25 is "too small" but the direct formula is already really bad there, ...
st.nM <- roundMpfr(stirlerr(nM, use.halves=FALSE, ## << on purpose;
verbose=TRUE), precBits = 128)
## NB: for x=xM <mpfr>; `cutoffs` are *not* used.
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0")
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*3e-15)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*1e-15)
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", abs=TRUE)
axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4)
abline(v = 2:6, lty=3, col=cola)
if(FALSE)## using exact values sferr_halves[] *instead* of MPFR ones: ==> confirmation they lay on top
lines((0:30)/2, abs(stirlerr((0:30)/2, cutoffs=cuts, verbose=TRUE)/DPQ:::sferr_halves - 1), type="o", col=2,lwd=2)
if(FALSE) ## nice (but unneeded) printout :
print(cbind(n = format(n, drop0trailing = TRUE),
stirlerr= format(st.,scientific=FALSE, digits=4),
relErr = signif(relE, 4))
, quote=FALSE)
showProc.time()
## ========== where should the cutoffs be ? ===================================================
.stirl.cutoffs <- function(scheme)
eval(do.call(substitute, list(formals(stirlerr)$cutoffs, list(scheme = scheme))))
drawCuts <- function(scheme, axis=NA, lty = 3, col = "skyblue", ...) {
abline(v = (ct <- .stirl.cutoffs(scheme)), lty=lty, col=col, ...)
if(is.finite(axis)) axisCuts(side = axis, at = ct, col=col, ...)
}
axisCuts <- function(scheme, side = 3, at = .stirl.cutoffs(scheme), col = "skyblue", line = -3/4, ...)
axis(side, at=at, labels=formatC(at), col.axis = col, col=NA, col.ticks=NA, line=line, ...)
mtextCuts <- function(cutoffs, scheme, ...) {
if(!missing(scheme)) cutoffs <- .stirl.cutoffs(scheme)
mtext(paste("cutoffs =", deparse1(cutoffs)), ...)
}
if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k.pdf") }
mK <- 20L # := max(k)
## order = k = 1:mK terms in series approx:
k <- 1:mK
n <- 2^seq(1, 28, by=1/16)
nM <- mpfr(n, 1024)
stnM <- stirlerr(nM) # the "true" values
stirlOrd <- sapply(k, function(k.) stirlerr(n, order = k.))
stirlO_lgcor <- cbind(stirlOrd, sapply(5:6, function(nal) lgammacor(n, nalgm = nal)))
relE <- asNumeric(stirlO_lgcor/stnM -1) # "true" relative error
## use a "smooth" but well visible polette :
palROBG <- colorRampPalette(c("red", "darkorange2", "blue", "seagreen"), space = "Lab")
palette(adjustcolor(palROBG(mK+2), 3/4))
## -- 2 lgamcor()'s
(tit.k <- substitute(list( stirlerr(n, order=k) ~~"error", k == 1:mK), list(mK = mK)))
(tit.kA <- substitute(list(abs(stirlerr(n, order=k) ~~"error"), k == 1:mK), list(mK = mK)))
lgammacorTit <- function(...) mtext("+ lgammacor(x, 5) [bad] + lgammacor(x, 6) [good]", col=1, ...)
matplotB(n, relE, cex=2/3, ylim = c(-1,1)*1e-13, col=k,
log = "x", xaxt="n", main = tit.k)
lgammacorTit()
eaxis(1, nintLog = 20)
drawCuts("R4.4_0")
## zoom in (ylim)
matplotB(n, relE, cex=2/3, ylim = c(-1,1)*5e-15, col=k,
log = "x", xaxt="n", main = tit.k)
lgammacorTit()
eaxis(1, nintLog = 20); abline(h = (-2:2)*2^-53, lty=3, lwd=1/2)
drawCuts("R4.4_0", axis = 3)
## log-log |rel.Err| -- "linear"
matplotB(n, abs19(relE), cex=2/3, col=k, ylim = c(8e-17, 1e-7), log = "xy", main=tit.kA)
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
lgammacorTit(line=-1)
drawCuts("R4.4_0", axis = 3)
## zoom in -- still "large" {no longer carry the lgammacor() .. }:
n2c <- 2^seq(2, 8, by=1/256)
nMc <- mpfr(n2c, 1024)
stnMc <- stirlerr(nMc) # the "true" values
stirlOrc <- sapply(k, function(k.) stirlerr(n2c, order = k.))
relEc <- asNumeric(stirlOrc/stnMc -1) # "true" relative error
matplotB(n2c, relEc, cex=2/3, ylim = c(-1,1)*1e-13, col=k,
log = "x", xaxt="n", main = tit.k)
eaxis(1, sub10 = 2)
drawCuts("R4.4_0", axis=3)
## log-log |rel.Err| -- "linear"
matplotB(n2c, abs19(relEc), cex=2/3, col=k, ylim = c(8e-17, 1e-3), log = "xy", main=tit.kA)
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
drawCuts("R4.4_0", axis = 3)
## zoom into the critical n region
nc <- seq(3.5, 11, by=1/128)
ncM <- mpfr(nc, 256)
stncM <- stirlerr(ncM) # the "true" values
stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k))
relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error
## log-log |rel.Err| -- "linear"
matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(2e-17, 1e-8),
log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2)
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6))
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4))
## lines(nc, abs19(asNumeric(stirlerr_simpl(nc,"MM2")/stncM - 1)), lwd=3, col=adjustcolor("purple", 2/3))
legend(10^par("usr")[1], 1e-9, legend=paste0("k=", k), bty="n", lwd=2,
col=k, lty=1:5, pch= c(1L:9L, 0L, letters)[seq_along(k)])
drawCuts("R4.4_0", axis=3)
## Zoom-in [only]
matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(4e-17, 1e-11), xlim = c(4.8, 6.5),
log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2)
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6))
lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4))
k. <- k[-(1:6)]
legend("bottomleft", legend=paste0("k=", k.), bty="n", lwd=2,
col=k., lty=1:5, pch= c(1L:9L, 0L, letters)[k.])
drawCuts("R4.4_0", axis=3)
showProc.time()
##--- Accuracy of "R4.4_0" -------------------------------------------------------
for(nc in list(seq(4.75, 28, by=1/512), # for a bigger pix
seq(4.75, 9, by=1/1024)))
{
ncM <- mpfr(nc, 1024)
stncM <- stirlerr(ncM) # the "true" values
stirl.440 <- stirlerr(nc, scheme = "R4.4_0")
stirl.3 <- stirlerr(nc, scheme = "R3")
relE440 <- asNumeric(relErrV(stncM, stirl.440))
relE3 <- asNumeric(relErrV(stncM, stirl.3 ))
##
plot(nc, abs19(relE440), xlab=quote(n), main = quote(abs(relErr(stirlerr(n, '"R4.4_0"')))),
type = "l", log = "xy", ylim = c(4e-17, 1e-13))
mtextCuts(scheme="R4.4_0", cex=4/5)
drawCuts("R4.4_0", lty=2, lwd=2, axis=4)
drawEps.h()
if(max(nc) <= 10) abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(4, 1/2))
if(TRUE) { # but just so ...
c3 <- adjustcolor("royalblue", 1/2)
lines(nc, pmax(abs(relE3), 1e-18), col=c3)
title(quote(abs(relErr(stirlerr(n, '"R3"')))), adj=1, col.main = c3)
drawCuts("R3", lty=4, col=c3); mtextCuts(scheme="R3", adj=1, col=c3)
}
addOrd <- TRUE
addOrd <- dev.interactive(orNone=TRUE)
if(addOrd) {
if(max(nc) >= 100) {
i <- (15 <= nc & nc <= 85) # (no-op in [4.75, 7] !)
ni <- nc[i]
} else { i <- TRUE; ni <- nc }
for(k in 8:17) lines(ni, abs19(asNumeric(relErrV(stncM[i], stirlerr(ni, order=k)))), col=adjustcolor(k, 1/3))
title(sub = "stirlerr(*, order k = 8:17)")
}
} ## for(nc ..)
if(FALSE)
lines(nc, abs19(relE440))# *re* draw!
showProc.time()
palette("Tableau")
## Focus more:
stirlerrPlot <- function(nc, k, res=NULL, legend.xy = "left", full=TRUE, precB = 1024,
ylim = c(3e-17, 2e-13), cex = 5/4) {
stopifnot(require("Rmpfr"), require("graphics"))
if(is.list(res) && all(c("nc", "k", "relEc","splarelE") %in% names(res))) { ## do *not* recompute
list2env(res, envir = environment())
} else { ## compute
stopifnot(is.finite(nc), nc > 0, length(nc) >= 100, k == as.integer(k), 0 <= k, k <= 20)
ncM <- mpfr(nc, precB)
stncM <- stirlerr(ncM) # the "true" values
stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k))
relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error
## log |rel.Err| -- "linear"
## smooth() on log-scale {and transform back}:
splarelE <- apply(log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y))
## the direct formulas (default "R3", "MM2"):
arelEs0 <- abs(asNumeric(stirlerr_simpl(nc )/stncM - 1))
arelEs2 <- abs(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1))
}
pch. <- c(1L:9L, 0L, letters)[k]
if(full)
matplotB(nc, abs19(relEc), col=k, pch = pch., cex=cex, ylim=ylim, log = "y",
xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
else ## smooth only
matplotB(nc, splarelE, col=adjustcolor(k,2/3), pch=pch., lwd=2, cex=cex, ylim=ylim, log = "y",
xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k)))))
mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2)
legend(legend.xy, legend=paste0("k=", k), bty="n", lwd=2, col=k, lty=1:5, pch = pch.)
abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(10, 1/2))
drawCuts("R4.4_0", axis=3)
if(full) {
matlines(nc, splarelE, col=adjustcolor(k,2/3), lwd=4)
lines(nc, pmax(arelEs0, 1e-19), lwd=1.5, col=adjustcolor( 2, 0.2))
lines(nc, pmax(arelEs2, 1e-19), lwd=1, col=adjustcolor(10, 0.2))
}
lines(nc, smooth.spline(arelEs0, df=12)$y, lwd=3, col= adjustcolor( 2, 1/2))
lines(nc, smooth.spline(arelEs2, df=12)$y, lwd=3, col= adjustcolor(10, 1/2))
invisible(list(nc=nc, k=k, relEc = relEc, splarelE = splarelE, arelEs0=arelEs0, arelEs2=arelEs2))
}
rr1 <- stirlerrPlot(nc = seq(4.75, 9.0, by=1/1024),
k = 7:20)
stirlerrPlot(res = rr1, full=FALSE, ylim = c(8e-17, 1e-13))
if(interactive())
stirlerrPlot(res = rr1)
rr <- stirlerrPlot(nc = seq(5, 6.25, by=1/2048), k = 9:18)
stirlerrPlot(res = rr, full=FALSE, ylim = c(8e-17, 1e-13))
showProc.time()
palette("default")
if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k-vs-k1.pdf") }
##' Find 'cuts', i.e., a region c(k) +/- s(k) i.e. intervals [c(k) - s(k), c(k) + s(k)]
##' where c(k) is such that relE(n=c(k), k) ~= eps)
##' 1. Find the c1(k) such that |relE(n, k)| ~= c1(k) * n^{-2k}
findC1 <- function(n, ks, e1 = 1e-15, e2 = 1e-5, res=NULL, precBits = 1024, do.plot = TRUE, ...)
{
if(is.list(res) && all(c("n", "ks", "arelE") %in% names(res))) { ## do *not* recompute, take from 'res':
list2env(res, envir = environment())
} else { ## compute
stopifnot(require("Rmpfr"),
is.numeric(ks), ks == (k. <- as.integer(ks)), length(ks <- k.) >= 1,
length(e1) == 1L, length(e2) == 1L, is.finite(c(e1,e2)), e1 >= 0, e2 >= e1,
0 <= ks, ks <= 20, is.numeric(n), n > 0, is.finite(n), length(n) >= 100)
nM <- mpfr(n, precBits)
stirM <- stirlerr(nM) # the "true" values
stirlOrd <- sapply(ks, function(k) stirlerr(n, order = k))
arelE <- abs(asNumeric(stirlOrd/stirM -1)) # "true" relative error
}
arelE19 <- pmax(arelE, 1e-19)
## log |rel.Err| -- "linear"
## on log-scale {and transform back; for linear fit, only use values inside [e1, e2]
if(do.plot) # experi
matplotB(n, arelE19, log="xy", ...)
## matplot(n, arelE19, type="l", log="xy", xlim = c(min(n), 20))
## re-compute these, as they *also* depend on (e1, e2)
c1 <- vapply(seq_along(ks), function(i) {
k <- ks[i]
y <- arelE19[,i]
iUse <- e1 <= y & y <= e2
if(sum(iUse) < 10) stop("only", sum(iUse), "values in [e1,e2]")
## .lm.fit(cbind(1, log(n[iUse])), log(y[iUse]))$coefficients
## rather, we *know* the error is c* n^{-2k} , i.e.,
## log |relE| = log(c) - 2k * log(n)
## <==> c = exp( log|relE| + 2k * log(n))
exp(mean(log(y[iUse]) + 2*k * log(n[iUse])))
}, numeric(1))
if(do.plot) {
drawEps.h()
for(i in seq_along(ks))
lines(n, c1[i] * n^(-2*ks[i]), col=adjustcolor(i, 1/3), lwd = 4, lty = 2)
}
invisible(list(n=n, ks=ks, arelE = arelE, c1 = c1))
} ## findC1()
c1.Res <- findC1(n = 2^seq(2, 26, by=1/128), ks = 1:18)
(s.c1.fil <- paste0("stirlerr-c1Res-", myPlatform(), ".rds"))
saveRDS(c1.Res, file = s.c1.fil)
if(!exists("c1.Res")) {
c1.Res <- readRDS(s.c1.fil)
## re-do the "default" plot of findC1():
findC1(res = c1.Res, xaxt="n"); eaxis(1, sub10=2)
}
## the same, zoomed in:
findC1(res = c1.Res, xlim = c(4, 40), ylim = c(2e-17, 1e-12))
ks <- c1.Res$ks; pch. <- c(1L:9L, 0L, letters)[ks]
legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.)
## smaller set : larger e1 :
c1.r2 <- findC1(res = c1.Res, xlim = c(4, 30), ylim = c(4e-17, 1e-13), e1 = 4e-15)
legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.)
print(digits = 4,
cbind(ks, c1. = c1.Res$c1, c1.2 = c1.r2$c1,
relD = round(relErrV(c1.r2$c1, c1.Res$c1), 4)))
c1.. <- c1.r2[c("ks", "c1")] # just the smallest part is needed here:
(s.c1.fil <- paste0("stirlerr-c1r2-", myPlatform(), ".rds"))
saveRDS(c1.., file = s.c1.fil) # was "stirlerr-c1.rds"
## 2. Now, find the n(k) +/- se.n(k) intervals
## Use these c1 from above
##' Given c1-results relErr(n,k); |relE(n,k) ~= c1 * n^{-2k} , find n such that
##' |relE(n,k)| line {in log-log scale} cuts y = eps, i.e., n* such that |relE(n,k)| <= eps for all n >= n*
n.ep <- function(eps, c1Res, ks = c1Res$ks, c1 = c1Res$c1, ...) {
stopifnot(is.finite(eps), length(eps) == 1L, eps > 0,
length(ks) == length(c1), is.numeric(c1), is.integer(ks), ks >= 1)
## n: given k, the location where the |relE(n,k)| line {in log-log} cuts y = eps
## |relE(n,k) ~= c1 * n^{-2k} <==>
## log|relE(n,k)| ~= log(c1) - 2k* log(n) <==>
## c := mean{ exp( log|relE(n,k)| + 2k* log(n) ) } ------- see findC1()
## now, solve for n :
## c1 * n^{-2k} == eps
## log(c1) - 2k* log(n) == log(eps)
## log(n) == (log(eps) - log(c1)) / (-2k) <==>
## n == exp((log(c1) - log(eps)) / 2k)
exp((log(c1) - log(eps))/(2*ks))
}
## get c1..
if(!exists("c1..")) c1.. <- readRDS("stirlerr-c1.rds")
ne2 <- n.ep(2^-51, c1Res = c1..) ## ok
ne1 <- n.ep(2^-52, c1Res = c1..)
ne. <- n.ep(2^-53, c1Res = c1..)
form <- function(n) format(signif(n, 3), scientific=FALSE)
data.frame(k = ks, ne2 = form(ne2), ne1 = form(ne1), ne. = form(ne.),
cutoffs = form(rev(.stirl.cutoffs("R4.4_0")[-1])))
## ------- Linux F 36/38 x86_64 (nb-mm5|v-lynne)
## k ne2 ne1 ne. cutoffs
## 1 8660000.00 12200000.00 17300000.00 17400000.00
## 2 2150.00 2560.00 3040.00 3700.00
## 3 159.00 178.00 200.00 200.00
## 4 46.60 50.80 55.40 81.00
## 5 23.40 25.10 26.90 36.00
## 6 15.20 16.10 17.10 25.00
## 7 11.50 12.10 12.70 19.00
## 8 9.43 9.85 10.30 14.00
## 9 8.20 8.53 8.86 11.00
## 10 7.42 7.68 7.95 9.50
## 11 6.89 7.11 7.34 8.80
## 12 6.53 6.72 6.92 8.25
## 13 6.27 6.44 6.62 7.60
## 14 6.10 6.25 6.41 7.10
## 15 5.98 6.12 6.26 6.50 * (not used)
## 16 5.90 6.03 6.16 6.50 * " "
## 17 5.85 5.97 6.10 6.50 * " "
## 18 5.83 5.95 6.06 6.50 << used all the way down to 5.25
## ok --- correct order of magnitude ! --- good!
## 2b. find *interval* around the 'n(eps)' values
## -- Try simply
d.k <- ne. - ne1
## interval
int.k <- cbind(ne1 - d.k,
ne1 + d.k)
## look at e.g.
data.frame(k=ks, `n(k)` = form(ne1), int = form(int.k))
## k n.k. int.1 int.2
## 1 12200000.00 7180000.00 17300000.00
## 2 2560.00 2070.00 3040.00
## 3 178.00 156.00 200.00
## 4 50.80 46.20 55.40
## 5 25.10 23.30 26.90
## 6 16.10 15.20 17.10
## 7 12.10 11.40 12.70
## 8 9.85 9.41 10.30
## 9 8.53 8.19 8.86
## 10 7.68 7.41 7.95
## 11 7.11 6.88 7.34
## 12 6.72 6.52 6.92
## 13 6.44 6.27 6.62
## 14 6.25 6.10 6.41
## 15 6.12 5.98 6.26
## 16 6.03 5.90 6.16
## 17 5.97 5.85 6.10
## 18 5.95 5.83 6.06
##' as function {well, *not* computing c1.k from scratch
nInt <- function(k, c1.k, ep12 = 2^-(52:53)) {
if(length(k) == 1L) { # special convention to call for *one* k, with c1.k vector
stopifnot(k == (k <- as.integer(k)), k >= 1, length(c1.k) >= k)
c1.k <- c1.k[k]
}
## see n.ep() above
n_ <- function(eps, k, c1) {
stopifnot(is.finite(eps), length(eps) == 1L, eps > 0,
length(k) == length(c1), is.numeric(c1), is.integer(k), k >= 1)
exp((log(c1) - log(eps))/(2*k))
}
ne1 <- n_(ep12[1], k, c1.k)
ne. <- n_(ep12[2], k, c1.k)
d.k <- ne. - ne1
stopifnot(d.k > 0)
## interval: {"fudge" 0.5 / 2.5} from results -- also 'noLdbl' gives *quite* different pic!:
odd <- k %% 2 == 1
cbind(ne1 - ifelse( odd & noLdbl, 1.5, 0.5) * d.k,
ne1 + ifelse(!odd , 8, 2.5) * d.k)
}
nInt(k= 1, c1..$c1)
nInt(k= 2, c1..$c1)
nInt(k=18, c1..$c1)
nints.k <- nInt(ks, c1..$c1)
## for printing
form(as.data.frame( nints.k ))
## (-.5 , +2.5) ## originally ( -1, +1)
## 1 9710000.00 24900000.00 # 7180000.00 17300000.00
## 2 2320.00 3770.00 # 2070.00 3040.00
## 3 167.00 232.00 # 156.00 200.00
## 4 48.50 62.30 # 46.20 55.40
## 5 24.20 29.60 # 23.30 26.90
## 6 15.70 18.50 # 15.20 17.10
## 7 11.70 13.60 # 11.40 12.70
## 8 9.63 10.90 # 9.41 10.30
## 9 8.36 9.36 # 8.19 8.86
## 10 7.54 8.36 # 7.41 7.95
## 11 7.00 7.68 # 6.88 7.34
## 12 6.62 7.21 # 6.52 6.92
## 13 6.36 6.88 # 6.27 6.62
## 14 6.17 6.64 # 6.10 6.41
## 15 6.05 6.48 # 5.98 6.26
## 16 5.96 6.36 # 5.90 6.16
## 17 5.91 6.28 # 5.85 6.10
## 18 5.89 6.24 # 5.83 6.06
## 3. Then for each of the intervals, compare order k vs k+1
## -- ==> optimal cutoff { how much platform dependency ?? }
### Here, compute only
find1cuts <- function(k, c1,
n = nInt(k, c1), # the *set* of n's or the 'range'
len.n = 1000,
precBits = 1024, nM = mpfr(n, precBits),
stnM = stirlerr(nM),
stirlOrd = sapply(k+(0:1), function(.k.) stirlerr(n, order = .k.)),
relE = asNumeric(stirlOrd/stnM -1), # "true" relative error for the {k, k+1}
do.spl=TRUE, df.spline = 9, # df = 5 gives 2 cutpoints for k==1
do.low=TRUE, f.lowess = 0.2,
do.cobs = require("cobs", quietly=TRUE), tau = 0.90)
{
## check relErrV( stirlerr(n, order=k ) vs
## stirlerr(n, order=k+1)
if(length(n) == 2L)
if(n[1] < n[2]) n <- seq(n[1], n[2], length.out = len.n)
else stop("'n' must be *increasing")
force(relE)
y <- abs19(relE)
## NB: all smoothing --- as in stirlerrPlot() above -- should happen in log-space
## (log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y))
ly <- log(y) # == log(abs19(relE)) == log(max(|r|, 1e-19))
if(do.spl) {##
s1 <- exp(smooth.spline(ly[,1], df=df.spline)$y)
s2 <- exp(smooth.spline(ly[,2], df=df.spline)$y)
}
if(do.low) { ## lowess
s1l <- exp(lowess(ly[,1], f=f.lowess)$y)
s2l <- exp(lowess(ly[,2], f=f.lowess)$y)
}
## also use cobs() splines for the 90% quantile !!
EE <- environment() # so can set do.cobs to FALSE in case of error
if(do.cobs) { ## <==> require("cobs") # yes, this is in tests/
cobsF <- function(Y) cobs(n, Y, tau=tau, nknots = 6, lambda = -1,
print.warn=FALSE, print.mesg=FALSE)
## sparseM::chol(<matrix.csr>) now gives error when it gave warning {about singularity}
cobsF <- function(Y) {
r <- tryCatch(cobs(n, Y, tau=tau, nknots = 6, lambda = -1,
print.warn=FALSE, print.mesg=FALSE),
error = identity)
if(inherits(r, "error")) {
assign("do.cobs", FALSE, envir = EE) # and return
list(fitted = FALSE)
}
else
r
}
cs1 <- exp(cobsF(ly[,1])$fitted)
cs2 <- exp(cobsF(ly[,2])$fitted)
}
smooths <- list(spl = if(do.spl ) cbind(s1, s2 ),
low = if(do.low ) cbind(s1l,s2l),
cobs= if(do.cobs) cbind(cs1,cs2))
## diffL <- list(spl = if(do.spl ) s2 -s1,
## low = if(do.low ) s2l-s1l,
## cobs= if(do.cobs) cs2-cs1)
### FIXME: simplification does not always *work* -- (i, n.) are not always ok
### ------ notably within R-devel-no-ldouble {hence probably macOS M1 ... ..}
sapply(smooths, function(s12) { d <- s12[,2] - s12[,1]
## typically a (almost or completely) montone increasing function, crossing zero *once*
## compute cutpoint:
## i := the first n[i] with d(n[i]) >= 0
i <- which(d >= 0)[1]
if(length(i) == 1L && !is.na(i) && i > 1L) {
i_ <- i - 1L # ==> d(n[i_]) < 0
## cutpoint must be in [n[i_], n[i]] --- do linear interpolation
n. <- n[i_] - (n[i] - n[i_])* d[i_] / (d[i] - d[i_])
} else {
if(length(i) != 1L) i <- -length(i)
n. <- NA_integer_
}
c(i=i, n.=n.)
}) -> n.L
list(k=k, n=n, relE = unname(relE), smooths=smooths, i.n = n.L)
} ## find1cuts()
k. <- 1:15
system.time(
## Failed on lynne [2024-06-04, R 4.4.1 beta] with
## Error in .local(x, ...) : insufficient space ---> SparseM :: chol(<..>)
## ==> now we catch this inside find1cuts():
resL <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1))
) ## -- warnings, notably from cobs() not converging
## needs 12 sec (!!) user system elapsed =
(s.find15.fil <- paste0("stirlerr-find1_1-15_", myPlatform(), ".rds"))
## now we catch cobs() errors {from SparseM::chol},
## okCuts <- !inherits(resL, "error")
## if(okCuts) {
saveRDS(resL, file = s.find15.fil) # was "stirlerr-find1_1-15.rds"
## } else traceback()
## 11: stop(mess)
## 10: .local(x, ...)
## 9: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax)
## 8: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax)
## 7: rq.fit.sfnc(Xeq, Yeq, Xieq, Yieq, tau = tau, rhs = rhs, control = rqCtrl)
## 6: drqssbc2(x, y, w, pw = pw, knots = knots, degree = degree, Tlambda = if (select.lambda) lambdaSet else lambda,
## constraint = constraint, ptConstr = ptConstr, maxiter = maxiter,
## trace = trace - 1, nrq, nl1, neqc, niqc, nvar, tau = tau,
## select.lambda = select.lambda, give.pseudo.x = keep.x.ps,
## rq.tol = rq.tol, tol.0res = tol.0res, print.warn = print.warn)
## 5: cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE,
## print.mesg = FALSE) at stirlerr-tst.R!udBzpT#32
## 4: cobsF(ly[, 1]) at stirlerr-tst.R!udBzpT#34
## 3: find1cuts(k = k, c1 = c1..$c1) at #1
## 2: FUN(X[[i]], ...)
## 1: lapply(setNames(, k.), function(k) find1cuts(k = k, c1 = c1..$c1))
ok1cutsLst <- function(res) {
stopifnot(is.list(res), sapply(res, is.list)) # must be list of lists
cobsL <- lapply(lapply(res, `[[`, "smooths"), `[[`, "cobs")
vapply(cobsL, is.array, NA)
}
(resLok <- ok1cutsLst(resL))
## 1 2 3 4 5 6 ...... 15
## FALSE TRUE TRUE TRUE TRUE TRUE ...... TRUE
if(FALSE) {
## e.g. in R-devel-no-ldouble:
(r1 <- find1cuts(k=1, c1=c1..$c1))$i.n # list --- no longer ok {SparseM::chol -> insufficient space}
(r2 <- find1cuts(k=2, c1=c1..$c1))$i.n # "good"
(r3 <- find1cuts(k=3, c1=c1..$c1))$i.n # 3-vector: only 3x 'i' == 1
}
## if(okCuts) {
mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1")$old.par -> opar
invisible(lapply(resL[resLok], plot1cuts))
## plus the "last" ones {also showing that k=15 is worse here anyway than k=17}
str(r17 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = 1500), c1=c1..$c1))
plot1cuts(r17) # no-ldouble is *very* different than normal: k *much better* than k+1
(s.find17.fil <- paste0("stirlerr-find1_17_", myPlatform(), ".rds"))
saveRDS(r17, file = s.find17.fil) # was "stirlerr-find1_17.rds"
## }
## ditto for tau = 0.8 {quantile for cobs()}:
system.time(
resL.8 <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1, tau = 0.80))
) # warnings from cobs {again 12.1 sec}
(resL8ok <- ok1cutsLst(resL.8))
mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1 -- tau = 0.8")
invisible(lapply(resL.8[resL8ok], plot1cuts))
## plus "last" one
str(r17.8 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = 1500), c1=c1..$c1, tau = 0.80))
plot1cuts(r17.8)
par(opar)
n.mn <- sapply(resL, `[[`, "i.n", simplify = "array")
n.mn8 <- sapply(resL.8, `[[`, "i.n", simplify = "array")
## of course, only the cobs part differs:
## ... when I later see errors here ... what's going on??
str(n.mn) # all NULL for "no-double" !!
dim(n.mn["n.",,])
str(n.mn8["n.", "cobs",])
sessionInfo()
form(data.frame(k = k., cutoffs = rev(.stirl.cutoffs("R4.4_0")[-1])[k.],
t(n.mn ["n.",,]), cobs.80 = n.mn8["n.", "cobs",]))
## k cutoffs spl low cobs cobs.80
## 1 1 17400000.00 15600000.00 15800000.00 15800000.00 15700000.00
## 2 2 3700.00 NA NA NA NA == from tau=0.90 I'd use ~ 5000
## 3 3 200.00 200.00 201.00 206.00 205.00 205
## 4 4 81.00 NA NA 86.60 86.20 86
## 5 5 36.00 27.20 27.10 27.00 27.10 27
## 6 6 25.00 23.50 NA NA NA 23.5
## 7 7 19.00 12.80 12.80 12.80 12.80 12.8
## 8 8 14.00 NA NA 12.30 12.80 12.3
## 9 9 11.00 8.91 8.95 8.89 8.86 8.9
## 10 10 9.50 9.69 NA 9.72 NA --skip--
## 11 11 8.80 7.26 7.31 7.32 7.29 7.3
## 12 12 8.25 8.16 NA 8.10 7.76 --skip--
## 13 13 7.60 6.51 6.53 6.51 6.52 6.52
## 14 14 7.10 7.43 NA 7.47 7.43 --skip--
## 15 15 6.50 6.10 6.10 6.08 6.09 6.10
## 16 --skip--
## 17 5.25 or something all the way dow
##---- In the no-ldouble case --- this is very different:
## k cutoffs spl low cobs cobs.80
## 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00
## 2 3700.00 NA NA 5890.00 6180.00
## 3 200.00 159.00 159.00 160.00 159.00
## 4 81.00 87.10 87.10 84.80 86.00
## 5 36.00 23.50 23.50 23.50 23.50
## 6 25.00 23.70 23.70 NA NA
## 7 19.00 11.50 11.50 11.50 11.50
## 8 14.00 NA NA 12.90 13.00
## 9 11.00 8.20 8.20 8.21 8.21
## 10 9.50 NA NA NA 9.69
## 11 8.80 6.84 6.84 6.85 6.83
## 12 8.25 NA NA 8.27 7.95
## 13 7.60 NA NA NA NA
## 14 7.10 NA NA 7.43 7.40
## 15 6.50 NA NA NA NA
## M1 macOS is similar to no-ldouble : "==" if equaln
## M1 mac | aarch64-apple-darwin20 | macOS Ventura 13.3.1 | R-devel (2024-01-29 r85841) ; LAPACK version 3.12.0
## k cutoffs spl low cobs cobs.80
## 1 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00 ==
## 2 2 3700.00 NA NA 5900.00 NA ~=
## 3 3 200.00 159.00 159.00 160.00 159.00 ==
## 4 4 81.00 87.10 87.10 84.80 86.00 ==
## 5 5 36.00 23.50 23.50 23.50 23.50
## 6 6 25.00 23.70 23.70 NA NA ==
## 7 7 19.00 11.50 11.50 11.50 11.50
## 8 8 14.00 NA NA 12.90 13.00 ==
## 9 9 11.00 8.20 8.20 8.21 8.21
## 10 10 9.50 NA NA NA 9.69 ==
## 11 11 8.80 6.84 6.84 6.85 6.83
## 12 12 8.25 NA NA 8.27 7.95 ==
## 13 13 7.60 NA NA NA NA ==
## 14 14 7.10 NA NA 7.43 7.40 ==
## 15 15 6.50 NA NA NA NA ==
## ========== Try a slightly better direct formula ======================================================
## after some trial error: gamma(n+1) = n*gamma(n) ==> lgamma(n+1) = lgamma(n) + log(n)
## stirlerrD2 <- function(n) lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - log(2*pi))/2
palette("default")
if(do.pdf) { dev.off(); pdf("stirlerr-tst_others.pdf") }
n <- n5m # (1/64 ... 5000) -- goes with st.nM
p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, abs=TRUE)
axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4)
abline(v = 2:6, lty=3, col=cola)
i.n <- 1 <= n & n <= 15 ; cr2 <- adjustcolor(2, 0.6)
lines(n[i.n], abs(relErrV(st.nM[i.n], stirlerr_simpl(n[i.n], "MM2" ))), col=cr2, lwd=2)
legend(20, 1e-13, legend = quote( relErr(stirlerr_simpl(n, '"MM2"'))), col=cr2, lwd=3, bty="n")
n <- seq(1,6, by= 1/200)
stM <- stirlerr(mpfr(n, 512), use.halves = FALSE)
relE <- asNumeric(relErrV(stM, cbind(stD = stirlerr_simpl(n), stD2 = stirlerr_simpl(n, "MM2"))))
signif(apply(abs(relE), 2, summary), 4)
## stD stD2
## Min. 3.093e-17 7.081e-18
## 1st Qu. 2.777e-15 1.936e-15
## Median 8.735e-15 5.238e-15
## Mean 1.670e-14 1.017e-14 .. well "67% better"
## 3rd Qu. 2.380e-14 1.408e-14
## Max. 1.284e-13 9.382e-14
c2 <- adjustcolor(1:2, 0.6)
matplotB(n, abs19(relE), type="o", cex=3/4, log="y", ylim = c(8e-17, 1.3e-13), yaxt="n", col=c2)
eaxis(2)
abline(h = 2^-53, lty=1, col="gray")
smrelE <- apply(abs(relE), 2, \(y) lowess(n, y, f = 0.1)$y)
matlines(n, smrelE, lwd=3, lty=1)
legend("topleft", legend = expression(stirlerr_simpl(n), stirlerr_simple(n, "MM2")),
bty='n', col=c2, lwd=3, lty=1)
drawEps.h(-(53:48))
## should we e.g., use interpolation spline through sfserr_halves[] for n <= 7.5
## -- doing the interpolation on the log(1 - 12*x*stirlerr(x)) vs log2(x) scale -- maybe ?
stirL <- curve(1-12*x*stirlerr(x, cutoffs = cuts, verbose=TRUE), 1/64, 8, log="xy", n=2048)
## just need "true" values for x = 2^-(6,5,4,3,2) in addition to those we already have at x = 1/2, 1.5, 2, 2.5, ..., 7.5, 8
xM <- mpfr(stirL$x, 2048)
stilEM <- roundMpfr(1 - 12*xM*stirlerr(xM, verbose=TRUE), 128)
relEsml <- relErrV(stilEM, stirL$y)
plot(stirL$x, relEsml) # again: stirlerr() is "limited.." for ~ [2, 6]
## The function to interpolate:
plot(asNumeric( (stilEM)) ~ x, data = stirL, type = "l", log="x")
plot(asNumeric(sqrt(stilEM)) ~ x, data = stirL, type = "l", log="x")
plot(asNumeric( log(stilEM)) ~ x, data = stirL, type = "l", log="x")
y <- asNumeric( log(stilEM))
x <- stirL$x
spl <- splinefun(log(x), y)
plot(asNumeric(log(stilEM)) ~ x, data = stirL, type = "l", log="x")
lines(x, spl(log(x)), col=2)
summary(rE <- relErrV(target = y, current = spl(log(x)))) # all 0 {of course: interpolation at *these* x}
ssmpl <- smooth.spline(log(x), y)
str(ssmpl$fit)# only 174 knots, 170 coefficients
methods(class="smooth.spline.fit") #-> ..
str(yp <- predict(ssmpl$fit, x=log(x)))
lines(x, yp$y, col=adjustcolor(3, 1/3), lwd=3) # looks fine
## but
summary(reSpl <- relErrV(target = y, current = yp$y))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.35e-10 -7.55e-11 1.10e-12 0.00e+00 7.32e-11 5.53e-10
## which is *NOT* good enough, of course ....
showProc.time()
Try the DPQ package in your browser
Any scripts or data that you put into this service are public.
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.