Example: Simulation Study of elm
vs. Bootstrap
## Define elm function and its service functions
f <- scan('http://www-stat.stanford.edu/~owen/empirical/el.S',
what=list(x=''), sep='\n')$x
## Author used _ as assignment operator, not a good idea
## Otherwise could have used source('http:...') directly
f <- gsub('_', '<-', f)
source(textConnection(f))
library(boot)
elcl <- function(x, pl=FALSE) {
guess <- quantile(x, c(.1,.9))
mu <- seq(guess[1],guess[2],length=100)
n <- length(mu)
w <- sapply(mu, function(X) elm(x, X)$logelr)
i <- which(w == max(w))
cls <- c(approx(x=w[1:i], y=mu[1:i], xout=max(w)-3.84/2, rule=2)$y,
approx(x=w[i:n], y=mu[i:n], xout=max(w)-3.84/2, rule=2)$y)
if(pl) {
plot(mu, w)
abline(v=mean(x))
abline(v=cls)
abline(h=max(w)-3.84/2, lty=2)
}
cls
}
m <- 1000
bnames <- c('normal','basic','percent','bca')
cl <- array(NA, dim=c(m,1+length(bnames),2),
dimnames=list(NULL,c('el',bnames),c('lower','upper')))
set.seed(1)
for(i in 1:m) {
cat(i,'')
x <- exp(rnorm(25))
cl[i,'el',] <- elcl(x)
## clb[i,] <- smean.cl.boot(x, na.rm=FALSE)[-1]
b <- boot(x, function(x,i) mean(x[i]), R=999)
cb <- boot.ci(b)
## For each type of bootstrap confidence interval, fetch last
## two values in vector - these are the CLs; put into matrix
cb <- t(sapply(cb[bnames],
function(w) {
l <- length(w)
c(w[l-1],w[l])
}))
cl[i,rownames(cb),] <- cb
}
mu <- exp(.5)
cat('\nMean CL lengths:\n\n')
apply(apply(cl,1:2,function(x)x[2]-x[1]), 2, mean)
cat('\nLeft Tail Areas:\n\n')
apply(cl[,,'lower'] > mu, 2, mean)
cat('\nRight Tail Areas:\n\n')
apply(cl[,,'upper'] < mu, 2, mean)
Mean CL lengths:
> apply(apply(cl,1:2,function(x)x[2]-x[1]), 2, mean)
el normal basic percent bca
1.385881 1.440112 1.416760 1.416760 1.640999
>
> cat('\nLeft Tail Areas:\n\n')
Left Tail Areas:
> apply(cl[,,'lower'] > mu, 2, mean)
el normal basic percent bca
0.024 0.003 0.002 0.008 0.026
>
> cat('\nRight Tail Areas:\n\n')
Right Tail Areas:
> apply(cl[,,'upper'] < mu, 2, mean)
el normal basic percent bca
0.097 0.148 0.172 0.131 0.090