# NOT RUN {
a <- 2^seq(-58,10, length = 256)
fExpr <- expression(
log(1 - exp(-a)),
log(-expm1(-a)),
log1p(-exp(-a)),
log1mexp(a))
names(fExpr) <- c("DEF", "expm1", "log1p", "F")
str(fa <- do.call(cbind, as.list(fExpr)))
head(fa)# expm1() works here
tail(fa)# log1p() works here
## graphically:
lwd <- 1.5*(5:2); col <- adjustcolor(1:4, 0.4)
op <- par(mfcol=c(1,2), mgp = c(1.25, .6, 0), mar = .1+c(3,2,1,1))
matplot(a, fa, type = "l", log = "x", col=col, lwd=lwd)
legend("topleft", fExpr, col=col, lwd=lwd, lty=1:4, bty="n")
# expm1() & log1mexp() work here
matplot(a, -fa, type = "l", log = "xy", col=col, lwd=lwd)
legend("left", paste("-",fExpr), col=col, lwd=lwd, lty=1:4, bty="n")
# log1p() & log1mexp() work here
par(op)
curve(log1pexp, -10, 10, asp=1)
abline(0,1, h=0,v=0, lty=3, col="gray")
## Cutoff c1 for log1pexp() -- not often "needed":
curve(log1p(exp(x)) - log1pexp(x), 16, 20, n=2049)
## need for *some* cutoff:
x <- seq(700, 720, by=2)
cbind(x, log1p(exp(x)), log1pexp(x))
## Cutoff c2 for log1pexp():
curve((x+exp(-x)) - x, 20, 40, n=1025)
curve((x+exp(-x)) - x, 33.1, 33.5, n=1025)
# }
Run the code above in your browser using DataLab