BCamedci<-function(x,conflevel=0.95,nboots=120){
#This function will calculate the "BCa" #
#  confidence interval for the population median. #
#nboots = B = number of bootstrap samples to take #
#tstar will contain the median calculated for each #
#  bootstrap sample            BH 02/29/04 version #
##
#This method begins like the others we have seen#
##
alpha<-1-conflevel
tstar<-rep(0,nboots)
for (b in 1:nboots){
  tstar[b]<-median(sample(x,replace=T))
}
tstar<-sort(tstar)
##
#We must then calculate z0hat and ahat.  
#The thhat hat thhatdot are as in the formulas on page # 
#  186 of Eforn and Tishirani.  The median on the z0hat #
#  line and thhat line are there because we are using #
#  making the CI about the median.  The mean on the thhatdot #
#  line is part of the definition of thhatdot for any statstic. # 
##
z0hat<-qnorm(sum(tstar<median(x))/nboots)
thhat<-rep(0,length(x))
for (i in 1:length(x)){
  thhat[i]<-median(x[-i])
}
thhatdot<-mean(thhat)
ahat<-sum((thhatdot-thhat)^3)/(6*(sum((thhatdot-thhat)^2))^3/2)
##
#Once we have those quantities we can calculate a1 and a2 #
#  and finish the method analogously to the others that used #           
#  percentiles, with two additions to avoid errors. #
#Take out the commonly occurring cases where ahat could cause #
#  a2 < a1 by making 1/(z0hat+(za/2))<ahat<1/(z0hat+z(1-a/2)). #
##
ahat<-max(ahat,.00001+1/(z0hat+qnorm(alpha/2)))
ahat<-min(ahat,-.00001+1/(z0hat+qnorm(1-alpha/2)))
a1<-pnorm(z0hat+(z0hat+qnorm(alpha/2))/(1-ahat*(z0hat+qnorm(alpha/2))))
a2<-pnorm(z0hat+(z0hat+qnorm(1-alpha/2))/(1-ahat*(z0hat+qnorm(1-alpha/2))
))
##
#Use the 100pct/(n+1) definition of percentiles and make sure #
#  we don't need the 0th or (n+1)st observation.  #
##
lowtstar<-0.5*(tstar[max(1,floor((nboots+1)*a1))]+
  tstar[max(1,ceiling((nboots+1)*a1))])
hitstar<-0.5*(tstar[min(nboots,floor((nboots+1)*a2))]+
  tstar[min(nboots,ceiling((nboots+1)*a2))])
c(lowtstar,hitstar)
}




#This sets up a simulation to compare the normal approximation #
#  basic bootstrap, and BCa bootstrap confidence intervals for #
#  the median at a 95% level.  B=500 and 100 simulation runs are#
#  used.#

success<-matrix(0,100,3)
ci<-array(0,dim=c(100,2,3))
for (i in 1:100){
 x<-rnorm(20,3,4)
 ci[i,,1]<-bootmedci(x,nboots=500)
 ci[i,,2]<-bbmedci(x,nboots=500)
 ci[i,,3]<-BCamedci(x,nboots=500)
 success[i,1]<-((ci[i,1,1]<3)&(ci[i,2,1]>3))
 success[i,2]<-((ci[i,1,2]<3)&(ci[i,2,2]>3))
 success[i,3]<-((ci[i,1,3]<3)&(ci[i,2,3]>3))
}
apply(success,2,mean)