01 - Probability Spaces and Kolmogorov's Axioms

A diagnostic test comes back positive: how likely is it that you are really ill? And when you say that drawing a card from a deck of 52 gives you "one chance in four" of getting a heart, what exactly are you assigning, and to what? Probability was born to answer questions like these, but answering them consistently takes a small set of rules, fixed once and for all.

The answer of this article is that three axioms are enough, formulated by Andrey Kolmogorov in 1933 in the Grundbegriffe der Wahrscheinlichkeitsrechnung (Springer, Berlin), and that everything else, from Bayes' formula to independence, is a theorem. To follow along you only need sets, numerical series and a little arithmetic with fractions; measure theory is not required, and where it would be needed I say so. Every result below is computed in two ways, with exact fractions and with a fixed-seed simulation, and the two agree within the error I state.

What a probability space is

A random experiment is described by three objects, which together form the triple (Ω, F, P). Ω is the set of possible outcomes: for a card drawn it has 52 elements, for the toss of two coins it has four (HH, HT, TH, TT). An event is a subset of Ω: "a heart comes up" is the set of the 13 hearts. F is the family of events to which you want to be able to assign a probability, and P is the function that measures them, with values between 0 and 1.

Why not simply say that every subset is an event? When Ω is finite you can, and that is what you do with the deck. But Ω can be an interval such as [0, 1], and there measure theory proves that some subsets exist to which no coherent assignment of "length" can give a value (Vitali's result; the proof requires the axiom of choice and lies outside this article).

This is why F must be a σ-algebra: a family of subsets that contains Ω, contains the complement of each of its elements, and contains the union of every countable family of its elements. In plain words, if you can evaluate "A" you can evaluate "not A", and if you can evaluate each event in a list (even an infinite but countable one) you can evaluate "at least one of them happens".

A σ-algebra is also a way of saying how much information you have. With two coins, the family {∅, {HH, HT}, {TH, TT}, Ω} contains only the events that depend on the first coin: it can say "the first is heads", it cannot say "the second is heads". The program below takes the definition literally and checks it by exhaustion: on a set of n points it lists all the 2²ⁿ families of subsets (256 for n = 3, 65,536 for n = 4), keeps those that contain Ω and are closed under complement and union (on a finite set countable unions reduce to finite ones), and counts them. Then it computes the σ-algebra generated by a partition, that is the smallest one that contains the given sets, by closing them repeatedly.

# How many sigma-algebras exist on a finite set? Exhaustive check.
from itertools import combinations

def sigma_algebre(n):
    """All families F of subsets of {0..n-1} (as bit masks)
    that contain Omega and are closed under complement and union."""
    tutti = (1 << n) - 1                  # Omega as a mask
    trovate = []
    for famiglia in range(1 << (1 << n)):  # each family = a subset of the 2^n sets
        membri = [s for s in range(1 << n) if famiglia >> s & 1]
        if tutti not in membri:
            continue
        insieme = set(membri)
        if any((tutti ^ s) not in insieme for s in membri):
            continue
        if any((a | b) not in insieme for a, b in combinations(membri, 2)):
            continue
        trovate.append(membri)
    return trovate

for n in (3, 4):
    print(f"n={n}: {len(sigma_algebre(n))} sigma-algebras (out of {2**(2**n)} candidate families)")

def generata(n, generatori):
    """Closure of the generators (and of Omega) under complement and union."""
    tutti = (1 << n) - 1
    f = {tutti, 0} | set(generatori)
    cambiato = True
    while cambiato:
        cambiato = False
        for s in list(f):
            if (tutti ^ s) not in f:
                f.add(tutti ^ s); cambiato = True
        for a, b in combinations(list(f), 2):
            if (a | b) not in f:
                f.add(a | b); cambiato = True
    return f

def maschera(elementi):
    return sum(1 << (e - 1) for e in elementi)

C = [maschera({1, 2}), maschera({3, 4, 5}), maschera({6, 7, 8})]
print("sigma generated by {1,2},{3,4,5},{6,7,8} on 8 points:", len(generata(8, C)), "elements")
prima_moneta = [maschera({1, 2})]            # outcomes 1=HH, 2=HT, 3=TH, 4=TT: 'first heads'
print("sigma generated by 'first coin heads' on 4 outcomes:", len(generata(4, [0b0011])), "elements")

# Sigma-algebras on n points are as many as the partitions of n points (Bell numbers)
def bell(n):
    riga = [1]
    for _ in range(n - 1):
        nuova = [riga[-1]]
        for x in riga:
            nuova.append(nuova[-1] + x)
        riga = nuova
    return riga[-1]
print("Bell numbers B(1..6):", [bell(k) for k in range(1, 7)])

# Output:
# n=3: 5 sigma-algebras (out of 256 candidate families)
# n=4: 15 sigma-algebras (out of 65536 candidate families)
# sigma generated by {1,2},{3,4,5},{6,7,8} on 8 points: 8 elements
# sigma generated by 'first coin heads' on 4 outcomes: 4 elements
# Bell numbers B(1..6): [1, 2, 5, 15, 52, 203]

On three points there are 5 σ-algebras, on four there are 15, and the row of Bell numbers confirms the pattern: there are as many as the partitions of the set into blocks, because a finite σ-algebra is determined by its "atoms", the minimal events that cannot be split.

You can see it in the last two lines too: the partition {1,2}, {3,4,5}, {6,7,8} generates 8 = 2³ events, one for each way of choosing which blocks to unite, and the first coin alone generates 4. The count is exact only because Ω is finite; on an infinite set no enumeration of this kind exists, which is why there one works with the Borel σ-algebra, generated by the intervals, without ever writing it out in full.

The three axioms and what follows from them

A function P defined on the events of F is a probability measure if it satisfies three conditions, and nothing else:

1. Non-negativity: P(A) ≥ 0 for every event A.
2. Normalization: P(Ω) = 1.
3. σ-additivity: if A₁, A₂, … are pairwise disjoint events, P(A₁ ∪ A₂ ∪ …) = P(A₁) + P(A₂) + …

The first two fix the scale (no negative probabilities, the certain event is worth 1); the third is the only rule of calculation, and says that the probability of "one among many cases that exclude each other" is the sum. The order here is not interchangeable: the third uses the word "probability" and makes sense only after the first two have established what a probability is. Everything else is derived.

For instance P(∅) = 0: applying the third to the sequence of infinitely many empty sets (disjoint because empty) gives P(∅) = P(∅) + P(∅) + …, and a non-negative number that equals the sum of infinitely many copies of itself can only be zero.

From P(A) + P(A^c) = P(Ω) = 1 it follows that P(A^c) = 1 − P(A), where A^c is the complement; if A ⊂ B it follows that P(A) ≤ P(B), because B splits into A and B − A, which are disjoint.

For any two events the inclusion-exclusion formula holds: P(A ∪ B) = P(A) + P(B) − P(A ∩ B), because by adding P(A) and P(B) you count the common part twice.

An example with real numbers: you draw a card from a deck of 52. Let A be "hearts" (13 cards) and B "face card", that is J, Q or K (12 cards). Three cards belong to both sets: the J, Q and K of hearts. With the uniform measure, in which each card is worth 1/52, P(A ∪ B) = 13/52 + 12/52 − 3/52 = 22/52.

The program builds the deck, counts the cards of the union directly with exact fractions, compares with the formula, and then repeats the draw a million times with a fixed-seed generator. The acceptance criterion is declared: the observed frequency must fall within three standard errors, 3√(p(1 − p)/N), of the exact value.

from fractions import Fraction
from itertools import product
import numpy as np

semi = range(4)            # 0 = hearts
valori = range(13)         # 0=ace, 1..9 = 2..10, 10=J, 11=Q, 12=K
mazzo = list(product(semi, valori))
cuori = {c for c in mazzo if c[0] == 0}
figure = {c for c in mazzo if c[1] >= 10}      # J, Q, K: 12 cards in all

P = lambda E: Fraction(len(E), len(mazzo))     # uniform measure
esatto = P(cuori | figure)
incl_escl = P(cuori) + P(figure) - P(cuori & figure)
print("|hearts|, |face cards|, |hearts and face cards| =", len(cuori), len(figure), len(cuori & figure))
print("P(A or B) exact       =", esatto, "=", float(esatto))
print("incl.-excl.            =", incl_escl, "equal:", esatto == incl_escl)

N = 10**6
rng = np.random.default_rng(42)
pesca = rng.integers(0, 52, size=N)
seme, valore = pesca // 13, pesca % 13
freq = np.mean((seme == 0) | (valore >= 10))
p = float(esatto)
tol = 3 * (p * (1 - p) / N) ** 0.5
print(f"Monte Carlo N={N}: {freq:.5f}; expected error (3 sigma) {tol:.5f}; deviation {abs(freq-p):.5f}")
assert abs(freq - p) < tol

# Output:
# |hearts|, |face cards|, |hearts and face cards| = 13 12 3
# P(A or B) exact       = 11/26 = 0.4230769230769231
# incl.-excl.            = 11/26 equal: True
# Monte Carlo N=1000000: 0.42294; expected error (3 sigma) 0.00148; deviation 0.00014

The exact value is 11/26, about 0.4231, and the formula gives exactly the same count. The simulation deviates by 0.00014 against a tolerance of 0.00148: in the typical case the deviation is about one standard error, so a value this small is normal and should not be read as stronger proof. What the simulation cannot tell you is why inclusion-exclusion is true: the axioms tell you.

With n events the complete formula has 2ⁿ − 1 intersections to compute, and already with n = 20 that is more than a million; for this reason in practice one often uses only Boole's inequality, P(A₁ ∪ … ∪ A_n) ≤ P(A₁) + … + P(A_n), which follows from the same additivity and is always true even when the events overlap.

Conditional probability and Bayes' formula

If you know that B has happened, the space of possible outcomes shrinks to B, and the probabilities must be rescaled so that they still sum to 1. This is conditional probability: P(A | B) = P(A ∩ B) / P(B), defined only if P(B) > 0. Rewriting it as P(A ∩ B) = P(A | B) P(B) gives you the product rule, which iterated gives the chain rule over several events.

If the events B₁, …, B_m form a partition of Ω (they exclude each other and cover everything), the law of total probability decomposes P(A) = Σ_j P(A | B_j) P(B_j).

Writing P(A ∩ B_i) in two ways and dividing by P(A) you arrive at Bayes' formula: P(B_i | A) = P(A | B_i) P(B_i) / Σ_j P(A | B_j) P(B_j). The two directions, P(A | B) and P(B | A), are different numbers, and swapping them is the most common mistake in the whole subject.

Worked exercise. A test for a disease has sensitivity 99% (if you are ill, it comes out positive in 99% of cases) and specificity 95% (if you are healthy, it comes out negative in 95% of cases, that is it gives a false positive in 5%). The disease has prevalence 1%, that is it affects one person in a hundred.

Your test is positive: what is the probability that you are ill? To answer, you compute the positive predictive value with Bayes' formula, rewrite it in natural frequencies (out of 10,000 people, how many are in the two groups) because that is how the mind grasps it, let the prevalence vary, and simulate a population of a million people, counting how many of the positives are really ill.

from fractions import Fraction as F
import numpy as np

def ppv(sens, spec, prev):
    return sens * prev / (sens * prev + (1 - spec) * (1 - prev))

sens, spec, prev = F(99, 100), F(95, 100), F(1, 100)
print("PPV exact =", ppv(sens, spec, prev), "=", round(float(ppv(sens, spec, prev)), 4))

# The same figures as natural frequencies, out of 10000 people
malati = 100; veri_pos = 99; sani = 9900; falsi_pos = 495
print("natural frequencies:", veri_pos, "true positives out of", veri_pos + falsi_pos, "positives =", F(veri_pos, veri_pos + falsi_pos))

for p in (F(1, 1000), F(1, 100), F(1, 10), F(1, 2)):
    print(f"prevalence {float(p):5.3f}: PPV = {float(ppv(sens, spec, p)):.4f}")

N = 10**6
rng = np.random.default_rng(7)
malato = rng.random(N) < 0.01
u = rng.random(N)
positivo = np.where(malato, u < 0.99, u < 0.05)     # sick: 99% positive; healthy: 5% positive
n_pos = int(positivo.sum())
freq = malato[positivo].mean()
p = float(ppv(sens, spec, prev))
tol = 3 * (p * (1 - p) / n_pos) ** 0.5
print(f"Monte Carlo: {n_pos} positives, share of sick {freq:.4f}; expected {p:.4f} +/- {tol:.4f}")
assert abs(freq - p) < tol

# Output:
# PPV exact = 1/6 = 0.1667
# natural frequencies: 99 true positives out of 594 positives = 1/6
# prevalence 0.001: PPV = 0.0194
# prevalence 0.010: PPV = 0.1667
# prevalence 0.100: PPV = 0.6875
# prevalence 0.500: PPV = 0.9519
# Monte Carlo: 59337 positives, share of sick 0.1660; expected 0.1667 +/- 0.0046

The answer is 1/6, about 16.7%, and not the 99% that the sensitivity figure suggests. The natural frequencies show why: out of 10,000 people the ill are 100 and 99 of them test positive, but the healthy are 9,900 and 5% of them, that is 495, test positive by mistake; false positives outnumber true ones five to one. The table of prevalences is the most instructive part: with a rare disease (one in a thousand) the same test, of identical quality, gives a predictive value of 1.9%, while with prevalence at 50% it rises to 95%.

Neither characteristic of the test is enough on its own; the population you apply it to matters as well. The simulation finds 59,337 positives, of which 16.60% are ill, against an expected 16.67% with a tolerance of 0.46 percentage points: the tolerance is wider than the one for the cards because the positives are fewer than the whole sample. This is arithmetic, not medical advice: no test result should be interpreted without a healthcare professional.

Independence: pairwise is not enough

Two events are independent if P(A ∩ B) = P(A) P(B); if P(B) > 0 this is the same as saying that P(A | B) = P(A), that is knowing B does not change the probability of A. It is an algebraic condition on the chosen P, not an impression about "real life": two events can be independent for one measure and dependent for another.

For more events, checking the pairs is not enough: you need mutual independence, that is the product rule for every subfamily of at least two events, which makes 2ⁿ − n − 1 equalities. With three events there are four: the three pairs and the triple. This definition is 4.2 in chapter 4 of Grinstead and Snell, who give the same kind of counterexample as an exercise.

The counterexample can be built with two fair coins. Let A be "the first is heads", B "the second is heads", C "exactly one head comes up". Each has probability 1/2, and every pair has probability 1/4 = 1/2 · 1/2: the pairs are independent. But if two heads come up, "exactly one head" cannot, so A ∩ B ∩ C is empty. The program computes it with fractions on the set of the four outcomes, and also checks a widespread error: two disjoint events are not independent, unless one of them has probability zero.

from fractions import Fraction as F
from itertools import product, combinations

omega = list(product((0, 1), repeat=2))              # (first, second); 1 = heads
P = lambda E: F(len(E), len(omega))
A = {w for w in omega if w[0] == 1}                  # first heads
B = {w for w in omega if w[1] == 1}                  # second heads
C = {w for w in omega if w[0] != w[1]}               # exactly one heads

for (na, X), (nb, Y) in combinations([("A", A), ("B", B), ("C", C)], 2):
    print(f"P({na} and {nb}) = {P(X & Y)}  product = {P(X) * P(Y)}  independent: {P(X & Y) == P(X) * P(Y)}")
print("P(A and B and C) =", P(A & B & C), " product =", P(A) * P(B) * P(C))

# Disjoint does not mean independent: A = 'first heads', D = 'first tails'
D = {w for w in omega if w[0] == 0}
print("A and D disjoint:", A & D == set(), "| P(A and D) =", P(A & D), "versus P(A)P(D) =", P(A) * P(D))

# Output:
# P(A and B) = 1/4  product = 1/4  independent: True
# P(A and C) = 1/4  product = 1/4  independent: True
# P(B and C) = 1/4  product = 1/4  independent: True
# P(A and B and C) = 0  product = 1/8
# A and D disjoint: True | P(A and D) = 0 versus P(A)P(D) = 1/4

The three pairs are independent, the triple is not: P(A ∩ B ∩ C) = 0, while the product of the three probabilities would be 1/8. Practical consequence: knowing that A and B have happened tells you everything about C (it has not happened), even though neither piece of news alone tells you anything about it.

The last line corrects the habit of using "disjoint" and "independent" as synonyms: "first heads" and "first tails" exclude each other, so their intersection has probability 0, while the product is 1/4. If two events cannot happen together, knowing that one has happened gives you maximum information about the other, the opposite of independence.

The Borel-Cantelli lemmas: does it happen infinitely often?

Given a sequence of events A₁, A₂, …, the event "A_n happens infinitely often" is written as the intersection over n of the unions from n onwards: no tail of the sequence is free of events. The first Borel-Cantelli lemma says that if the sum of the probabilities converges, Σ_n P(A_n) < ∞, then almost surely only finitely many events happen.

The proof is one line: the probability that at least one event happens from n onwards is at most the tail of the series, Σ_{k ≥ n} P(A_k), which tends to 0. The second lemma requires independence: if the events are independent and the sum diverges, almost surely infinitely many of them happen. Without independence the second lemma is false: just take A_n = A for every n, with P(A) = 1/2.

A simulation cannot show "infinitely" because it runs for a finite time, and it is worth saying so before reading its output. What you can observe are the quantities the reasoning uses, which can be computed exactly. The program generates 20,000 independent sequences, with P(A_n) = 1/n² (convergent series) or 1/n (divergent series), and counts how many events happen for n from 21 to 2000, the "tail" after the first twenty.

It compares the mean with the exact sum of the probabilities, and the fraction of sequences with no event at all in the tail with the exact product of the complements, which by independence is the probability of seeing none.

import numpy as np
from fractions import Fraction

R, N, SOGLIA = 20000, 2000, 20     # replicates, horizon, "tail" = events with n > SOGLIA
n = np.arange(1, N + 1)
rng = np.random.default_rng(2024)
risultati = {}
for nome, p in (("1/n^2", 1 / n**2), ("1/n", 1 / n)):
    nella_coda = np.zeros(R, dtype=int)
    for inizio in range(0, R, 500):                  # in blocks, to save memory
        occ = rng.random((500, N)) < p               # occ[r, k] = A_{k+1} occurs in replicate r
        nella_coda[inizio:inizio + 500] = occ[:, SOGLIA:].sum(axis=1)
    risultati[nome] = nella_coda

coda = n[SOGLIA:]
for nome, p in (("1/n^2", 1 / coda**2), ("1/n", 1 / coda)):
    atteso = p.sum()                                  # expected number of events in the tail
    prob_nessuno = np.prod(1 - p)                     # independence: product
    r = risultati[nome]
    print(f"P(A_n)={nome:6}: events with n in ({SOGLIA},{N}]: mean {r.mean():.3f} (exact {atteso:.3f}); "
          f"P(none) simulated {np.mean(r == 0):.4f} (exact {prob_nessuno:.4f})")

    sd = (p * (1 - p)).sum() ** 0.5                   # standard deviation of the count
    assert abs(r.mean() - atteso) < 4 * sd / R ** 0.5, nome
print("means within 4 standard errors of the exact value")

# Output:
# P(A_n)=1/n^2 : events with n in (20,2000]: mean 0.049 (exact 0.048); P(none) simulated 0.9522 (exact 0.9529)
# P(A_n)=1/n   : events with n in (20,2000]: mean 4.603 (exact 4.581); P(none) simulated 0.0092 (exact 0.0100)
# means within 4 standard errors of the exact value

In the convergent case the tail contains 0.049 events on average and 95% of the sequences have none: after the first twenty, almost nothing happens any more, which is what the first lemma predicts (the sum of the tail, 0.048, is also an upper bound on the probability of seeing at least one).

In the divergent case the tail contains 4.6 on average, and the probability of seeing none up to 2000 is 1%: as the horizon is lengthened the mean count grows like the logarithm of N without stopping, which is the behaviour the second lemma carries to the limit. The two means are within four standard errors of the exact value, as the program checks with an assertion. Remember that these are coherence checks on a finite horizon, not a proof of the lemmas: the proof is the previous paragraph.

What is left out

Three things go beyond what has been proved here, and it is worth knowing where they begin. Conditioning on an event of probability zero, such as "X = 3" for a continuous variable, is not covered by the definition above and requires regular conditional probability. The extension of a probability from a simple family of events to the whole generated σ-algebra (Carathéodory's theorem) is what lets you construct Lebesgue measure and infinite stochastic processes.

Conditional expectation with respect to a σ-algebra generalizes total probability. For these, and for the complete proofs, the classic reference is a textbook on measure-theoretic probability, such as those by Billingsley or Kallenberg; I did not open them for this article, and I name them only to point to where to go next.

Frequently asked questions

Why does Kolmogorov use measure theory?

Because it lets you treat discrete and continuous cases together with the same rules, and construct probabilities on infinite spaces, such as the trajectories of a process. The axioms are those of a positive measure of total mass 1, and this connects probability to an apparatus already developed for lengths and volumes.

Is an event of probability zero impossible?

No. If you pick a number at random in [0, 1] with uniform distribution, every single value has probability zero, and yet one of them will come up. Zero means "almost surely does not happen", not "cannot happen"; the impossible event is only the empty set.

Why does a positive test not mean I am 99% likely to be ill?

Because 99% is the probability of a positive result if you are ill, not the probability of being ill if the result is positive. The two probabilities are linked by Bayes' formula through the prevalence: with prevalence at 1% and specificity at 95% the predictive value is 1/6.

What you can do now

You can define a probability space and check whether a family of events is a σ-algebra, derive from the three rules the properties you will use every day, and compute a conditional probability both with the formula and with natural frequencies. You can tell pairwise independence from mutual independence, and independence from incompatibility, and you can read the Borel-Cantelli lemmas as a statement about tails of series.

The next step is to give a number to every outcome: discrete random variables, which use the law of total probability to break a count down into independent trials. With a variable that can take any value in an interval the way probabilities are assigned changes, and they become areas: you see it in continuous variables. Further on, multivariate distributions take up conditional probability for several quantities together, and the central limit theorem uses the independence introduced here. The index of the whole series is on the Probability and Statistics page.

The same series also includes Exploratory Data Analysis (EDA) with Python.