require(boot)
mean.boot <- function(x, ind) {mean(x[ind])}
sd.boot <- function(x, ind) {sd(x[ind])}
var.boot <- function(x, ind) {var(x[ind])}

x <- rexp(10)
boot(x, mean.boot, 1000)
boot(x, sd.boot, 1000)

x <- rexp(25)
out <- boot(x, var.boot, 1000)
hist(out$t, breaks=30, col="gray60", border="white", main="n=25", xlab=expression(theta[b]^'*'))

## A function to construct a CI for the variance using functional delta method
nciv <- function(x,alpha=.05) {
  n <- length(x)
  xbar <- mean(x)
  sig <- mean((x-xbar)^2)
  L <- (x-xbar)^2-sig
  tau2 <- mean(L^2)
  SE <- sqrt(tau2/n)
  return(sig + c(-1,1)*qnorm(1-alpha/2)*SE)
}
sqrt(nciv(out$t))