Level 6 · Statistics & Probability: Modelling

The earlier levels used ready-made recipes: the mean estimates the centre, n - 1 fixes the variance, a t-test compares two groups. This level asks where such recipes come from, and how to build your own when none exists. The answer is a model: a description of how the data could have been generated, with some unknown parameters. Once you have a model, one general principle, maximum likelihood, tells you how to estimate its parameters, and Bayes' rule tells you how to combine them with prior knowledge. Most of the formulas you have used turn out to be special cases.

You will fit models by maximum likelihood, measure how good an estimator is with its bias and variance, use conjugate priors to do Bayesian updating in closed form, model sequences with Markov chains, discover hidden groups with a mixture model fitted by the EM algorithm, and fit regressions with several inputs by solving the normal equations. As always, everything is computed by hand in Python, and simulations check that the methods do what they claim.

This is the level where the Statistics and Machine Learning paths meet. The bias–variance trade-off is the heart of ML's generalisation level; MAP estimation is regularisation (ridge regression) in disguise; EM for a mixture is a soft version of k-means clustering; and multiple regression is the linear model that ML fits by gradient descent. Here you will see the statistical reasons behind those methods.

Maximum likelihood

A model gives the probability of the data for each value of its parameters. Read the other way round, as a function of the parameter with the data held fixed, that probability is the likelihood. The maximum likelihood estimate (MLE) is the parameter value that makes the observed data most probable.

A player makes 7 of 10 free throws. Model: each shot succeeds independently with probability θ. The likelihood is L(\theta) = \theta^7 (1 - \theta)^3. In practice you work with its logarithm, the log-likelihood, because a product of many probabilities becomes a sum of logs, and because the product itself is too small for floating point. For independent observations x_1, \dots, x_n:

In Python the product is math.prod(...) or a loop that multiplies, and the log-likelihood is sum(math.log(...) for ...). The code below first finds the peak on a grid, then computes both for 2000 coin flips:

import math
import random

k, n = 7, 10
grid = [i / 1000 for i in range(1, 1000)]
loglik = [k * math.log(t) + (n - k) * math.log(1 - t) for t in grid]
print("MLE on a grid:", grid[loglik.index(max(loglik))])

rng = random.Random(90)
flips = [rng.random() < 0.3 for _ in range(2000)]
product = 1.0
for f in flips:
    product *= 0.3 if f else 0.7
print("likelihood as a product:", product)
print("log-likelihood:", round(sum(math.log(0.3 if f else 0.7) for f in flips), 1))

MLE on a grid: 0.7
likelihood as a product: 0.0
log-likelihood: -1203.9

The grid search finds θ = 0.7, the observed proportion, as you would hope. Calculus confirms it: setting the derivative of k \log\theta + (n - k)\log(1 - \theta) to zero gives \hat{\theta} = k/n (the hat marks an estimate). Press Run to see the log-likelihood curve and its peak:

import math

grid = [i / 1000 for i in range(1, 1000)]
plot(grid, [7 * math.log(t) + 3 * math.log(1 - t) for t in grid], xlabel="θ", ylabel="log-likelihood of 7 out of 10")

The coin flips in the first program show why logs are not optional: the likelihood of 2000 coin flips is about 10^{-523}, far below the smallest positive float, so the product underflows to exactly 0.0, while its logarithm is an ordinary number.

The same principle recovers other familiar estimators. For counts modelled as Poisson with rate λ, the MLE is the sample mean; for normal data, it is the sample mean and the variance divided by n. For counts c_1, \dots, c_n the Poisson log-likelihood is the sum of the logs of the Poisson probabilities from Level 4,

and poisson_loglik below is that sum, term for term:

import math

counts = [2, 4, 3, 0, 5, 3, 2, 4, 1, 3]

def poisson_loglik(lam):
    return sum(c * math.log(lam) - lam - math.lgamma(c + 1) for c in counts)

grid = [i / 100 for i in range(1, 1001)]
values = [poisson_loglik(g) for g in grid]
print("Poisson MLE:", grid[values.index(max(values))], " sample mean:", sum(counts) / len(counts))

Poisson MLE: 2.7  sample mean: 2.7

(math.lgamma(c + 1) is log(c!), which does not depend on λ but keeps the value a true log-likelihood.) So maximum likelihood explains why the mean is the right estimate of a Poisson rate. It also raises a question: for the normal variance it gives the ÷ n version, which Level 2 showed to be biased. Is the MLE then wrong? The next section answers that.

Maximum likelihood works for any model you can write a probability for, including ones with no textbook formula; then you maximise the log-likelihood numerically (a grid, or gradient ascent as in the Machine Learning path). Logistic regression in that path is exactly this: the "log loss" it minimises is the negative log-likelihood of a Bernoulli model.

Common mistakes: multiplying probabilities instead of adding logs. Taking log(0) at the edges of a grid (start the grid just above 0 and end just below 1). Maximising over a grid too coarse to find the peak. And forgetting that the MLE is only as good as the model: a Poisson fit to bursty data gives a correct mean and a badly wrong picture of the extremes.

Practice: The Likeliest Success Rate, Cyclists per Minute, Lamps Still Burning

Bias and variance of estimators

An estimator is a random quantity: it depends on which sample you happened to draw. Two numbers describe how good it is. Its bias is how far its average (over all possible samples) is from the truth. Its variance is how much it scatters around that average. They combine into the mean squared error, the average squared distance from the truth:

where \hat{\theta} is the estimate and \theta the true value.

A small bias is not automatically good, and neither is zero bias. Compare three ways to estimate the variance σ² = 100 of a normal population from samples of size 5: divide the sum of squared deviations by n - 1 (unbiased), by n (the MLE), or by n + 1.

import random

rng = random.Random(91)
n, trials, true_var = 5, 40000, 100.0
estimates = {n - 1: [], n: [], n + 1: []}
for _ in range(trials):
    s = [rng.gauss(0, 10) for _ in range(n)]
    m = sum(s) / n
    ss = sum((x - m) ** 2 for x in s)
    for d in estimates:
        estimates[d].append(ss / d)

print("divide by   bias  variance   MSE")
for d, v in estimates.items():
    avg = sum(v) / trials
    var = sum((x - avg) ** 2 for x in v) / trials
    mse = sum((x - true_var) ** 2 for x in v) / trials
    print(f"{d:9d} {avg - true_var:7.2f} {var:9.0f} {mse:6.0f}")

divide by   bias  variance   MSE
        4   -0.04      5046   5046
        5  -20.03      3229   3631
        6  -33.36      2243   3355

Dividing by n - 1 has no bias, but the most variance. Dividing by n is biased low by about 20, yet its errors are smaller on average. Dividing by n + 1 is even more biased and has the smallest mean squared error of the three. Shrinking an estimate towards zero trades a little bias for a larger cut in variance.

That is the bias–variance trade-off, and it is everywhere. In the Machine Learning path it explains overfitting: a very flexible model has low bias but high variance (it changes a lot with the training sample), and regularisation deliberately adds bias to reduce variance. In statistics it explains why the "unbiased" estimator is a convention rather than a law: it is the right choice when you will average many estimates (their biases would add up), and not necessarily when you need one estimate that is as close as possible.

Common mistakes: treating "unbiased" as "best". Judging an estimator from one sample instead of from its behaviour over many. And forgetting that the trade-off depends on n: with 500 values the three estimators above are practically identical.

Practice: Report Card for a Rate Rule

MAP estimates and conjugate priors

Level 5 updated a prior on a grid. For many models there is a shortcut. The Beta(a, b) distribution on θ between 0 and 1 has density proportional to \theta^{a-1}(1 - \theta)^{b-1}. Multiply it by a binomial likelihood \theta^k (1 - \theta)^{n-k} and the exponents simply add: the posterior is again a Beta distribution,

A prior whose posterior belongs to the same family is called conjugate. The update is two additions, and a and b read as "pseudo-counts": a Beta(11, 11) prior behaves like having already seen 10 successes and 10 failures (the "league" prior of Level 5). Two summaries of a Beta(a, b) posterior are standard: its mean \frac{a}{a + b}, and its peak, the maximum a posteriori (MAP) estimate \frac{a - 1}{a + b - 2}.

def update(a, b, k, n):
    return a + k, b + n - k

def describe(a, b):
    return f"Beta({a}, {b}): mean {a / (a + b):.3f}, MAP {(a - 1) / (a + b - 2):.3f}"

print(describe(*update(1, 1, 7, 10)))      # flat prior
print(describe(*update(11, 11, 7, 10)))    # league prior, 10 shots
print(describe(*update(11, 11, 70, 100)))  # league prior, 100 shots

Beta(8, 4): mean 0.667, MAP 0.700
Beta(18, 14): mean 0.562, MAP 0.567
Beta(81, 41): mean 0.664, MAP 0.667

These match the grid of Level 5 (posterior means 0.667, 0.562 and 0.664) without any grid. With a flat prior, Beta(1, 1), the MAP estimate equals the MLE 0.7: MAP is maximum likelihood with the log-prior added to the log-likelihood, so a flat prior adds nothing. An informative prior pulls the estimate towards its centre, strongly when data is scarce and hardly at all when it is plentiful.

That pull is useful. A spam filter has seen 200 spam messages, and none contained the word "invoice". The MLE of P("invoice" | spam) is 0/200 = 0, and a Naive Bayes classifier multiplying probabilities would then declare any message containing "invoice" impossible as spam, however spammy the rest of it. The posterior mean with a flat prior gives \frac{k + 1}{n + 2} instead: small, but not zero.

k, n = 0, 200
print("MLE:", k / n)
print("posterior mean, Beta(1, 1) prior:", round((k + 1) / (n + 2), 4))

MLE: 0.0
posterior mean, Beta(1, 1) prior: 0.005

This is Laplace (add-one) smoothing, used in the Naive Bayes classifier of the Machine Learning path; now you know it is a Bayesian estimate. In the same way, ridge regression in that path is the MAP estimate of a regression with a normal prior centred on zero for each coefficient: the penalty on large coefficients is the log of the prior.

Conjugate pairs exist for other models too: a Gamma prior for a Poisson rate, a normal prior for a normal mean with known variance. When no conjugate prior fits, the grid (for one or two parameters) or simulation methods take over.

Common mistakes: confusing the MAP (the peak) with the posterior mean; they differ for skewed posteriors. Using the MAP formula when a + b - 2 is zero or when a or b is below 1 (the peak is then at an edge). And forgetting that a strong prior is a strong assumption; report which prior you used.

Practice: Where the Posterior Peaks, A Prior Is a Penalty

Markov chains

Every model so far assumed independent observations. Many sequences are not like that: tomorrow's weather depends on today's, the next word depends on the last one, a board-game token's next square depends on where it is. A Markov chain is the simplest model with memory. The system is in one of a few states, and the probability of the next state depends only on the current one, not on the path that led there (the Markov property). The model is a table of transition probabilities, each row adding up to 1:

import random

P = {"sun":   {"sun": 0.7, "cloud": 0.2, "rain": 0.1},
     "cloud": {"sun": 0.3, "cloud": 0.4, "rain": 0.3},
     "rain":  {"sun": 0.2, "cloud": 0.4, "rain": 0.4}}

rng = random.Random(100)
state = "rain"
week = []
for _ in range(14):
    options = list(P[state])
    state = rng.choices(options, weights=[P[state][o] for o in options])[0]
    week.append(state)
print(" ".join(week))

sun sun cloud rain rain cloud rain cloud sun sun sun rain rain cloud

Simulated weather comes in spells, because each state tends to repeat itself. What is the chance of sun in k days' time, if it is raining today? Keep a probability distribution over the states and push it through the table one day at a time: the probability of each state tomorrow is the sum, over today's states, of P(today) · P(today → tomorrow). With \pi_t(i) the probability of state i on day t and P_{ij} the probability of moving from i to j:

In step below the sum over i is the generator expression, and the dict comprehension around it does this once for every state j.

P = {"sun":   {"sun": 0.7, "cloud": 0.2, "rain": 0.1},
     "cloud": {"sun": 0.3, "cloud": 0.4, "rain": 0.3},
     "rain":  {"sun": 0.2, "cloud": 0.4, "rain": 0.4}}
states = list(P)

def step(dist):
    return {t: sum(dist[s] * P[s][t] for s in states) for t in states}

dist = {"sun": 0.0, "cloud": 0.0, "rain": 1.0}
for day in range(1, 11):
    dist = step(dist)
    if day in (1, 2, 3, 5, 10):
        print(day, {s: round(p, 3) for s, p in dist.items()})

1 {'sun': 0.2, 'cloud': 0.4, 'rain': 0.4}
2 {'sun': 0.34, 'cloud': 0.36, 'rain': 0.3}
3 {'sun': 0.406, 'cloud': 0.332, 'rain': 0.262}
5 {'sun': 0.45, 'cloud': 0.313, 'rain': 0.237}
10 {'sun': 0.461, 'cloud': 0.308, 'rain': 0.231}

Press Run to watch the three probabilities settle:

P = {"sun":   {"sun": 0.7, "cloud": 0.2, "rain": 0.1},
     "cloud": {"sun": 0.3, "cloud": 0.4, "rain": 0.3},
     "rain":  {"sun": 0.2, "cloud": 0.4, "rain": 0.4}}
history = [{"sun": 0.0, "cloud": 0.0, "rain": 1.0}]
for day in range(15):
    history.append({t: sum(history[-1][s] * P[s][t] for s in P) for t in P})
for s in P:
    plot([h[s] for h in history], label=s, xlabel="day", ylabel="probability")

After a week or so the starting weather is forgotten: the distribution settles at about 46% sun, 31% cloud and 23% rain, and starting from sun instead would lead to the same place. That limit is the stationary distribution: the distribution that one more step leaves unchanged, \pi(j) = \sum_i \pi(i)\, P_{ij} for every state j. For this chain it is exactly 6/13, 4/13 and 3/13, as exact arithmetic confirms:

from fractions import Fraction as F

P = {"sun":   {"sun": F(7, 10), "cloud": F(2, 10), "rain": F(1, 10)},
     "cloud": {"sun": F(3, 10), "cloud": F(4, 10), "rain": F(3, 10)},
     "rain":  {"sun": F(2, 10), "cloud": F(4, 10), "rain": F(4, 10)}}
pi = {"sun": F(6, 13), "cloud": F(4, 13), "rain": F(3, 13)}
after = {t: sum(pi[s] * P[s][t] for s in P) for t in P}
print(after == pi, after)

True {'sun': Fraction(6, 13), 'cloud': Fraction(4, 13), 'rain': Fraction(3, 13)}

The stationary distribution is also the long-run fraction of time in each state. Count the states over a long simulation:

import random

P = {"sun":   {"sun": 0.7, "cloud": 0.2, "rain": 0.1},
     "cloud": {"sun": 0.3, "cloud": 0.4, "rain": 0.3},
     "rain":  {"sun": 0.2, "cloud": 0.4, "rain": 0.4}}
rng = random.Random(101)
state = "rain"
counts = {s: 0 for s in P}
days = 100000
for _ in range(days):
    options = list(P[state])
    state = rng.choices(options, weights=[P[state][o] for o in options])[0]
    counts[state] += 1
print({s: round(c / days, 4) for s, c in counts.items()})

{'sun': 0.4622, 'cloud': 0.3062, 'rain': 0.2316}

This is the law of large numbers again, now for dependent observations. The idea behind Google's original PageRank was a Markov chain of a random web surfer, whose stationary distribution ranks pages. And Markov chain Monte Carlo methods reverse the logic: to sample from a posterior distribution too complicated for a grid, they design a chain whose stationary distribution is that posterior, then simulate it.

Estimating the transition table from data is maximum likelihood again: the MLE of P(sun → rain) is the number of sun-to-rain transitions divided by the number of sunny days that had a next day.

Common mistakes: rows that do not add up to 1. Assuming the Markov property when the process has longer memory (a three-day heatwave may be more likely to continue than a one-day one); the fix is to put more history into the state. And reading the stationary distribution as the forecast for tomorrow; it is the forecast for the distant future.

Practice: Coffee Brand Switching, Learning the Table from a Log, Where the Robot Spends Its Time

Mixture models and the EM algorithm

Some datasets are two populations mixed together without labels. Eruption lengths of a geyser, for example, come in short and long kinds; the heights of a mixed group of children and adults have two humps. A mixture model says: each value was produced by first picking a hidden group at random (group 1 with probability w, group 2 with probability 1 - w) and then drawing from that group's own normal distribution. Here are 300 simulated eruption lengths in minutes:

import random

rng = random.Random(110)
data = [rng.gauss(2.0, 0.4) for _ in range(120)] + [rng.gauss(4.3, 0.5) for _ in range(180)]
rng.shuffle(data)
counts = {}
for x in data:
    b = int(x * 2) / 2                 # bins half a minute wide
    counts[b] = counts.get(b, 0) + 1
for b in sorted(counts):
    print(f"{b:4.1f} {'#' * (counts[b] // 2)}")

 1.0 ####
 1.5 ######################
 2.0 ############################
 2.5 #####
 3.0 ########
 3.5 ##############
 4.0 #######################################
 4.5 ########################
 5.0 ####

Two humps are visible, but no value says which group it came from. If you knew the groups, fitting each normal would be easy (the mean and ÷ n standard deviation of each group, by maximum likelihood). If you knew the normals, you could compute for each value the probability that it came from each group (Bayes' rule). The EM algorithm (expectation–maximisation) solves the chicken-and-egg problem by alternating the two:

• E step: with the current parameters, compute each value's responsibility r, the posterior probability that it belongs to group 1. With f_1 and f_2 the two normal densities, Bayes' rule gives r_i = \frac{w f_1(x_i)}{w f_1(x_i) + (1 - w) f_2(x_i)}.
• M step: refit w, the means and the standard deviations, counting each value as r of a group-1 value and 1 - r of a group-2 value. The new group-1 mean, for example, is the weighted mean \mu_1 = \frac{\sum_i r_i x_i}{\sum_i r_i}, which the code writes as sum(r * x for r, x in zip(resp, data)) / n1.

Each round is guaranteed not to decrease the log-likelihood, so the fit improves until it settles.

import math
import random

rng = random.Random(110)
data = [rng.gauss(2.0, 0.4) for _ in range(120)] + [rng.gauss(4.3, 0.5) for _ in range(180)]
rng.shuffle(data)

def normal_pdf(x, mu, sd):
    return math.exp(-((x - mu) / sd) ** 2 / 2) / (sd * math.sqrt(2 * math.pi))

w, mu1, sd1, mu2, sd2 = 0.5, 1.0, 1.0, 5.0, 1.0      # a rough first guess
print("round  loglik     w    mu1   sd1   mu2   sd2")
for rnd in range(1, 21):
    resp, loglik = [], 0.0
    for x in data:                                    # E step
        a = w * normal_pdf(x, mu1, sd1)
        b = (1 - w) * normal_pdf(x, mu2, sd2)
        resp.append(a / (a + b))
        loglik += math.log(a + b)
    n1 = sum(resp)                                    # M step
    n2 = len(data) - n1
    mu1 = sum(r * x for r, x in zip(resp, data)) / n1
    mu2 = sum((1 - r) * x for r, x in zip(resp, data)) / n2
    sd1 = math.sqrt(sum(r * (x - mu1) ** 2 for r, x in zip(resp, data)) / n1)
    sd2 = math.sqrt(sum((1 - r) * (x - mu2) ** 2 for r, x in zip(resp, data)) / n2)
    w = n1 / len(data)
    if rnd in (1, 2, 3, 5, 10, 20):
        print(f"{rnd:5d} {loglik:8.2f} {w:5.2f} {mu1:6.2f} {sd1:5.2f} {mu2:5.2f} {sd2:5.2f}")

for x in (2.0, 2.8, 2.9, 3.0, 4.5):
    a = w * normal_pdf(x, mu1, sd1)
    b = (1 - w) * normal_pdf(x, mu2, sd2)
    print(f"P(short group | {x}) = {a / (a + b):.3f}")

round  loglik     w    mu1   sd1   mu2   sd2
    1  -611.76  0.40   2.07  0.47  4.24  0.57
    2  -387.55  0.40   2.04  0.39  4.26  0.51
    3  -379.54  0.40   2.02  0.37  4.26  0.51
    5  -378.76  0.39   2.01  0.36  4.25  0.52
   10  -378.64  0.39   2.01  0.35  4.25  0.52
   20  -378.64  0.39   2.01  0.35  4.25  0.52
P(short group | 2.0) = 1.000
P(short group | 2.8) = 0.773
P(short group | 2.9) = 0.507
P(short group | 3.0) = 0.229
P(short group | 4.5) = 0.000

(Each printed log-likelihood belongs to the parameters at the start of that round.) The log-likelihood rises and levels off within a few rounds, and the fitted mixture is close to the simulated truth: 40% short eruptions around 2.0 minutes with standard deviation 0.4, and 60% long ones around 4.3 with 0.5. The last lines show what a mixture offers over a hard split: a value of 2.0 is almost certainly short, 4.5 almost certainly long, and between 2.8 and 3.0 the probability swings from about 3 in 4 to about 1 in 4: those values are genuinely ambiguous, and the model says so.

The k-means clustering of the Machine Learning path is EM with hard assignments: every point is given entirely to its nearest centre (responsibility 0 or 1) and the spreads are ignored. Both share EM's weaknesses. They can converge to a local optimum: a bad start (both means on the same hump) can settle on a poor fit, so run from several starts and keep the one with the highest log-likelihood. And the number of groups must be chosen; adding groups always raises the likelihood, so compare fits on held-out data or with a penalty for extra parameters.

Common mistakes: a standard deviation collapsing to 0 when one group captures a single point (the likelihood then goes to infinity; put a floor under the standard deviations). Starting both groups with identical parameters (they then stay identical forever). And reading the groups as real kinds of thing: a mixture of two normals also fits one skewed population well, and that alone is no evidence of two causes.

Practice: Two Kinds of Fish, How Many Operating Modes?

Multiple regression and the normal equations

Level 3 fitted a line with one input. Most outcomes depend on several: a house's price depends on its size and its age. Multiple regression fits

by least squares, choosing the coefficients that minimise the sum of squared residuals. Write each house as a row [1, size, age] of a table X (the leading 1 multiplies the intercept b0) and the prices as a column y. Setting the derivative of the squared error to zero with respect to each coefficient gives a small linear system, the normal equations:

X^\top X is a 3 × 3 table of sums of products of the columns, and X^\top \mathbf{y} is a list of three sums:

where r runs over the rows (the houses). In fit below each entry is one generator expression over the rows, and a list comprehension over i and j builds the table; with one input these equations reduce to the covariance-over-variance formula of Level 3. Solving the system takes Gaussian elimination, which is a few lines:

def solve(A, b):
    """Solve A x = b by Gaussian elimination with partial pivoting."""
    n = len(A)
    M = [row[:] + [bi] for row, bi in zip(A, b)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(M[r][col]))
        M[col], M[pivot] = M[pivot], M[col]
        for r in range(n):
            if r != col:
                f = M[r][col] / M[col][col]
                M[r] = [x - f * y for x, y in zip(M[r], M[col])]
    return [M[i][n] / M[i][i] for i in range(n)]

def fit(X, y):
    k = len(X[0])
    XtX = [[sum(row[i] * row[j] for row in X) for j in range(k)] for i in range(k)]
    Xty = [sum(row[i] * yi for row, yi in zip(X, y)) for i in range(k)]
    return solve(XtX, Xty)

size  = [55, 70, 80, 95, 100, 120, 130, 150]        # square metres
age   = [8, 5, 20, 15, 35, 22, 40, 38]               # years
price = [135, 172, 177, 205, 177, 241, 247, 288]     # thousands of pounds

b0, b1, b2 = fit([[1, s, a] for s, a in zip(size, age)], price)
print(f"size and age: intercept {b0:.1f}, per m2 {b1:.2f}, per year {b2:.2f}")
c0, c1 = fit([[1, s] for s in size], price)
print(f"size only:    intercept {c0:.1f}, per m2 {c1:.2f}")

size and age: intercept 34.4, per m2 2.05, per year -1.51
size only:    intercept 53.8, per m2 1.51

(The prices were simulated from 40 + 2.0 · size - 1.5 · age plus noise, so you can see that the fit recovers the true slopes well.) Read the coefficients in words: among houses of the same age, each extra square metre goes with about £2,050 more; among houses of the same size, each extra year goes with about £1,500 less.

Now compare the size-only fit: it says each square metre is worth only about £1,500. Both fits are correct least-squares answers to different questions. In this data the bigger houses also tend to be older, so the size-only slope mixes the value of size with the discount for age. Multiple regression separates them by "holding the other inputs fixed". It is the regression form of the confounding you met in Level 3: leaving out a variable that is related to both the input and the outcome biases the slope. It is not magic, though: it can only adjust for confounders you measured and included.

The same machinery fits curves. For the curved data of Level 3, add a column x²; the model is still linear in its coefficients, so the normal equations apply unchanged:

def solve(A, b):
    n = len(A)
    M = [row[:] + [bi] for row, bi in zip(A, b)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(M[r][col]))
        M[col], M[pivot] = M[pivot], M[col]
        for r in range(n):
            if r != col:
                f = M[r][col] / M[col][col]
                M[r] = [x - f * y for x, y in zip(M[r], M[col])]
    return [M[i][n] / M[i][i] for i in range(n)]

xs = [0, 1, 2, 3, 4, 5, 6]
ys = [x * x for x in xs]
X = [[1, x, x * x] for x in xs]
XtX = [[sum(r[i] * r[j] for r in X) for j in range(3)] for i in range(3)]
Xty = [sum(r[i] * y for r, y in zip(X, ys)) for i in range(3)]
coef = solve(XtX, Xty)
print([round(c, 6) + 0.0 for c in coef])
print([round(y - sum(c * v for c, v in zip(coef, r)), 6) + 0.0 for r, y in zip(X, ys)])

[0.0, 0.0, 1.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

The fit finds y = 0 + 0·x + 1·x² and the residual pattern of Level 3 disappears. (Adding 0.0 turns a rounded -0.0 into 0.0 for printing.)

In the Machine Learning path the same model is fitted by gradient descent, which avoids forming Xᵀ X and scales to many inputs; the normal equations give the exact answer that descent approaches. Two problems show up there as well. If two inputs are nearly copies of each other (size in m² and size in square feet), Xᵀ X is nearly singular and the coefficients become huge and unstable. Ridge regression adds a small λ to the diagonal of Xᵀ X, which is exactly the MAP estimate with a normal prior from earlier in this level, and makes the system well-behaved. And with as many coefficients as data points the fit is perfect and useless: the bias–variance trade-off again.

Common mistakes: forgetting the column of ones (the fit is then forced through the origin). Reading a coefficient as the effect of changing one input while ignoring that the others are held fixed. Including two nearly identical inputs. Dividing by a zero pivot (partial pivoting, choosing the largest available entry, avoids that when the system is solvable). And fitting many inputs to few data points, which fits the noise.

Practice: Mileage or Age?

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.

• The Likeliest Success Rate
• Cyclists per Minute
• Lamps Still Burning
• Coffee Brand Switching
• Learning the Table from a Log
• Where the Robot Spends Its Time

Checklist

You can now:

• Write down a likelihood, maximise the log-likelihood on a grid, and explain why logs are needed.
• Show that maximum likelihood gives the sample proportion, the Poisson mean and the ÷ n normal variance.
• Measure an estimator's bias, variance and mean squared error by simulation, and explain the bias–variance trade-off.
• Update a Beta prior with binomial data in closed form, compute the posterior mean and the MAP estimate, and derive add-one smoothing.
• Simulate a Markov chain, propagate its distribution step by step, find its stationary distribution, and check it against long-run frequencies.
• Fit a two-component normal mixture with EM, read the responsibilities, and guard against collapse and bad starts.
• Fit a multiple regression by solving the normal equations, interpret each coefficient with the others held fixed, and use extra columns to fit curves.
• Name the links to the Machine Learning path: log loss is a negative log-likelihood, ridge is MAP, k-means is hard EM, and gradient descent approaches the normal-equation solution.