alpha <- .10 normval <- qnorm(1-alpha/2) numsamp <- 50 numsim <- 25 numboot <- 500 bootmat <- matrix(0,nrow=numsim,ncol=2); normmat <- bootmat y <- 1:numsim; ymat <- rbind(y,y) for (i in 1:numsim) { cat(i," ") samp <- rexp(numsamp) # generate random exponentials sampmean <- mean(samp); sampse <- sqrt(var(samp)/numsamp) normmat[i,] <- c(sampmean - normval*sampse,sampmean + normval*sampse) bootmean <- 0 # initialize this to hold our bootstrap estimates for (j in 1:numboot) { bootsamp <- sample(samp,numsamp,replace=T) bootmean[j] <- mean(bootsamp) } bootmat[i,] <- quantile(bootmean,c(alpha/2,1-(alpha/2))) } matplot(t(bootmat),ymat,pch=" ",yaxt="n",ylab="", xlab=expression(paste(lambda))) # empty plot matlines(t(normmat),ymat,lty=rep(1,numsim),col=1) matlines(t(bootmat),ymat+.3,lty=rep(2,numsim),col=1) abline(v=1) legend(1.4,numsim,legend=c("bootstrap","normal approx"),lty=c(2,1))