#
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# fix the random.seed, so I'll always get the same answer:
set.seed(10131985)

y = c(4.313, 4.513, 5.489, 4.265, 3.641, 5.106, 8.006, 5.087)

# Part (a):
#
# Generate bootstrapped sample means \bar{y}: 
nBoot = 1000
B=array(0,dim=c(nBoot,1))
for(i in 1:nBoot){
  ystar = sample(y,size=8,replace=T)
  B[i,] = mean(ystar)
}

# Plot all of these means:
# 
postscript("../../WriteUp/Graphics/Chapter1/ex_1_7.eps", onefile=FALSE, horizontal=FALSE)
hist(B,freq=T)
dev.off()

# The 95% bootstrapped estimate of \bar{y} is then given by
#
print( sprintf("the bootstrapped estimate of q_{0.95} is %10.6f",quantile(B, probs = 0.95)) )


# Part (b):
#
nBoot1=1000
nBoot2=1000
B1=array(0,dim=c(nBoot1,1))
B2=array(0,dim=c(nBoot2,1)) # initialize once and then overwrite elements as needed
for( ii in 1:nBoot1 ){
  # generate a single bootstrap sample from the original data set y:
  ystar = sample(y,size=8,replace=T)
  # using *this* bootstrap sample we compute an estimate of q_{0.95}(\bar{y}) via an inner bootstrap
  for( jj in 1:nBoot2 ){
    B2[jj] = mean(sample(ystar,size=8,replace=T))
  }
  B1[ii] = quantile(B2,probs=0.95)
}

# Lets display the histogram of the bootstrapped q_{0.95}:
#
postscript("../../WriteUp/Graphics/Chapter1/ex_1_7_boot_q95.eps", onefile=FALSE, horizontal=FALSE)
hist(B1,freq=T)
dev.off()

print("the bootstrapped estimated confidence interval for q_{0.95} is")
print( sprintf("(%10.6f,%10.6f)", quantile(B1, probs = 0.025), quantile(B1, probs = 1-0.025)) )