# This is a figure of the theoretical minimum uncertainty regions
# sensu Ellner and Holmes figure 1
# Code written by Steven Ellner and Eli Holmes
# The figure function
CSEGtmufigure <- function(N = 20, u = -0.1, s2p = 0.01, make.legend = TRUE) {
mu <- as.vector(u) # function requires scalar not 1x1 matrix
sigma2.b <- as.vector(s2p)
if (s2p == 0) stop("Stopped in CSEGtmufigure(): function does not work with s2p=0.\n", call. = FALSE)
# Set up some figure parameters
ngrid <- 100
Tvals <- seq(1, 110, length = ngrid)
# win.graph(6,6); par(mfrow=c(1,1),cex.axis=1.35,cex.lab=1.35,yaxs="i",xaxs="i");
par(yaxs = "i", xaxs = "i")
# Define a number of functions
Pext <- function(U, V) {
pnorm(U - V) + exp(2 * U * V + pnorm(-(U + V), log.p = T))
}
TMU <- function(T, RelNe, side = "one") {
qt <- switch(EXPR = side, one = qnorm(0.95), two = qnorm(0.975))
a <- log(1 / RelNe)
U <- -mu * sqrt(T / sigma2.b)
V <- a / (sqrt(sigma2.b * T))
upci <- Pext(U + qt * sqrt(T / N), V)
lowci <- Pext(U - qt * sqrt(T / N), V)
return(list(upper = upci, lower = lowci))
}
TMUpper <- function(T, RelNe, side = "one") TMU(T, RelNe, side)$upper
TMLower <- function(T, RelNe, side = "one") TMU(T, RelNe, side)$lower
CIWidth <- function(T, RelNe) {
qt <- qnorm(0.975)
a <- log(1 / RelNe)
U <- -mu * sqrt(T / sigma2.b)
V <- a / (sqrt(sigma2.b * T))
upci <- Pext(U + qt * sqrt(T / N), V)
lowci <- Pext(U - qt * sqrt(T / N), V)
return(upci - lowci)
}
xlabs <- expression(paste("Projection interval ", italic(T), " time steps"))
ylabs <- "xe = log10(N0/Ne)"
safe.limits <- numeric(ngrid)
dead.limits <- numeric(ngrid)
minval <- 1e-16
maxval <- 1 - minval
for (j in 1:ngrid) {
T <- Tvals[j]
safe.limits[j] <- uniroot(function(x) {
TMUpper(T, x) - 0.05
}, lower = minval, upper = maxval)$root
dead.limits[j] <- uniroot(function(x) {
TMLower(T, x) - 0.95
}, lower = minval, upper = maxval)$root
}
graphics::matplot(Tvals, log10(cbind(safe.limits, dead.limits)),
ylim = c(-2.2, 0), type = "l", lty = 1, col = "white",
xlim = c(1, 100), xlab = xlabs, ylab = ylabs
)
polygon(c(Tvals, rev(Tvals)), log10(c(safe.limits, rev(dead.limits))), col = "grey85")
polygon(c(Tvals, rev(Tvals)), log10(c(dead.limits, rep(1, ngrid))), col = "black")
polygon(c(Tvals, rev(Tvals)), c(log10(safe.limits), rep(-3, ngrid)), col = "white")
logRelNe <- seq(-3.125, -0.001, length = ngrid)
RelNe <- 10^(logRelNe)
CIW <- outer(Tvals, RelNe, CIWidth)
z <- contourLines(Tvals, log10(RelNe), CIW, levels = 0.8)
if (length(z) != 0) {
xvals <- z[[1]]$x
yvals <- z[[1]]$y
polygon(c(xvals, 100), c(yvals, yvals[1]), col = "grey45")
}
abline(h = -2.2)
abline(h = 0)
abline(v = 100)
abline(v = 1)
offset <- -0.05
abline(h = -.3)
text(5, -.3 + offset, "50%")
abline(h = -1)
text(5, -1 + offset, "90%")
abline(h = -2)
text(5, -2 + offset, "99%")
# title(expression(paste("nyrs = ", N, ", hat(mu), = ", mu, sigma[b]," = ",sigma2.b)));
title(paste("time steps = ", N, "\n mu = ", format(mu, digits = 2), "s2.p = ", format(sigma2.b, digits = 2)))
if (make.legend) legend("topright", inset = 0.02, bg = "white", c("high certainty P<0.05", "high certainty P>0.95", "uncertain", "highly uncertain"), fill = c("white", "black", "grey85", "grey45"))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.