Mdthbs

Section 1.7.2 of Discovering Statistics Using R by Andy Fields, et all, while listing a 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 median is relatively unaffected by extreme scores, I'd have thought it to be more stable across samples. So I was puzzled by authors's 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.

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.

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 jump: 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 jumps instantly to that point and again doesn’t move as you continue moving your point downwards. This is very digital-like and obviously the median is not as smooth as the mean. In fact it has a discontinuity where it instantaneously changes value. Unstable.

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.

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.