Search code examples
rmedianmontecarloconfidence-intervalchi-squared

R code check for bootsrapping ChiSq Median Estimation


I have a homework problem to write a function that bootstraps its way to finding the 95% CI of a df = 2 Chisq distribution median. My function seems to work, but from wiki I get the formula for the median as k(1 - 2/9k)^3 which yields .343 for k = 2. My function's CI is estimated at (1.28, 1.51) with a plenty large distribution, sample size, and number of simulations (100000). So the theoretical answer isn't anywhere in that interval. Can someone tell me where my code is failing please?

ChisqMedian_CI <- function(chiN, nsim, sampleN){
y <- rchisq(n=chiN, df=2)
medy_resample <- NULL
for (i in 1:nsim) {
  y_resample <- sample(y, replace=TRUE, size=sampleN)
  medy_resample[i] <- median(y_resample)
}
  LB <- quantile(medy_resample, probs=c(0.025))
  UB <- quantile(medy_resample, probs=c(0.975))
  return(c(LB, UB))
}

Solution

  • There are two issues here:

    1. The median of a chi-squared distributed random variable is approximately k * (1 - 2 / (9k))^3. For k = 2, this is about 1.4. This is in line with your results.

    2. For (standard) bootstrap, the bootstrap samples are always as large as the original sample. So chiN is not necessary and would be set to sampleN.