Probability

Gamblers in Galileo's time played a game with three dice and bet on the total. Experienced players had noticed that a total of 10 turns up a little more often than a total of 9. That seemed strange, because each total can be made in six ways:

• $9 = 1 + 2 + 6 = 1 + 3 + 5 = 1 + 4 + 4 = 2 + 2 + 5 = 2 + 3 + 4 = 3 + 3 + 3$;
• $10 = 1 + 3 + 6 = 1 + 4 + 5 = 2 + 2 + 6 = 2 + 3 + 5 = 2 + 4 + 4 = 3 + 3 + 4$.

Were the players imagining it? Galileo explained the answer around 1620, and by the end of this lesson you will be able to explain it too — and even to check it with a computer and with Lean.

Probability is the mathematics of uncertainty. It helps a weather service say “70% chance of rain”, a doctor judge a test result, and an engineer design a network that rarely fails. It begins with a simple idea: count the equally likely possibilities.

We measure chance on a scale from $0$ to $1$:

• $0$ means impossible (rolling a $7$ on an ordinary die);
• $1$ means certain (rolling a number less than $7$);
• $\tfrac12$ means as likely as not (a fair coin landing heads).

When an experiment has a few outcomes that are all equally likely, the probability of an event is the fraction of outcomes that make it happen. A fair die has six equally likely faces, and three are even, so the probability of an even number is $\tfrac36 = \tfrac12$.

But “equally likely” is an assumption that must be checked. The spinner below has three outcomes — but they are not equally likely.

For one die, $S = \{1, 2, 3, 4, 5, 6\}$, and “an even number” is the event $A = \{2, 4, 6\}$.

Here $|A|$ is the number of elements of $A$, as in the lesson on sets. The event “not $A$”, written $A^c$ (the complement), is the set of outcomes of $S$ that are not in $A$.

Probability is a statement about the experiment; relative frequency is a record of what happened. The two are linked, but they are not the same thing.

Throughout, $S$ is a finite, non-empty sample space of equally likely outcomes.

The complement rule is often the easy way in. The probability of at least one six in four rolls of a die is $1 - P(\text{no six}) = 1 - \tfrac{5^4}{6^4} = 1 - \tfrac{625}{1296} = \tfrac{671}{1296} \approx 0.518$.

For one die, $P(\text{even or more than } 4) = \tfrac36 + \tfrac26 - \tfrac16 = \tfrac46$: the outcome $6$ is both even and more than $4$, and must not be counted twice.

This is why we multiply along the branches of a tree diagram. For a double six, $P = \tfrac16 \times \tfrac16 = \tfrac1{36}$. (When the second stage does depend on the first — drawing two cards without putting the first back — the rule needs changing; that is the subject of conditional probability.)

We do not prove this here; a precise version was first proved by Jacob Bernoulli (published 1713), and you will meet proofs in probability theory. Notice what it does not say: it does not say the relative frequency ever equals the probability exactly, and it does not say that the coin “corrects” a run of tails. Earlier results are not undone; they are swamped by the many trials that follow.

In Lean, a finite sample space of equally likely outcomes is a list S, and an event is a test A that answers true or false for each outcome. The probability is computed exactly as a rational number.

/-- The number of outcomes in the list `S` that belong to the event `A`. -/
def count (S : List α) (A : α → Bool) : Nat := (S.filter A).length

/-- **Probability with equally likely outcomes:** `P(A) = |A| / |S|`, the number of favourable
outcomes divided by the number of possible outcomes. -/
def prob (S : List α) (A : α → Bool) : Rat := (count S A : Rat) / (S.length : Rat)

Theorem 1. The proofs of count_compl and count_union_inter go through the list one outcome at a time (induction): each outcome is counted on both sides of the equation in the same way.

/-- Probabilities are never negative: a count divided by a count. -/
theorem prob_nonneg (S : List α) (A : α → Bool) : 0 ≤ prob S A := by
  unfold prob
  rw [Rat.div_def]
  apply Rat.mul_nonneg (Rat.natCast_nonneg)
  by_cases h : (S.length : Rat) = 0
  · rw [h, Rat.inv_zero]
    exact Rat.le_refl
  · have hpos : (0 : Rat) < S.length := Rat.lt_of_le_of_ne Rat.natCast_nonneg (Ne.symm h)
    exact Rat.le_of_lt (Rat.inv_pos.mpr hpos)

/-- Probabilities are at most `1`: the event has at most as many outcomes as `S`. -/
theorem prob_le_one (S : List α) (A : α → Bool) : prob S A ≤ 1 := by
  unfold prob
  have hle : (count S A : Rat) ≤ S.length := Rat.natCast_le_natCast.mpr (count_le_length S A)
  by_cases h : S.length = 0
  · rw [h]
    simp [Rat.div_def]
    decide
  · have hpos : (0 : Rat) < S.length := Rat.natCast_pos.mpr (Nat.pos_of_ne_zero h)
    -- `c / n ≤ 1` means "not `1 < c / n`", i.e. not `1 · n < c`.
    rw [← Rat.not_lt, Rat.lt_div_iff hpos]
    grind

/-- The certain event (every outcome) has probability `1`, when there is at least one outcome. -/
theorem prob_certain (S : List α) (hS : S ≠ []) : prob S (fun _ => true) = 1 := by
  unfold prob count
  simp only [List.filter_eq_self.mpr (fun _ _ => rfl)]
  have hpos : (S.length : Rat) ≠ 0 := by
    have : S.length ≠ 0 := by
      intro h0
      exact hS (List.eq_nil_of_length_eq_zero h0)
    simpa using this
  rw [Rat.div_def]
  exact Rat.mul_inv_cancel _ hpos

/-- Every outcome is either in `A` or in "not `A`", never both:
`|A| + |not A| = |S|`. Proof by induction on the list of outcomes. -/
theorem count_compl (S : List α) (A : α → Bool) :
    count S A + count S (fun x => !A x) = S.length := by
  induction S with
  | nil => rfl
  | cons a S ih =>
    unfold count at ih ⊢
    cases h : A a <;> simp [h] <;> omega

/-- **The complement rule:** `P(not A) = 1 - P(A)` (for a non-empty sample space).
Divide `|A| + |not A| = |S|` by `|S|`. -/
theorem prob_compl (S : List α) (hS : S ≠ []) (A : α → Bool) :
    prob S (fun x => !A x) = 1 - prob S A := by
  have hlen : S.length ≠ 0 := by
    intro h0
    exact hS (List.eq_nil_of_length_eq_zero h0)
  have hne : (S.length : Rat) ≠ 0 := by simpa using hlen
  have hsum : (count S A : Rat) + (count S (fun x => !A x) : Rat) = S.length := by
    have := count_compl S A
    have hc := congrArg (fun n : Nat => (n : Rat)) this
    simpa using hc
  unfold prob
  grind

Theorem 2, the addition rule:

/-- **Counting a union.** Outcomes in `A` or `B`, plus outcomes in both, equals outcomes in `A`
plus outcomes in `B`: `|A ∪ B| + |A ∩ B| = |A| + |B|`. Induction on the list, checking the
four cases for whether the first outcome is in `A` and in `B`. -/
theorem count_union_inter (S : List α) (A B : α → Bool) :
    count S (fun x => A x || B x) + count S (fun x => A x && B x) = count S A + count S B := by
  induction S with
  | nil => rfl
  | cons a S ih =>
    unfold count at ih ⊢
    cases hA : A a <;> cases hB : B a <;> simp [hA, hB] <;> omega

/-- **The addition rule:** `P(A or B) = P(A) + P(B) - P(A and B)`. Divide `count_union_inter`
by `|S|`. When `A` and `B` cannot happen together, this is `P(A) + P(B)`. -/
theorem prob_union (S : List α) (A B : α → Bool) :
    prob S (fun x => A x || B x) = prob S A + prob S B - prob S (fun x => A x && B x) := by
  have h := congrArg (fun n : Nat => (n : Rat)) (count_union_inter S A B)
  push_cast at h
  unfold prob
  rw [Rat.div_def, Rat.div_def, Rat.div_def, Rat.div_def]
  grind

Theorem 3, the product rule. The sample space of two stages is the list of all pairs.

/-- All ordered pairs `(a, b)` with `a` from `S` and `b` from `T`: the sample space of a
two-stage experiment, as in a table or a tree diagram. -/
def pairs (S : List α) (T : List β) : List (α × β) := S.flatMap (fun a => T.map (fun b => (a, b)))

/-- A two-stage experiment with `|S|` and `|T|` outcomes has `|S| · |T|` outcomes. -/
theorem length_pairs (S : List α) (T : List β) : (pairs S T).length = S.length * T.length := by
  induction S with
  | nil => simp [pairs]
  | cons a S ih =>
    simp only [pairs, List.flatMap_cons, List.length_append, List.length_map] at ih ⊢
    rw [ih, List.length_cons, Nat.succ_mul, Nat.add_comm]

/-- The pairs with first outcome in `A` and second outcome in `B` number `|A| · |B|`. -/
theorem count_pairs (S : List α) (T : List β) (A : α → Bool) (B : β → Bool) :
    count (pairs S T) (fun p => A p.1 && B p.2) = count S A * count T B := by
  induction S with
  | nil => simp [pairs, count]
  | cons a S ih =>
    unfold count pairs at ih ⊢
    simp only [List.flatMap_cons, List.filter_append, List.length_append] at ih ⊢
    rw [ih]
    cases hA : A a
    · simp [hA, List.filter_map, Function.comp_def]
    · simp [hA, List.filter_map, Function.comp_def]
      rw [Nat.succ_mul, Nat.add_comm]

/-- **The product rule for independent stages:** if the second stage does not depend on the
first, then `P(A at stage 1 and B at stage 2) = P(A) · P(B)` — multiply along the branches of
the tree. -/
theorem prob_pairs (S : List α) (T : List β) (A : α → Bool) (B : β → Bool) :
    prob (pairs S T) (fun p => A p.1 && B p.2) = prob S A * prob T B := by
  unfold prob
  rw [count_pairs, length_pairs]
  push_cast
  by_cases hS : (S.length : Rat) = 0
  · rw [hS]
    simp [Rat.div_def]
  · by_cases hT : (T.length : Rat) = 0
    · rw [hT]
      simp [Rat.div_def]
    · grind

Finally, Lean computes the dice examples and Galileo's counts by checking all the outcomes:

/-- The six equally likely outcomes of a fair die. -/
def die : List Nat := [1, 2, 3, 4, 5, 6]

/-- The 36 equally likely outcomes of rolling two dice (a red one and a blue one). -/
def twoDice : List (Nat × Nat) := pairs die die

/-- With two dice, a sum of `7` has probability `6/36 = 1/6`. Computed by Lean. -/
theorem prob_sum_seven : prob twoDice (fun p => p.1 + p.2 == 7) = 1 / 6 := by
  decide +kernel

/-- The eleven possible sums are **not** equally likely: a sum of `2` has probability `1/36`,
but a sum of `7` has probability `1/6`. -/
theorem sums_not_equally_likely :
    prob twoDice (fun p => p.1 + p.2 == 2) = 1 / 36 ∧
    prob twoDice (fun p => p.1 + p.2 == 2) ≠ prob twoDice (fun p => p.1 + p.2 == 7) := by
  decide +kernel

/-- The 216 equally likely outcomes of rolling three dice. -/
def threeDice : List (Nat × Nat × Nat) := pairs die (pairs die die)

/-- **Galileo's dice problem.** With three dice, a sum of `10` can be made in `27` ways but a
sum of `9` in only `25`, although each sum has six partitions into three numbers from 1 to 6.
Counted by Lean over all 216 outcomes. -/
theorem galileo :
    count threeDice (fun t => t.1 + t.2.1 + t.2.2 == 9) = 25 ∧
    count threeDice (fun t => t.1 + t.2.1 + t.2.2 == 10) = 27 := by
  decide

/-- **d'Alembert's question.** Toss two coins (`true` = heads). The outcomes HH, HT, TH, TT
are equally likely, so "at least one head" has probability `3/4`, not `2/3`. -/
theorem at_least_one_head : prob (pairs [true, false] [true, false]) (fun p => p.1 || p.2) = 3 / 4 := by
  decide +kernel

What is not formalised: the law of large numbers, and the modelling assumption that real dice are fair. Lean proves facts about the mathematical model; whether a real die fits the model is a question for experiment.

"""Famous dice and coin questions, settled by listing every outcome.

When all outcomes are equally likely, P(A) = |A| / |S|. This experiment
lists sample spaces with itertools.product and counts exactly with Fraction:
the sums of two dice; Galileo's question about three dice (why does 10 come
up more often than 9?); d'Alembert's mistake about two coins; and the
Chevalier de Mere's two bets. Where listing is impossible (36**24 outcomes)
it uses the complement rule P(not A) = 1 - P(A) instead, and says so.
"""

from collections import Counter
from fractions import Fraction
from itertools import product

DIE = range(1, 7)

def sum_counts(dice):
    """Counter: sum -> number of ordered outcomes, for `dice` fair dice."""
    return Counter(sum(roll) for roll in product(DIE, repeat=dice))

def partitions(total, parts=3):
    """Unordered ways to write total as a sum of `parts` numbers from 1-6."""
    return sorted({tuple(sorted(roll)) for roll in product(DIE, repeat=parts)
                   if sum(roll) == total})

def probability(event, outcomes):
    """Exact P(event) when every outcome in the list is equally likely."""
    outcomes = list(outcomes)
    return Fraction(sum(1 for o in outcomes if event(o)), len(outcomes))

def at_least_one_six(rolls):
    """P(at least one 6 in `rolls` rolls), by listing all 6**rolls outcomes."""
    return probability(lambda r: 6 in r, product(DIE, repeat=rolls))

def at_least_one_double_six(throws):
    """P(at least one double six in `throws` throws of two dice).

    Uses the complement rule: no double six in one throw has probability
    35/36, and independent throws multiply.
    """
    return 1 - Fraction(35, 36) ** throws

def main():
    two = sum_counts(2)
    print("Two dice: ordered outcomes giving each sum (out of 36)")
    print("  sum: " + " ".join(f"{s:2d}" for s in range(2, 13)))
    print("  ways:" + " ".join(f"{two[s]:2d}" for s in range(2, 13)))
    print()
    three = sum_counts(3)
    print("Galileo: three dice, 216 equally likely ordered outcomes")
    for total in (9, 10):
        p = Fraction(three[total], 216)
        print(f"  sum {total:2d}: {len(partitions(total))} partitions, "
              f"{three[total]} outcomes, P = {p} = {float(p):.4f}")
    print()
    coins = list(product("HT", repeat=2))
    p = probability(lambda c: "H" in c, coins)
    print("d'Alembert: two coins, outcomes " + " ".join("".join(c)
                                                       for c in coins))
    print(f"  P(at least one head) = {p}, not 2/3: of his three cases")
    print("  H, TH, TT, the first (H) covers two outcomes, HH and HT.")
    print()
    four = at_least_one_six(4)
    print("de Mere's bets:")
    print(f"  at least one 6 in 4 rolls: listing 1296 outcomes gives {four}"
          f" = {float(four):.4f}")
    print(f"  complement rule 1 - (5/6)^4 = {1 - Fraction(5, 6) ** 4}"
          " (the same)")
    p24 = at_least_one_double_six(24)
    print(f"  at least one double 6 in 24 throws: 1 - (35/36)^24 = "
          f"{float(p24):.4f}")
    print("  (36^24 outcomes are far too many to list: complement rule only)")
    print()
    print("Epistemic status: exact counts over complete lists of outcomes,")
    print("assuming fair dice and coins with independent throws. Whether real")
    print("dice behave like this is a question about the world, tested by")
    print("experiment (see the simulation), not settled by counting.")

if __name__ == "__main__":
    main()

The Chevalier de Méré, a French gambler, is said to have asked Blaise Pascal about these two bets in 1654: the first wins slightly more than half the time, the second slightly less.

"""Simulated dice and coins against exact probabilities.

A seeded pseudo-random generator rolls two dice thousands of times. The
relative frequency of each sum is compared with its exact probability
(ways / 36). Repeating with more and more rolls shows the frequencies
settling down - the law of large numbers - with typical errors shrinking
roughly like 1/sqrt(n). Two misconceptions are then tested: that a coin is
"due" heads after three tails (the gambler's fallacy), and that a
probability of 1/6 means exactly one six in every six rolls.
"""

import math
import random
from collections import Counter

EXACT = {s: (6 - abs(s - 7)) / 36 for s in range(2, 13)}

def roll_two_dice(rng, n):
    """Counter of sums for n rolls of two dice."""
    return Counter(rng.randint(1, 6) + rng.randint(1, 6) for _ in range(n))

def largest_error(counts, n):
    """Largest |relative frequency - exact probability| over all sums."""
    return max(abs(counts[s] / n - EXACT[s]) for s in range(2, 13))

def heads_after_tails(flips, run):
    """Among flips that follow `run` tails in a row, the fraction of heads.

    Returns (fraction, number of such flips). Flips are 'H' or 'T'.
    """
    after = [flips[i] for i in range(run, len(flips))
             if all(f == "T" for f in flips[i - run:i])]
    if not after:
        return None, 0
    return after.count("H") / len(after), len(after)

def exactly_one_six_rate(rng, blocks):
    """Fraction of blocks of six rolls that contain exactly one six."""
    hits = 0
    for _ in range(blocks):
        sixes = sum(1 for _ in range(6) if rng.randint(1, 6) == 6)
        hits += sixes == 1
    return hits / blocks

def main():
    rng = random.Random(2024)
    counts = roll_two_dice(rng, 3600)
    print("3,600 rolls of two dice: relative frequency vs exact probability")
    print("  sum   count   frequency   exact (ways/36)")
    for s in range(2, 13):
        print(f"  {s:3d}   {counts[s]:5d}   {counts[s] / 3600:9.4f}   "
              f"{EXACT[s]:.4f}")
    print()
    print("More rolls, smaller errors (largest error over the 11 sums):")
    for n in (60, 600, 6000, 60000):
        err = largest_error(roll_two_dice(rng, n), n)
        print(f"  n = {n:6d}: largest error {err:.4f}   1/sqrt(n) = "
              f"{1 / math.sqrt(n):.4f}")
    print()
    flips = ["H" if rng.random() < 0.5 else "T" for _ in range(100000)]
    frac, cases = heads_after_tails(flips, 3)
    print(f"Coin: after three tails in a row ({cases} times in 100,000 "
          f"flips),")
    print(f"  the next flip was heads {frac:.3f} of the time. Not 'due'.")
    rate = exactly_one_six_rate(rng, 20000)
    print(f"Blocks of six rolls with exactly one six: {rate:.3f} "
          f"(exact (5/6)^5 = {(5 / 6) ** 5:.3f})")
    print()
    print("Epistemic status: statistical evidence from one seeded run of a")
    print("pseudo-random generator, compared with exact values. Frequencies")
    print("vary from run to run; the law of large numbers (a theorem) says")
    print("they settle near the probabilities as n grows, not that they hit")
    print("them exactly.")

if __name__ == "__main__":
    main()

The simulated frequencies are close to the exact probabilities but never exactly equal, and the largest error shrinks as $n$ grows — roughly in step with $\tfrac{1}{\sqrt n}$, so a hundred times as many rolls makes the error about ten times smaller. This is statistical evidence: another run would give slightly different numbers. The exact counts are a different kind of knowledge, true for the model of fair dice.

On the two-dice tab, compare the histogram with the exact triangle of ways out of $36$. Increase the number of rolls: the bars hug the triangle more closely, but never perfectly. On the coin tab, watch the running proportion of heads wander and then settle near $\tfrac12$, mostly inside the band $\tfrac12 \pm \tfrac{1}{\sqrt n}$. Press “New run” several times: every run is different, and every run settles.

Probability speaks the language of sets. Outcomes are elements, the sample space is the universe, and events are subsets. “Or” is union, “and” is intersection, “not” is complement — so a Venn diagram is a picture of events. The addition rule is the counting formula $|A \cup B| = |A| + |B| - |A \cap B|$, and the product rule is $|A \times B| = |A| \cdot |B|$, both from the sets lesson.

Optional, and never graded: a question that leads somewhere further.

• India. In the Mahābhārata, King Ṛtuparṇa estimates the number of leaves and fruits on a whole tree by counting those on one branch, and connects this skill with knowledge of dice. The philosopher Ian Hacking has pointed to this story as an early link between sampling and chance.
• The Islamic world. In the 9th century, al-Kindī described how to break secret codes by comparing the relative frequencies of letters in a message with their frequencies in ordinary Arabic text.
• Italy. Gerolamo Cardano's Book on Games of Chance (written in the 16th century, published in 1663) compared the number of favourable outcomes with the number of all equally likely outcomes. Galileo's short note on three dice (around 1620) explained the 9 versus 10 puzzle by counting ordered outcomes.
• France and the Netherlands. In 1654 Blaise Pascal and Pierre de Fermat exchanged letters on how to divide the stakes of an unfinished game, often seen as the start of mathematical probability. Christiaan Huygens wrote the first printed textbook on probability (1657). In 1754 Jean le Rond d'Alembert argued that the chance of at least one head in two tosses is $\tfrac23$ — a famous mistake.
• Switzerland and Russia. Jacob Bernoulli's Ars Conjectandi (1713) proved the first law of large numbers. Andrey Kolmogorov's book of 1933 gave probability the axioms used today.

• Ian Hacking, An Introduction to Probability and Inductive Logic (2001).
• Charles M. Grinstead and J. Laurie Snell, Introduction to Probability (1997), with many simulations.
• William Feller, An Introduction to Probability Theory and Its Applications, volume 1 (1950) — a classic, for later.
• Ian Hacking, The Emergence of Probability (1975), on the history of the subject.