set.seed(454323)

nboot<-10000
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


nvec<-c(50,100,1000,10000)
nreps<-5

pdf('GammaBootTest.pdf',paper='USr',width=11,height=9)
par(mfrow=c(1,2),oma=c(0,0,1,0))
for(i in 1:length(nvec)){
	n<-nvec[i]
	print(i)
	for(irep in 1:nreps){
		X<-rgamma(n,al,1)

		bmean<-boot(X,meancalc,R=nboot)
		tvec<-bmean$t[,1]

		mc.vec<-replicate(nboot,mean(rgamma(n,al,1)))

		minval<-min(c(tvec,mc.vec))
		maxval<-max(c(tvec,mc.vec))
		lo<-minval-0.05*minval/sqrt(n)
		hi<-maxval+0.05*maxval/sqrt(n)
		del<-(hi-lo)/50

		hist(tvec,breaks=seq(lo,hi,length=51),col='gray',main='Bootstrap',ylim=range(0,nboot/10));box()
		xv<-seq(lo,hi,by=0.001)
		yv<-dgamma(xv,n*al,n)
		lines(xv,yv*del*nboot,col='red')
		hist(mc.vec,breaks=seq(lo,hi,length=51),col='gray',main='Truth',ylim=range(0,nboot/10));box()
		lines(xv,yv*del*nboot,col='red')
		title(substitute(n==nval, list(nval=n)), outer=TRUE)

	}
}
dev.off()


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

ci.mat<-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[irep,1,]<-bci$normal[2:3]
	ci.mat[irep,2,]<-bci$basic[4:5]
	ci.mat[irep,3,]<-bci$student[4:5]
	ci.mat[irep,4,]<-bci$percent[4:5]
	ci.mat[irep,5,]<-bci$bca[4:5]
	
}

#quit(save = "yes")

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

apply(cover.mat,2,mean)