Variance of sample median

Question

Given a sample of n units out of a population of N, population median can be estimated by the sample median.

How can we get the variance of this estimator?

Answer 1

The variance of the sample median depends on the distribution you are sampling from. If you know the sampling distribution you can use the distribution of order statistics to find the distribution of the median and thence its variance.

If you don't know and don't want to make assumptions about the distribution, then you can do something like bootstrapping to estimate the variance.

Answer 2

To add to Jonathan's answer and your follow-up question.

Everything depends on the actual population distribution. You know it, but what is it? In the example below (using R) I know the population distribution but it is somewhat unusual - the product of a gamma distribution multiplied by a normal distribution. The sampling distribution of the median could be calculated but is unlikely to be worth the effort.

Bootstrapping is a good practical alternative. My object results contains 1000 sample medians for samples of size 10 drawn from that population and is a nice way to illustrate its sampling distribution.

> population <- rgamma(1000,1,1)*rnorm(1000,2,1)
> plot(density(population), main="Distribution of population")

> n <-10
> samp <- sample(population, n)
> c(median(population), median(samp))
[1] 1.136544 1.738544
> 
> reps <- 1000
> results <- numeric(reps)
> for (i in 1:1000){
+ samp <- sample(population, n)
+ results[i] <- median(samp)
+ }
> 
> plot(density(results), main="Distribution of sample median")

> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1134  0.8320  1.1670  1.2530  1.6090  5.8890 
> quantile(results, c(0.025,0.975))
    2.5%    97.5% 
0.373812 2.573625  

Answer 3

Depending on how serious you are about the finite population nature of your problem, you might want to look at papers like Woodruff (1952) and Francisco and Fuller (1991). Technically, the variance of the median in the case of an SRS/i.i.d. sample depends on the density at the population median; however, the concept of density is not defined for finite populations, so Woodruff was approximating the median with a primitive kernel-type estimator.

Answer 4

My opinion is that an efficient and simple solution in practice is perhaps possible for small sample sizes. First to quote Wikipedia on the topic of Median:

"For univariate distributions that are symmetric about one median, the Hodges–Lehmann estimator is a robust and highly efficient estimator of the population median.[21]"

The HL median estimate is especially simple for small samples of size n, just compute all possible two point (including repeats) averages. From these n(n+1)/2 new constructs, compute the HL Median Estimator as the usual sample median.

Now, per the same Wikipedia article on the median, the cited variance of the median 1/(4*n*f(median)*f(median)). However, for a discrete sample of size n, I would argue that a conservative estimate to assume for the value of the density function at the median point is 1/n, as we are dividing by this term. As a consequence, the variance of the median is expected to be n/4 or lower. For large n, this would be poor, so yes a more complex (and some would suggest subjective) exercise involving re-sampling could be employed to construct bins of the optimal width so as provide a greater probability mass for f(median).

Now if the purpose of the variance estimate is to gain a precision estimate on the median, may I suggest employing the following bounds due to Mallow assuming the Median is greater than Mean, namely: Median - Mean is less than or equal to Sigma (or, -Sigma when the Median is less than the Mean). Equivalently, the Median lies between the Mean plus Sigma and the Mean minus Sigma.

So, inserting population estimators for the mean and sigma, possibly robust, one can establish a bound for the median that would be consisent with the provided mean and sigma estimates based on the sample population.

Answer 5

Building on the aswer by StasK and assuming that we don't know the distribution:

Sort the values ($x_0,..,x_i,..,x_{N-1}$). Then Woodruff (1952) suggests to take the difference between the values at half the square root of the sample size above and below the sample median $\tilde{\mu}_s$: \begin{align} \sigma_{\tilde{\mu}_s} \approx \dfrac{x_{\lceil \frac{1}{2}\left(N+\sqrt{N}\right) \rceil} - x_{\lfloor \frac{1}{2}\left(N-\sqrt{N}\right) \rfloor}}{2} \end{align} This can intuitively be understood, because the median value deviates from the middle position in a sorted list of random samples by $\frac{\sqrt{N}}{2}$ on average. Expanding this idea, you can also calculate: \begin{align} \sigma^2_{\tilde{\mu}_s} \approx \sum_{i=0}^{N-1} {N-1 \choose i} \left(\frac{1}{2}\right)^{1-N} \left(x_i - \tilde{\mu}_s\right)^2  \end{align} Analogously to the variance of the sample mean, and weighting the quadratic difference of each sample value by its probability of being the closest value to the true median.

The first estimator is easier, faster and more robust, the second may be a little more accurate.