Bayes' theorem

Intuition

A test that catches 99% of the people who have a disease comes back positive. Most people conclude they are probably sick. They are wrong, and not by a little: the answer here is 7.4%.

The mistake is reading a number about the test as a number about the patient. "Sadly, p(A|B) is not the same as p(B|A), but there is a way to get from one to the other" [1, §1.3], and that way is Bayes' theorem. What it adds to the test's accuracy is the thing the accuracy leaves out — how many people had the disease before anyone was tested.

The theorem is a piece of arithmetic nobody disputes. What makes it worth a page is that its answer contradicts a strong intuition, in situations that matter, and the contradiction is not a subtlety. It is a factor of more than thirteen.

Formal definition

The derivation is three lines. A conjunction can be written two ways — p(A and B) = p(A) p(B|A) for any A and B [1, §1.2], and the same with the letters swapped — so the two products are equal, and dividing through by p(B) gives the theorem [1, §1.4].

The version worth memorising is the one in the frontmatter above, with H for a hypothesis and D for data. This is the diachronic interpretation, so called because "something is happening over time; in this case the probability of the hypotheses changes, over time, as we see new data" [1, §1.5]. Each term has a name:

• p(H) — the prior, the probability of the hypothesis before the data.
• p(D|H) — the likelihood, the probability of the data if the hypothesis is true. "The likelihood is usually the easiest part to compute."
• p(D) — the normalizing constant, "the probability of the data under any hypothesis at all", which "in the most general case ... is hard to nail down".
• p(H|D) — the posterior, what you wanted.

All four are quoted from [1, §1.5].

The normalizing constant stops being hard the moment the hypotheses are mutually exclusive — "at most one hypothesis in the set can be true" — and collectively exhaustive — "there are no other possibilities" [1, §1.5]. Then the law of total probability gives p(D) as a sum over the hypotheses, and it is just the column total.

Often you can skip it entirely. Comparing hypotheses against the same data, "the denominator is the normalizing term p(x) which does not depend on the value of the label y. As a result, we only need to worry about comparing the numerator across different values of y" [2, §22.9.3]. Every classifier that scores prior × likelihood and takes the largest is doing exactly this.

How it works

The figure is the screening test drawn out. Four people in a thousand have the disease and the test finds all but a hundredth of them; fifty of the remaining 996 are flagged anyway, because 5% of a large number is larger than 99% of a small one. Fifty-four positives, four of them real.

Nothing about the test is defective. Sensitivity and specificity are properties of the instrument; the prior is a property of the population, and the posterior is the product. Change the population — test people with symptoms rather than everyone — and the same instrument gives a completely different answer [3, §2.6.5].

Where it's used

• Diagnostics and screening, where the base rate is usually the dominant term and is usually the one omitted from the report.
• Naive Bayes classifiers, which apply the theorem per class and then assume the features are conditionally independent given the label. Without some such assumption the chain rule leaves "roughly 2^d parameters" to estimate [2, §22.9.3] — and a model that needs "to see every single possible example in order to predict the corresponding label" is "not really learning a pattern but just memorizing the dataset" [2, §22.9.2].
• Spam filtering, A/B tests, sensor fusion — anywhere evidence arrives in pieces and yesterday's posterior becomes today's prior.
• Bayesian regression and hierarchical models, where the same update runs over a continuous parameter instead of a handful of hypotheses. See Bayesian linear regression.
• Generative classifiers such as GDA, which model p(x|y) and invert it.

Why & when

Use it whenever you have the conditional the wrong way round. That is the whole trick: "if p(A|B) is hard to compute, or hard to measure experimentally, check whether it might be easier to compute the other terms in Bayes's theorem" [1, §1.8].

Do not fear the prior on the grounds that it is subjective. With enough data it stops mattering — "with enough data, people who start with different priors will tend to converge on the same posterior" [1, §4.3]. The worked example measures how quickly.

Do fear a prior of zero. "In a Bayesian update, we multiply each prior probability by a likelihood, so if p(H) is 0, p(H|D) is also 0, regardless of D" [1, §4.6]. This is Cromwell's rule — "avoid giving a prior probability of 0 to any hypothesis that is even remotely possible" [1, §4.6].

Reach for a conjugate prior when the update has to be cheap. For a binomial likelihood the beta distribution is conjugate, so "if the prior is a beta distribution with parameters alpha and beta, and we see data with h heads and t tails, the posterior is a beta distribution with parameters alpha+h and beta+t. In other words, we can do an update with two additions" [1, §4.5].

Worked example

Two questions, both answered by counting rather than by asserting. First: how much of a positive result is the test and how much is the population? Second: how long does a strongly held prior survive contact with data?

import numpy as np

rng = np.random.default_rng(0)

# --- 1. a screening test, simulated one person at a time -------------------
PREV, SENS, SPEC = 0.004, 0.99, 0.95
N = 2_000_000

sick = rng.random(N) < PREV
pos1 = np.where(sick, rng.random(N) < SENS, rng.random(N) >= SPEC)

print(f"  {N:,} people, prevalence {PREV:.1%}, sensitivity {SENS:.0%}, specificity {SPEC:.0%}")
print(f"    tested positive          {pos1.sum():>8,}")
print(f"    of those, actually sick  {(pos1 & sick).sum():>8,}")
print(f"    P(sick | positive)         {(pos1 & sick).sum() / pos1.sum():>8.4f}   counted")
bayes = PREV * SENS / (PREV * SENS + (1 - PREV) * (1 - SPEC))
print(f"                               {bayes:>8.4f}   Bayes's theorem")

# --- 2. the same people take the test a second time ------------------------
pos2 = np.where(sick, rng.random(N) < SENS, rng.random(N) >= SPEC)
both = pos1 & pos2
print(f"\n    positive twice           {both.sum():>8,}")
print(f"    of those, actually sick  {(both & sick).sum():>8,}")
print(f"    P(sick | positive twice)   {(both & sick).sum() / both.sum():>8.4f}   counted")
b2 = bayes * SENS / (bayes * SENS + (1 - bayes) * (1 - SPEC))
print(f"                               {b2:>8.4f}   yesterday's posterior as today's prior")

# --- 3. three analysts who disagree about the coin -------------------------
TRUE_P = 0.35
priors = {
    "uniform      Beta(1, 1)": (1.0, 1.0),
    "sure it's ~0.9  B(30, 3)": (30.0, 3.0),
    "sure it's ~0.1  B(3, 30)": (3.0, 30.0),
}
flips = rng.random(10_000) < TRUE_P
print(f"\n  a coin with p = {TRUE_P}, three priors, posterior mean after n flips")
print("    " + "prior".ljust(26) + "".join(f"{n:>9}" for n in (0, 10, 100, 1_000, 10_000)))
for label, (a, b) in priors.items():
    row = []
    for n in (0, 10, 100, 1_000, 10_000):
        h = int(flips[:n].sum())
        row.append(f"{(a + h) / (a + b + n):>9.3f}")
    print("    " + label.ljust(26) + "".join(row))

# --- 4. an analyst who ruled the truth out in advance ----------------------
grid = np.linspace(0.0005, 0.9995, 1000)
flat = np.ones_like(grid)
closed = np.where(grid >= 0.5, 1.0, 0.0)  # "I am certain the coin favours heads"
print()
for prior, name in ((flat, "open mind    "), (closed, "p >= 0.5 only")):
    for n in (100, 10_000):
        h = int(flips[:n].sum())
        # in log space: theta**3500 underflows to zero in float64
        loglik = h * np.log(grid) + (n - h) * np.log(1 - grid)
        post = prior * np.exp(loglik - loglik.max())
        post = post / post.sum()
        print(f"    {name}  after {n:>6,} flips   posterior mean {(grid * post).sum():.3f}")

  2,000,000 people, prevalence 0.4%, sensitivity 99%, specificity 95%
    tested positive           107,359
    of those, actually sick     7,940
    P(sick | positive)           0.0740   counted
                                 0.0737   Bayes's theorem

    positive twice             12,776
    of those, actually sick     7,851
    P(sick | positive twice)     0.6145   counted
                                 0.6116   yesterday's posterior as today's prior

  a coin with p = 0.35, three priors, posterior mean after n flips
    prior                             0       10      100     1000    10000
    uniform      Beta(1, 1)       0.500    0.417    0.314    0.352    0.351
    sure it's ~0.9  B(30, 3)      0.909    0.791    0.459    0.370    0.353
    sure it's ~0.1  B(3, 30)      0.091    0.163    0.256    0.344    0.351

    open mind      after    100 flips   posterior mean 0.314
    open mind      after 10,000 flips   posterior mean 0.351
    p >= 0.5 only  after    100 flips   posterior mean 0.512
    p >= 0.5 only  after 10,000 flips   posterior mean 0.501

The arithmetic and the counting agree, which is the point of running it. 0.0737 against 0.0740 over two million people: Bayes' theorem is not an approximation or a modelling choice here, it is what counting gives you.

A second positive changes everything. The posterior goes from 0.0740 to 0.6145, because the update ran with 0.0737 as the prior rather than 0.004. This is the diachronic reading made concrete [1, §1.5], and it is why repeat testing is the standard response to a rare-disease screen.

Note that the twice-positive figure is slightly above the chained arithmetic (0.6145 against 0.6116). The two tests were simulated as independent given the true status, which is the assumption the chaining makes; the gap is sampling noise over 12,776 people, not a flaw in either.

Two analysts who could not disagree more end up in the same place. One starts at 0.909, the other at 0.091 — a nine-to-one belief in opposite directions. After ten thousand flips they are at 0.353 and 0.351, against a truth of 0.35. Neither was argued out of their position; they were outvoted by data. That is swamping the priors [1, §4.3], measured.

The one prior that never recovers is the one that ruled the answer out. The analyst certain the coin favours heads is still at 0.501 after ten thousand flips of a coin that lands heads 35% of the time — pinned against the edge of what they were willing to consider, and no quantity of evidence moves them, because every hypothesis that could explain the data was multiplied by zero [1, §4.6].

Pitfalls

Confusing p(D|H) with p(H|D). The test is 99% sensitive and a positive result is 7.4% informative, at the same time, with no contradiction. This is the single most consequential arithmetic error in applied statistics.

Dropping the base rate. It is the term nobody measures and everybody assumes. If the prevalence is unknown, the posterior is unknown — the accuracy of the test cannot supply it.

Chaining updates that are not conditionally independent. The worked example runs two tests that are independent given the true status. Two runs of the same assay on the same badly drawn sample are not, and treating them as if they were inflates the posterior: run the same test twice and "no additional insight is gained" [3, §2.6.5].

A prior of zero, or a hypothesis set that is not exhaustive. Both are the same error: the truth is outside the suite, and "if p(H) is 0, p(H|D) is also 0, regardless of D" [1, §4.6]. So the posterior piles up on the least-wrong hypothesis it was allowed and reports high confidence — the worked example's analyst who allowed only p ≥ 0.5, pinned at 0.501.

Reading a posterior as if the model were certain. "All Bayesian analysis is based on modeling decisions. If you and I do not choose the same model, we might interpret data differently" [1, §4.6]. Priors get swamped; a wrong likelihood function does not.

Multiplying raw probabilities over a long dataset. The last block above needs its log-space rewrite: theta ** 3500 is zero in float64, and the posterior comes out as nan. See maximum likelihood.