bootmedci<-function(x,conflevel=0.95,nboots=100){

#This function will calculate a naive confidence interval#
#  for the population median. #
#nboots = B = number of bootstrap samples to take#
#tstar will contain the median calculated for each bootstrap#
#  sample #
alpha<-1-conflevel

tstar<-rep(0,nboots)

for (b in 1:nboots){

   tstar[b]<-median(sample(x,replace=T))

}

mtstar<-mean(tstar)

vtstar<-var(tstar)

biasest<-mtstar-median(x)

lcl<-(median(x)-biasest)-qnorm(1-alpha/2)*sqrt(vtstar)

ucl<-(median(x)-biasest)+qnorm(1-alpha/2)*sqrt(vtstar)

c(lcl,ucl)

}





#This will generate one sample of size 20 that is N(3,4)#
#  and calculate the confidence interval for that sample#

x<-rnorm(20,3,4)

bootmedci(x)





#This will perform a brief simulation study (n=1000) and#
#  give the observed coverage probability (nominal=0.95)#

success<-rep(0,1000)

ci<-matrix(0,1000,2)

for (i in 1:1000){

  x<-rnorm(20,3,4)

  ci[i,]<-bootmedci(x)

  success[i]<-((ci[i,1]<3)&(ci[i,2]>3))

}

sum(success/1000)




#This will draw a graphical display of the effectiveness#
#  of the confidence interval for the 1,000 samples#

low<-min(ci[,1])

hi<-max(ci[,2])

plot(0,0,xlim=c(1,1000),ylim=c(low,hi),type="n")

for (i in 1:1000){

  lines(c(i,i),c(low,ci[i,1]))

  lines(c(i,i),c(ci[i,2],hi))

}

lines(c(0,1001),c(3,3))





#These lines compare the width for the classic large sample#
#  CI for the median in this case to the observed widths#
#  and also gives the summary of the observed centers#

2*1.96*sqrt(1/(4*dnorm(3,mean=3,sd=4)^2))/sqrt(20)

summary(ci[,2]-ci[,1])



summary((ci[,2]+ci[,1])/2)