Chapter 3
Probability

Lectures 5–8 · the language of uncertainty, how to reason with it, and how it goes wrong

You use probability every day without calling it that. "70% chance of rain." "One in a billion chance of a random DNA match." "This email is probably spam." "If I take the express I'll probably make it." Each is a statement about something uncertain, and each guides a decision: take an umbrella, convict or acquit, delete the email, run for the train.

Chapter 1 drew the map of the course: probability reasons forwards, from a mechanism to the data it produces. Before we can reason backwards from data to mechanisms (inference), we need a precise language for uncertainty. This chapter builds it. The rules turn out to be few and simple. The trouble is that human intuition about probability is famously bad, and some of the mistakes it leads to have sent innocent people to prison. So alongside the rules, we'll keep meeting situations where your gut says one thing and the mathematics says another, and you'll learn to trust the mathematics.

Probability serves both of the course's goals, in opposite directions. Reasoning forwards, from a known mechanism to the data it will produce, is prediction: if a test is 90% sensitive, gives false alarms for 9% of healthy people, and the disease affects 1% of people, how many positive tests should a clinic expect? Reasoning backwards, from observed data to the mechanism that produced them, is inference: given that this patient tested positive, how likely is it that she has the disease? Most of this chapter builds the forward machinery; Section 3.6 then shows how Bayes' theorem runs it in reverse, and that reversal is the seed from which all of statistical inference grows.

Here is the route through the chapter.

• Sections 3.1 to 3.3 set up the language: what a probability means, how to describe events as sets, and the three rules every probability must obey.
• Section 3.4 turns probability into counting, for situations where every outcome is equally likely.
• Sections 3.5 and 3.6 are the heart of the chapter: how probabilities change when we learn something (conditional probability), how to translate a word problem into symbols, and how to reason backwards from evidence to causes (Bayes' theorem).
• Section 3.7 shows how to let a computer do the work by simulation.
• Section 3.8 attaches numbers to outcomes, which sets up the next two chapters.

Lecture 5 — What probability means and its rules

3.1 What does probability mean?

Three ways to read a probability

What does it mean to say a coin lands heads with probability $\frac{1}{2}$? People have given three different answers.

1. Equally likely outcomes (classical). There are two faces and nothing to favour either, so each gets half. This works for coins, dice and lotteries, where symmetry tells us the outcomes are equally likely. It doesn't tell us the probability that a drug works or that it rains.
2. Long-run frequency (frequentist). Flip the coin many, many times: the fraction of heads settles down to $\frac{1}{2}$. A probability is the fraction of the time something happens in a long run of repetitions. This is how a casino, an insurer, or a lab estimating a test's error rate thinks.
3. Degree of belief (Bayesian). A probability measures how strongly a reasonable person should believe something, given what they know. This makes sense even for one-off events: "the probability that this suspect left this fingerprint", "the probability that this vaccine works".

The three interpretations disagree about what probability is, but they agree on the rules probability must follow. That's what makes probability so useful: the same rules serve the casino, the scientist and the juror.

Here is the long-run idea in action. Toss a coin over and over, and track the fraction of heads.

Early on the fraction swings wildly; after hundreds of tosses it barely moves. This settling-down is called the Law of Large Numbers. We'll state it precisely and prove it in Chapter 7; for now, it justifies estimating a probability by a long-run fraction.

set.seed(3)
for (n in c(10, 100, 1000, 100000)) {
  rolls <- sample(1:6, n, replace = TRUE)
  cat(format(n, big.mark = ",", scientific = FALSE), "rolls: fraction of sixes =", mean(rolls == 6), "\n")
}

10 rolls: fraction of sixes = 0 
100 rolls: fraction of sixes = 0.15 
1,000 rolls: fraction of sixes = 0.174 
100,000 rolls: fraction of sixes = 0.16699 

set.seed(6)
rolls <- sample(1:6, 10000, replace = TRUE)
running <- cumsum(rolls == 6) / seq_along(rolls)
ggplot(data.frame(n = seq_along(rolls), running), aes(n, running)) +
  geom_line() +
  geom_hline(yintercept = 1/6, linetype = "dashed") +
  scale_x_log10() +
  labs(x = "number of rolls (log scale)", y = "fraction of sixes")

Empirical probability: probabilities from data

For a die we know the probability in advance. Most of the time we don't. What is the probability that a randomly chosen eruption of Old Faithful lasts more than 3 minutes? That a hospital patient with these symptoms needs intensive care? That a new email is spam? Nobody can work these out from symmetry. Instead, we look at what happened in the past and count.

The Law of Large Numbers is what makes empirical probabilities worth computing: if the observations are independent repetitions of the same experiment, then as $n$ grows, $\widehat{\mathrm{\text{Pr}}}(A)$ settles down to the true $\mathrm{\text{Pr}}(A)$. The empirical probability is the statistic; the true probability is the parameter (Section 1.3).

e <- faithful$eruptions
sum(e > 3)              # n_A: eruptions longer than 3 minutes
length(e)               # n
mean(e > 3)             # empirical probability n_A / n
mean(e[1:20] > 3)       # using only the first 20 eruptions

[1] 175
[1] 272
[1] 0.6433824
[1] 0.55

running <- cumsum(e > 3) / seq_along(e)
ggplot(data.frame(k = seq_along(e), running), aes(k, running)) +
  geom_line() +
  geom_hline(yintercept = mean(e > 3), linetype = "dashed") +
  labs(x = "number of eruptions observed", y = "fraction over 3 min")

3.2 Experiments, outcomes and events

To reason precisely we need a vocabulary. Any process with an uncertain result, from rolling a die to running a DNA test, is called an experiment.

Some examples:

• A drug screen returns negative, inconclusive or positive: $S$$=\{\mathrm{\text{neg}},$$\mathrm{\text{inc}},$$\mathrm{\text{pos}}\}$. The event "not negative" is $\{\mathrm{\text{inc}},\mathrm{\text{pos}}\}$.
• We count failed login attempts on a server in one hour: $S$$=\{0,$$1,$$2,$$\dots \}$. The event "more than 5 attempts" is $\{6,7,8,\dots \}$.
• We measure the time, in hours, until a lab centrifuge fails: $S=[0,\infty )$. The event "it lasts at least 1,000 hours" is the interval $[1000,\infty )$.

Why bother with sets? Because it turns vague questions into precise ones. "What's the chance the test isn't negative?" becomes "what's the probability of the set $\{\mathrm{\text{inc}},\mathrm{\text{pos}}\}$?", and sets can be combined, compared and counted.

Combining events

Events are sets, so we combine them with set operations. In words and in symbols:

In mathematics, "or" always means "one or the other or both." Venn diagrams make these operations visible: draw $S$ as a rectangle and events as regions inside it. The shaded outcomes belong to the event.

Some combinations are built from these: "$A$ but not $B$" is $A\cap B^c$, and "exactly one of $A$ and $B$" is $(A\cap B^c)\cup (A^c\cap B)$. The best way to get fluent is to practise in both directions: from a picture to symbols, and from symbols to a picture.

You can check De Morgan's laws in the widget: shade "neither" in the two-event diagram, and the readout offers both ways of writing it.

3.3 The rules of probability

What rules should any sensible assignment of probabilities follow? Think about long-run fractions, or empirical probabilities. A fraction of the time is never negative and never more than 1. Something certain to happen happens 100% of the time. And if two events can't happen together, the fraction of times that one or the other happens is the sum of their fractions. Remarkably, these three properties are all we need.

Everything else follows from these. Here are the consequences you'll use constantly.

The addition rule is easiest to see in a picture. Move the sliders, and watch which areas get counted twice.

tt <- as.data.frame(Titanic)                 # one row per combination, with a count
xtabs(Freq ~ Sex + Age, data = tt)           # cross-tabulate sex by age

        Age
Sex      Child Adult
  Male      64  1667
  Female    45   425

Lecture 6 — Counting and equally likely outcomes

3.4 Equally likely outcomes and counting

When all outcomes in a finite sample space are equally likely (fair dice, shuffled cards, random samples, random passwords), the axioms give a simple formula.

So probability becomes counting. The trouble is that the counts are often astronomically large: nobody can list all the possible poker hands, let alone all the possible 8-character passwords. We need shortcuts, and they all grow from one idea.

The multiplication rule

For example, a meal deal with 3 mains, 4 sides and 2 drinks gives $3\times 4\times 2=24$ different meals.

The four ways to choose

Most counting problems ask: in how many ways can we choose $k$ items from $n$ distinct items? The answer depends on two questions, which you should ask every time:

1. Does order matter? Is "Alice then Bob" different from "Bob then Alice"? For a race (gold, silver), yes; for a committee, no.
2. Is replacement allowed? Can the same item be chosen more than once? For the digits of a PIN, yes; for people on a committee, no.

Two yes/no questions give four cases. We take them one at a time, each with a motivating example, a formula and a proof.

Here $n!$$=n$$\times (n$$-1)$$\times \cdots$$\times 2$$\times 1$ ("$n$ factorial"), with the convention $0!=1$. In R: factorial(n) and choose(n, k).

Case 1: order matters, with replacement

A PIN has 4 digits, each from 0 to 9, and digits can repeat. There are 10 choices for the first digit, 10 for the second, and so on: ${10}^4=\mathrm{10,000}$ PINs.

m <- c(26, 62, 94)
n_passwords <- m^8
seconds <- n_passwords / 1e10
data.frame(symbols = m, passwords = signif(n_passwords, 3),
           hours = round(seconds / 3600, 3), days = round(seconds / 86400, 2))

  symbols passwords   hours days
1      26  2.09e+11   0.006 0.00
2      62  2.18e+14   6.065 0.25
3      94  6.10e+15 169.325 7.06

len <- 4:14
grid <- expand.grid(length = len, symbols = factor(c(26, 62, 94)))
grid$seconds <- as.numeric(as.character(grid$symbols))^grid$length / 1e10
ggplot(grid, aes(length, seconds, linetype = symbols)) +
  geom_line() +
  scale_y_log10(breaks = 10^c(-6, 0, 6, 12, 18), labels = scales::math_format(10^.x)(c(-6, 0, 6, 12, 18))) +
  scale_linetype_manual(values = c("solid", "dashed", "dotted")) +
  geom_hline(yintercept = c(86400, 3.15e9), colour = "grey60") +
  annotate("text", x = 4, y = c(86400, 3.15e9) * 6, label = c("1 day", "1 century"), hjust = 0, vjust = -0.4, size = 3.8, colour = "grey40") +
  labs(x = "password length", y = "seconds (log scale)",
       linetype = "symbols")

Case 2: order matters, without replacement

Eight runners compete for gold, silver and bronze. Any of the 8 can win gold; once gold is decided, any of the remaining 7 can win silver; then any of the remaining 6, bronze: $8$$\times 7$$\times 6$$=336$ possible podiums. Each choice uses up one runner, so the number of options falls by one at each stage.

Factorials grow astonishingly fast. A deck of 52 cards can be shuffled into $52!\approx 8\times {10}^{67}$ orders, more than the number of atoms in the Earth, so every thorough shuffle you make has almost certainly never happened before.

factorial(6)
factorial(5) / factorial(6)

[1] 720
[1] 0.1666667

1 - prod(365:343) / 365^23      # by the formula
pbirthday(23)                   # R has a built-in function

[1] 0.5072972
[1] 0.5072972

people <- 1:80
p_match <- sapply(people, pbirthday)
ggplot(data.frame(people, p_match), aes(people, p_match)) +
  geom_line() +
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  geom_vline(xintercept = 23, linetype = "dotted") +
  labs(x = "number of people", y = "Pr(shared birthday)")

Case 3: order doesn't matter, without replacement

A lab must choose 2 of its 5 analysts to review a case. Ordered pairs (first reviewer, second reviewer) number $5\times 4=20$ by case 2. But a review panel of Ana and Ben is the same panel as Ben and Ana: each panel was counted twice, once in each order. So there are $20/2=10$ panels. The same idea, dividing out the orderings, works in general.

Two quick checks: $\left(\genfrac{}{}{0ex}{}{n}{0}\right)$$=\left(\genfrac{}{}{0ex}{}{n}{n}\right)$$=1$ (one way to choose nothing, one way to choose everything), and $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$$=\left(\genfrac{}{}{0ex}{}{n}{n-k}\right)$ (choosing which $k$ to take is the same as choosing which $n-k$ to leave). Combinations are the most important counting tool in the course: Chapter 4's binomial distribution is built on them.

choose(69, 5)
choose(69, 5) * 26

[1] 11238513
[1] 292201338

choose(1900, 20) / choose(2000, 20)

[1] 0.356686

k <- 1:100
p_miss <- choose(1900, k) / choose(2000, k)
ggplot(data.frame(k, p_miss), aes(k, p_miss)) +
  geom_line() +
  geom_hline(yintercept = 0.05, linetype = "dashed") +
  labs(x = "number of bags tested", y = "Pr(no fakes found)")

Case 4: order doesn't matter, with replacement

The last case is the least obvious. A lab orders 10 reagent kits, and each kit can be one of 4 types (A, B, C, D). How many different orders are possible? Order doesn't matter (an order is just "how many of each type"), and types can repeat. An order might be 3 of A, 0 of B, 5 of C and 2 of D.

The trick, called stars and bars, is to draw each order as a row of 10 stars (the kits) separated by $4-1=3$ bars into 4 compartments (the types). The order above is $$\star \star \star ||\star \star \star \star \star |\star \star$$ Three stars before the first bar (type A), none between the first and second bars (B), five between the second and third (C), and two after the third (D). Every order gives exactly one such row, and every row of 10 stars and 3 bars gives exactly one order. So we just count rows: 13 positions, of which we choose 3 for the bars, giving $\left(\genfrac{}{}{0ex}{}{13}{3}\right)=286$ orders.

The widget lists every selection for small cases, so you can see where all four formulas come from.

Lecture 7 — Conditional probability and independence

3.5 Conditional probability and independence

Zooming in

Of the 2,201 people aboard the Titanic, 711 survived: $\mathrm{\text{Pr}}(\mathrm{\text{survived}})$$=711/2201$$\approx 0.32$. But that number hides enormous differences. Suppose you learn that the person you picked was a woman. Your probability should change, because women were far more likely to be given places in the lifeboats. Learning something narrows the sample space to the outcomes consistent with what you learned. Try it.

Among the 470 women, 344 survived, so the probability of survival given that the person was female is $344/470$$\approx 0.73$. Notice how this was computed: restrict attention to the women, then take the fraction who survived. In a Venn diagram, the new information shrinks the sample space from the whole rectangle to the circle "female", and we ask what share of that circle lies inside "survived". Dividing both counts by 2,201 turns this into probabilities, and gives the general definition.

Check it on the Titanic: $\mathrm{\text{Pr}}(\mathrm{\text{survived}}\mid \mathrm{\text{female}})$$=\frac{\mathrm{\text{Pr}}(\mathrm{\text{survived}}\cap \mathrm{\text{female}})}{\mathrm{\text{Pr}}(\mathrm{\text{female}})}$$=\frac{344/2201}{470/2201}$$=\frac{344}{470}$$\approx 0.73.$ The 2,201s cancel: conditioning is just counting within the smaller group.

Rearranging the definition gives a rule for "and" that we'll use all the time.

In words: the chance that both happen is the chance the first happens, times the chance the second happens given the first. Drawing two cards from a deck without replacement, the chance both are aces is $\frac{4}{52}\times \frac{3}{51}$: after the first ace, only 3 aces remain among 51 cards.

tt <- as.data.frame(Titanic)
by_class <- xtabs(Freq ~ Class + Survived, data = tt)
by_class
round(prop.table(by_class, 1), 3)     # Pr(survived | class): divide by row totals
by_sex <- xtabs(Freq ~ Sex + Survived, data = tt)
round(prop.table(by_sex, 1), 3)       # Pr(survived | sex)
round(prop.table(by_sex, 2), 3)       # Pr(sex | survived): divide by column totals

      Survived
Class   No Yes
  1st  122 203
  2nd  167 118
  3rd  528 178
  Crew 673 212
      Survived
Class     No   Yes
  1st  0.375 0.625
  2nd  0.586 0.414
  3rd  0.748 0.252
  Crew 0.760 0.240
        Survived
Sex         No   Yes
  Male   0.788 0.212
  Female 0.268 0.732
        Survived
Sex         No   Yes
  Male   0.915 0.516
  Female 0.085 0.484

rate <- prop.table(by_class, 1)[, "Yes"]
ggplot(data.frame(class = names(rate), rate), aes(class, rate)) +
  geom_col(width = 0.6) +
  geom_hline(yintercept = 711 / 2201, linetype = "dashed") +
  labs(x = NULL, y = "Pr(survived | class)")

Reading a probability word problem

Most mistakes in probability happen before any calculation: in reading the question. A word problem hides its probabilities in everyday phrases, and the hardest part is deciding which phrase is a plain probability, which is an "and", and which is a "given". This algorithm works for almost every problem in the course.

The phrases below are clues to the translation. The most important one is that "of", "among", "if", "given that" and "for those who" all introduce the group after the bar.

ggplot(data.frame(source = c("genuine victims", "non-victims"), n = c(420, 88)),
       aes(x = "reports", y = n, fill = source)) +
  geom_col(width = 0.5) + coord_flip() +
  scale_fill_manual(values = c("grey40", "orange")) +
  labs(x = NULL, y = "reports per 10,000 adults", fill = NULL)

Independence

Sometimes learning that $B$ happened tells us nothing about $A$: the probability stays the same.

Why are the two forms the same? If $\mathrm{\text{Pr}}(B)>0$, dividing $\mathrm{\text{Pr}}(A\cap B)$$=\mathrm{\text{Pr}}(A)\mathrm{\text{Pr}}(B)$ by $\mathrm{\text{Pr}}(B)$ gives $\mathrm{\text{Pr}}(A\cap B)/\mathrm{\text{Pr}}(B)$$=\mathrm{\text{Pr}}(A)$, which says $\mathrm{\text{Pr}}(A\mid B)$$=\mathrm{\text{Pr}}(A)$; the steps reverse too.

For the Titanic, $\mathrm{\text{Pr}}(\mathrm{\text{female}})\mathrm{\text{Pr}}(\mathrm{\text{survived}})$$=\frac{470}{2201}$$\times \frac{711}{2201}$$\approx 0.069$, but $\mathrm{\text{Pr}}(\mathrm{\text{female}}\cap \mathrm{\text{survived}})$$=\frac{344}{2201}$$\approx 0.156$. Not even close: sex and survival were strongly dependent.

In practice we usually assume independence because of how an experiment works (separate coin flips, separate patients, separate random draws with replacement), and then use the product rule to compute probabilities that would otherwise be impossible. Assuming independence is a modelling decision, and it deserves the same scepticism as any other assumption.

1 - (5/6)^4
1 - (35/36)^24

[1] 0.5177469
[1] 0.4914039

set.seed(1654)
n <- 20000
bet1 <- replicate(n, any(sample(1:6, 4, replace = TRUE) == 6))
bet2 <- replicate(n, any(sample(1:6, 24, replace = TRUE) == 6 & sample(1:6, 24, replace = TRUE) == 6))
d <- rbind(data.frame(play = 1:n, wins = cumsum(bet1) / (1:n), bet = "bet 1: one six in 4 rolls"),
           data.frame(play = 1:n, wins = cumsum(bet2) / (1:n), bet = "bet 2: double six in 24 rolls"))
ggplot(d, aes(play, wins, linetype = bet)) +
  geom_line() +
  geom_hline(yintercept = 0.5, colour = "grey60") +
  scale_x_log10() + coord_cartesian(ylim = c(0.3, 0.7)) +
  labs(x = "number of plays (log scale)", y = "fraction won so far", linetype = NULL)

set.seed(46)
n <- 100000
A <- rbinom(n, 1, 0.10)                  # disease A: 10% of people
B <- rbinom(n, 1, 0.10)                  # disease B: 10%, independent of A
# admitted if ill and the illness is serious enough (50%), or for other reasons (2%)
hosp <- ((A == 1 | B == 1) & runif(n) < 0.5) | (runif(n) < 0.02)
# Pr(B | A) and Pr(B | not A), in everyone and among hospital patients
c(with_A = mean(B[A == 1]), without_A = mean(B[A == 0]))
c(with_A = mean(B[A == 1 & hosp]), without_A = mean(B[A == 0 & hosp]))

    with_A  without_A 
0.09671642 0.10058922 
    with_A  without_A 
0.09903661 0.73431199 

Lecture 8 — Bayes' theorem, simulation and random variables

3.6 Thinking backwards: total probability and Bayes' theorem

From clues to causes

A detective arrives at a crime scene. She doesn't see the crime; she sees its effects: a footprint, a broken lock, a missing file. Her job is to reason backwards from the effects to the cause. A doctor does the same with symptoms, a mechanic with a strange noise, a scientist with an experiment's results.

Probability, so far, has gone forwards: if this person has the disease, how likely is a positive test? If this suspect is the thief, how likely is this footprint? Those forward probabilities are usually the ones we know, because we can study them: run the test on people known to be ill, measure the suspect's shoes. But the question we need answered goes backwards: given the positive test, how likely is the disease? Given the footprint, how likely is this suspect?

Two results connect the directions. The law of total probability says how likely a clue is overall, adding up over all the possible causes. Bayes' theorem then says what share of that total each cause is responsible for. Both are easiest to understand, and much harder to get wrong, if you start by counting.

Where do the clues come from? Total probability