set.seed(454323)

library(boot)
meancalc<-function(x,i){
     m <- mean(x[i])
     n <- length(i)
     v <- (n-1)*var(x[i])/n^2
     return(c(m, v))
}

al<-2



n<-20
nboot<-10000
nreps<-1000

ci.mat.20<-array(0,c(nreps,5,2))

for(irep in 1:nreps){

	print(irep)
	X<-rgamma(n,al,1)

	bmean<-boot(X,meancalc,R=nboot)
	tvec<-bmean$t[,1]
	bci<-boot.ci(bmean)
	ci.mat.20[irep,1,]<-bci$normal[2:3]
	ci.mat.20[irep,2,]<-bci$basic[4:5]
	ci.mat.20[irep,3,]<-bci$student[4:5]
	ci.mat.20[irep,4,]<-bci$percent[4:5]
	ci.mat.20[irep,5,]<-bci$bca[4:5]
	
}

#quit(save = "yes")

cover.mat.20<-apply(ci.mat.20,1:2,function(x){return(as.numeric(x[1]<=al & x[2] >= al))})

apply(cover.mat.20,2,mean)