#################
## Rat example ##
#################
x <- c(52, 10, 40, 104, 51, 27, 146, 30, 46)
SE <- sqrt(var(x)/length(x))

## Bootstrap (Slide 7)
require(boot)
mean.boot <- function(x, ind) {
  c(mean(x[ind]), var(x[ind])/length(x))
}
out <- boot(x, mean.boot, 999)
boot.ci(out)

## Slide 6
z <- (out$t[,1]-out$t0[1])/sqrt(out$t[,2])
That <- ecdf(z)
quantile(z, p=c(.025,.975))
mean(x) + c(-1,1)*qnorm(.975)*SE
mean(x) + c(-1,1)*qt(.975,8)*SE
mean(x) - quantile(z, p=c(.975,.025))*SE ## Same as studentized interval

## Slide 12
n <- 10
y <- rnorm(n)
exp(mean(y)+c(-1,1)*qnorm(.975)*sqrt(var(y)/n))
exp.mean.boot <- function(x, ind) {
  exp(mean(x[ind]))
}
Out <- boot(y, exp.mean.boot, 999)
boot.ci(Out, type=c("norm", "perc"))

## Slide 13
quantile(out$t[,1], p=c(.025, .975)) ## Same as percentile interval

## Slide 17
z0 <- qnorm(sum(out$t[,1] < out$t0[1])/out$R)
L <- empinf(out)
a <- sum(L^3)/(6 * sum(L^2)^1.5)
za <- qnorm(c(.025, .975))
zz <- pnorm(z0 + (z0 + za)/(1 - a * (z0 + za)))
quantile(out$t[,1], p=zz)

## SE of quantiles
x <- seq(-3, 3, len=99)
f <- function(x){
  p <- pnorm(x)
  sqrt(p*((1-p)/(dnorm(x)))^2 + (1-p)*(p/(dnorm(x)))^2)
}
plot(x, f(x), type="l", ylab="SE")