How can I create a $1-\alpha$ nonparametric confidence interval for the median using order statistics?

Question

I need to create a 95% confidence interval for the median based on 9 ordered statistics. I know how to determine the confidence level for a given interval but I admit I'm stuck when it comes to creating the interval for a given level. Here is my work so far:

$0.95=P(Y_i<m<Y_j)$

$\implies0.05=P(i-0.5<W<j+0.5)$ adjusting for correction, with $W\sim\mathrm{bin}(9,0.5)$ with $\mu=4.5$ and $\sigma^2=2.25$

By the central limit theorem, we have that $\frac{W-\mu}{\sigma}\approx Z$, where $Z\sim N(0,1)$.

Therefore, $0.95=P(\frac{i-5}{1.5}<Z<\frac{j-4}{1.5})$

There must clearly be a more efficient way of finding $i$ and $j$ then by looking in my tables for the normal distribution in order to find a match. Am I doing something wrong? Also, is there a way to solve a problem like this in R?

Answer 1

Exact binomial. If $W \sim Binom(9, .05)$ it is easy to make a CDF table of the distribution of $w$ using R.

 > n = 9;  p = .5;  x = 0:9;  pdf=dbinom(x, n, p)
 > cbind(x, pdf)
      x         pdf
 [1,] 0 0.001953125
 [2,] 1 0.017578125
 [3,] 2 0.070312500
 [4,] 3 0.164062500
 [5,] 4 0.246093750
 [6,] 5 0.246093750
 [7,] 6 0.164062500
 [8,] 7 0.070312500
 [9,] 8 0.017578125
[10,] 9 0.001953125

Because binomial is discrete, you probably can't get $exactly$ probability $.95,$ but I suppose you want to get an interval as close to $.95$ as possible, without going under.

Looking at the above it seems as if including probabilities from 2 through 7 will go a little over .95; verified below.

$P(2 \le W \le 7) = P(W \le 7) - P(W \le 1) = 0.9609.$

 > diff(pbinom(c(1,7), n, p))
 [1]  0.9609375
 > sum(dbinom(2:7, n, p))
 [1] 0.9609375

Normal approximation. Using the normal approximation $Norm(\mu = 4.6, \sigma = 1.5)$, you could get quantiles .025 and .975, which are 1.56 and 7.43 respectively. So again you can include integers from 2 through 7 to get as close to .95 without going under.

 qnorm(c(.025, .975), n*p, sqrt(n*p*(1-p)))
 [1] 1.560054 7.439946

[Of course, you could revert to the 20th century, standardize, and use printed standard normal tables to get the same result: $(j-4.5)/1.5 = 1.96$ implies $j \approx 7.44.$]

Either way, in your formulation with $<$ instead of $\leq$ and continuity correction, I guess you would say $i = 1$ and $j = 8.$

plot(x, pdf, type="h", ylim=c(0,.3), main="BINOM(9, 1/2) and NORM(4.5, 1.5)")
abline(h=-.002, col="green2")
lines(2:7, pdf[3:8], type="h", lwd=3, col="blue")
curve(dnorm(x, 4.5, 1.5), col="maroon", add=T)