Level 3 · Statistics & Probability: Relationships

So far every dataset has been a single list: one measurement per person or per day. The interesting questions in statistics are usually about two things at once. Do students who study longer score higher? Does a positive test mean you have the disease? This level is about relationships: how two quantities move together, and how knowing one thing changes the probability of another.

The first half treats paired numbers. You will measure how two variables vary together (covariance), turn that into a scale-free number between −1 and 1 (correlation), and fit the straight line that predicts one from the other with the smallest squared error (least-squares regression). The residuals, what the line gets wrong, turn out to be as informative as the line itself. The second half does the same for events: conditional probability, independence, and Bayes' rule, all computed from tables of counts, where they are easiest to see and hardest to get wrong.

Both halves reappear in the Machine Learning path. The regression line here is the "fitting a line" model there, found by a formula here and by gradient descent at Level 4 of that path. Bayes' rule on tables is the whole idea behind the Naive Bayes classifier. And the warning of this level, that a relationship in data is not proof of cause and effect, applies to every model that learns from examples.

Covariance: moving together

Here are hours of revision and test scores for eight students, as two lists where position i belongs to student i:

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
for h, s in zip(hours, scores):
    print(h, s, "#" * ((s - 50) // 2))

1 52 #
2 55 ##
2 61 #####
3 60 #####
4 68 #########
5 70 ##########
6 77 #############
8 85 #################

The bars grow as the hours grow: the two variables move together. To measure that, look at each student's deviation from the mean in both variables and multiply them. A student above average in both gives a positive product, as does one below average in both. A student above in one and below in the other gives a negative product. The mean of these products is the covariance:

Here x_i and y_i are student i's hours and score, and \bar{x}, \bar{y} the two means. A sum over pairs is a loop over zip: each pair (x_i, y_i) arrives together as x, y.

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
n = len(hours)
mx, my = sum(hours) / n, sum(scores) / n
products = [(x - mx) * (y - my) for x, y in zip(hours, scores)]
print("means:", mx, my)
print([round(p, 1) for p in products])
print("covariance:", sum(products) / n)

means: 3.875 66.0
[40.2, 20.6, 9.4, 5.2, 0.2, 4.5, 23.4, 78.4]
covariance: 22.75

Every product is positive here, so the covariance is clearly positive: more hours go with higher scores. A negative covariance would mean that one goes up as the other goes down, and a covariance near zero means no linear tendency either way. As with the variance, divide by n - 1 instead of n when the data is a sample and you want to estimate the population covariance (that is what statistics.covariance does).

The covariance has a serious weakness: its size depends on the units. Measure revision in minutes instead of hours and the covariance becomes sixty times larger, although the relationship is exactly the same. So you can read its sign, but not its size. The next section fixes that.

Common mistakes: pairing the wrong values (sorting one list but not the other destroys the pairing; sort pairs, never columns separately). Mixing the n and n - 1 versions between the covariance and the variances they are compared with. Reading a large covariance as a strong relationship when it only reflects large units.

Practice: Wind and Kites

Correlation

Dividing the covariance by both standard deviations removes the units. The result is the correlation coefficient r (Pearson's correlation), which is always between −1 and 1:

where s_x and s_y are the standard deviations of the two lists (all three dividing by n, as in the code).

import math

def correlation(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    cov = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / n
    sx = math.sqrt(sum((x - mx) ** 2 for x in xs) / n)
    sy = math.sqrt(sum((y - my) ** 2 for y in ys) / n)
    return cov / (sx * sy)

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
minutes = [60 * h for h in hours]
print(round(correlation(hours, scores), 4))
print(round(correlation(minutes, scores), 4))

0.984
0.984

Changing hours to minutes leaves r unchanged, as it should. Another way to see it: r is the covariance of the z-scores. Standardise both lists (Level 2), multiply pairwise and average.

How to read r: +1 means the points lie exactly on a rising straight line, −1 exactly on a falling one, and 0 means no linear trend. Here r ≈ 0.98, a very strong positive linear relationship. Three warnings go with every correlation you compute.

Correlation only measures straight-line relationships. A perfect curve can have r = 0:

import math

def correlation(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    cov = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / n
    sx = math.sqrt(sum((x - mx) ** 2 for x in xs) / n)
    sy = math.sqrt(sum((y - my) ** 2 for y in ys) / n)
    return cov / (sx * sy)

xs = [-3, -2, -1, 0, 1, 2, 3]
ys = [x * x for x in xs]
print(ys, round(correlation(xs, ys), 4))

xs2 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 30]
ys2 = [5, 3, 6, 2, 4, 5, 3, 6, 4, 30]
print(round(correlation(xs2[:-1], ys2[:-1]), 2), round(correlation(xs2, ys2), 2))

[9, 4, 1, 0, 1, 4, 9] 0.0
0.07 0.94

y is completely determined by x in the first example, yet r is 0, because the left half falls and the right half rises and the two cancel.

One point can make a correlation. In the second example nine points show essentially no relationship (r ≈ 0.07), and adding one extreme point (30, 30) pushes r above 0.9. Always look at the data, even as a quick text plot, before trusting r.

Correlation is not causation. Towns with more ice-cream sales have more drownings, because hot weather drives both. A third variable that influences both (here, temperature) is a confounder. The revision data is consistent with "revision raises scores", but also with "students who already understand the topic enjoy revising more". Only an experiment that assigns the revision time can separate the two.

Common mistakes: concluding "no relationship" from r ≈ 0 without looking for curves. Reporting r from a handful of points (with three points, a large r happens easily by chance). Averaging correlations from different groups as if they were ordinary numbers. And treating r = 0.5 as "half as related" as r = 1; the useful scale is r², explained with residuals below.

Practice: What Goes With Hot Chocolate?, Two Series That Just Drift, Big Feet, Good Readers?

The least-squares line

A correlation says the points lie close to a line; regression finds that line, y = a + b·x, so you can predict a score from the hours. Which line? Each student has a residual, the actual score minus the predicted score. The least-squares line is the one that makes the sum of squared residuals as small as possible. (Squares again, for the same reasons as the variance.) Solving that minimisation gives two short formulas, the slope b and the intercept a:

The slope is the covariance divided by the variance of x, and each sum is one generator expression in fit_line below. The second formula says the line always passes through the point of means (\bar{x}, \bar{y}).

def fit_line(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    cov = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / n
    var = sum((x - mx) ** 2 for x in xs) / n
    b = cov / var
    return my - b * mx, b

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
a, b = fit_line(hours, scores)
print(f"score = {a:.2f} + {b:.2f} * hours")
for h in (0, 4, 7):
    print(h, "hours ->", round(a + b * h, 1))

score = 47.86 + 4.68 * hours
0 hours -> 47.9
4 hours -> 66.6
7 hours -> 80.6

Read the numbers in words. The slope says: each extra hour of revision goes with about 4.7 more points, on average, in this data. The intercept is the predicted score for zero hours, about 48. Whether that means anything depends on whether zero hours is within the data; here nobody revised less than an hour, so it is a mild extrapolation.

Press Run to see the students and the fitted line together:

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
a, b = 47.86, 4.68
plot(hours, scores, kind="scatter", label="students", xlabel="hours of revision", ylabel="score")
plot([0, 9], [a, a + 9 * b], label="least-squares line")

You can check that the formula really gives the minimum by nudging the line and watching the squared error \sum_i \big(y_i - (a + b x_i)\big)^2 grow (the function sse is that sum, word for word):

def sse(a, b, xs, ys):
    return sum((y - (a + b * x)) ** 2 for x, y in zip(xs, ys))

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
a, b = 47.8585, 4.6817
print("best:", round(sse(a, b, hours, scores), 2))
for da, db in ((1, 0), (-1, 0), (0, 0.2), (0, -0.2)):
    print("nudged", (da, db), round(sse(a + da, b + db, hours, scores), 2))

best: 27.94
nudged (1, 0) 35.94
nudged (-1, 0) 35.93
nudged (0, 0.2) 34.3
nudged (0, -0.2) 34.29

Every nudge makes it worse. In the Machine Learning path you will find the same line a different way: start with any a and b and repeatedly step downhill on this squared error (gradient descent). The formula is the exact answer that gradient descent approaches.

Two cautions. The line predicting y from x is not the same as the line predicting x from y (each minimises errors in a different direction), so decide which variable you want to predict. And a line is a summary of this data: predicting 20 hours of revision gives about 141 points on a test out of 100. Straight lines rarely stay straight beyond the data they were fitted to.

Common mistakes: using the n - 1 covariance with the n variance (the factors must match; then they cancel). Extrapolating far beyond the data. Reading the slope as a causal effect ("revising one more hour will add 4.7 points"); it describes a pattern among these students. And fitting a line to data that is clearly curved, which the residuals below will reveal.

Practice: Minutes per Kilometre, Is the Line Better Than a Guess?, The Flat That Moves the Line

Residuals and R²

The residuals are what the line fails to explain. They are worth printing every time you fit a line, for three reasons.

First, they add up to zero (up to rounding): the least-squares line is a balance point in the same sense as the mean. Second, their size says how good the predictions are. Third, and most useful, their pattern says whether a straight line was the right shape at all.

def fit_line(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    b = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    return my - b * mx, b

hours  = [1, 2, 2, 3, 4, 5, 6, 8]
scores = [52, 55, 61, 60, 68, 70, 77, 85]
a, b = fit_line(hours, scores)
res = [y - (a + b * x) for x, y in zip(hours, scores)]
print([round(r, 2) for r in res])
print("residuals sum to zero:", abs(sum(res)) < 1e-9)

my = sum(scores) / len(scores)
sse = sum(r * r for r in res)
sst = sum((y - my) ** 2 for y in scores)
print("R^2 =", round(1 - sse / sst, 4))

[-0.54, -2.22, 3.78, -1.9, 1.41, -1.27, 1.05, -0.31]
residuals sum to zero: True
R^2 = 0.9683

(The fit divides the two sums directly: the / n in the covariance and the variance cancel.) The residuals go up and down with no pattern, and none is larger than about 4 points. That is what a good straight-line fit looks like.

The last number is R², the fraction of the variation in the scores that the line accounts for:

where \hat{y}_i = a + b x_i (read "y hat") is the line's prediction for student i, so the top sum is sse and the bottom one sst. sst is the total squared variation of the scores around their mean (what you would get predicting the mean for everyone); sse is what remains after using the line. Here the line removes about 97% of it. For a straight line, R² is exactly r²: 0.984² ≈ 0.968. That is why r = 0.5 is weaker than it sounds: it explains only a quarter of the variation.

Now fit a line to curved data and look at the residuals' signs:

def fit_line(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    b = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    return my - b * mx, b

xs = [0, 1, 2, 3, 4, 5, 6]
ys = [x * x for x in xs]
a, b = fit_line(xs, ys)
res = [y - (a + b * x) for x, y in zip(xs, ys)]
print([round(r, 1) for r in res])
print(" ".join("+" if r > 0.5 else "-" if r < -0.5 else "0" for r in res))

[5.0, 0.0, -3.0, -4.0, -3.0, 0.0, 5.0]
+ 0 - - - 0 +

Positive at both ends and negative in the middle: the line is too high in the middle and too low at the edges, which is the signature of a curve. The correlation here is high (about 0.96), so r alone would never have told you the line is the wrong shape. A pattern in the residuals is the signal to try something else, such as fitting against x² (multiple regression in Level 6 makes that routine).

Common mistakes: judging a fit by r or R² alone without looking at the residuals. Comparing R² between datasets with very different spreads (a narrow range of x makes R² small even for a perfect trend with noise). Forgetting that residuals are actual minus predicted; with the sign reversed, "above the line" and "below the line" swap.

Practice: Is a Straight Line the Right Shape?, Straighten the Curve

Conditional probability from a table

Now switch from numbers to events. A school records, for 200 students, whether they handed in the homework and whether they passed the test. The natural summary is a two-way table of counts:

The conditional probability of A given B, written P(A | B), is the probability of A when you only look at the cases where B happened. On a table that is literally what you do: restrict to one row (or column) and divide within it.

Build the 200 records and compute it by filtering:

records = ([("hw", "pass")] * 90 + [("hw", "fail")] * 30 +
           [("none", "pass")] * 40 + [("none", "fail")] * 40)

def prob(event, given=lambda r: True):
    pool = [r for r in records if given(r)]
    return sum(1 for r in pool if event(r)) / len(pool)

passed = lambda r: r[1] == "pass"
did_hw = lambda r: r[0] == "hw"
no_hw  = lambda r: r[0] == "none"
print("P(pass)             =", prob(passed))
print("P(pass | homework)  =", prob(passed, did_hw))
print("P(pass | none)      =", prob(passed, no_hw))
print("P(homework | pass)  =", round(prob(did_hw, passed), 3))

P(pass)             = 0.65
P(pass | homework)  = 0.75
P(pass | none)      = 0.5
P(homework | pass)  = 0.692

Overall 65% passed; among students who did the homework, 75%; among those who did not, 50%. Conditioning changed the probability, so the two events are related.

The last line is a different question with a different answer. P(pass | homework) is 0.75, but P(homework | pass) is about 0.69: the first divides by the 120 students who did the homework, the second by the 130 who passed. Confusing P(A | B) with P(B | A) is the most common error in probability, and it has sent people to prison ("the chance of this evidence if he is innocent is tiny" is not "the chance he is innocent given this evidence"). Bayes' rule, below, is the correct way to turn one into the other.

Conditioning can also reverse a comparison. In the made-up hospital data below, treatment A has the better success rate for small cases and for large cases, yet the worse rate overall:

data = {("A", "small"): (81, 87),  ("A", "large"): (192, 263),
        ("B", "small"): (234, 270), ("B", "large"): (55, 80)}
for t in "AB":
    for size in ("small", "large"):
        ok, n = data[(t, size)]
        print(t, size, f"{ok / n:.0%}")
    ok = sum(data[(t, s)][0] for s in ("small", "large"))
    n = sum(data[(t, s)][1] for s in ("small", "large"))
    print(t, "overall", f"{ok / n:.0%}")

A small 93%
A large 73%
A overall 78%
B small 87%
B large 69%
B overall 83%

This is Simpson's paradox. Treatment A was mostly given to large, harder cases, and B mostly to small, easy ones, so the overall figures compare the case mix rather than the treatments. Case size is a confounder, and the fair comparison is within each size group.

Common mistakes: dividing by the grand total instead of the row or column you conditioned on. Swapping P(A | B) and P(B | A). Trusting a combined table when a variable like case size affects both who gets which treatment and how it turns out.

Practice: Which Way Round?, Given at Least One Six

Independence

Two events are independent when knowing that one happened does not change the probability of the other: P(A | B) = P(A). Equivalently, P(A and B) = P(A) · P(B), which is the multiplication rule you used for separate dice rolls in Level 2. Independence is a property you check, not one you assume because it sounds plausible.

With two dice, is "the total is 7" independent of "the first die is even"? And "the total is 8"? Count:

from itertools import product
from fractions import Fraction

rolls = list(product(range(1, 7), repeat=2))

def p(event, given=lambda r: True):
    pool = [r for r in rolls if given(r)]
    return Fraction(sum(1 for r in pool if event(r)), len(pool))

first_even = lambda r: r[0] % 2 == 0
for target in (7, 8):
    total = lambda r: r[0] + r[1] == target
    print(target, "P(total) =", p(total), " P(total | first even) =", p(total, first_even))

7 P(total) = 1/6  P(total | first even) = 1/6
8 P(total) = 5/36  P(total | first even) = 1/6

Total 7 is independent of the first die being even: whatever the first die shows, exactly one value of the second die makes 7, so knowing the first die tells you nothing. Total 8 is not: an even first die (2, 4 or 6) can always be completed to 8, while a first die of 1 never can. Independence is about the numbers, not about whether the events "seem connected".

With real data you estimate the probabilities, and they will never match exactly. Here is a seeded simulation of 6000 pairs of rolls:

import random
rng = random.Random(11)
pairs = [(rng.randint(1, 6), rng.randint(1, 6)) for _ in range(6000)]
even = [r for r in pairs if r[0] % 2 == 0]
for target in (7, 8):
    overall = sum(a + b == target for a, b in pairs) / len(pairs)
    cond = sum(a + b == target for a, b in even) / len(even)
    print(target, round(overall, 3), round(cond, 3))

7 0.159 0.158
8 0.138 0.16

For total 7 the two estimates are close but not equal: that gap is sampling noise. For total 8 the gap is larger, in line with the exact 5/36 ≈ 0.139 against 1/6 ≈ 0.167. Deciding whether an observed gap is bigger than noise can explain is what the hypothesis tests of Level 5 do (the chi-square test works on exactly this kind of table).

Independence is also a modelling assumption with consequences. Multiplying probabilities as if events were independent, when they are not, gives numbers that are wildly too small: two smoke detectors on the same faulty circuit do not fail independently. The Naive Bayes classifier in the Machine Learning path deliberately assumes that words in a message are independent given the class; the assumption is false, and the method still works well, but it is worth knowing that it is an assumption.

Common mistakes: confusing independent with disjoint. Disjoint events (cannot happen together) are strongly dependent: if one happens, the other's probability drops to zero. Assuming independence for repeated measurements of the same person or the same day. Concluding dependence from a small gap in estimated probabilities.

Practice: If They Were Independent, Does Knowing One Tell You the Other?

Bayes' rule on a table

A disease affects 1% of people. A test for it detects 90% of people who have it (its sensitivity) and wrongly flags 5% of people who do not (its false positive rate). You test positive. How likely is it that you have the disease?

Most people answer about 90%. The honest way to find out is to imagine 10,000 people and fill in the table of counts:

population = 10000
sick = population * 1 // 100
healthy = population - sick
true_pos = sick * 90 // 100
false_pos = healthy * 5 // 100
print("sick:", sick, "of whom test positive:", true_pos)
print("healthy:", healthy, "of whom test positive:", false_pos)
positives = true_pos + false_pos
print("all positives:", positives)
print("P(sick | positive) =", round(true_pos / positives, 3))

sick: 100 of whom test positive: 90
healthy: 9900 of whom test positive: 495
all positives: 585
P(sick | positive) = 0.154

Only about 15% of positives are sick. The test is good, but the disease is rare: 9,900 healthy people times a 5% error rate produce 495 false alarms, which swamp the 90 true cases. The 1% you started from is the prior probability (the base rate), and ignoring it is called base-rate neglect.

The table is Bayes' rule. Written as a formula, it is:

The second line adds up over the two kinds of people. With more possible causes c it becomes \sum_c P(\text{pos} \mid c)\, P(c), which in Python is one generator expression over the causes.

def bayes(prior, sensitivity, false_pos_rate):
    p_pos = sensitivity * prior + false_pos_rate * (1 - prior)
    return sensitivity * prior / p_pos

print(round(bayes(0.01, 0.90, 0.05), 3))
print(round(bayes(0.20, 0.90, 0.05), 3))    # a patient with symptoms: higher prior
print(round(bayes(bayes(0.01, 0.90, 0.05), 0.90, 0.05), 3))   # a second positive test

0.154
0.818
0.766

The same test means very different things for different people. For someone with symptoms (a prior of 20%) a positive result makes the disease about 82% likely. And a second, independent positive test uses the first result as its new prior and raises 15% to about 77%. Updating a prior with evidence, again and again, is the core of Bayesian statistics, which Levels 5 and 6 develop.

And, as always on this path, check by simulation:

import random
rng = random.Random(5)
pos = sick_and_pos = 0
for _ in range(200000):
    sick = rng.random() < 0.01
    positive = rng.random() < (0.90 if sick else 0.05)
    if positive:
        pos += 1
        sick_and_pos += sick
print(pos, sick_and_pos, round(sick_and_pos / pos, 3))

11804 1826 0.155

Naive Bayes in the Machine Learning path is this calculation with more evidence: the prior is how common each class is (spam or not), and the evidence is each word in a message, each with its own P(word | class) estimated from counts in a table.

Common mistakes: reporting the sensitivity P(pos | sick) when the question asks for P(sick | pos). Leaving out the base rate. Using a prior from the wrong population (a hospital's patients are not the general public). And treating repeated tests as independent evidence when they share the same cause of error.

Practice: Whose Cracked Egg?, Listening for Leaks, Was It Really a Hoopoe?

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• Which Way Round?
• If They Were Independent
• Listening for Leaks
• Wind and Kites
• What Goes With Hot Chocolate?
• Is a Straight Line the Right Shape?

Checklist

You can now:

• Compute the covariance of paired data by hand and interpret its sign.
• Compute Pearson's correlation, explain why it has no units, and show how a curve or a single point can mislead it.
• Explain confounding and why correlation does not establish cause and effect.
• Fit a least-squares line from the covariance and variance, interpret its slope and intercept in words, and avoid extrapolating.
• Compute residuals and R², and read a pattern in the residuals as a sign that a line is the wrong shape.
• Compute conditional probabilities from a two-way table, tell P(A | B) from P(B | A), and recognise Simpson's paradox.
• Check independence exactly by counting and approximately from simulated or real data, and tell independent from disjoint.
• Apply Bayes' rule with a table of natural frequencies and with the formula, and explain base-rate neglect.