log1pmx-etc.R
In DPQ: Density, Probability, Quantile ('DPQ') Computations

### R code from vignette source 'log1pmx-etc.Rnw'

###################################################
### code chunk number 1: require
###################################################
require(DPQ)


###################################################
### code chunk number 2: plot_D0tilde
###################################################
par(mar = 0.1 + c(2.5, 3, 0, 0), mgp = c(1.5, 0.6, 0), las=1)
curve(x*log(x)+1-x, 1e-7, 6, n=1001, col=2, lwd=2,
      panel.first=grid(), xlab=quote(u), ylab="")
mtext(quote(tilde(D[0])(u) == ~ u %.%~ log(u)+1-~u), line=-2, col=2)
abline(a = 1-exp(1), b=1, col=adjustcolor(4, 3/4), lwd=1.5, lty=2)
text(5, 2.5, quote(1%.%u - e+1), col=4)
axis(1, at=exp(1), quote(e), tck=0.2, col=4, col.axis=4, lty=3)


###################################################
### code chunk number 3: p.l1p1-def
###################################################
p.p1l1 <- function(from, to, ylim=NULL, cS = adjustcolor(6, 1/2),
                   n=1024, do.leg=TRUE) {
    stopifnot(is.numeric(from), is.numeric(to),
              is.character(cS), length(cS) == 1)
    cols <- palette()[c(2,4, 6, 3,5)]; cols[3] <- cS
    c1 <- curve(x*log1p(x), from=from, to=to, n=n,
                col=2, ylab="", ylim=ylim,
                panel.first = abline(h=0, v=0, lty=3, lwd=1/2))
    c2 <- curve(log1pmx(x), add=TRUE, n=n, col=4)
    with(c1, {
        lines(x, y+c2$y, col=cS, lwd=3)
        lines(x, x^2/2          , col=3, lty=2)
        lines(x, x^2/2*(1 - x/3), col=5, lty=4)
    })
    if(do.leg) {
        labexpr <- expression(
            x %.% log1p(x), log1pmx(x),
            p1l1(x) == log1pmx(x) + x %.% log1p(x),
            x^2 / 2, # frac(x^2, 2)  # is too large
            x^2 / 2 %.% (1 - x/3))
        legend("top", labexpr, col = cols,
               lwd=c(1,1,3,1,1), lty=c(1,1,1:2,4), bty="n")
    }
}


###################################################
### code chunk number 4: l1p1-curves
###################################################
par(mfcol=1:2, mar = 0.1 + c(2.5, 3, 1, 2), mgp = c(1.5, 0.6, 0), las=1)
p.p1l1(-7/8, 2, ylim = c(-1,2))
zoomTo <- function(x,y=x, tx,ty){ arrows(x,-y, tx, ty)
                                    text(x,-y, "zoom in", adj=c(1/3,9/8)) }
zoomTo0 <- function(x,y=x) zoomTo(x,y, 0,0)
zoomTo0(.3)
p.p1l1(-1e-4, 1.5e-4, ylim=1e-8*c(-.6, 1), do.leg=FALSE)


###################################################
### code chunk number 5: log1pmx-curves
###################################################
lcurve <- function(Fn, a,b, ylab = "", lwd = 1.5, ...)
    plot(Fn, a,b, n=1001, col=2, ylab=ylab, lwd=lwd, ...,
         panel.last = abline(h=0, v=-1:0, lty=3))

par(mfrow=c(2,2), mar = 0.1 + c(2.5, 3, 1, 2), mgp = c(1.5, 0.6, 0), las=1)
lcurve(log1pmx, -.9999, 7, main=quote(log1pmx(x) == log(1+x)-x))
                            rect(-.1,  log1pmx(-.1  ), .1  , 0); zoomTo0(1/2, 1)
lcurve(log1pmx, -.1,  .1 ); rect(-.01, log1pmx(-.01 ), .01 , 0); zoomTo0(.02, .001)
lcurve(log1pmx, -.01, .01); rect(-.002,log1pmx(-.002), .002, 0); zoomTo0(2e-3,1e-5)
lcurve(function(x) -log1pmx(x), -.002, .002, log="y", yaxt="n") -> l1r
sfsmisc::eaxis(2); abline(v=0, lty=3)
d1r <- cbind(as.data.frame(l1r), y.naive = with(l1r, -(log(1+x)-x)))
c4 <- adjustcolor(4, 1/3)
lines(y.naive ~ x, data=d1r, col=c4, lwd=3, lty=2)
legend("left", legend=expression(- log1pmx(x), -(log(1+x)-x)),
       col=c(palette()[2],c4), lwd=c(1,3), lty=1:2, bty="n")


###################################################
### code chunk number 6: log1pmx-naive-error
###################################################
par(mfrow=1:2, mar = 0.1 + c(2.5, 3, 1, 2), mgp = c(1.5, 0.6, 0), las=1)
d1r[, "relE.naive"] <- with(d1r, sfsmisc::relErrV(y, y.naive))
plot(relE.naive ~ x, data=d1r, type="l", ylim = c(-1,1)*1e-6)
y2 <- 1e-8
rect(-.002, -y2, .002, y2, col=adjustcolor("gray",1/2), border="transparent")
zoomTo(15e-4, 9*y2, 13e-4, -y2)
plot(relE.naive ~ x, data=d1r, type="l", ylim = c(-1,1)*y2); abline(h=0,lty=3)