The median of many small averages

Chebyshev’s inequality says that a quantity with standard deviation $\sigma$ lands more than $k\sigma$ from its mean with chance at most $1/k^2$, knowing nothing else about it. Applied to the average of $n$ independent samples, whose standard deviation is $\sigma/\sqrt n$, it gives a promise: with confidence $1 - \delta$, the average is within $\sigma/\sqrt{n\delta}$ of the true mean. That bound is the answer to a search — there are distributions that meet it — so if the only thing known is the spread, the average cannot be promised to do better.

The trouble is the $1/\sqrt\delta$. Asking for ten times more confidence costs a factor of about three in accuracy, and asking for confidence of one in a million costs a factor of a thousand. For well-behaved data the average is far more accurate than that, because its errors follow the bell curve and the bell curve’s tails fall like $e^{-x^2/2}$. But the bell-curve promise needs assumptions about the tails, and the Chebyshev promise is all that survives without them.

This essay is about a way of combining the same $n$ samples that keeps the tail assumption out and gets the confidence cost down from $1/\sqrt\delta$ to $\sqrt{\log(1/\delta)}$. It is almost embarrassingly simple: cut the samples into blocks, average each block, and report the median of the block averages. The median of means was used implicitly by several people and made explicit, with its guarantee, by Arkadi Nemirovski and David Yudin in 1983, by Mark Jerrum, Leslie Valiant and Vijay Vazirani in 1986, and by Noga Alon, Yossi Matias and Mario Szegedy in 1996.

One sample, one jump #

The failure it is designed against is easy to draw. Take draws that are ordinary normal numbers, except that each has a small chance of being replaced by a very large number, positive or negative — the kind of distribution that makes Chebyshev’s bound nearly tight, because almost all of its spread is concentrated in rare events.

One draw in this sample is $-21.9$, and it drags the ordinary average of all sixty to $-0.445$ — about three standard errors from the true centre, from a single value. Cut the same sixty draws into six blocks of ten. Five of the block means are ordinary; one, the block holding the jump, is near $-2$. The median of six numbers is the average of the middle two, and the jump’s block is at the extreme, so the median is $-0.080$. The jump has spoiled one block, and the median does not look at one block.

That is the whole idea. A single wild value can move an average arbitrarily far; it can move a median only by one position in the ordering. By averaging first, within blocks, the method keeps the averaging that reduces ordinary noise; by taking the median across blocks, it caps the damage any one block can do.

Why the median of blocks is safe #

The guarantee has two steps, each elementary. Split $n$ samples into $k$ blocks of $n/k$. Each block mean has standard deviation $\sigma\sqrt{k/n}$, so by Chebyshev each one misses the true mean by more than twice that with chance at most $1/4$. The blocks are independent, so the number of blocks that miss is at most a binomial count with success chance $1/4$.

The median misses by more than that margin only if at least half the blocks miss, because otherwise more than half the block means are within the margin, and the median is among them. So the median’s failure chance is at most the chance that a coin landing heads a quarter of the time comes up heads at least half the time in $k$ throws.

That binomial tail is small and falls exponentially: Hoeffding’s inequality bounds it by $e^{-k/8}$, and the figure computes it exactly, from $0.156$ at three blocks to $1.8 \times 10^{-5}$ at fifty-nine, always under the bound. So taking $k = 8\ln(1/\delta)$ blocks makes the median fail with chance at most $\delta$, and when it does not fail it is within $2\sigma\sqrt{k/n} = 2\sigma\sqrt{8\ln(1/\delta)/n}$ of the true mean. The confidence cost is now logarithmic: asking for a million times more confidence multiplies the error by about four, not a thousand.

Put numbers in to see the size of the gain. With $n = 10{,}000$ samples and a confidence of one failure in a million, $\delta = 10^{-6}$, the plain mean’s Chebyshev promise is $\sigma/\sqrt{n\delta} = \sigma/\sqrt{0.01}$, which is $10\sigma$ — a guarantee that says the average of ten thousand samples could be ten standard deviations off, and is useless. The median of means with $k = 8\ln(10^6) \approx 111$ blocks of about ninety samples promises $2\sigma\sqrt{111/10{,}000} \approx 0.21\sigma$. The same data, the same single assumption of finite variance, and a promise nearly fifty times sharper, bought by nothing but the order in which the averaging and the voting are done.

The two steps use two different kinds of concentration, which is what makes the result feel like a trick. Chebyshev, the weakest bound there is, handles each block, and needs only a finite variance. Hoeffding, a strong exponential bound, handles the count of failed blocks — and it applies because the failure indicators are bounded, even though the data are not. The median turns an unbounded problem into a bounded one.

On data that behave, the mean is better #

None of this makes the median of means a better estimator in general, and it is worth seeing first where it loses.

On normal data the mean is the best estimator there is, and the figure shows it winning at every confidence level: at $99.9\%$ it is within $0.518$ where the median of means is within $0.626$. The ratio, about $1.2$, is close to $\sqrt{\pi/2} = 1.25$, which is the price of using a median instead of a mean on normal data — a median throws away the information in how far each value is from the centre, and normal data have no outliers for that to protect against. The median of means pays a fixed tax of about a fifth for insurance.

On heavy tails, a draw #

Heavy tails alone are not enough to make the insurance pay.

Student’s $t$ with three degrees of freedom has tails that fall off only like a power — the chance of exceeding $x$ falls like $x^{-3}$ — which is heavy enough to make its fourth moment infinite. Yet with $120$ samples the two estimators are nearly equal, the mean slightly better at ordinary confidence and the median of means slightly better at $99.9\%$. The heavy tail here is spread evenly across many moderately large values, and the central limit theorem, though slower, is still doing most of the averaging. The median of means’ advantage is a guarantee about the worst case, and Student’s $t$ is not the worst case.

Against rare jumps, a rout #

The worst case is the one that makes Chebyshev tight: almost all of the variance concentrated in rare, large events.

Here the jump happens in only about one sample of $120$ in twenty. Most of the time there is no jump and both estimators are as good as on normal data, which is why they agree at $80\%$ confidence. But one repetition in twenty contains a jump, and in those the mean is off by roughly the jump divided by $120$ — about $0.6$ — while the median of means, whose twelve blocks include at most one spoiled one, barely moves. At $99.9\%$ confidence the mean’s error is $0.783$ and the median of means’ is $0.361$. The mean’s error curve rises towards its Chebyshev guarantee, the thin line, which is what “Chebyshev is tight” looks like in practice; the median of means’ curve rises only slowly.

The point of the method is not that it is better on this distribution. It is that its guarantee does not depend on which distribution it faces. On normal data it costs a fifth; on rare jumps it saves more than half; and without knowing in advance which kind of data will arrive, its worst case is enormously better than the mean’s.

Why the plain mean cannot be rescued #

It might seem that a cleverer analysis of the ordinary mean would find a better guarantee. It cannot, and the reason is the distribution in the last figure. For a distribution with a rare large jump chosen so that the jump occurs with chance about $\delta$ among $n$ samples, the mean is off by the jump divided by $n$ whenever it occurs — and with the jump’s size set to keep the variance at $\sigma^2$, that error is of order $\sigma/\sqrt{n\delta}$. So the mean’s $1/\sqrt\delta$ is not a weakness of Chebyshev’s proof; it is a property of the mean, and it is attained.

The median of means escapes because it is not an average. No estimator that is a fixed weighted average of the samples can do better than the mean does here, since a large enough jump in any one sample moves it proportionally. The estimator has to be able to ignore a sample, and the median is the simplest operation that can. There is a matching lower bound: no estimator of any kind can beat $\sigma\sqrt{\log(1/\delta)/n}$ by more than a constant factor, even on normal data, so the median of means is optimal up to its constant — a constant that Gábor Lugosi and Shahar Mendelson, and later Samuel Hopkins, reduced with more elaborate constructions.

Why not take the median of everything #

If a median ignores wild values, why average within blocks at all — why not simply take the median of all $n$ samples? For symmetric data that is a fine estimator, and for Student’s $t$ it would even beat both of the estimators drawn. But in general it estimates the wrong thing. The median of the samples estimates the median of the distribution, and for a skewed distribution that is not the mean: incomes, file sizes, insurance claims and waiting times all have means well above their medians, and a median of samples converges confidently to the wrong answer.

The blocks are what fix that. Each block mean is an average, so it is centred exactly on the true mean, and by the central limit theorem its distribution becomes more symmetric as the block grows. The median of block means therefore sits near the mean rather than near the median of the raw distribution, and the Chebyshev step in the proof does not care about symmetry at all — it guarantees each block is within its margin with chance three in four, whatever the skew. That is the division of labour: averaging within blocks aims the estimator at the mean, and the median across blocks protects it from blocks that went wrong. The number of blocks sets the balance, and $8\ln(1/\delta)$ is the fewest that make the protection strong enough for the confidence asked.

What the method costs #

It costs the constant. The guarantee $2\sigma\sqrt{8\ln(1/\delta)/n}$ has a constant of about $5.7$ in front of $\sigma\sqrt{\ln(1/\delta)/n}$, where the normal distribution’s own tail would give about $1.4$. On well-behaved data the median of means is noticeably worse than the mean, as the first comparison showed, and a practitioner who knows the data are close to normal should use the mean.

It costs a choice. The number of blocks is set by the confidence wanted, not by the data, and a different confidence means a different estimator. Asking for 99% and 99.99% gives two different numbers from the same samples, which is disconcerting until it is seen as the honest price of a guarantee.

And it needs independence between blocks. The Hoeffding step counts failed blocks as independent coin flips. If the samples arrive in an order that correlates nearby values — a drift, a burst — the blocks must be formed at random, or the argument fails. The concentration of the count of failures is exactly the kind of statement that breaks when inputs can move together.

Where it is used #

The median-of-means construction appears wherever a tail guarantee is needed and the data cannot be trusted to be normal. Alon, Matias and Szegedy used it in 1996 to estimate the second moment of a data stream — how concentrated the frequencies of items passing by are — from a tiny memory: each block is a crude randomised estimate, and the median of many crude estimates is a reliable one. The same device amplifies the confidence of any randomised estimate whose single-run success chance is merely better than one half, by repeating it and taking the median; it is the standard way to turn “right three times in four” into “right with confidence $1 - \delta$” at the cost of $\log(1/\delta)$ repetitions.

The same problem arises inside Monte Carlo computation. An estimate built from random samples has an error that falls like one over the square root of the sample size only if the quantity being averaged has finite variance and reasonable tails, and sampling where the answer lives can produce weights with very heavy tails when the sampling distribution is poorly matched to the target. A median of block averages is one of the standard safeguards: it cannot repair an estimator whose variance is infinite, but it turns a merely finite variance into a reliable confidence statement.

In statistics it is one of a family of robust estimators — the trimmed mean, Catoni’s estimator, the Huber estimator — designed to lose little on normal data and a great deal less than the mean on contaminated data. What distinguishes the median of means is that its guarantee is proved for every distribution with finite variance, with no assumption about contamination, from nothing more than Chebyshev and a coin count.

What the curves cannot show #

They cannot show the worst case. Each comparison is for one distribution, and the guarantee is a statement about all of them at once. The rare-jump distribution shows the mean near its worst; no finite set of figures shows the median of means at its worst, because its worst depends on the confidence asked for and the distribution chosen to attack it.

They cannot show confidences beyond one in a thousand. Sixty thousand repetitions measure the $99.9\%$ error with only sixty repetitions beyond it, and the $99.99\%$ error would need ten times more. The logarithmic growth is the theorem, and the figures show its beginning.

And the one-sample figure is one sample. It was chosen by a stated rule — the first sample in a fixed order containing exactly one jump — to show the mechanism. On a single sample either estimator can be the closer; the claim is about how often each is far off.

Still open: the best constant #

The median of means is optimal up to a constant, and the constant is not settled. For estimating a single mean with confidence $1 - \delta$ from samples of finite variance, Jasper Lee and Paul Valiant showed in 2021 that an estimator exists whose error is within a factor $1 + o(1)$ of what the normal distribution’s own tail would give — as good, in the limit, as if the data were normal — but the constructions are intricate, and for estimating a mean in many dimensions, where the medians of vectors are themselves subtle, the question of how close a computationally simple estimator can come to the optimal constant is active research.

Averages that cannot be hijacked #

Given only a finite variance, the average of $n$ samples can be promised to within $\sigma/\sqrt{n\delta}$ with confidence $1-\delta$, and that is tight — rare large jumps attain it. Cutting the samples into $8\ln(1/\delta)$ blocks and taking the median of the block averages improves the promise to about $\sigma\sqrt{\log(1/\delta)/n}$: each block is safe with chance three in four by Chebyshev, and the chance that half of them fail falls like $e^{-k/8}$ by Hoeffding.

On normal data the method costs about a fifth in accuracy, on Student’s $t$ the two estimators are close, and against rare jumps the median of means is more than twice as accurate at $99.9\%$. Its advantage is not on any one distribution but in a guarantee that holds on all of them.

When the enemy is a rare catastrophe, do not average over everything — average in small groups, and let a vote decide which groups to trust.