#### slides for WSUG presentation, also KTN

## This code will write out figures as pdf files, so choose an
## appropriate directory ...

setwd('~/working/Presentations/')

## Use Pareto distribution for illustration (this not available in R)

ppareto <- function(x, mode, shape = 1)
  ifelse(x <= mode, 0, 1 - (x / mode)^(-shape))

qpareto <- function(u, mode, shape = 1)
  ifelse(u <= 0, mode, mode / (1 - u)^(1/shape))

rpareto <- function(n, mode, shape = 1)
  qpareto(runif(n), mode = mode, shape = shape)

mode <- 0.5
shape <- 2
v <- rpareto(1e3, mode = mode, shape = shape); hist(v) # seems to work OK

## first picture, shows PE curve and then the area under it

par(mar = c(3, 2, 3, 1), mgp = c(1.5, 0.5, 0), cex.axis = 0.8,
    cex.lab = 0.8, cex.main = 1, las = 1, tcl = -0.25)

xlim <- c(0, 3)
x <- seq(from = xlim[1], to = xlim[2], length = 201)
y <- 1 - ppareto(x, mode = mode, shape = shape)

plot(x, y,
     xlim = xlim,
     type = 'l', lwd = 2,
     main = 'Prob. of exceedence (PE)',
     xlab = 'Maximum water depth (m)',
     ylab = '', bty = 'n')

dev.print(pdf, file = 'PEcurve1.pdf', width = 3)

polygon(c(x, xlim[2], xlim[1]), c(y, 0, 0), col = 'grey', border = NA)
lines(x, y, lwd = 2)
text(1.25, 0.6, 'Shaded area is\nexpected maximum\nwater depth',
     pos = 4, cex = 0.8)

dev.print(pdf, file = 'PEcurve2.pdf', width = 3)

## second picture shows multiple PE curves and their combination

shapes <- c(shape, 1.5, 2.5, 3, 4)
pis <- c(0.3, 0.2, 0.2, 0.2, 0.1)

par(mar = c(3, 2, 3, 1), mgp = c(1.5, 0.5, 0), cex.axis = 0.8,
    cex.lab = 0.8, cex.main = 0.8, las = 1, tcl = -0.25, xpd = NA)

xlim <- c(0, 5)
x <- seq(from = xlim[1], to = xlim[2], length = 201)
Y <- do.call('cbind', lapply(shapes, function(ss)
                             1 - ppareto(x, mode = mode, shape = ss)))

col <- palette()[1 + seq(along = shapes)]
plot(x, Y[, 1],
     xlim = xlim,
     type = 'l',
     main = 'Prob. of exceedence (PE)',
     xlab = 'Maximum water depth (m)',
     ylab = '', bty = 'n', col = col[1],
     axes = FALSE)
axis(1, pos = 0)
axis(2)

dev.print(pdf, file = 'shapes1.pdf', height = 3.75, width = 5)

for (i in seq(along = shapes)[-1])
  lines(x, Y[, i], col = col[i])

dev.print(pdf, file = 'shapes2.pdf', height = 3.75, width = 5)

legend('topright', inset = c(0.05, 0.05),
       legend = c(formatC(100 * pis, format = 'f', width = 2, digits = 0), 'Combined'),
       lty = 1, pch = NA,
       lwd = rep(c(1, 2), c(length(shapes), 1)),
       col = c(col, 'black'), bty = 'n',
       title = 'Probability (%)',
       cex = 0.75, y.intersp = 3)

y <- drop(Y %*% pis)
lines(x, y, lwd = 2)

text(2.5, 0.6, 'This step\nis quite\ncontentious!', cex = 1)

dev.print(pdf, file = 'shapes3.pdf', height = 3.75, width = 5)

abline(h = 0.1, lty = 2)

for (i in 1:ncol(Y)) {
  foo <- approx(Y[, i], x, xout = 0.1)$y
  lines(c(foo, foo), c(0.1, 0), col = col[i], lty = 2)
  points(foo, 0, col = col[i], pch = 16)
}

dev.print(pdf, file = 'shapes4.pdf', height = 3.75, width = 5)

lines(x, y, lwd = 2)
foo <- approx(y, x, xout = 0.1)$y
lines(c(foo, foo), c(0.1, 0), col = 'black', lty = 2, lwd = 2)
points(foo, 0, col = 'black', pch = 16)

axis(1, at = foo, formatC(foo, format = 'f', digits = 1),
     tick = FALSE, pos = 0)

dev.print(pdf, file = 'shapes5.pdf', height = 3.75, width = 5)

## additional analysis for KTN slides, showing confidence band for
## different numbers of simulations.  The confidence band material can
## be found in L. Wasserman, All of Statistics, Springer.

N <- c(100, 500, 1000, 1500)
xlim <- c(0, 3)
x <- seq(from = xlim[1], to = xlim[2], length.out = 301)
alpha <- 0.05

par(mfrow = c(2, 2), mar = c(3, 2, 2, 1), oma = c(0, 0, 0, 0),
    cex.axis = 0.8, xpd = FALSE)

for (n in N) {

  epsn <- sqrt((1/(2*n)) * log(2 / alpha))
  
  X <- rpareto(n, mode = mode, shape = shape)
  F <- function(xval) sapply(xval, function(xv) mean(X <= xv))
  plot(x, 1 - F(x), type = 'l', lwd = 2,
       xlim = xlim, ylim = c(0, 1),
       xlab = 'Maximum water depth (m)',
       ylab = '', main = sprintf("N = %i", n))
  abline(h = 0.1, lty = 2)
  text(0.1, 0.1, ' 10%', pos = 3, offset = 0.25, cex = 0.8)
#   rug(X)

  if (n == N[1])
    dev.print(pdf, file = sprintf('cband%iinit.pdf', n), height = 4.25,
              width = 5)
    

  Lx <- 1 - pmax(F(x) - epsn, 0)
  Ux <- 1 - pmin(F(x) + epsn, 1)
  lines(x, Lx, lty = 2, lwd = 2) # upper bound
  lines(x, Ux, lty = 2, lwd = 2) # lower bound

  dev.print(pdf, file = sprintf('cband%i.pdf', n), height = 4.25,
            width = 5)
}