Probability theory

A seed cooperative in western Kenya tests a new maize variety: each seed germinates with probability about $0.9$, and nobody can say which seeds will fail. Yet of $10\,000$ seeds, the cooperative can promise that between roughly $8\,940$ and $9\,060$ will germinate, and be right about $95$ times in $100$. A call centre in Manila cannot predict when the next call will come, but it can staff its desks because the number of calls per minute follows a stable pattern. A detector in Chile's Atacama desert counts cosmic-ray muons arriving at random moments; the counts per hour cluster tightly around their average.

Individual outcomes are unpredictable; aggregates are predictable, within quantifiable limits. Probability theory explains this. It answers four questions about a random quantity:

• What is its long-run average? (Expectation.)
• How far does it typically stray from that average? (Variance.)
• Why do averages of many independent repetitions settle down, and how fast? (The law of large numbers.)
• Why do their fluctuations so often take the shape of the same bell curve? (The central limit theorem.)

This lesson assumes the rules of probability — conditional probability, independence of events, Bayes' theorem — from probability rules, and sums, integrals and series from single-variable calculus.

Drop a ball through $n$ rows of pins. At each pin it bounces left or right with probability $\tfrac12$, independently. Its final slot is the number $S_n$ of rightward bounces, a sum of $n$ independent coin flips. One ball's path is unpredictable; but drop a thousand and a shape appears.

Three things are visible already, and the rest of the lesson makes each precise:

1. The balls pile up around the middle, $n/2$: the expectation of $S_n$ is $n/2$.
2. The width of the pile grows, but only like $\sqrt n$: the standard deviation of $S_n$ is $\sqrt n/2$. So the proportion $S_n/n$ of right bounces has spread $1/(2\sqrt n)$, which shrinks: this is the law of large numbers.
3. The pile has a characteristic bell shape, whatever the details: this is the central limit theorem.

When $\Omega$ is countable we may take every subset to be an event. For continuous models (a uniformly random point of $[0, 1]$, say) this is impossible, and making the definitions fully rigorous requires measure theory; we use its results where needed and say so.

Two kinds of distribution cover most examples:

• $X$ is discrete if it takes values in a countable set $\{x_1, x_2, \dots\}$; its probability mass function is $p(x) = P(X = x)$, with $\sum_i p(x_i) = 1$.
• $X$ is continuous with density $f$ if $f \ge 0$ and $P(a \le X \le b) = \int_a^b f(x)\,dx$ for all $a \le b$; then $\int_{-\infty}^\infty f = 1$, $F_X' = f$ wherever $f$ is continuous, and $P(X = x) = 0$ for every single $x$.

Expanding the square and using linearity (proved below) gives the computing formula $\operatorname{Var} X = E[X^2] - (E[X])^2$, and also $\operatorname{Var}(aX + b) = a^2 \operatorname{Var} X$: shifting does not change spread, and scaling by $a$ scales the standard deviation by $|a|$.

Seven distributions do most of the work in applications. Write $q = 1 - p$.

The binomial row is derived with indicators in the next section. Here are the others.

So when there are many independent opportunities, each with a small chance, the number of successes is approximately Poisson. Ladislaus Bortkiewicz's 1898 study of deaths from horse kicks in Prussian cavalry corps — rare, independent accidents — is the classic fit; typing errors per page, mutations per genome per generation and goals per football match are often modelled the same way.

The indicator of an event $A$ is the random variable $\mathbf 1_A$ that equals $1$ on $A$ and $0$ off it. It is Bernoulli$(P(A))$, so $E[\mathbf 1_A] = P(A)$. Combining this with linearity gives one of the most useful tricks in probability: to find the expected number of things that happen, add up their probabilities — even when the events are dependent.

Two consequences we use constantly. The binomial variance: the indicators of independent trials are independent, each with variance $pq$, so $\operatorname{Var} S = npq$. The variance of a sample mean: if $X_1, \dots, X_n$ are pairwise uncorrelated with common variance $\sigma^2$, then $\bar X_n = \frac1n\sum X_i$ has $\operatorname{Var} \bar X_n = \frac{1}{n^2} \cdot n\sigma^2 = \frac{\sigma^2}{n}$. Averaging $n$ measurements divides the standard deviation by $\sqrt n$ — the $1/\sqrt n$ law behind every margin of error.

Chebyshev's inequality cannot be improved using the variance alone: for $k \ge 1$, the variable taking the values $\pm k$ with probability $\frac{1}{2k^2}$ each and $0$ otherwise has mean $0$, variance $1$ and $P(|X| \ge k) = \frac1{k^2}$ exactly.

The law of large numbers works by swamping, not compensation. After ten heads in a row the coin is not "due" a tail; the next thousand flips simply dilute the early surplus. The number of excess heads $S_n - n/2$ typically grows like $\sqrt n$; only the proportion settles.

How sharp is the bound? For $10\,000$ fair coin flips and $\varepsilon = 0.02$, Chebyshev gives $P(|\bar X_n - \tfrac12| \ge 0.02) \le \frac{1/4}{10\,000 \cdot 0.0004} = 0.0625$, while the true value is about $6 \times 10^{-5}$. Chebyshev uses only the variance, so it must cover the worst distribution with that variance; for bounded variables, Hoeffding's inequality (1963) gives the exponentially small bound $2e^{-2n\varepsilon^2}$. The second experiment below compares all of these with exact binomial tails.

Two stronger theorems, stated without proof:

• Khinchin's weak law (1929): for independent, identically distributed (i.i.d.) $X_i$, finite mean alone suffices.
• The strong law of large numbers (Kolmogorov, 1933; Etemadi's 1981 proof needs only pairwise independence): if the $X_i$ are pairwise independent, identically distributed, with $E|X_1| < \infty$, then $P(\bar X_n \to \mu) = 1$. The strong law says almost every individual sequence of averages converges; the weak law only says that, at each large $n$, a deviation is unlikely. It is proved in Mathlib (see below).

The hypothesis $E|X_1| < \infty$ matters: for i.i.d. Cauchy variables, $\bar X_n$ has exactly the same Cauchy distribution as $X_1$ for every $n$ — averaging does not help at all. The law of large numbers is also the principle behind Monte Carlo integration: to estimate $\int_0^1 g$, average $g(U_1), \dots, g(U_n)$ for independent uniform $U_i$; the error is of order $\sigma/\sqrt n$, whatever the dimension.

The law of large numbers says $\bar X_n - \mu \to 0$. The fluctuations have size $\sigma/\sqrt n$, so we magnify them by $\sqrt n/\sigma$ and ask about the shape of what we see.

How fast? The Berry–Esseen theorem (1941–42) quantifies the convergence: if $\rho = E|X_1 - \mu|^3 < \infty$, then

where $C$ is a universal constant (known to satisfy $C \le 0.4748$, Shevtsova 2011). Skewed distributions, with large $\rho/\sigma^3$, need larger $n$; the third experiment measures this for a die and for a coin with $P(\text{heads}) = 0.1$.

What the theorem does not say.

1. It does not say that data become normal. A sample of $10\,000$ household incomes is as skewed as the population it came from; what is approximately normal is the distribution of the average (or sum) of many independent draws, across hypothetical repetitions.
2. It needs finite variance. Averages of i.i.d. Cauchy variables are Cauchy for every $n$, and heavy-tailed sums converge (when they converge at all) to other "stable" laws, with a different scaling than $\sqrt n$.
3. It needs independence, or at least weak dependence. Sums of strongly dependent variables — the same measurement error repeated, contagion in financial losses — need not look normal.
4. It controls absolute, not relative, error. Far in the tails, where $\Phi$ is tiny, the normal approximation may be off by orders of magnitude; such "large deviation" probabilities need other tools (Cramér, 1938).

The theorem also explains why the normal distribution is everywhere: a measurement influenced by many small, roughly independent effects (genes and environment on height, many small errors in an astronomical observation) is approximately normal — a heuristic, valid exactly to the extent that its hypotheses are.

Mathlib builds probability on measure theory: a probability space is a measure μ with IsProbabilityMeasure μ, a random variable is a function X : Ω → ℝ, its expectation is the integral ∫ ω, X ω ∂μ (written μ[X]), its variance is Var[X; μ], and X ⟂ᵢ[μ] Y says that X and Y are independent. "Has finite variance" is MemLp X 2 μ. Mathlib's library has grown to include the strong law of large numbers and, recently, the central limit theorem, so this lesson can connect its whole arc to formal proofs.

Linearity and indicators. The expected number of events that occur is the sum of their probabilities, with no independence assumed; and the "fixed points of a shuffle" example, proved through Burnside's lemma (the total number of fixed points of a group action is the number of orbits times the size of the group — here there is one orbit):

/-- **Linearity of expectation for indicators.** For any events `A i` (`i ∈ s`), the expected
number of them that occur is the sum of their probabilities: `E[∑ 1_{A i}] = ∑ P(A i)`.
No independence is assumed. Proof: the integral of a finite sum is the sum of the integrals
(`integral_finsetSum`), and the integral of an indicator is the measure of the set
(`integral_indicator_one`). -/
theorem expectation_sum_indicators [IsProbabilityMeasure μ] {ι : Type*} (s : Finset ι)
    (A : ι → Set Ω) (hA : ∀ i, MeasurableSet (A i)) :
    ∫ ω, (∑ i ∈ s, (A i).indicator (1 : Ω → ℝ) ω) ∂μ = ∑ i ∈ s, μ.real (A i) := by
  rw [integral_finsetSum]
  · exact sum_congr rfl fun i _ => integral_indicator_one (hA i)
  · intro i _
    exact (integrable_const (1 : ℝ)).indicator (hA i)

/-- **A random permutation has, on average, exactly one fixed point.** Summing the number of
fixed points over all permutations of a non-empty finite set gives the number of permutations.
The lesson proves this with indicators and linearity of expectation; here it follows from
Mathlib's **Burnside lemma** (`MulAction.sum_card_fixedBy_eq_card_orbits_mul_card_group`:
the total number of fixed points is (number of orbits) × (size of the group)), because the
permutations of `α` move any point to any other, so there is exactly one orbit. -/
theorem perm_fixed_points_average (α : Type*) [Fintype α] [DecidableEq α] [Nonempty α] :
    (∑ σ : Equiv.Perm α, Fintype.card (MulAction.fixedBy α σ)) =
      Fintype.card (Equiv.Perm α) := by
  classical
  rw [MulAction.sum_card_fixedBy_eq_card_orbits_mul_card_group]
  have h1 : Fintype.card (MulAction.orbitRel.Quotient (Equiv.Perm α) α) = 1 := by
    rw [Fintype.card_eq_one_iff]
    obtain ⟨a⟩ := ‹Nonempty α›
    refine ⟨Quotient.mk _ a, fun q => ?_⟩
    induction q using Quotient.inductionOn with
    | h b =>
      apply Quotient.sound
      exact MulAction.exists_smul_eq (Equiv.Perm α) a b |>.imp fun g hg => hg
  rw [h1, one_mul]

Distributions. Mathlib has the binomial, geometric, Poisson, exponential and Gaussian distributions. We cite the binomial mean and the Gaussian mean and variance, prove the geometric mean from the differentiated geometric series exactly as in the lesson, prove the memoryless property from Mathlib's distribution function of the exponential law, and cite the Poisson limit theorem:

/-- **The binomial mean.** If `X` has the binomial distribution `Bin(n, p)`, then `E[X] = np`.
This is Mathlib's `ProbabilityTheory.integral_of_hasLaw_binomial`, proved there with the
identity `k·C(n, k) = n·C(n − 1, k − 1)` and the binomial theorem; the lesson's proof writes
`X` as a sum of `n` indicators instead. -/
theorem binomial_mean {n : ℕ} {p : unitInterval} {X : Ω → ℝ}
    (hX : HasLaw X (Bin(ℝ, n, p)) μ) : μ[X] = (p : ℝ) * n :=
  integral_of_hasLaw_binomial hX

/-- **The geometric mean.** If each trial succeeds with probability `p ∈ (0, 1]`, the number
of trials up to and including the first success takes the value `k + 1` with probability
`p (1 − p)ᵏ`, and its expectation is `∑ (k + 1) p (1 − p)ᵏ = 1/p`. Proof, as in the lesson:
with `r = 1 − p`, `∑ k rᵏ = r/(1 − r)²` (Mathlib's `hasSum_coe_mul_geometric_of_norm_lt_one`,
"differentiate the geometric series") and `∑ rᵏ = 1/(1 − r)`, so
`∑ (k + 1) rᵏ = 1/(1 − r)² = 1/p²`; multiply by `p`. -/
theorem geometric_mean {p : ℝ} (hp0 : 0 < p) (hp1 : p ≤ 1) :
    HasSum (fun k : ℕ => ((k : ℝ) + 1) * (p * (1 - p) ^ k)) (1 / p) := by
  have hr : ‖1 - p‖ < 1 := by
    rw [Real.norm_eq_abs, abs_lt]
    constructor <;> linarith
  have h1 := hasSum_coe_mul_geometric_of_norm_lt_one hr
  have h2 := hasSum_geometric_of_norm_lt_one hr
  have h3 := (h1.add h2).mul_left p
  have hp : p ≠ 0 := hp0.ne'
  convert h3 using 1
  · funext k
    ring
  · simp only [sub_sub_cancel]
    field_simp
    ring

/-- **The exponential distribution is memoryless.** For Mathlib's exponential distribution
`expMeasure r` with rate `r > 0`, the survival function `P(T > x) = 1 − F(x)` satisfies
`P(T > s + t) = P(T > s) · P(T > t)` for `s, t ≥ 0`; equivalently
`P(T > s + t | T > s) = P(T > t)`. Proof: Mathlib's `cdf_expMeasure_eq` gives
`1 − F(x) = e^{−rx}` for `x ≥ 0`, and `e^{−r(s + t)} = e^{−rs} e^{−rt}`. -/
theorem exponential_memoryless {r : ℝ} (hr : 0 < r) {s t : ℝ} (hs : 0 ≤ s) (ht : 0 ≤ t) :
    1 - cdf (expMeasure r) (s + t) = (1 - cdf (expMeasure r) s) * (1 - cdf (expMeasure r) t) := by
  simp only [cdf_expMeasure_eq hr, add_nonneg hs ht, hs, ht, ↓reduceIte, sub_sub_cancel]
  rw [← Real.exp_add]
  ring_nf

/-- The normal distribution `N(m, v)` (Mathlib's `gaussianReal m v`) has mean `m`. This is
Mathlib's `ProbabilityTheory.integral_id_gaussianReal`. -/
theorem gaussian_mean (m : ℝ) (v : NNReal) : ∫ x, x ∂(gaussianReal m v) = m :=
  integral_id_gaussianReal

/-- The normal distribution `N(m, v)` has variance `v`. This is Mathlib's
`ProbabilityTheory.variance_id_gaussianReal`. -/
theorem gaussian_variance (m : ℝ) (v : NNReal) : Var[id; gaussianReal m v] = v :=
  variance_id_gaussianReal

/-- **Poisson limit theorem.** If `n · pₙ → λ`, then the binomial probabilities
`C(n, k) pₙᵏ (1 − pₙ)ⁿ⁻ᵏ` tend to the Poisson probability `e^{−λ} λᵏ / k!` for each fixed `k`.
This is Mathlib's `ProbabilityTheory.tendsto_choose_mul_pow_of_tendsto_mul_atTop`; the lesson
proves the case `pₙ = λ/n`. -/
theorem poisson_limit {p : ℕ → ℝ} {lam : ℝ} (k : ℕ)
    (hp : Tendsto (fun n : ℕ => (n : ℝ) * p n) atTop (𝓝 lam)) :
    Tendsto (fun n : ℕ => ((n.choose k : ℕ) : ℝ) * p n ^ k * (1 - p n) ^ (n - k)) atTop
      (𝓝 (Real.exp (-lam) * lam ^ k / (k.factorial : ℝ))) :=
  tendsto_choose_mul_pow_of_tendsto_mul_atTop k hp

Independence. The product rule and the additivity of variance, from Mathlib, and the uncorrelated-but-dependent example, computed in exact rational arithmetic (a three-point distribution needs no measure theory, so we average over Fin 3 directly):

/-- **Independent variables multiply in expectation:** if `X` and `Y` are independent then
`E[XY] = E[X] E[Y]`. This is Mathlib's `IndepFun.integral_fun_mul_eq_mul_integral`.
Without independence (or at least zero covariance) the formula fails: see
`uncorrelated_not_independent` and the lesson's example `X = Y` a fair coin. -/
theorem expectation_mul_indep {X Y : Ω → ℝ} (hXY : X ⟂ᵢ[μ] Y)
    (hX : AEStronglyMeasurable X μ) (hY : AEStronglyMeasurable Y μ) :
    ∫ ω, X ω * Y ω ∂μ = (∫ ω, X ω ∂μ) * ∫ ω, Y ω ∂μ :=
  hXY.integral_fun_mul_eq_mul_integral hX hY

/-- **Variances of independent variables add:** `Var(X + Y) = Var X + Var Y` when `X ⟂ Y` and
both have finite variance. This is Mathlib's `ProbabilityTheory.IndepFun.variance_add`.
(Mathlib also has the general formula with the covariance term.) -/
theorem variance_add_indep {X Y : Ω → ℝ} (hX : MemLp X 2 μ) (hY : MemLp Y 2 μ)
    (hXY : X ⟂ᵢ[μ] Y) : Var[X + Y; μ] = Var[X; μ] + Var[Y; μ] :=
  hXY.variance_add hX hY

/-- The average of a function on three equally likely outcomes `0, 1, 2`: the expectation for
the uniform distribution on `Fin 3`, computed exactly in `ℚ`. -/
def avg3 (f : Fin 3 → ℚ) : ℚ := (∑ i, f i) / 3

/-- **Uncorrelated does not imply independent.** Let `X` be uniform on `{−1, 0, 1}` and
`Y = X²`. Then `E[XY] = E[X] E[Y]` (both sides are `0`), yet `X` and `Y` are dependent:
`P(X = 1 and Y = 0) = 0` while `P(X = 1) · P(Y = 0) = 1/3 · 1/3`. The probabilities of events
are averages of indicator functions. Checked by exact rational arithmetic. -/
theorem uncorrelated_not_independent :
    let X : Fin 3 → ℚ := ![-1, 0, 1]
    let Y : Fin 3 → ℚ := fun i => X i ^ 2
    avg3 (fun i => X i * Y i) = avg3 X * avg3 Y ∧
      avg3 (fun i => if X i = 1 ∧ Y i = 0 then 1 else 0) ≠
        avg3 (fun i => if X i = 1 then 1 else 0) * avg3 (fun i => if Y i = 0 then 1 else 0) := by
  simp only [avg3, Fin.sum_univ_three]
  norm_num

Markov, Chebyshev and the weak law. Markov's inequality is cited from Mathlib. Chebyshev's inequality is then proved from it, line by line as in the human proof, and compared with Mathlib's own version (which states the probability as an extended non-negative real, hence the ENNReal.ofReal):

/-- **Markov's inequality:** for a non-negative integrable `f` and any `a`,
`a · P(f ≥ a) ≤ E[f]`. This is Mathlib's `mul_meas_ge_le_integral_of_nonneg`. (The measure
theory lesson proves the version for `[0, ∞]`-valued functions.) -/
theorem markov {f : Ω → ℝ} (hf : 0 ≤ᵐ[μ] f) (hint : Integrable f μ) (a : ℝ) :
    a * μ.real {ω | a ≤ f ω} ≤ ∫ ω, f ω ∂μ :=
  mul_meas_ge_le_integral_of_nonneg hf hint a

/-- **Chebyshev's inequality**, proved from Markov's as in the lesson: for `X` with finite
variance and `c > 0`, `c² · P(|X − E X| ≥ c) ≤ Var X`. Proof: `|X − EX| ≥ c` happens exactly
when `(X − EX)² ≥ c²`; apply Markov's inequality to the non-negative variable `(X − EX)²`,
whose expectation is `Var X`. -/
theorem chebyshev [IsProbabilityMeasure μ] {X : Ω → ℝ} (hX : MemLp X 2 μ) {c : ℝ}
    (hc : 0 < c) : c ^ 2 * μ.real {ω | c ≤ |X ω - μ[X]|} ≤ Var[X; μ] := by
  set m := μ[X]
  -- The events `{c ≤ |X − m|}` and `{c² ≤ (X − m)²}` are equal.
  have hset : {ω | c ≤ |X ω - m|} = {ω | c ^ 2 ≤ (X ω - m) ^ 2} := by
    ext ω
    simp only [Set.mem_ofPred_eq]
    rw [sq_le_sq, abs_of_pos hc]
  -- `(X − m)²` is non-negative and integrable.
  have hint : Integrable (fun ω => (X ω - m) ^ 2) μ :=
    (hX.sub (memLp_const m)).integrable_sq
  have hnn : 0 ≤ᵐ[μ] fun ω => (X ω - m) ^ 2 := Eventually.of_forall fun ω => sq_nonneg _
  rw [hset, variance_eq_integral hX.aestronglyMeasurable.aemeasurable]
  exact mul_meas_ge_le_integral_of_nonneg hnn hint (c ^ 2)