Why does mean tend be more stable in different samples than median?

Question

Section 1.7.2 of Discovering Statistics Using R by Andy Field, et al., while listing virtues of mean vs median, states:

... the mean tends to be stable in different samples.

This after explaining median's many virtues, e.g.

... The median is relatively unaffected by extreme scores at either end of the distribution ...

Given that the median is relatively unaffected by extreme scores, I'd have thought it to be more stable across samples. So I was puzzled by the authors' assertion. To confirm I ran a simulation — I generated 1M random numbers and sampled 100 numbers 1000 times and computed mean and median of each sample and then computed the sd of those sample means and medians.

nums = rnorm(n = 10**6, mean = 0, sd = 1)
hist(nums)
length(nums)
means=vector(mode = "numeric")
medians=vector(mode = "numeric")
for (i in 1:10**3) { b = sample(x=nums, 10**2); medians[i]= median(b); means[i]=mean(b) }
sd(means)
>> [1] 0.0984519
sd(medians)
>> [1] 0.1266079
p1 <- hist(means, col=rgb(0, 0, 1, 1/4))
p2 <- hist(medians, col=rgb(1, 0, 0, 1/4), add=T)

As you can see the means are more tightly distributed than medians.

In the attached image the red histogram is for medians — as you can see it is less tall and has fatter tail which also confirms the assertion of the author.

I’m flabbergasted by this, though! How can median which is more stable tends to ultimately vary more across samples? It seems paradoxical! Any insights would be appreciated.

Answer 1

The median is maximally robust to outliers, but highly susceptible to noise. If you introduce a small amount of noise to each point, it will enter the median undampened as long as the noise is small enough to not change the relative order of the points. For the mean it's the other way around. Noise is averaged out, but a single outlier can change the mean arbitrarily.

                                        mean  median   
original   [1.0, 2.0, 3.0, 4.0, 5.0]       3       3
noise      [1.1, 1.9, 3.1, 4.1, 4.9]    3.02     3.1
outlier    [100, 2.0, 3.0, 4.0, 5.0]    22.8       4

Your test mostly measures robustness to noise, but you can easily create a sample where the median performs better. If you want an estimator that is robust to both outliers and noise, just throw away the top and bottom third and average the remainder.

Answer 2

As @whuber and others have said, the statement is not true in general. And if you’re willing to be more intuitive — I can’t keep up with the deep math geeks around here — you might look at other ways mean and median are stable or not. For these examples, assume an odd number of points so I can keep my descriptions consistent and simple.

1. Imagine you have spread of points on a number line. Now imagine you take all of the points above the middle and move them up to 10x their values. The median is unchanged, the mean moved significantly. So the median seems more stable.

2. Now imagine these points are fairly spread out. Move the center point up and down. A one-unit move changes the median by one, but barely moved the mean. The median now seems less stable and more sensitive to small movements of a single point.

3. Now imagine taking the highest point and moving it smoothly from the highest to the lowest point. The mean will also smoothly move. But the median will not move continuously: it won’t move at all until your high point becomes lower than the previous median, then it starts following the point until it goes below the next point, then the median sticks to that point and again doesn’t move as you continue moving your point downwards. [Edited per comment]

So different transformations of your points cause either mean or median to look less smooth or stable in some sense. The math heavy-hitters here have shown you distributions from which you can sample, which more closely matches your experiment, but hopefully this intuition helps as well.

Answer 3

Suppose you have $n$ data points from some underlying continuous distribution with mean $\mu$ and variance $\sigma^2 < \infty$. Let $f$ be the density function for this distribution and let $m$ be its median. To simplify this result further, let $\tilde{f}$ be the corresponding standardised density function, given by $\tilde{f}(z) = \sigma \cdot f(\mu+\sigma z)$ for all $z \in \mathbb{R}$. The asymptotic variance of the sample mean and sample median are given respectively by:

$$\mathbb{V}(\bar{X}_n) = \frac{\sigma^2}{n} \quad \quad \quad \quad \quad \mathbb{V}(\tilde{X}_n) \rightarrow \frac{\sigma^2}{n} \cdot \frac{1}{4} \cdot \tilde{f}\Big( \frac{m-\mu}{\sigma} \Big)^{-2}.$$

We therefore have:

$$\frac{\mathbb{V}(\bar{X}_n)}{\mathbb{V}(\tilde{X}_n)} \rightarrow 4 \cdot \tilde{f}\Big( \frac{m-\mu}{\sigma} \Big)^2.$$

As you can see, the relative size of the variance of the sample mean and sample median is determined (asymptotically) by the standardised density value at the true median. Thus, for large $n$ we have the asymptotic correspondence:

$$\mathbb{V}(\bar{X}_n) < \mathbb{V}(\tilde{X}_n) \quad \quad \iff \quad \quad f_* \equiv \tilde{f} \Big( \frac{m-\mu}{\sigma} \Big) < \frac{1}{2}.$$

That is, for large $n$, and speaking asymptotically, the variance of the sample mean will be lower than the variance of the sample median if and only if the standardised density at the standardised median value is less than one-half. The data you used in your simulation example was generated from a normal distribution, so you have $f_* = 1 / \sqrt{2 \pi} = 0.3989423 < 1/2$. Thus, it is unsurprising that you found a higher variance for the sample median in that example.

Answer 4

Comment: Just to echo back your simulation, using a distribution for which SDs of means and medians have the opposite result:

Specifically, nums are now from a Laplace distribution (also called 'double exponential'), which can be simulated as the difference of two exponential distribution with the same rate (here the default rate 1). [Perhaps see Wikipedia on Laplace distributions.]

set.seed(2019)
nums = rexp(10^6) - rexp(10^6)
means=vector(mode = "numeric")
medians=vector(mode = "numeric")
for (i in 1:10^3) { b = sample(x=nums, 10^2); 
  medians[i]= median(b); means[i]=mean(b) }
sd(means)
[1] 0.1442126
sd(medians)
[1] 0.1095946   # <-- smaller

hist(nums, prob=T, br=70, ylim=c(0,.5),  col="skyblue2")
 curve(.5*exp(-abs(x)), add=T, col="red")

Note: Another easy possibility, explicitly mentioned in @whuber's link, is Cauchy, which can be simulated as Student's t distribution with one degree of freedom, rt(10^6, 1). However, its tails are so heavy that making a nice histogram is problematic.