18 - Probability: Reasoning About Uncertainty

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 3

Goal. Give the word "probability" a definition that can be computed, derive its rules from three axioms, and use them for two tasks that machine learning asks for constantly: updating a belief when new information arrives (Bayes' theorem) and understanding when two events are really independent. At the end you will be able to explain why a test "90% accurate" that comes out positive can mean only a 15% probability of being ill.

Myth to dispel. "After five heads in a row, tails is more likely: the coin must make up for it." The coin has no memory. This belief (the gambler's fallacy) is the most widespread and the most expensive; we will take it apart with a simulation.

Introduction

A machine learning model almost never gives a certain answer: it gives a probability ("this email is spam with probability 0.93"). To read it, to trust it and to understand when it errs, you need the language of probability. We build it from the ground up, with analogies first and the formulas derived in full. We take up sets (lesson 1) and counts (lesson 2): they are the bricks.

The path goes from the random experiment to Kolmogorov's axioms, then to conditional probability, independence and Bayes' theorem, and ends with a toy classifier, the Monty Hall problem and the most common reasoning errors. For the counts I refer to Combinatorics, for the language of events to Sets and logic.

Experiment, sample space, events

The analogy: the roulette with numbered slots

A random experiment is a process of which you know all the possible results but not which will come out: I roll a die, I draw a card, a customer opens an email or not. Each single result is an outcome. The set of all outcomes is the sample space, denoted Ω ("omega"). An event is a subset of Ω: a set of outcomes we care about. It is exactly the language of the sets lesson: "an even number comes up" is the set {2, 4, 6} ⊆ Ω = {1, 2, 3, 4, 5, 6}; "and" becomes intersection, "or" union, "not" complement.

Probability with equally likely outcomes

If the |Ω| outcomes are equally possible (a fair die), the probability of an event A is the proportion of favourable outcomes:

P(A) = |A| / |Ω|.

Here combinatorics enters: to compute |A| and |Ω| we count.

from itertools import product
from fractions import Fraction

omega = set(product(range(1, 7), repeat=2))     # sample space of two dice: 36 pairs
def P(evento):
    return Fraction(len(evento), len(omega))

somma_7 = {(a, b) for (a, b) in omega if a + b == 7}
doppio = {(a, b) for (a, b) in omega if a == b}
print("|Ω| =", len(omega))
print("P(sum = 7) =", P(somma_7))
print("P(double)  =", P(doppio))
print("P(sum = 7 or double) =", P(somma_7 | doppio), "| P(sum = 7 and double) =", P(somma_7 & doppio))

Output obtained:

|Ω| = 36
P(sum = 7) = 1/6
P(double)  = 1/6
P(sum = 7 or double) = 1/3 | P(sum = 7 and double) = 0

What the code does. product(range(1, 7), repeat=2) generates the 36 pairs (die1, die2); the events are sets of pairs built with a comprehension; the function P divides the cardinalities. We use Fraction to obtain exact fractions (6/36 is printed as 1/6) instead of approximate decimals. The last line shows the union ("or") and the intersection ("and") between events: between sum = 7 and doubles there is no common outcome (the sum 7 is odd, a double has an even sum), so the intersection is empty.

Kolmogorov's axioms

If the outcomes are not equally likely (a loaded coin) the rule |A|/|Ω| no longer holds. Andrey Kolmogorov (1933) fixed three minimal rules that any probability must respect, and from which everything else is derived:

1. Non-negativity: P(A) ≥ 0 for every event.
2. Normalisation: P(Ω) = 1 (something must happen).
3. Additivity: if A and B are disjoint (A ∩ B = ∅), then P(A ∪ B) = P(A) + P(B).

These are very few rules, and breaking just one gives a "probability" that does not behave like one: with a negative value one could buy bets at a negative price, with P(Ω) ≠ 1 the percentages would not add up to a hundred, without additivity two exclusive events could be worth more than their sum. For infinite spaces axiom 3 extends to countably infinitely many disjoint events; in our finite examples the two-event version is enough.

What is derived, step by step

The complement. A and Aᶜ are disjoint and their union is Ω. By axiom 3: P(A) + P(Aᶜ) = P(Ω) = 1 (axiom 2). Hence P(Aᶜ) = 1 − P(A). It is the formula we used for the birthdays.

The impossible event. Ω and ∅ are disjoint and Ω ∪ ∅ = Ω: P(Ω) + P(∅) = P(Ω), hence P(∅) = 0.

Monotonicity and upper bound. If A ⊆ B, write B = A ∪ (B ∖ A) (disjoint): P(B) = P(A) + P(B∖A) ≥ P(A) by axiom 1. In particular P(A) ≤ P(Ω) = 1: no probability exceeds 1.

The union rule. For any two events (even non-disjoint): A ∪ B breaks into the three disjoint parts A∖B, A∩B, B∖A. With P(A) = P(A∖B) + P(A∩B) and P(B) = P(B∖A) + P(A∩B) and summing the three parts one obtains

P(A ∪ B) = P(A) + P(B) − P(A ∩ B)

(the inclusion-exclusion of lesson 1, now for probabilities). The same proof, with the words "number of elements" in place of "probability".

# check of the consequences on the two dice
A, B = somma_7, {(a, b) for (a, b) in omega if a == 3 or b == 3}      # B: at least one 3
print("P(A) =", P(A), "| P(B) =", P(B), "| P(A∩B) =", P(A & B))
print("union rule:", P(A | B) == P(A) + P(B) - P(A & B))
print("complement:", P(omega - A) == 1 - P(A))
print("monotonicity (A∩B ⊆ A):", P(A & B) <= P(A))

Output obtained:

P(A) = 1/6 | P(B) = 11/36 | P(A∩B) = 1/18
union rule: True
complement: True
monotonicity (A∩B ⊆ A): True

What the code does. It takes two overlapping events (sum 7; at least one 3) and verifies with exact numbers three of the consequences just proved. If one of the three came out False, we would have made a mistake in a step: the computer acts as a check, not as a proof.

Probability as frequency: the simulation

What does "the probability of heads is 0.5" mean? An operational reading: if you repeat the experiment many times, the fraction of heads approaches 0.5. It is the law of large numbers, which we will prove in elementary form in the lesson on random variables; here we observe it.

import numpy as np

rng = np.random.default_rng(42)
lanci = rng.integers(0, 2, size=1_000_000)       # 0 = tails, 1 = heads

for n in (10, 100, 1_000, 10_000, 100_000, 1_000_000):
    freq = lanci[:n].mean()
    print(f"n = {n:>9,}: frequency of heads = {freq:.4f}   gap from 0.5 = {abs(freq - 0.5):.4f}")

Output obtained:

n =        10: frequency of heads = 0.4000   gap from 0.5 = 0.1000
n =       100: frequency of heads = 0.5200   gap from 0.5 = 0.0200
n =     1,000: frequency of heads = 0.5090   gap from 0.5 = 0.0090
n =    10,000: frequency of heads = 0.4939   gap from 0.5 = 0.0061
n =   100,000: frequency of heads = 0.5009   gap from 0.5 = 0.0009
n = 1,000,000: frequency of heads = 0.4998   gap from 0.5 = 0.0002

What the code does. It generates a million simulated tosses (zeros and ones), then looks at the first n for increasing n and prints the fraction of ones. The mean of an array of 0s and 1s is exactly the relative frequency of the ones. The gap from 0.5 tends to shrink, even if not regularly: at 10 tosses it can be large, at a million it is hundreds of times smaller. One important thing: the absolute gap in the number of heads does not vanish (it can grow), it is the fraction that vanishes. This is why the gambler's fallacy is a fallacy: the future does not "repay" the past, it dilutes it.

Conditional probability

The analogy: narrowing the room

You are in a room with 100 people: 40 wear glasses, 20 are students, 10 are students with glasses. What is the probability that a person chosen at random wears glasses? 40/100. Now they tell you: "the person chosen is a student". The room has narrowed to the 20 students: in that room, 10 have glasses. The probability becomes 10/20 = 0.5. You have changed the reference universe.

Definition. The probability of A knowing that B has happened is

P(A | B) = P(A ∩ B) / P(B), with P(B) > 0.

(The numerator is the share of B also occupied by A; the denominator is the size of the new room.) Rewritten, it is the product rule: P(A ∩ B) = P(A | B) · P(B): to have both B and A, first B is needed, then A given B.

def P_cond(A, B):
    """P(A | B) from the definition: P(A ∩ B) / P(B)."""
    return P(A & B) / P(B)

somma_8 = {(a, b) for (a, b) in omega if a + b == 8}
primo_pari = {(a, b) for (a, b) in omega if a % 2 == 0}
print("P(sum = 8)                    =", P(somma_8))
print("P(sum = 8 | first die even)   =", P_cond(somma_8, primo_pari))
print("check by counting in the narrowed room:",
      Fraction(len(somma_8 & primo_pari), len(primo_pari)))

Output obtained:

P(sum = 8)                    = 5/36
P(sum = 8 | first die even)   = 1/6
check by counting in the narrowed room: 1/6

What the code does. P_cond is the definition translated literally. The last line counts directly inside the narrowed "room" (primo_pari) and obtains the same number: it is the analogy turned into code. Knowing that the first die is even changes the probability of the sum 8 (from 5/36 to 1/6, because some pairs with sum 8 leave the room).

Mind the direction. P(A | B) and P(B | A) are in general different. "The probability that a culprit has a certain blood group" is not "the probability that whoever has that blood group is guilty". Confusing them is the prosecutor's fallacy and has sent wrong conclusions to court.

Independence

Two events are independent if knowing one does not change the probability of the other: P(A | B) = P(A). With the definition of conditional it is equivalent to

P(A ∩ B) = P(A) · P(B).

Two coin tosses are independent: the coin does not remember. But beware of a very common confusion: disjoint does not mean independent. Two disjoint events with positive probabilities are in fact maximally dependent: if one happens, the other cannot happen (P(A | B) = 0 ≠ P(A)).

primo_6 = {(a, b) for (a, b) in omega if a == 6}
somma_pari = {(a, b) for (a, b) in omega if (a + b) % 2 == 0}
secondo_6 = {(a, b) for (a, b) in omega if b == 6}

def indipendenti(A, B):
    return P(A & B) == P(A) * P(B)

print("first = 6  vs second = 6   independent:", indipendenti(primo_6, secondo_6))
print("first = 6  vs even sum     independent:", indipendenti(primo_6, somma_pari))
print("sum = 7  vs sum = 8 (disjoint) independent:", indipendenti(somma_7, somma_8))
print("sum = 7  vs first = 6          independent:", indipendenti(somma_7, primo_6))

Output obtained:

first = 6  vs second = 6   independent: True
first = 6  vs even sum     independent: True
sum = 7  vs sum = 8 (disjoint) independent: False
sum = 7  vs first = 6          independent: True

What the code does. It verifies the definition P(A∩B) = P(A)P(B) with exact fractions, hence without rounding errors (floats lesson). The third case shows the trap: sum = 7 and sum = 8 are disjoint but not independent. The fourth is an instructive case: knowing that the first die is 6 does not change the probability of obtaining sum 7 (whatever the first die, there is exactly one value of the second that gives 7), and indeed they are independent even if they "seem" linked.

In machine learning. The assumption of independence between the features (given the outcome) is the simplification of the Naive Bayes classifier; the assumption of independence between the examples (independent samples) is the basis of almost all statistics. When it is false (time series, several rows of the same customer) the results become too optimistic.

Total probability and Bayes' theorem

The law of total probability

If the events A₁, A₂, …, A_k form a partition of Ω (disjoint, with union Ω), then every event B breaks into the disjoint pieces B ∩ Aᵢ, and by the product rule

P(B) = Σᵢ P(B | Aᵢ) · P(Aᵢ).

It is a weighted mean: the probability of B is the mean of its probabilities "in each case", weighted by how likely each case is.

Bayes' theorem

From the product rule written in two ways, P(A ∩ B) = P(B | A) P(A) = P(A | B) P(B), dividing by P(B):

P(A | B) = P(B | A) · P(A) / P(B)

It is the theorem of Thomas Bayes (published posthumously in 1763) in the form made famous by Laplace. It reads like this: the new belief in A after seeing B (posterior) is the starting belief (prior, P(A)) corrected by the ratio between how much B is expected if A is true and how much B is expected in general (the likelihood P(B | A) normalised by P(B), computed with total probability).

The example that fools almost everyone: the diagnostic test

Invented and declared parameters: a disease affects 1% of the population; a test recognises it in 90% of the ill (sensitivity); it gives a false alarm in 5% of the healthy (false positive rate). A patient tests positive. How likely is it that he is ill?

Intuition answers "about 90%". Let us compute.

prevalenza = 0.01           # P(ill)
sensibilita = 0.90          # P(positive | ill)
falsi_positivi = 0.05       # P(positive | healthy)

p_pos = sensibilita * prevalenza + falsi_positivi * (1 - prevalenza)      # total probability
p_malato_dato_pos = sensibilita * prevalenza / p_pos                      # Bayes
print(f"P(positive)          = {p_pos:.4f}")
print(f"P(ill | positive)    = {p_malato_dato_pos:.4f}")

# the same thing counting people, out of 10,000
n = 10_000
malati = n * prevalenza
veri_pos = malati * sensibilita
falsi_pos = (n - malati) * falsi_positivi
print(f"out of {n:,}: ill = {malati:.0f}, true positives = {veri_pos:.0f}, false positives = {falsi_pos:.0f}")
print(f"among the positives, the ill are {veri_pos:.0f} out of {veri_pos + falsi_pos:.0f} = {veri_pos / (veri_pos + falsi_pos):.4f}")

Output obtained:

P(positive)          = 0.0585
P(ill | positive)    = 0.1538
out of 10,000: ill = 100, true positives = 90, false positives = 495
among the positives, the ill are 90 out of 585 = 0.1538

What the code does. The first part applies the two formulas (total probability in the denominator, then Bayes). The second repeats the calculation counting people out of 10,000: it is the method of natural frequencies, which the studies of Gigerenzer and Hoffrage (1995) showed makes Bayesian reasoning much easier to follow. Result: about 15.4%.

Why so low? Because the disease is rare: out of 10,000 people the ill are 100, the healthy 9,900. The 5% of false alarms on the healthy produces 495 wrong positives, many more than the 90 true positives. The base rate (the prior) weighs as much as the precision of the test. Neglecting it is the base-rate fallacy.

Updating again. After a second positive test (assuming, declared assumption, that the two tests are independent given the health state), the posterior of the first becomes the prior of the second:

def aggiorna(prior, sens=0.90, fp=0.05):
    pos = sens * prior + fp * (1 - prior)
    return sens * prior / pos

p = prevalenza
for i in range(1, 4):
    p = aggiorna(p)
    print(f"after {i} positive test{'s' if i > 1 else ''}: P(ill) = {p:.4f}")

Output obtained:

after 1 positive test: P(ill) = 0.1538
after 2 positive tests: P(ill) = 0.7660
after 3 positive tests: P(ill) = 0.9833

What the code does. It applies aggiorna repeatedly, reusing the output as input: the posterior becomes the next prior. With three positive tests the probability goes from 1% to about 98%. It is the idea at the basis of Bayesian filters and incremental learning, and shows that evidence accumulates.

The link with the metrics (lesson 3 of module 00). The probability of being ill given a positive is the precision of a classifier, and the sensitivity is the recall. The same test, applied to a population with a different prevalence, has a different precision: the quality of a model is not expressed with a single number, but with the number and the population in which it is used.

A toy Bayesian classifier

We apply Bayes to a text problem, with an invented mini-corpus: 8 messages classified as spam or normal, and the words they contain.

corpus = [
    ({"free", "win"}, "spam"), ({"free", "offer"}, "spam"), ({"win", "offer", "click"}, "spam"),
    ({"click", "free"}, "spam"), ({"meeting", "tomorrow"}, "normal"), ({"tomorrow", "offer"}, "normal"),
    ({"meeting", "project"}, "normal"), ({"project", "free"}, "normal"),
]
vocabolario = {p for parole, _ in corpus for p in parole}
classi = ("spam", "normal")

prior = {c: sum(1 for _, y in corpus if y == c) / len(corpus) for c in classi}

def p_parola_dato_classe(parola, classe, alpha=1.0):
    """P(word in the message | class) with Laplace smoothing."""
    messaggi = [parole for parole, y in corpus if y == classe]
    con = sum(1 for parole in messaggi if parola in parole)
    return (con + alpha) / (len(messaggi) + 2 * alpha)

def posteriore(parole_osservate):
    punteggio = {}
    for c in classi:
        p = prior[c]
        for w in vocabolario:                         # naive: words independent given the class
            pw = p_parola_dato_classe(w, c)
            p *= pw if w in parole_osservate else (1 - pw)
        punteggio[c] = p
    tot = sum(punteggio.values())
    return {c: v / tot for c, v in punteggio.items()}

print("prior:", prior)
print("{'free','click'}      ->", {c: round(v, 3) for c, v in posteriore({'free', 'click'}).items()})
print("{'meeting','tomorrow'} ->", {c: round(v, 3) for c, v in posteriore({'meeting', 'tomorrow'}).items()})

Output obtained:

prior: {'spam': 0.5, 'normal': 0.5}
{'free','click'}      -> {'spam': 0.926, 'normal': 0.074}
{'meeting','tomorrow'} -> {'spam': 0.024, 'normal': 0.976}

What the code does, line by line. prior is the fraction of messages per class. p_parola_dato_classe estimates P(word present | class) as the fraction of messages of that class that contain the word, with Laplace smoothing (+alpha): without it a word never seen in a class would give probability 0 and annul the whole product (a single zero cancels every other piece of evidence). posteriore applies Bayes: it multiplies the prior by the likelihood of each word (present or absent), assuming conditional independence of the words given the class (the "naive" assumption), and then normalises by dividing by the sum (the denominator P(B) of total probability). The messages with "free" and "click" turn out much more probably spam; those with "meeting" and "tomorrow" much more probably normal.

Trade-off. The independence assumption is almost always false (the words of a text are correlated), and yet Naive Bayes works surprisingly well because to classify what is needed is the order of the posteriors, not their exact value. But the probabilities it returns are often poorly calibrated: they must not be read as true percentages.

Monty Hall

Three doors: behind one there is a car, behind the other two a goat. You choose door 1. The host, who knows where the car is, opens one of the other two doors that hides a goat and asks whether you want to switch. Is it worth switching?

Intuition says "50 and 50". The analysis by cases (with the initial choice fixed, door 1) says otherwise:

• The car is behind 1 (probability 1/3): by switching you lose.
• The car is behind 2 (1/3): the host must open 3, by switching you win.
• The car is behind 3 (1/3): the host must open 2, by switching you win.

Switching wins in 2 cases out of 3. The key is that the opening of the door is not random: the host cannot open the one with the car, so his gesture carries information. (The problem was discussed by Selvin in 1975; the version with explicit assumptions on the host's behaviour is essential.)

rng = np.random.default_rng(1)

def gioca(cambia, prove=100_000):
    auto = rng.integers(0, 3, prove)
    scelta = np.zeros(prove, dtype=int)                 # you always choose door 0
    # the host opens a door different from `scelta` and without the car
    apre = np.empty(prove, dtype=int)
    for i in range(prove):
        candidate = [d for d in range(3) if d != scelta[i] and d != auto[i]]
        apre[i] = candidate[rng.integers(len(candidate))]
    if cambia:
        scelta_finale = 3 - scelta - apre               # the only door left (0+1+2 = 3)
    else:
        scelta_finale = scelta
    return float((scelta_finale == auto).mean())

print("staying   :", round(gioca(False), 4))
print("switching :", round(gioca(True), 4))

Output obtained:

staying   : 0.334
switching : 0.6649

What the code does. It simulates 100,000 games. For each it draws the door with the car, fixes the initial choice on door 0 (without loss of generality) and has the host open one of the doors different from yours and without the car (if he has two available, he picks one at random). If you switch, the final door is the only one left: since the doors are 0, 1, 2 and the sum is 3, 3 − scelta − apre identifies it. The result confirms 1/3 against 2/3. A simulation does not replace the analysis, but gives you the courage to trust a counter-intuitive result.

Reasoning errors (symptom, cause, remedy)

• Gambler's fallacy. "After 5 heads, tails is more likely." We show it with the simulation below: the frequency of heads after five heads in a row stays 0.5.
• Prosecutor's fallacy. Confusing P(A | B) with P(B | A). Remedy: write explicitly which is the data and which the hypothesis.
• Base-rate fallacy. Ignoring the prior. Remedy: natural frequencies (tree with 10,000 people).
• Disjoint = independent. False: disjoint with positive probability are dependent.
• Conjunction fallacy. Judging "A and B" more probable than "A". Impossible, because A ∩ B ⊆ A and so P(A ∩ B) ≤ P(A) (monotonicity): documented by Tversky and Kahneman (1983).
• Multiplying non-independent probabilities. Remedy: use the product rule with the conditionals.
• Selected samples. If you condition on a chosen subset (only those who answer the survey) you are computing a conditional different from the one you wanted.

Six mistakes out of seven are the same distraction: answering a question different from the one asked. Writing out in full "probability of what, knowing what" before computing forces you to see which datum is the hypothesis and which the observation, and is enough on its own to avoid the prosecutor's, base-rate and selected-sample fallacies. The simulation that follows addresses instead the first one, the gambler's.

lanci = rng.integers(0, 2, size=2_000_000)
# positions where the last 5 tosses are all heads (1)
cinque_teste = (lanci[:-5] + lanci[1:-4] + lanci[2:-3] + lanci[3:-2] + lanci[4:-1]) == 5
successivo = lanci[5:][cinque_teste]
print("cases with 5 heads in a row:", int(cinque_teste.sum()))
print("frequency of heads on the next toss:", round(float(successivo.mean()), 4))

Output obtained:

cases with 5 heads in a row: 63267
frequency of heads on the next toss: 0.5022

What the code does. lanci[:-5] + lanci[1:-4] + … sums five sliding windows of tosses: the sum is 5 only if all five are heads. lanci[5:][cinque_teste] selects the toss immediately after each streak. The mean of these is the frequency of heads after five heads: it is 0.5, as without the condition.

How to work in practice

Always declare the sample space and the event before computing. For small discrete problems use exact fractions (Fraction), for large ones simulate and compare. For Bayes start from the tree with the frequencies and only afterwards write the formula; verify the consequences of the axioms (sum equal to 1, complement, monotonicity) with an automatic test; and declare the independence assumptions by asking yourself whether they hold in your data.

The essential vocabulary: the sample space Ω is the set of possible outcomes and an event is a subset of it. Two events are disjoint if they have no common outcomes and independent if knowing one does not change the probability of the other. The conditional P(A|B) is the probability of A in the universe narrowed to B; prior and posterior are the belief before and after a datum, and the likelihood measures how expected the datum was if the hypothesis is true. For a classifier, the sensitivity is the recall and the precision is the share of true positives among the positives; the base rate is how common the phenomenon is in the population.

Mini-project

1. With Fraction, compute for two dice P(sum ≥ 10), P(sum ≥ 10 | at least one 6) and prove, with code, that the two events are not independent.
2. Take up the diagnostic test again and produce a table with the precision (P(ill | positive)) for prevalences of 0.1%, 1%, 10%, 30%, 50%, with sensitivity 90% and false positives 5%. Comment on the progression.
3. Extend the toy Naive Bayes by adding 4 messages of your own: does the posterior of the message {'offer'} change? Why?

The first point trains you to recognise dependence even when the two events "seem" unconnected, the second anticipates the worked challenge below, the third lets you observe, by measuring it yourself, how much Laplace smoothing weighs when the corpus is as small as this one (eight messages).

A worked challenge: the precision of the test as the prevalence varies

The mini-project asked for a table of the precision P(ill | positive) for different prevalences. Here it is, with Bayes' formula in a function and with the numbers side by side.

def precisione(prevalenza, sensibilita=0.90, falsi_positivi=0.05):
    positivi_veri = sensibilita * prevalenza
    positivi_falsi = falsi_positivi * (1 - prevalenza)
    return positivi_veri / (positivi_veri + positivi_falsi)

print(f"{'prevalence':>11} {'precision':>11} {'true ill per 10,000 positives':>38}")
for prev in (0.001, 0.01, 0.10, 0.30, 0.50, 0.80):
    p = precisione(prev)
    print(f"{prev:>11.3f} {p:>11.3f} {10_000 * p:>38.0f}")

Output obtained:

 prevalence   precision          true ill per 10,000 positives
      0.001       0.018                                    177
      0.010       0.154                                   1538
      0.100       0.667                                   6667
      0.300       0.885                                   8852
      0.500       0.947                                   9474
      0.800       0.986                                   9863

How to read the code. The function precisione is Bayes' theorem in the form of a count: in the numerator the positives coming from the ill (sensitivity × prevalence), in the denominator all the positives (true plus false). Varying only the prevalence, keeping sensitivity and false positives fixed, the precision goes from about 1.8% (out of 1000 people only one is ill, and the false alarms of the healthy swamp the true positives) to over 90% (with a prevalence of 80% almost every positive is ill: 98.6% in the output above). Here is the lesson: a test does not have "the" precision, it has one precision for each population in which it is used. A screening test applied to the whole population (low prevalence) produces many false alarms; the same test applied to people already suspected (high prevalence) is much more reliable. It is also the reason why a classifier evaluated in the lab on balanced data may disappoint in production, where the class of interest is rare.

Challenge

• Easy. Prove P(A ∪ B ∪ C) for three events, by analogy with inclusion-exclusion, and verify it on dice.
• Medium. Prove that, if A and B are independent, so are A and Bᶜ. (Hint: P(A ∩ Bᶜ) = P(A) − P(A ∩ B).)
• Hard. An urn contains two coins: one fair and one with two heads. You draw one at random and toss it three times, obtaining three heads. Compute with Bayes the probability that it is the loaded one and verify by simulation.
• Expert. In Monty Hall, suppose the host opens a door at random among the two remaining (and luckily finds a goat). Show with Bayes' theorem that now switching and staying are equivalent and simulate both cases conditioning on the scenarios with a goat.
• Enterprise. For a real classifier of your organisation, estimate sensitivity and false positive rate, the base rate in real traffic, and compute the expected precision in production. Compare it with the one measured in the lab.

The first two challenges are short proofs, the third and the fourth apply Bayes to a case with hidden data and must be checked by simulation as in the lesson. If you have time for only one, pick Enterprise: it is the only one that produces a number useful for your work, and if the positive class is rare in production, the expected precision will be lower than the one measured on balanced data.

To practise on the simulated gate you will find the skeleton to complete, with the commented solutions at the bottom, in esercizi.py: it reads conditionals and independence from a contingency table with Fraction, computes Bayes as the prevalence varies and estimates the same quantity by simulation with a fixed seed. It requires only Python 3.12 or later and is launched with python3 esercizi.py; the file was not executed for this article.

What you can do now

You can translate a problem into a sample space and events, derive the rules for the complement and the union from the three axioms, compute a conditional by narrowing the room and tell disjoint from independent. You can use Bayes' theorem with natural frequencies, and explain why the precision of a test depends on the population in which you use it. The next step is random variables, which turn outcomes into numbers; statistics will take the opposite road, estimating unknown probabilities from data, and the Bayesian classifier will be taken up again in Naive Bayes.

The same series also includes: 28 - Bayesian Statistics: Priors, Posteriors, Conjugates and Why Regularisation Is a Belief (which refers back to this lesson).