###
### BOOTSTRAP SE'S FOR SAMPLE MEDIAN AS ESTIMATE OF POPULATION MEDIAN
###

set.seed(13)
n <- 100
y <- rexp(n)  ### NOTE: population median is  log(2)=0.69

est <- median(y)
print("sample median")
print(est)
# [1] 0.5797601

rep <- 500; est.rep <- numeric(rep)

for (i in 1:rep) {
  ### sample with replacement
  ind <- sample(1:n, replace=T)
  est.rep[i] <- median(y[ind]) }

print("Bootstrap (and other) SEs")
print(sqrt(var(est.rep)))
# [1] 0.1108093

### what about using asymptotic theory plus density estimate?
dens.est <-  density(y, n=1, from=est)$y
print(1/(sqrt(2)*sqrt(n)*dens.est))
[1] 0.1225824

### or if we happened to "know" the sample is from an exponential dist
print(1/( sqrt(2) * sqrt(n) * (1/mean(y))*exp(-(1/mean(y))*est) ))
# [1] 0.1212503


### note that the bootstrap distribution is centred near the original
### estimate

print(mean(est.rep))
#[1] 0.6117738

print("Various Boostrap 95% CIs")

print(c(est-1.96*sqrt(var(est.rep)), est+1.96*sqrt(var(est.rep))))
# [1] 0.3625738 0.7969464

print(quantile(est.rep, c(.025,.975)))
# [1] 0.4562955 0.8905167 


print(2*est - quantile(est.rep, c(.975,.025)))
# [1] 0.2690036 0.7032247