Level 5 · Machine Learning Fundamentals: Generalisation

Level 4 made training precise: choose the parameters that minimise the loss on the training data. This level is about the catch. The training data is only a sample, and a model with enough freedom fits the noise in its sample as faithfully as the pattern, and then predicts badly. Level 2 glimpsed this in 1-NN's perfect training score; here you will measure it, control it and explain it.

You will watch polynomials of rising degree overfit, choose settings honestly with a validation set and k-fold cross-validation, change the loss itself with regularisation, and split error into bias and variance, the decomposition Statistics level 6 uses for estimators. Two new models complete the level: decision trees, whose depth is a complexity dial, and k-means clustering, which learns without labels by minimising a loss of its own.

The optimisation thread continues: regularisation changes the loss being minimised, choosing a hyperparameter is an optimisation one level up, a tree is grown by greedy search, and k-means is a descent method on a clustering loss that gets stuck in local optima just as the non-convex losses of Level 4 did.

Overfitting and underfitting

The data for most of this level comes from a smooth curve, sin(3x) on the interval from -1 to 1, with random noise of standard deviation 0.3 added to every observation:

import math, random

def make_wave(n, seed):
    rng = random.Random(seed)
    xs = [round(rng.uniform(-1, 1), 2) for _ in range(n)]
    ys = [round(math.sin(3 * x) + rng.gauss(0, 0.3), 2) for x in xs]
    return xs, ys

x_train, y_train = make_wave(12, seed=3)
x_new, y_new = make_wave(500, seed=4)            # stands in for future data
print(list(zip(x_train, y_train))[:4])

[(-0.52, -0.82), (0.09, -0.03), (-0.26, -0.81), (0.21, 0.46)]

The models are polynomials of degree d:

\hat{y} = w_0 + w_1 x + w_2 x^2 + \dots + w_d x^d = \sum_{j=0}^{d} w_j x^j

The sum is predict below, sum(wj * x ** j for j, wj in enumerate(w)), with enumerate supplying the power j next to the weight w_j. A polynomial is still a linear model, in the sense of Levels 1 and 3: it is a dot product of the weights with the feature vector [1, x, x^2, \dots, x^d]. So the least-squares weights come from the normal equations (X^\top X)\mathbf{w} = X^\top \mathbf{y}, where X has one row of features per example and X^\top is X turned on its side (the multiple-regression formula of Statistics level 6), which are a small system of linear equations. Gaussian elimination solves it:

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 c in range(n):
        p = max(range(c, n), key=lambda r: abs(M[r][c]))    # largest pivot, for stability
        M[c], M[p] = M[p], M[c]
        for r in range(c + 1, n):
            f = M[r][c] / M[c][c]
            for k in range(c, n + 1):
                M[r][k] -= f * M[c][k]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        x[r] = (M[r][n] - sum(M[r][k] * x[k] for k in range(r + 1, n))) / M[r][r]
    return x

def fit_poly(xs, ys, degree, lam=0.0):
    F = [[x ** j for j in range(degree + 1)] for x in xs]
    A = [[sum(f[i] * f[j] for f in F) + (lam if i == j and i > 0 else 0)
          for j in range(degree + 1)] for i in range(degree + 1)]
    b = [sum(f[i] * y for f, y in zip(F, ys)) for i in range(degree + 1)]
    return solve(A, b)

def predict(w, x):
    return sum(wj * x ** j for j, wj in enumerate(w))

def mse(w, xs, ys):
    return sum((predict(w, x) - y) ** 2 for x, y in zip(xs, ys)) / len(xs)

print(solve([[2, 1], [1, 3]], [3, 5]))          # 2a + b = 3, a + 3b = 5
print([round(v, 3) for v in fit_poly([0, 1, 2], [1, 3, 5], 1)])   # the line y = 1 + 2x

[0.8, 1.4]
[1.0, 2.0]

(The lam parameter is for regularisation, later in this level; ignore it for now.) Now fit every degree from 0 to 9 and compare the error on the 12 training points with the error on 500 new ones:

for d in range(10):
    w = fit_poly(x_train, y_train, d)
    print(d, f"train {mse(w, x_train, y_train):.4f}", f"new {mse(w, x_new, y_new):.4f}")

0 train 0.4060 new 0.7898
1 train 0.1198 new 0.3283
2 train 0.1198 new 0.3283
3 train 0.0167 new 0.1530
4 train 0.0125 new 0.1462
5 train 0.0117 new 0.1392
6 train 0.0084 new 0.1397
7 train 0.0082 new 0.1416
8 train 0.0064 new 2.4404
9 train 0.0044 new 60.1854

Read the two columns separately. The training error falls with every extra degree, because a higher-degree polynomial can do everything a lower one can and more; with enough degrees it would pass through every training point exactly. The error on new data tells a different story. Degrees 0 to 2 are too stiff to follow the wave: they underfit, and are wrong on training and new data alike. Degrees 3 to 7 capture the shape, with new-data error around 0.14, not far above the 0.09 that the noise alone guarantees (0.3^2). From degree 8 the curve bends wildly between and beyond the training points to hit them all, and the error on new data explodes: the model overfits.

Plotting the two columns of the table against the degree shows this best. Press Run to see training error fall while the error on new data turns up again (the vertical axis is logarithmic, because degree 9 is so far off):

degrees = list(range(10))
train = [0.406, 0.1198, 0.1198, 0.0167, 0.0125, 0.0117, 0.0084, 0.0082, 0.0064, 0.0044]
new = [0.7898, 0.3283, 0.3283, 0.153, 0.1462, 0.1392, 0.1397, 0.1416, 2.4404, 60.1854]
plot(degrees, train, label="training error", xlabel="degree", ylabel="MSE", logy=True)
plot(degrees, new, label="error on new data")

The general picture: as a model's flexibility grows, training error keeps falling, while error on new data falls, bottoms out and rises again. The goal is the bottom of that U, and the training error cannot tell you where it is.

Common mistakes: choosing the model with the lowest training error, which always picks the most flexible one. Thinking overfitting needs a complicated model; it depends on flexibility relative to the amount of data, and a degree-8 polynomial that is wild with 12 points behaves well with 1,000. Blaming noise for all the error; compare with the irreducible level when you know it.

Practice: Twelve Points, Ten Curves, Bending the Line

Train, validation and test

The degree of the polynomial is not learned by fit_poly; you choose it before training. Settings like that are called hyperparameters: the degree here, k in kNN, the learning rate in gradient descent, the depth of a tree below. Choosing them is itself a search, one level above training: for each candidate setting, train, then score. The question is what to score on.

Not the training data, as the previous section showed. Not the test data either: if you try ten degrees and keep the one with the best test score, that score is optimistic, because the winner was picked partly for its luck on those test examples. The standard answer is a third set: the validation set compares settings, and the test set is used once, at the very end, to report the chosen model.

x_val, y_val = make_wave(12, seed=5)             # 12 more examples, used only for choosing
scores = {}
for d in range(10):
    w = fit_poly(x_train, y_train, d)
    scores[d] = mse(w, x_val, y_val)
    print(d, round(scores[d], 4))
best = min(scores, key=scores.get)
print("chosen degree:", best)

0 0.9301
1 0.2887
2 0.2891
3 0.1067
4 0.1117
5 0.1072
6 0.12
7 0.1155
8 6.3285
9 168.8768
chosen degree: 3

Validation picks degree 3, from the good range, and firmly rejects degrees 8 and 9. Only now is the test set opened, once:

x_test, y_test = make_wave(12, seed=6)
w = fit_poly(x_train, y_train, best)
print("test MSE of the chosen model:", round(mse(w, x_test, y_test), 4))

test MSE of the chosen model: 0.0946

With 12 test points this figure is itself noisy (here it happens to land slightly below the noise level of 0.09), which is why real test sets are much larger. A common refinement: after choosing the degree, retrain on training and validation data together (more data gives better weights), then test.

In practice you split one dataset, perhaps 60% training, 20% validation and 20% test, shuffled first as in Level 2. The weakness shows here: with 12 validation points the choice rests on a few noisy numbers, and every example set aside is one the model cannot learn from. Cross-validation addresses both.

Common mistakes: tuning on the test set, then reporting the test score as if it were unseen (the most common way that published results turn out too good). Opening the test set "just to check" during development. Preprocessing (such as computing the means for standardisation) on all the data before splitting, which leaks information from the validation and test sets into training.

Practice: Stop While It Still Generalises

k-fold cross-validation

k-fold cross-validation reuses every example for both training and validation. Shuffle the data and cut it into k equal folds. For each fold in turn, train on the other k - 1 folds and validate on that one. The average of the k validation errors is the cross-validation score, and every example has been used for validation exactly once.

Here the 24 examples of the training and validation sets are pooled and scored with 4 folds:

def k_fold_indices(n, k, seed):
    idx = list(range(n))
    random.Random(seed).shuffle(idx)
    return [idx[i::k] for i in range(k)]          # k folds of (nearly) equal size

def cross_val_mse(xs, ys, degree, k=4, lam=0.0, seed=0):
    folds = k_fold_indices(len(xs), k, seed)
    total = 0.0
    for fold in folds:
        held_out = set(fold)
        train = [i for i in range(len(xs)) if i not in held_out]
        w = fit_poly([xs[i] for i in train], [ys[i] for i in train], degree, lam)
        total += mse(w, [xs[i] for i in fold], [ys[i] for i in fold])
    return total / k

x_all, y_all = x_train + x_val, y_train + y_val
print([len(f) for f in k_fold_indices(len(x_all), 4, 0)])
for d in range(10):
    print(d, round(cross_val_mse(x_all, y_all, d), 4))

[6, 6, 6, 6]
0 0.5967
1 0.1963
2 0.1931
3 0.0634
4 0.0691
5 0.0838
6 0.0789
7 0.0828
8 0.2444
9 0.3458

Cross-validation also chooses degree 3, and its curve has the U shape more clearly than a single validation split, because each score averages four different splits. The price is k trainings per candidate instead of one. Five or ten folds are common choices; with k equal to the number of examples it is called leave-one-out cross-validation, which is cheap for kNN (no training) and expensive for most other models.

Cross-validation estimates how well a procedure ("fit a degree-3 polynomial") performs on data like this. Once the procedure is chosen, fit it once on all the available training data to get the final model.

Common mistakes: shuffling differently for each candidate setting, so the settings are compared on different splits (fix the seed, as cross_val_mse does). Doing feature scaling or feature selection on all the data before the folds are cut. Treating the cross-validation score of the winning setting as an unbiased estimate; like any score used for choosing, it is slightly optimistic, and the untouched test set is still needed for the final report.

Practice: Every Day Takes a Turn, Four Folds Choose the Degree

Regularisation

Choosing a lower degree is one way to stop overfitting. Regularisation is another: keep the flexible model, but change the loss so that wild solutions cost extra. Look at the weights of the overfitted degree-9 polynomial:

w9 = fit_poly(x_train, y_train, 9)
print([round(v, 1) for v in w9])

[-0.4, 4.0, 12.2, -25.0, -178.3, -88.2, 490.6, 434.2, -327.2, -328.9]

Large positive and negative weights, in the hundreds, cancel each other almost exactly on the training points and swing enormously in between. Ridge regression discourages this by adding a penalty on the size of the weights to the loss:

L(\mathbf{w}) = \sum_{i=1}^{n} (\hat{y}_i - y_i)^2 + \lambda \sum_{j=1}^{d} w_j^2

The regularisation strength \lambda (lambda) sets the trade-off: \lambda = 0 is ordinary least squares, and a large \lambda pushes all weights towards zero. The penalty's sum starts at j = 1: the intercept w_0 is usually left unpenalised, because shifting the whole curve up or down is not a sign of overfitting. For squared error the minimum still has a formula, the normal equations with \lambda added to the diagonal, (X^\top X + \lambda I)\mathbf{w} = X^\top \mathbf{y} (I is the identity matrix, 1 on the diagonal and 0 elsewhere), and that is exactly what the lam argument of fit_poly does. (With gradient descent you would add 2\lambda w_j to each partial derivative instead.)

for lam in (0, 0.001, 0.01, 0.1, 1):
    w = fit_poly(x_train, y_train, 9, lam)
    print(f"λ={lam:<6} train {mse(w, x_train, y_train):.4f}  new {mse(w, x_new, y_new):.4f}"
          f"  largest weight {max(abs(v) for v in w):.1f}")

λ=0      train 0.0044  new 60.1854  largest weight 490.6
λ=0.001  train 0.0092  new 0.1389  largest weight 3.0
λ=0.01   train 0.0140  new 0.1555  largest weight 2.1
λ=0.1    train 0.0232  new 0.1897  largest weight 1.6
λ=1      train 0.0955  new 0.3135  largest weight 0.9

A tiny penalty changes everything. At λ = 0.001 the largest weight drops from 490 to 3, the training error rises a little (the model is no longer free to chase every point), and the error on new data falls from 60 to 0.14, as good as the best low-degree polynomial. Too much penalty underfits again: at λ = 1 the curve is flattened towards a line.

So λ is a hyperparameter, a dial from overfit to underfit, and it is chosen by cross-validation in exactly the same way as the degree:

for lam in (0.0001, 0.001, 0.01, 0.1, 1):
    print(lam, round(cross_val_mse(x_all, y_all, 9, lam=lam), 4))

0.0001 0.0977
0.001 0.079
0.01 0.0843
0.1 0.0965
1 0.1452

Cross-validation picks λ = 0.001, the value that also did best on new data.

This is the optimisation thread at work. Nothing about the model changed; the objective did, from "fit the data" to "fit the data, but prefer small weights". Much of the craft of machine learning is choosing what to minimise. The lasso penalises |wⱼ| instead and tends to set some weights exactly to zero; early stopping simply halts gradient descent before it reaches the overfitted minimum.

Common mistakes: regularising features on different scales, where the penalty falls unevenly (standardise first, as in Level 3). Penalising the intercept. Choosing λ on the test set.

Practice: A Cautious Slope, Open the Test Set Once, Bikes, Rain and a Fair Penalty

Bias and variance

Why does error on new data trace a U? Imagine repeating the whole experiment many times: draw a fresh training set, train the model, record its predictions. Two different things can go wrong.

Bias: even averaged over all those training sets, the predictions are off, because the model family cannot represent the truth. A straight line through a wave is wrong in the same way every time.
Variance: the predictions change a lot from one training set to the next, because the model follows the noise of whichever sample it got. A degree-9 polynomial is wild in a different way every time.

For squared error the expected error at a point splits exactly into three parts:

\text{expected error} = \text{bias}^2 + \text{variance} + \text{noise}

It is the same decomposition that Statistics level 6 uses for estimators, where an estimator's mean squared error is its bias squared plus its variance. With generated data you can measure both parts directly. Keep the 12 inputs fixed and evenly spaced, draw the noise 200 times, and look at the predictions at a few test points:

fixed_x = [round(-1 + 2 * i / 11, 3) for i in range(12)]
grid = [-0.9, -0.6, -0.3, 0.0, 0.3, 0.6, 0.9]

def bias_variance(degree, runs=200):
    preds = {g: [] for g in grid}
    for seed in range(runs):
        rng = random.Random(seed)
        ys = [math.sin(3 * x) + rng.gauss(0, 0.3) for x in fixed_x]
        w = fit_poly(fixed_x, ys, degree)
        for g in grid:
            preds[g].append(predict(w, g))
    bias2 = variance = 0.0
    for g, p in preds.items():
        mean_p = sum(p) / len(p)
        bias2 += (mean_p - math.sin(3 * g)) ** 2
        variance += sum((q - mean_p) ** 2 for q in p) / len(p)
    return bias2 / len(grid), variance / len(grid)

for d in (1, 3, 5, 7, 9):
    b2, v = bias_variance(d)
    print(d, f"bias² {b2:.4f}  variance {v:.4f}  expected error {b2 + v + 0.09:.4f}")

1 bias² 0.1770  variance 0.0156  expected error 0.2826
3 bias² 0.0043  variance 0.0270  expected error 0.1212
5 bias² 0.0003  variance 0.0383  expected error 0.1285
7 bias² 0.0004  variance 0.0659  expected error 0.1563
9 bias² 0.0026  variance 0.1924  expected error 0.2850

The table is the U curve taken apart. Bias drops sharply from degree 1 to degree 3 and stays near zero after that (from degree 3 the family is rich enough to follow the wave), while variance rises steadily: each extra degree gives the noise more ways to move the curve. Their sum, plus the 0.09 of noise that no model can remove, is smallest at degree 3. Underfitting is high bias; overfitting is high variance; regularisation and more data both reduce variance, which is why they fix overfitting.

The same trade-off runs through the other models: kNN with small k has low bias and high variance, large k the reverse; so do deep and shallow decision trees. Averaging many high-variance models trained on resampled data (bagging, the idea behind random forests) reduces variance without adding bias, the bootstrap of Statistics level 4 put to work.

Common mistakes: thinking bias and variance can be read off one training run; they are properties of the procedure over many possible training sets. Forgetting the noise term, and expecting a perfect model to reach zero error.

Practice: A Thousand Imaginary Days

Decision trees

The threshold rule of Level 2 asked one question. A decision tree asks a sequence of them: "is sweetness above 4.9? If so, is crunchiness above 7.1? …", and each answer leads either to another question or to a leaf with a prediction. Trees handle features on any scale without standardisation, capture interactions between features, and can be printed and read.

A tree is grown greedily from the top. At each node, try every feature and every threshold (the midpoints between sorted values, as in Level 2) and pick the split that makes the two sides as pure as possible, meaning each side is dominated by one class. Purity is measured by the Gini impurity: the probability that two examples drawn at random from the node (with replacement) have different labels,

G = 1 - \sum_{i} p_i^2

over the class proportions p_i. It is 0 for a pure node and 0.5 for a 50/50 node of two classes. (Entropy, H = -\sum_i p_i \log_2 p_i, is the other common choice and gives very similar trees.) In gini below, c / n is p_i and the for runs over the classes:

from collections import Counter

def gini(labels):
    n = len(labels)
    return 1 - sum((c / n) ** 2 for c in Counter(labels).values())

print(gini(["a", "a", "a", "a"]), gini(["a", "a", "b", "b"]), round(gini(["a", "a", "a", "b"]), 3))

0.0 0.5 0.375

A split is scored by the Gini impurity of its two sides, weighted by their sizes, and the best split is the one with the lowest score. The tree then recurses on each side, and stops when a node is pure, too small to split, or at the maximum depth. Here a tree is a nested dictionary and a leaf is just a label:

def best_split(X, y):
    best = None
    for f in range(len(X[0])):
        values = sorted(set(x[f] for x in X))
        for a, b in zip(values, values[1:]):
            t = (a + b) / 2
            left = [lab for x, lab in zip(X, y) if x[f] <= t]
            right = [lab for x, lab in zip(X, y) if x[f] > t]
            score = (len(left) * gini(left) + len(right) * gini(right)) / len(y)
            if best is None or score < best[0]:
                best = (score, f, round(t, 2))
    return best

def grow(X, y, depth, max_depth, min_size=2):
    majority = Counter(y).most_common(1)[0][0]
    if depth == max_depth or len(y) < min_size or gini(y) == 0:
        return majority
    split = best_split(X, y)
    if split is None:                          # every example identical
        return majority
    _, f, t = split
    L = [i for i, x in enumerate(X) if x[f] <= t]
    R = [i for i, x in enumerate(X) if x[f] > t]
    return {"feature": f, "threshold": t,
            "left": grow([X[i] for i in L], [y[i] for i in L], depth + 1, max_depth, min_size),
            "right": grow([X[i] for i in R], [y[i] for i in R], depth + 1, max_depth, min_size)}

def classify(tree, x):
    while isinstance(tree, dict):
        tree = tree["left"] if x[tree["feature"]] <= tree["threshold"] else tree["right"]
    return tree

The data: spots on a 10 km by 10 km map, with a phone mast at (5, 5) that reaches 3.5 km. A spot has "signal" if it is within range, except that 10% of the survey readings are wrong. The true boundary is a circle, which a tree must approximate with horizontal and vertical cuts.

def make_survey(n, seed):
    rng = random.Random(seed)
    X, y = [], []
    for _ in range(n):
        spot = [round(rng.uniform(0, 10), 1), round(rng.uniform(0, 10), 1)]   # km east, km north
        label = "signal" if math.dist(spot, (5, 5)) < 3.5 else "none"
        if rng.random() < 0.1:                                                # a wrong reading
            label = "none" if label == "signal" else "signal"
        X.append(spot); y.append(label)
    return X, y

def accuracy(tree, X, y):
    return sum(classify(tree, x) == lab for x, lab in zip(X, y)) / len(y)

X_tr, y_tr = make_survey(150, seed=3)
X_nw, y_nw = make_survey(2000, seed=99)
print(grow(X_tr, y_tr, 0, max_depth=2))
for depth in (1, 2, 3, 4, 5, 6, 8, 12):
    tree = grow(X_tr, y_tr, 0, depth)
    print(depth, f"train {accuracy(tree, X_tr, y_tr):.3f}  new {accuracy(tree, X_nw, y_nw):.3f}")

{'feature': 1, 'threshold': 2.15, 'left': {'feature': 0, 'threshold': 5.35, 'left': 'none', 'right': 'none'}, 'right': {'feature': 1, 'threshold': 8.25, 'left': 'signal', 'right': 'none'}}
1 train 0.620  new 0.557
2 train 0.740  new 0.692
3 train 0.827  new 0.770
4 train 0.907  new 0.845
5 train 0.913  new 0.835
6 train 0.940  new 0.819
8 train 0.960  new 0.797
12 train 1.000  new 0.794

The depth-2 tree is readable: it first cuts off the southern strip (north below 2.15 km, out of range), then the northern strip. With only four rectangles it cannot follow a circle, so it underfits. Depth is the complexity dial. Accuracy on new data rises until depth 4 (0.845, close to the 0.9 that the 10% of wrong readings allows), then falls, while training accuracy climbs to a perfect 1.000 at depth 12: the deepest tree has carved out a tiny rectangle for every wrong reading. It is the polynomial story again, with depth in place of degree, and the same tools apply: choose the depth (or the minimum leaf size) by cross-validation. The split inside the southern strip ends in two leaves that both say "none": a split can reduce impurity without changing any prediction, which is harmless.

Growing a tree is a greedy optimisation, like the constructive heuristics of Algorithms level 6: each split is the best available at that moment, with no look-ahead, because finding the best tree overall is NP-hard. Deep trees have high variance, which is why practical systems average many trees (random forests, gradient-boosted trees), among the most successful models for tabular data.

Common mistakes: growing trees without a depth or leaf-size limit and trusting their perfect training accuracy. Recomputing every candidate split from scratch on large data (sorting each feature once and sweeping the threshold is the efficient way, as for the stump in Level 2). Forgetting a stopping rule for nodes whose examples are all identical.

Practice: How Mixed Is the Basket?, Is Rain the Better Question?, The Question That Sorts Best

k-means clustering

Every model so far learned from labels. Unsupervised learning has no labels: the task is to find structure in the feature vectors alone. The classic example is clustering, grouping points so that points in the same group are close together. k-means does this for a chosen number of clusters k, and it is an optimisation problem from the start. Its loss (often called the inertia) is the total squared distance from each point to the centre of its cluster, \sum_{i} \lVert \mathbf{x}_i - \mathbf{c}_{a_i} \rVert^2, where a_i is the cluster of point \mathbf{x}_i and \mathbf{c}_j the centre of cluster j, and the algorithm alternates two steps, each of which can only lower that loss:

Assign each point to its nearest centre (the closest-vector search of Level 1).
Update each centre to the mean of the points assigned to it (the centroid of Level 1, which by the Level 2 argument is the point with the smallest total squared distance to them).

Repeat until the assignments stop changing. Here are 30 points drawn around three centres:

def make_blobs(seed):
    rng = random.Random(seed)
    points = []
    for cx, cy in ((2, 2), (8, 3), (5, 8)):
        for _ in range(10):
            points.append([round(rng.gauss(cx, 1), 1), round(rng.gauss(cy, 1), 1)])
    return points

def sq_dist(a, b):
    return sum((x - y) ** 2 for x, y in zip(a, b))

def k_means(points, k, seed, max_rounds=100):
    centres = [p[:] for p in random.Random(seed).sample(points, k)]    # k random points to start
    for _ in range(max_rounds):
        assign = [min(range(k), key=lambda j: sq_dist(p, centres[j])) for p in points]
        new = []
        for j in range(k):
            members = [p for p, a in zip(points, assign) if a == j]
            new.append([sum(c) / len(members) for c in zip(*members)] if members else centres[j])
        if new == centres:
            break
        centres = new
    inertia = sum(sq_dist(p, centres[a]) for p, a in zip(points, assign))
    return centres, assign, inertia

points = make_blobs(1)
for seed in (0, 4, 9):
    centres, assign, inertia = k_means(points, 3, seed)
    print(seed, round(inertia, 1), [[round(v, 1) for v in c] for c in centres])

0 42.2 [[1.9, 2.1], [8.2, 3.2], [5.5, 7.6]]
4 174.6 [[1.9, 2.6], [6.9, 5.4], [1.9, 1.0]]
9 255.1 [[3.7, 4.9], [7.8, 3.5], [8.9, 2.7]]

From seed 0 the centres land close to the true ones, (2, 2), (8, 3) and (5, 8). From seeds 4 and 9 k-means converges too, but to a poor answer: two centres share one blob while a single centre sits between the other two, and neither step can fix that: every point is already assigned to its nearest centre, and every centre is already the mean of its points. These are local minima of a non-convex loss, as in Level 4, and the same remedy works: run from several starting points and keep the lowest inertia. (The k-means++ initialisation, which spreads the starting centres apart, makes bad starts much rarer.)

How many clusters? The inertia always falls as k grows (with one centre per point it is zero), so it cannot choose k by itself. Plot it against k and look for the elbow, where adding clusters stops helping much:

for k in range(1, 7):
    best = min(k_means(points, k, seed)[2] for seed in range(10))
    print(k, round(best, 1))

The drop from 3 clusters to 4 is small compared with the drops before it, so 3 is the natural answer here. Real data is rarely this clear-cut, and the right k often depends on what the clusters are for.

k-means inherits the assumptions of its distance: features must be on comparable scales (standardise first, as in Level 3), and clusters are assumed to be round blobs of similar size. Statistics level 6 relaxes that with Gaussian mixture models fitted by EM, which is k-means with soft assignments and clusters of different shapes.

Common mistakes: running k-means once and trusting the result. Forgetting that an empty cluster has no mean (keep the old centre, or restart it at a far-away point). Using unscaled features. Comparing inertia across different k as if lower were always better.

Practice: Move the Depots Once, Rounds Until the Depots Settle, How Many Depots?

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.

Checklist

Fit polynomials with the normal equations, and recognise underfitting and overfitting from training and new-data error.
Explain hyperparameters, and choose them with a validation set while keeping the test set for one final report.
Implement k-fold cross-validation with fixed folds, and use it to choose a degree or a regularisation strength.
Add a ridge penalty to the loss, solve the regularised normal equations, and explain why small weights generalise better.
Measure bias² and variance by retraining on many simulated datasets, and relate them to underfitting and overfitting.
Grow a decision tree with Gini impurity and greedy splits, and control its complexity with the maximum depth.
Implement k-means as alternating assignment and update steps, restart it against local minima, and choose k with the elbow of the inertia.
Describe regularisation, hyperparameter search, tree growing and k-means each as an optimisation problem.