19. Order statistics

Take a sample X_1, \dots, X_n and sort it. The sorted values are the order statistics:

X_{(1)} \le X_{(2)} \le \cdots \le X_{(n)}

X_{(1)} is the minimum, X_{(n)} the maximum, and X_{(\lceil n/2 \rceil)} is the median. These are random variables in their own right, with their own distributions — and they're what you need whenever the question is about extremes rather than averages.

Note immediately: even if the X_i are independent, the order statistics are not. Sorting creates dependence — X_{(1)} \le X_{(2)} always, by construction.

The maximum

The trick is that a maximum is small only if everything is small:

X_{(n)} \le x \iff X_1 \le x \text{ and } X_2 \le x \text{ and } \cdots \text{ and } X_n \le x

For independent, identically distributed X_i, that's a product:

F_{X_{(n)}}(x) = P(X_{(n)} \le x) = \big[F(x)\big]^n

Differentiate for the density:

f_{X_{(n)}}(x) = n\big[F(x)\big]^{n-1}f(x)

The minimum

Symmetrically, a minimum is large only if everything is large — so work with the survival function:

P(X_{(1)} > x) = \big[1 - F(x)\big]^n

F_{X_{(1)}}(x) = 1 - \big[1 - F(x)\big]^n, \qquad f_{X_{(1)}}(x) = n\big[1 - F(x)\big]^{n-1}f(x)

Both formulas come straight from the CDF method of the last lesson. Extremes are one of the places where going through the CDF isn't just safer — it's the only sensible route.

The general k-th order statistic

f_{X_{(k)}}(x) = \frac{n!}{(k-1)!\,(n-k)!}\big[F(x)\big]^{k-1}\big[1 - F(x)\big]^{n-k}f(x)

The structure is a multinomial count. For X_{(k)} to sit at x, you need:

• k-1 observations below x, each with probability F(x)
• exactly one observation at x, contributing the density f(x)
• n-k observations above x, each with probability 1 - F(x)

and \frac{n!}{(k-1)!\,1!\,(n-k)!} ways to choose which observations play which role. Setting k = n or k = 1 recovers the two special cases above.

The uniform case, and why it matters

For X_i \sim \text{Uniform}(0,1) we have F(x) = x on [0,1], so:

f_{X_{(k)}}(x) = \frac{n!}{(k-1)!(n-k)!}x^{k-1}(1-x)^{n-k}

That is exactly a Beta(k,\; n-k+1) density. So

E[X_{(k)}] = \frac{k}{n+1}

The n order statistics of a uniform sample divide [0,1] into n+1 gaps of equal expected width. The minimum of 9 uniforms sits at 0.1 on average, the median at 0.5, the maximum at 0.9.

This isn't only about uniforms. Because F(X) \sim \text{Uniform}(0,1) for any continuous F (the probability integral transform, which we'll use again in §11), the uniform case is the universal one in disguise.

Worked example

n = 5 independent Uniform(0,1) draws. Find the density of the maximum, its mean, and P(X_{(5)} > 0.9).

With F(x) = x:

F_{X_{(5)}}(x) = x^5, \qquad f_{X_{(5)}}(x) = 5x^4 \quad (0 \le x \le 1)

Mean:

E[X_{(5)}] = \int_0^1 x \cdot 5x^4\,dx = \int_0^1 5x^5 dx = \left[\frac{5x^6}{6}\right]_0^1 = \frac{5}{6} \approx 0.833

which matches \frac{k}{n+1} = \frac{5}{6}. ✓

Tail probability:

P(X_{(5)} > 0.9) = 1 - (0.9)^5 = 1 - 0.59049 = 0.41

Read that last number carefully. Any single draw exceeds 0.9 with probability 0.1. But with five draws, there's a 41% chance that at least one does. Maxima drift upward fast, and the more samples you take the further out they go.

This is the mathematics behind a practical trap: if you run 20 experiments and report only the best result, the best of 20 is systematically extreme even when nothing real is happening. It's the same computation that makes multiple-testing corrections necessary in §5.

Doing it in Python

Simulate by sorting, and compare against the formulas:

import numpy as np

rng = np.random.default_rng(0)
n, trials = 5, 200_000

samples = rng.uniform(0, 1, size=(trials, n))
ordered = np.sort(samples, axis=1)          # each row sorted ascending

print(f"{'k':>3} {'simulated mean':>16} {'k/(n+1)':>10}")
for k in range(1, n + 1):
    print(f"{k:>3} {ordered[:, k-1].mean():>16.4f} {k/(n+1):>10.4f}")

mx = ordered[:, -1]
print("\nP(max > 0.9) simulated:", round((mx > 0.9).mean(), 4))
print("P(max > 0.9) theory   :", round(1 - 0.9**n, 4))

The Beta connection, confirmed directly:

import numpy as np
from scipy.stats import beta

n, k = 5, 3
rng = np.random.default_rng(1)
ordered = np.sort(rng.uniform(0, 1, size=(200_000, n)), axis=1)
kth = ordered[:, k - 1]

dist = beta(k, n - k + 1)
print(f"{k}th of {n} order statistics vs Beta({k}, {n-k+1}):")
print("  mean:", round(kth.mean(), 4), "vs", round(dist.mean(), 4))
print("  var :", round(kth.var(), 4), "vs", round(dist.var(), 4))
for q in (0.1, 0.5, 0.9):
    print(f"  q={q}: {np.quantile(kth, q):.4f} vs {dist.ppf(q):.4f}")

How the maximum grows with sample size — the point about extremes drifting:

import numpy as np

rng = np.random.default_rng(2)
print(f"{'n':>8} {'E[max]':>10} {'P(max > 0.99)':>15}")
for n in (1, 5, 10, 100, 1_000):
    mx = rng.uniform(0, 1, size=(50_000, n)).max(axis=1)
    print(f"{n:>8} {mx.mean():>10.4f} {(mx > 0.99).mean():>15.4f}")

By n = 1000, a value in the top 1% is essentially guaranteed to appear — which is exactly why "the best run out of a thousand" tells you very little.

Your turn

1. n independent Exponential(\lambda) lifetimes. Show the minimum is Exponential(n\lambda).

2. For 5 Uniform(0,1) draws, find E[X_{(1)}] and P(X_{(1)} < 0.1).

3. A system needs all 4 components working. Components fail independently with lifetimes Uniform(0, 10) years. What's the expected time until the system fails?

Solutions

1. For an Exponential, 1 - F(x) = e^{-\lambda x}. Apply the minimum formula:

P(X_{(1)} > x) = \big[e^{-\lambda x}\big]^n = e^{-n\lambda x}

That is the survival function of an Exponential with rate n\lambda. So

X_{(1)} \sim \text{Exponential}(n\lambda), \qquad E[X_{(1)}] = \frac{1}{n\lambda}

Ten components each lasting 100 hours on average means the first failure arrives after 10 hours on average. Combining independent memoryless processes adds their rates — the same fact that makes superposed Poisson processes Poisson (§8).

2. Minimum with F(x) = x:

P(X_{(1)} > x) = (1-x)^5 \implies f_{X_{(1)}}(x) = 5(1-x)^4

E[X_{(1)}] = \frac{1}{n+1} = \frac{1}{6} \approx 0.167

P(X_{(1)} < 0.1) = 1 - (0.9)^5 = 0.41

Symmetric with the maximum example, as it must be — X_{(1)} for uniforms behaves like 1 - X_{(n)}.

3. "All 4 must work" means the system dies at the first component failure — so it's the minimum, not the maximum.

P(X_{(1)} > x) = \left(1 - \frac{x}{10}\right)^4, \qquad 0 \le x \le 10

For a non-negative variable, integrate the survival function:

E[X_{(1)}] = \int_0^{10}\left(1 - \frac{x}{10}\right)^4 dx = 10\int_0^1 u^4\,du = \frac{10}{5} = 2 \text{ years}

(Or use E[X_{(1)}] = \frac{10}{n+1} = \frac{10}{5} = 2 directly.)

Each component lasts 5 years on average, but the system lasts 2. Series reliability is governed by the weakest link, and adding components makes it worse — the same structural point as the redundancy example in §0, running in the opposite direction.

Check yourself in code

Simulate 5 uniform draws many times and confirm the expected values of all five order statistics match k/(n+1).

Print exactly this:

k=1 0.1667
k=2 0.3333
k=3 0.5
k=4 0.6667
k=5 0.8333
max tail matches theory: True

Print each theoretical mean k/(n+1) rounded to 4 decimal places. For the last line, compare the simulated P(X_{(5)} > 0.9) with 1 - 0.9^5 to 2 decimal places, using numpy.random.default_rng(0) and 200000 trials.

import numpy as np

rng = np.random.default_rng(0)
n, trials = 5, 200_000
ordered = np.sort(rng.uniform(0, 1, size=(trials, n)), axis=1)

for k in range(1, n + 1):
    print(f"k={k}", round(k / (n + 1), 4))

# Compare the simulated P(max > 0.9) against 1 - 0.9**n, to 2 decimal places.

import numpy as np

rng = np.random.default_rng(0)
n, trials = 5, 200_000
ordered = np.sort(rng.uniform(0, 1, size=(trials, n)), axis=1)

for k in range(1, n + 1):
    print(f"k={k}", round(k / (n + 1), 4))

simulated = (ordered[:, -1] > 0.9).mean()
theory = 1 - 0.9**n
print("max tail matches theory:", round(simulated, 2) == round(theory, 2))

---

Order statistics are what you sort a sample into. Maxima come from [F(x)]^n, minima from [1-F(x)]^n, and the general case is a multinomial count that turns out to be a Beta distribution for uniform samples. They're the right tool whenever the question is about the best, the worst, or the middle — and they explain why extremes get more extreme the more you look.

Next: the multivariate Normal, the one joint distribution you'll meet everywhere.