library(boot)

mouse.c <- scan("http://www.rohan.sdsu.edu/~babailey/stat672/mouse.c.data")
mouse.t <- scan("http://www.rohan.sdsu.edu/~babailey/stat672/mouse.t.data")


y<-mouse.c
sd(y)

n.c<-length(mouse.c)

B<-100000

thetahat.c<-mean(mouse.c)
sehat.c<-sd(mouse.c)*sqrt((n.c-1)/n.c^2)

theta.boot<-rep(0,B)

Zb<-rep(0,B)
Zbt<-rep(0,B)
for(b in 1:B){

	yboot<-sample(mouse.c,size=n.c,rep=T)
	thetahat.boot<-mean(yboot)
	sehat.boot<-sd(yboot)*sqrt((n.c-1)/n.c^2)

	Zb[b]<-(thetahat.boot-thetahat.c)/sehat.c

	Zbt[b]<-(thetahat.boot-thetahat.c)/sehat.boot

	theta.boot[b]<-thetahat.boot

}
hist(theta.boot,xlab=expression(theta),br=seq(0,130,by=5),col='gray',main=expression(paste('Bootstrap resamples of ',hat(theta))));box()
abline(v=thetahat.c,col='red',lwd=2)

abline(v=quantile(theta.boot,prob=c(0.025,0.975)),col='red',lty=2,lwd=2)

alpha<-0.05
qs<-quantile(Zb,prob=c(1-alpha/2,alpha/2))
qst<-quantile(Zbt,prob=c(1-alpha/2,alpha/2))
abline(v=thetahat.c-qst*sehat.c,col='green',lty=2,lwd=2)
abline(v=thetahat.c-qs*sehat.c,col='blue',lty=2,lwd=2)

Bvec<-c(100,200,500,1000,2000,5000,10000,20000,50000,100000)
Qmat<-matrix(0,nrow=length(Bvec),ncol=2)
for(i in 1:length(Bvec)){
	Qmat[i,]<-quantile(Zbt[1:Bvec[i]],c(1-alpha/2,alpha/2))
}

#Studentized bootstrap intervals for different B
thetahat.c-Qmat*sehat.c


#Normal Interval
thetahat.c+c(-1,1)*qnorm(1-alpha/2)*sehat.c

#Student Interval
thetahat.c+c(-1,1)*qt(1-alpha/2,n.c-1)*sehat.c


##################################################################
#Using the boot function

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

mouse.ci<-boot.ci(bmean)

hist(bmean$t[,1],xlab=expression(theta),br=seq(-10,125,by=5),col='gray',main=expression(paste('Bootstrap resamples of ',hat(theta))));box()
abline(v=thetahat.c,col='red',lwd=2)

abline(v=mouse.ci$normal[2:3],col='green')
abline(v=mouse.ci$basic[4:5],col='blue')
abline(v=mouse.ci$student[4:5],col='green',lty=2)
abline(v=mouse.ci$percent[4:5],col='blue',lty=2)
abline(v=mouse.ci$bca[4:5],col='orange')

legend(-10,14000,c('Est.','Normal','Basic','Student','Percent','BCa'),
col=c('red','green','blue','green','blue','orange'),lty=c(1,1,1,2,2,1))

dev.print(pdf,'MouseBootCI.pdf',width=11,height=9)