A line that predicts

Given an input x, predict a continuous output y with a straight line: ŷ = w·x + b. The slope w says how much y changes per unit of x; the intercept b is the value when x = 0. "Learning" means choosing the single w and b that make the line pass as close as possible to all the data at once — and "close" has to be defined precisely before we can optimise it.

A tiny real dataset

Everything below uses these five points — small enough to compute entirely by hand, so nothing is hidden. Think of x as hours studied and y as a quiz score:

There is an upward trend but no perfect line through all five — real data has noise. The job is the line that is wrong by the least, totalled across every point.

Measuring "close"

For each point the line makes an error, the residual ŷ_i − y_i. We square it (so positive and negative misses both count, and big misses count disproportionately) and average over the dataset. That is the mean squared error:

L is a smooth bowl-shaped function of the two unknowns w and b — and a bowl has exactly one lowest point. Finding it is the whole game, and there are two routes there: solve for it directly, or roll downhill to it.

Solving for the bottom of the bowl

At the minimum the bowl is flat in every direction — nudging w a hair does not lower the cost, and neither does nudging b. "How fast the cost changes as you nudge w" is exactly the partial derivative ∂L/∂w, so the bottom of the bowl is simply where both partial derivatives are zero. Working out those two conditions turns the calculus into two ordinary equations in w and b (the normal equations), which we then solve like any pair of simultaneous equations. Every step, with the arithmetic written out on our five points:

Gradient descent finds the same answer

The closed form is a luxury of linear regression; most models have no such formula and must descend the loss surface step by step. So watch the other route here: start from a flat line and repeatedly nudge w and b downhill along the gradient. The dashed grey line is the exact least-squares answer from above (ŷ = 0.6x + 2.2); the solid line is gradient descent chasing it, with the vertical sticks showing the current residuals:

The bowl we slide down

The loss is a surface — a bowl floating above the plane of all possible parameter pairs, its height the error that choice of line makes. The lowest point is the best line, and gradient descent is a ball that rolls straight downhill: each step it reads the slope beneath it (the gradient) and moves a little way down, the step length set by the learning rate α. One set-up choice makes the picture clean — we plot the slope w against the centred intercept c, the line's height at the average input x̄ = 3. Centring the inputs leaves the line unchanged (the slope is still 0.6) but turns the loss surface from a long, skewed valley into a near-round bowl, so the ball's descent is easy to follow all the way down. Step it by hand toward the bottom at (w, c) = (0.6, 4.0):

Linear regression in Python, step by step

The math is the algorithm — so here it is as a short Python program you can paste into any console and run. We build it up one piece at a time, in the same order as the ideas above, using nothing but plain Python for the core (no libraries).

Step 1 — the data and the model

First the same five points, and the model itself: a function that, given a slope w and intercept b, predicts y for any x.

# the data: hours studied (x) and quiz score (y)
xs = [1, 2, 3, 4, 5]
ys = [2, 4, 5, 4, 5]
n  = len(xs)

# the model: a straight line. given w and b, predict y for an x
def predict(x, w, b):
    return w * x + b

Step 2 — the cost function (what we optimise)

This is the single number the whole algorithm tries to shrink: the mean squared error. "Optimising" the model means nothing more than finding the w and b that make cost return its smallest value.

# the cost: mean squared error -- the single number we want to make small
def cost(w, b):
    total = 0.0
    for i in range(n):
        error = predict(xs[i], w, b) - ys[i]   # residual for point i
        total += error ** 2
    return total / n

Step 3 — the gradients

To roll downhill we need the slope of the cost in the w and b directions — the exact partial derivatives we derived by hand, transcribed line for line.

# the gradients: how the cost shifts when we nudge w and b
# (these are the two partial derivatives we worked out above)
def gradients(w, b):
    dw = db = 0.0
    for i in range(n):
        error = predict(xs[i], w, b) - ys[i]
        dw += (2 / n) * error * xs[i]   # dCost/dw
        db += (2 / n) * error           # dCost/db
    return dw, db

Step 4 — gradient descent

Start with a flat line at zero, then repeatedly step each parameter a little way against its gradient (downhill). The learning rate lr sets the step size.

# gradient descent: start anywhere, step downhill, repeat
w, b = 0.0, 0.0          # start with a flat line at zero
lr   = 0.02              # learning rate: the size of each step

for step in range(2001):
    dw, db = gradients(w, b)
    w -= lr * dw         # move opposite the gradient
    b -= lr * db
    if step % 500 == 0:
        print(f"step {step:4d}   cost {cost(w,b):.4f}   w {w:.3f}   b {b:.3f}")

print(f"learned line:  y = {w:.3f} x + {b:.3f}")

Running it prints the cost falling and the parameters settling — onto the exact values we found by hand:

The cost bottoms out at 0.4800 — which is exactly the SSE of 2.40 from the R² section divided by n = 5 — and the line is ŷ = 0.600x + 2.200, matching the by-hand answer to three decimals.

Step 5 — the shortcut (closed form)

Because linear regression has a formula, you rarely descend in practice — you solve it directly. Here is the same answer with no loop, plus the one-line library call that does it all internally.

import numpy as np
x = np.array(xs); y = np.array(ys)

# closed form for a line -- no iteration
slope     = np.sum((x - x.mean()) * (y - y.mean())) / np.sum((x - x.mean())**2)
intercept = y.mean() - slope * x.mean()
print(f"closed form:  y = {slope:.3f} x + {intercept:.3f}")

# or let numpy's least-squares fit do it in one call
w_np, b_np = np.polyfit(x, y, 1)
print(f"np.polyfit:   y = {w_np:.3f} x + {b_np:.3f}")

Why it is called "least squares"

There is a beautiful picture behind the algebra. Stack the targets into a vector y, and consider every line the model can produce as a flat plane of possibilities. Least squares finds the point on that plane closest to y — its projection — and the line from y to that closest point is perpendicular to the plane. That perpendicularity is the normal equations: at the optimum the residual vector is orthogonal to the inputs, Σ(ŷ−y) = 0 and Σ(ŷ−y)x = 0, which is exactly why those two sums vanished above.

Check it on our own fit. The residuals are y − ŷ = [−0.8, 0.6, 1.0, −0.6, −0.2]. Their sum is −0.8 + 0.6 + 1.0 − 0.6 − 0.2 = 0 (the first normal equation), and their dot product with x is (−0.8)(1) + (0.6)(2) + (1.0)(3) + (−0.6)(4) + (−0.2)(5) = −0.8 + 1.2 + 3.0 − 2.4 − 1.0 = 0 (the second). The residuals are exactly perpendicular to the data — geometrically, that is what makes this line the best one.

R², the variance explained

A loss of "2.4" means nothing on its own. The standard yardstick compares the line against the laziest possible model — always guessing the mean ȳ — and reports the fraction of the variance the line removes:

From a line to a hyperplane

Real problems have many inputs — size, bedrooms, age. Each gets its own weight, ŷ = b + w₁x₁ + w₂x₂ + …, and the whole dataset stacks into a design matrix X (a column of 1s for the intercept, then one column per feature). The same "set the gradient to zero" step becomes a single matrix equation:

Polynomial regression

"Linear" refers to the weights, not the shape. Feed the model x, x², x³ … as extra features and the very same least-squares machinery fits a curve — it is still linear in the coefficients. But more flexibility is a double-edged sword: too low a degree and the curve cannot follow the data (underfitting); too high and it contorts through every noisy point, memorising rather than generalising (overfitting). Slide the degree and watch the trade-off:

Watch the training error fall to near zero at high degree while the curve becomes visibly absurd between points. That gap — fitting the data you have versus the data you will see — is the central drama of machine learning, and the subject of the bias–variance dive.

When the line lies

Linear regression is fast, interpretable and exact, but it leans on assumptions. It expects a roughly linear relationship (curved data needs features like x² or a different model); it is sensitive to outliers, because squaring a large residual gives it enormous pull on the line; and when features are strongly correlated (multicollinearity) the matrix XᵀX becomes nearly singular, so the weights turn unstable and uninterpretable even when predictions stay fine.

It also assumes the noise is roughly even across the range (homoscedasticity) and that you scale features before comparing weights. None of these are fatal — they are signposts pointing to regularisation, robust losses, and feature engineering, each a tool in later dives.

Why it is everywhere

For all the sophistication that came after, linear regression remains a default first model — and often a hard one to beat. It trains instantly, needs little data, and every coefficient is a plain-English statement ("each extra bedroom adds $20k"). It is the baseline every other model must justify beating, the building block inside logistic regression and neural networks (a single neuron is a linear regression followed by a squashing function), and the cleanest place to meet the loss-and-gradient pattern that runs through all of machine learning.

What you saw

Linear regression fits ŷ = wx + b by minimising mean squared error — a convex bowl with a single minimum. Setting both partial derivatives to zero gave the normal equations, which we solved by hand to the exact w = 0.60, b = 2.20; gradient descent then rediscovered the identical numbers by rolling downhill. R² rated the fit at 0.60 of the variance explained. The design matrix carried the whole method to many features via XᵀX w = Xᵀy, and polynomial features extended it to curves — straight into the overfitting that motivates regularisation next.

Charge for complexity

Take the familiar least-squares loss and bolt a penalty onto it — a term that grows with the size of the weights. Now the model minimises error plus weight-size, so it only grows a weight when the error reduction is worth the price. A single knob λ sets how steep that price is: λ = 0 is ordinary regression; large λ forces the weights toward zero. The two classic penalties differ only in how they measure "size."

When weights blow up

Recall where linear regression left off: extra flexibility invites overfitting, and when features are correlated the matrix XᵀX turns nearly singular, so the closed-form weights explode. Two features that carry almost the same information can be handed enormous opposing weights — one feature +800, its near-twin −790 — which cancel on the training data but make the model violently unstable on anything new.

Both failures share a symptom: weights that are far larger than the signal warrants. So the cure is direct — make large weights costly. The penalty term does exactly that, and as a bonus it makes the unstable matrix invertible again, as we are about to see.

Here is that instability in actual numbers, on two nearly-identical features:

Shrinking, in closed form

Ridge adds the squared length of the weight vector to the loss. Because the penalty is smooth, the whole thing still has an exact closed-form solution — the OLS formula with one small change:

There is a deeper way to see the shrinkage. Rewrite the data along its principal directions (the natural axes of how it varies); Ridge then scales each direction by a factor d²∕(d² + λ), where d measures how much the data spreads along that direction. Directions the data barely explores — small d, which are exactly the unstable, near-duplicate ones — get shrunk the hardest, while the strong, well-measured directions are left almost untouched. That is precisely how Ridge calms the instability without erasing the genuine signal.

One change, a different world

Lasso changes a single thing: it penalises the absolute value of the weights instead of their square.

That swap has a dramatic consequence. Ridge shrinks weights toward zero but never quite reaches it; Lasso drives many weights to exactly zero, switching features off entirely. A model that selects its own features is enormously useful — it is interpretable and it ignores noise. But why should replacing w² with |w| turn gentle shrinkage into outright deletion? The answer is geometric, and worth seeing.

Circles shrink, diamonds delete

Read the penalty as a budget: minimise the error while keeping the weights inside a region of fixed size. L2 makes that region a circle (w₁² + w₂² ≤ t); L1 makes it a diamond (|w₁| + |w₂| ≤ t). The unconstrained least-squares answer sits outside the budget, so the solution is where the expanding rings of equal loss first touch the region. A circle is smooth, so they meet at a slanted point with both weights nonzero. A diamond has sharp corners sitting right on the axes — and an expanding contour almost always strikes a corner first, where one weight is exactly zero. Toggle the penalty and slide λ to watch the solution travel inward:

(The “budget of size t” picture and the “add λ·penalty” formula are two views of one problem: every budget t corresponds to some λ, the link being a Lagrange multiplier. A bigger λ is a tighter budget. The slider below uses the λ form.)

The cutoff that deletes

"Touches a corner" can be made exact. Consider one weight at a time (the clean case where features don't overlap); the Lasso solution for that weight is a simple rule with a hard cutoff built in:

Because |w| has no derivative at zero, Lasso has no one-shot matrix formula the way Ridge does. But this single-weight rule is all we need: cycle through the weights, applying soft-thresholding to each in turn (against the residual the other weights leave behind), and repeat until nothing moves. That loop is coordinate descent — exactly the Lasso code in the “In code” chapter.

Watch every coefficient move

Sweep λ from zero upward and track all the weights at once. This dataset has five standardised features, but only two carry real signal: x₀ (true weight +3) and x₁ (true weight −2). x₂ and x₃ are pure noise, and x₄ is a near-copy of x₀ (90% correlated). Toggle Ridge versus Lasso and slide λ:

Both penalties at once

The path above exposes Lasso's one weakness: faced with a group of correlated features (our x₀ and x₄), it picks among them somewhat arbitrarily rather than sharing weight sensibly. Elastic Net simply adds both penalties, inheriting Lasso's sparsity and Ridge's stable "grouping" of correlated features, which it shrinks together:

With λ₂ = 0 it is pure Lasso; with λ₁ = 0 it is pure Ridge; in between it is the pragmatic default when you have many predictors and suspect some are both correlated and irrelevant.

The dial between under- and overfitting

λ is the bias–variance dial. At λ = 0 you have plain OLS: low bias, high variance — it overfits. Crank λ enormous and every weight is crushed to zero: high bias, no variance — it underfits. The useful value lives in between and is found by cross-validation — try a grid of λ values, measure error on held-out data, and keep whichever generalises best. Here is the overfitting degree-8 polynomial from the linear-regression dive, now with a Ridge penalty; slide λ and watch the frantic curve settle down:

Ridge and Lasso in Python, step by step

Both methods are short to write. Ridge is the OLS solve with λ added to the diagonal; Lasso is a loop that applies the soft-thresholding rule one weight at a time until the weights stop moving (coordinate descent). Here they are on a five-feature dataset where only the first two features matter.

Step 1 — build a dataset with junk features

import numpy as np
np.random.seed(0)

# 60 samples, 5 features. only x0 (+3) and x1 (-2) truly matter.
# x2, x3 are pure noise;  x4 is a near-copy of x0 (correlated, redundant)
n  = 60
x0 = np.random.randn(n)
x1 = np.random.randn(n)
x2 = np.random.randn(n)
x3 = np.random.randn(n)
x4 = 0.9 * x0 + 0.3 * np.random.randn(n)
X  = np.column_stack([x0, x1, x2, x3, x4])
y  = 3*x0 - 2*x1 + 0.4*np.random.randn(n)

# standardize features and center y so the penalty is fair to every weight
X = (X - X.mean(0)) / X.std(0)
y = y - y.mean()

Step 2 — Ridge: just add λ to the diagonal

The entire method is the normal equations with + lam*I. That single change both shrinks the weights and rescues the matrix when features are correlated.

# Ridge closed form:  w = (X'X + lam*I)^-1 X'y
def ridge(lam):
    p = X.shape[1]
    return np.linalg.solve(X.T @ X + lam*np.eye(p), X.T @ y)

Step 3 — Lasso: soft-threshold one weight at a time

Lasso has no closed form, so we cycle through the weights, each time fitting the leftover residual and applying the cutoff rule we derived. This is coordinate descent.

# the soft-thresholding rule:  shrink toward 0, and snap to 0 if small
def soft(z, t):
    return np.sign(z) * max(abs(z) - t, 0.0)

def lasso(lam, iters=500):
    p    = X.shape[1]
    w    = np.zeros(p)
    col2 = (X**2).sum(0)
    for _ in range(iters):
        for j in range(p):
            rho  = X[:, j] @ (y - X @ w + X[:, j]*w[j])   # fit ignoring feature j
            w[j] = soft(rho, lam) / col2[j]
    return w

Step 4 — compare them

print("feature:            x0     x1     x2     x3     x4")
print("OLS  (no penalty): ", ridge(0.001))
print("Ridge (lam = 20):  ", ridge(20))
print("Lasso (lam = 8):   ", lasso(8))

Read the rows. OLS spreads a little weight onto every feature, noise included. Ridge shrinks all five toward zero and, unable to delete anything, splits the shared signal across the twins x₀ and x₄ (note x₄ climbs to 1.21). Lasso sets the two pure-noise features to exactly 0.00 — it has selected features — and trims the redundant x₄ to a sliver.

Step 5 — scikit-learn computes exactly this

The library wraps the same two formulas. Its penalty constant is scaled a little differently — Ridge uses α = λ, while Lasso divides by n (α = λ/n) because its objective carries a 1/2n factor — but the coefficients come out identical to the digit:

from sklearn.linear_model import Ridge, Lasso
print(Ridge(alpha=20,    fit_intercept=False).fit(X, y).coef_)
print(Lasso(alpha=8/n,   fit_intercept=False).fit(X, y).coef_)

Almost always on

Regularization is the default, not the exception. Ridge — under the name weight decay — sits inside virtually every neural network as the standard stabiliser. Lasso is the tool of choice when you suspect most features are useless and want the model to tell you which ones to keep. Elastic Net hedges between them when predictors are both numerous and correlated. The one non-negotiable step is to standardise features first: the penalty compares weights directly, so a feature measured in larger units would otherwise be punished merely for its scale.

What you saw

Regularization adds a price for large weights, curing both overfitting and the weight explosions caused by correlated features. Ridge's L2 penalty keeps a clean closed form, w = (XᵀX + λI)⁻¹Xᵀy, and shrinks every weight smoothly — the +λI literally rescuing the matrix. Lasso's L1 penalty has no closed form but, through the soft-thresholding rule w* = sign(z)·max(|z|−λ, 0), drives small weights to exactly zero — feature selection for free. The constraint geometry explained the gap: a smooth circle versus a cornered diamond. You watched a live path shrink under Ridge and snap to zero under Lasso, met Elastic Net's blend, and wrote both in a few lines of Python whose Lasso matched scikit-learn to the digit.

The most confident boundary

Suppose two classes are cleanly separated. Infinitely many straight lines split them — so which is best? A boundary jammed against the points will misread the next slightly-shifted example. The most robust choice sits in the middle of the widest gap, leaving the most room for error on both sides. That widest gap is the margin, and maximising it is the whole goal.

Four points, two classes

We'll work with four points. Two of them — one per class — will turn out to be the only ones that matter.

Score, then take the sign

A linear classifier scores a point with f(x) = w·x + b and predicts by the sign: positive → class +1, negative → class −1. The boundary is where the score is zero. The weight vector w points perpendicular to that boundary — it is the direction "across the street."

How wide is the street?

The distance from a point x to the line is |w·x + b| / ‖w‖. The same line can be written with many different scales of w and b, so we pin the scale down by insisting the closest points score exactly ±1:

Each support vector then lies a distance 1/‖w‖ from the boundary, so the full street — edge to edge — has width

A smaller ‖w‖ means a wider street. That is the lever we will pull.

The optimisation problem

Maximising 2/‖w‖ is the same as minimising ‖w‖, and minimising ‖w‖ is cleanest written as ½‖w‖². Every point must stay off the street — on the correct side and at least up to its margin line:

That's the entire hard-margin SVM: a bowl-shaped (convex quadratic) objective with one linear constraint per point. There is exactly one solution.

The four points, worked out

Because the data is symmetric we can solve this with pencil and paper. The two support vectors sit on the edges of the street, where the constraint is tight (= 1).

Only the edge points matter

The solution touched just two points — those with y(w·x+b) = 1, sitting on the margin lines. These are the support vectors. Move or delete any other point and the boundary is unchanged; nudge a support vector and the whole street shifts. Sweep the boundary's angle below and watch the street width: it is widest at exactly one orientation — the max-margin solution (45°, width 2√2 ≈ 2.83). The support vectors are the points the street bumps into.

The same loop, with a kinked loss

An SVM can be trained exactly like every other machine in this lab. Replace the hard constraint with a penalty that is zero while a point sits safely outside the street and grows linearly once it crosses the margin line — the hinge loss:

Minimising this by (sub-)gradient descent is the familiar loop: predict f = w·x + b, measure the hinge loss, take its gradient (for a violating point, nudge w by +C·y·x; the ½‖w‖² term gently shrinks w), and let the optimiser step — the very same SGD slot as linear regression. The hinge is just a convex, kink-at-1 stand-in for the un-trainable 0/1 mistake count.

When the classes overlap

Real data overlaps, and then no line separates it perfectly. The soft margin lets points slip inside the street, or even onto the wrong side, each paying a hinge penalty. The constant C sets the price of a violation: a large C insists on classifying everyone correctly and accepts a narrow street; a small C tolerates some mistakes in exchange for a wider, simpler boundary. Slide C and watch the street breathe.

Bending the boundary

Some data no straight line can split — like one class ringed by another. The kernel trick lifts the points into a higher dimension where a flat boundary does work. Add one feature, z = x₁² + x₂²: the inner ring now sits low, the outer ring high, and a flat plane slices cleanly between them. Projected back to the plane, that cut is a circle. Slide the boundary radius below — the strip on the right shows the lifted feature z, where the split is just a horizontal threshold (a straight cut). SVMs do this implicitly through a kernel and never compute the lift; on these rings a linear kernel scores 70%, an RBF kernel 100%.

SVMs in scikit-learn

The same problem in a few lines — the linear solver recovers our hand solution exactly, a from-scratch hinge-loss SGD lands in the same place, and on the rings the kernel earns its keep.

# the tiny worked example
X = np.array([[2,2],[3,2],[0,0],[-1,0]]); y = np.array([1,1,-1,-1])
clf = SVC(kernel="linear", C=1e6).fit(X, y)   # big C -> hard margin
print(clf.coef_[0], clf.intercept_[0])
print("margin =", 2/np.linalg.norm(clf.coef_[0]))
print("support vectors =", clf.support_vectors_)

# from scratch: sub-gradient descent on the hinge loss
def fit_hinge(X, y, C=50, epochs=3000):
    w = np.zeros(2); b = 0.0
    for t in range(1, epochs+1):
        lr = 1/(0.5*t)
        m = y*(X@w + b) < 1                  # margin-violating points
        w -= lr*(w - C*(y[m,None]*X[m]).sum(0)/len(X))
        b -= lr*(   - C*(y[m]).sum()/len(X))
    return w, b

# kernels on concentric circles
print(SVC(kernel="linear").fit(Xc, yc).score(Xc, yc),
      SVC(kernel="rbf", gamma=0.5).fit(Xc, yc).score(Xc, yc))

output
w = [0.5 0.5]   b = -1.0
margin width = 2/||w|| = 2.828
support vectors = [[0.0, 0.0], [2.0, 2.0]]
from-scratch hinge SGD: w = [0.51 0.51]  b = -1.01
circles  linear-kernel acc = 0.7   RBF-kernel acc = 1.0

Using SVMs well

SVMs are distance-based, so always scale your features first. Tune C and the kernel width γ together by cross-validation: too large and the model memorises, too small and it underfits. Reach for an RBF kernel on nonlinear data and a linear kernel for high-dimensional sparse data such as text. SVMs are excellent on small-to-medium clean datasets and produce sparse models defined by a handful of support vectors, but they scale poorly to very large datasets, where the SGD/logistic and tree families take over.

A flowchart that learns

Picture a game of twenty questions. Each internal node of the tree asks about one feature — "is x ≤ 0.5?" — and sends you left or right. You follow the answers down until you reach a leaf, which holds the prediction. Training a tree means choosing which questions to ask, and in what order, so that each one separates the classes as cleanly as possible.

A staircase of two classes

Here are 32 points in two features. Class A sits on the left and in the top-right; class B fills the bottom-right. No single straight line splits them — but a short sequence of axis-aligned cuts will.

Cutting on a single feature

A split picks one feature and one threshold and divides the points into two groups: those at or below the threshold go left, the rest go right. Every node makes exactly one such axis-aligned cut. The only question is which feature and which threshold — and to choose, we need a way to score how "clean" the two resulting groups are.

How mixed is a group?

Impurity measures how mixed the labels in a node are: zero when the node holds a single class, largest when it is an even split. Two common measures, where p is the fraction of one class:

Both peak at a 50/50 mix and fall to zero for a pure node. Gini is slightly cheaper to compute and is scikit-learn's default; the two almost always pick the same splits.

Reward the cleanest cut

A good split lowers impurity. The information gain is the parent's impurity minus the size-weighted average impurity of the two children:

The tree tries every feature and every threshold and keeps the split with the highest gain. That single greedy rule, applied over and over, is the whole training algorithm.

Try it: pick a feature, slide the threshold, and watch each side’s impurity and the resulting gain. The amber dashed line marks the best split this feature allows.

Searching the thresholds

Take six points on one feature with classes A, A, A, B, A, B. We'll score every threshold and keep the best.

Split, then split the pieces

After the first cut, each side is just a smaller dataset — so we recurse, splitting each child the same way, until a node is pure (or we decide to stop). The result partitions the plane into axis-aligned rectangles, one per leaf. Slide the depth: at depth 1 a single cut already gets 81%, and at depth 2 a second cut on the other feature reaches 100% — after which the tree stops, since every leaf is already pure.

When the tree memorises

A tree can always reach 100% on the training set — just keep splitting until every leaf holds one point. On clean data that's fine, but on noisy, overlapping data it carves out tiny islands around individual points: it has memorised the noise, not learned the pattern. Slide the depth on this messier dataset and watch the regions fragment as accuracy creeps to 100%.

The cures are all forms of "ask fewer questions": cap the depth, require a minimum number of points per leaf, or grow the full tree and then prune back the branches that don't help on held-out data.

One tree is twitchy

A single tree is high-variance: move a few points and the whole structure can change. The fix is to grow many trees on bootstrapped samples and random feature subsets, then average their votes — a random forest. Boosting takes a different route, growing small trees in sequence where each fixes the last one's mistakes. Both turn the humble, twitchy tree into one of the strongest off-the-shelf models there is.

From scratch, then scikit-learn

The heart of a tree is one routine: try every split and keep the one with the lowest weighted Gini, then recurse. The whole thing is about fifteen lines of NumPy — and it chooses the same root split and reaches the same accuracy as scikit-learn:

# --- decision tree from scratch ---
def gini(y):
    p = np.bincount(y, minlength=2) / len(y)
    return 1 - (p**2).sum()

def best_split(X, y):
    best = (None, None, gini(y))
    for f in range(X.shape[1]):
        v = np.unique(X[:, f]); thrs = (v[:-1] + v[1:]) / 2     # midpoints
        for thr in thrs:
            L = X[:, f] <= thr
            w = (L.sum()*gini(y[L]) + (~L).sum()*gini(y[~L])) / len(y)
            if w < best[2]: best = (f, thr, w)                  # lowest impurity wins
    return best

def build(X, y, depth, max_depth):
    if gini(y) == 0 or depth == max_depth:
        return int(np.bincount(y, minlength=2).argmax())        # leaf = majority class
    f, thr, _ = best_split(X, y); L = X[:, f] <= thr
    return (f, thr, build(X[L], y[L], depth+1, max_depth),
                    build(X[~L], y[~L], depth+1, max_depth))     # recurse both sides

print("root split:", best_split(X, y)[:2])

from scratch
root split: (1, -0.44)      # feature 1 = y, threshold -0.44
depth-1 acc: 0.81    depth-2 acc: 1.00

scikit-learn builds the very same tree, just faster:

clf = DecisionTreeClassifier(max_depth=2).fit(X, y)
print(round(clf.score(X, y), 2))    # 1.00

Using trees well

Trees are wonderfully practical: they need no feature scaling, handle numeric and categorical features together, are robust to outliers, and a shallow one is fully interpretable — you can read the rules off it. The catch is variance: deep single trees overfit and their axis-aligned cuts struggle with diagonal boundaries. In practice, keep individual trees shallow for interpretability, and reach for random forests or gradient boosting when you want accuracy.

Centres and their followers

k-means is unsupervised: there are no labels, just points, and we want to discover structure. We propose k cluster centres (centroids) and declare that every point belongs to whichever centre is closest. A good clustering is one where points sit tightly around their centre. The trick is that we don't know the centres or the memberships in advance — so we bootstrap them from each other.

Points with no labels

Here are 60 points with no class information at all. Our eyes see three blobs immediately; the job is to make an algorithm find them.

Tight clusters = low inertia

k-means measures a clustering by its inertia — the total squared distance from each point to its assigned centre, also called the within-cluster sum of squares:

Lower is tighter. Two facts make the whole algorithm fall out: for fixed centres, each point lowers J most by joining its nearest centre; and for fixed memberships, the centre that minimises J for a group is exactly its mean. k-means just alternates those two optimal moves.

Assign, then update

This back-and-forth is Lloyd's algorithm:

Each step can only lower J or leave it unchanged, and there are finitely many ways to group the points — so the loop must stop. It halts when an assign step changes nothing.

Six points on a line

Take six points on a number line — 1, 2, 3, 8, 9, 10 — and ask for k = 2 clusters, starting from a deliberately poor guess.

The loop, running

Three centres start on random points. Click Assign to colour each point by its nearest centre, then Update to slide the centres to their cluster means — or just Step through both, or Run to the finish. Watch the inertia fall each round until the centres stop moving. Hit New init to drop fresh centres and run it again.

It can settle for less

Each step lowers inertia, but only to the nearest valley — not necessarily the deepest one. A poor starting placement can leave two centres splitting one true blob while a third covers two. Press New init a few times above: most runs find the three blobs, but some settle on a worse arrangement with higher inertia. The standard fixes are to run k-means several times from different starts and keep the lowest-inertia result, and to seed cleverly with k-means++, which spreads the initial centres out instead of dropping them at random.

Below is the same data and the same algorithm from two different starts. Almost every start finds the three blobs (left); an unlucky one leaves one centre stranded between blobs while two centres share the rest (right) — a valid stopping point with far higher inertia.

The elbow

More clusters always lower inertia — at k = n every point is its own centre and inertia is zero — so we can't just minimise it. Instead we look for the elbow: the k after which adding clusters barely helps. Slide k and watch both the curve and the clustering; the bend at k = 3 is where the real structure is captured and extra centres only start splitting blobs that were already tight.

Straight borders only

Because a point joins its nearest centre, the border between two clusters is always a straight line — the perpendicular bisector of their centres. That makes k-means assume clusters are roughly round, similar in size, and convex. Give it a core wrapped in a ring and it can't comply: no straight cut separates inside from outside, so it slices through both. For shapes like these, density- or connectivity-based methods (DBSCAN, spectral clustering) are the right tools.

From scratch, then scikit-learn

Lloyd’s algorithm is a dozen lines of NumPy — assign each point to the nearest centre, move each centre to its mean, repeat — wrapped in a few random restarts to dodge bad local minima:

# --- k-means from scratch ---
def kmeans(X, k, seed):
    rng = np.random.RandomState(seed)
    C = X[rng.choice(len(X), k, replace=False)]            # init centres from points
    for _ in range(30):
        d = ((X[:, None, :] - C[None, :, :])**2).sum(2)     # distances to centres
        a = d.argmin(1)                                      # (1) assign nearest
        C = np.array([X[a == j].mean(0) for j in range(k)])  # (2) move to means
    return C, a, ((X - C[a])**2).sum()                       # centres, labels, inertia

best = min((kmeans(X, 3, s) for s in range(10)), key=lambda t: t[2])
print("best inertia:", round(best[2], 1))

from scratch
best inertia: 44.6

scikit-learn’s KMeans, with its own 10 restarts, finds the identical optimum:

km = KMeans(n_clusters=3, n_init=10).fit(X)
print(round(km.inertia_, 1))     # 44.6

Using k-means well

Always scale your features — k-means is built on Euclidean distance, so a large-range feature will dominate. Run it with several initialisations (scikit-learn's n_init) and let k-means++ seed them. Use the elbow or the silhouette score to choose k, and remember the model's assumptions: round, comparable, convex clusters. When those hold, k-means is fast, simple, and scales to huge datasets; when they don't, reach for a density- or model-based method.

A better set of axes

Your features come with arbitrary axes — x and y, or a hundred sensor channels. PCA looks for a new, rotated set of axes lined up with the data itself: the first points along the direction of greatest variance, the next along the greatest remaining variance perpendicular to it, and so on. Re-expressed in these axes, most of the interesting variation lives in the first few, so the rest can be discarded with little loss.

A tilted cloud

Here are 80 points in two features. They're correlated — the cloud is an ellipse tilted at roughly 35°. Neither x nor y alone captures that tilt; the most informative direction runs diagonally along the cloud.

Subtract the mean

PCA is about spread, not location, so the first step is always to centre the data — subtract the mean of each feature so the cloud sits at the origin. Everything that follows measures how points vary around that centre.

How spread out, which way?

Pick a unit direction u. Project every centred point onto it — the projection of x is the scalar u·x — and measure the variance of those projections. That variance has a compact form involving the covariance matrix Σ:

Different directions give different spreads. PCA asks: which unit direction makes this as large as possible?

Spread in every direction at once

For centred data X (one row per point), the covariance matrix packs all the variances and co-variances together:

The diagonal holds each feature's own variance; the off-diagonal holds how they move together. A tilted cloud has a large off-diagonal term — that's exactly the correlation PCA will straighten out.

The maximum-variance directions

Maximising uᵀΣu subject to ‖u‖ = 1 is a classic constrained problem, and its solution is clean: the maximising directions are the eigenvectors of Σ, and the variance captured by each is its eigenvalue.

The eigenvector with the largest eigenvalue is the first principal component (PC1); the next-largest, perpendicular to it, is PC2; and so on. Diagonalising Σ is PCA.

Four points, two components

Take four centred points and find their principal components with pencil and paper.

Hunting the maximum

Sweep a direction through the cloud below. The grey spokes drop each point onto the line, and the strip underneath plots the variance of those projections as the angle turns. It rises and falls in a smooth wave with one clear peak — and that peak direction, capturing the most variance, is exactly PC1.

Dropping a dimension

To reduce to one dimension, keep only each point's PC1 coordinate and throw PC2 away. Reconstructing from that single number lands every point on the PC1 line; the leftover gaps are the reconstruction error, and they add up to exactly the variance you discarded (PC2's eigenvalue). Toggle how many components you keep:

Explained variance

Each eigenvalue, divided by their total, is the fraction of variance that component explains. Plotting them gives a scree chart; you keep components until the cumulative total crosses a threshold (often 90–95%) or the bars flatten into an "elbow." Here PC1 alone already carries most of the spread.

From scratch, then scikit-learn

Under the hood, PCA is only centre → covariance → eigen-decomposition. Here it is by hand in NumPy, with the one-line scikit-learn version below confirming it:

# --- PCA from scratch ---
Xc   = X - X.mean(0)                  # 1. centre the data
cov  = Xc.T @ Xc / len(Xc)            # 2. covariance matrix
vals, vecs = np.linalg.eigh(cov)      # 3. eigen-decompose it
order = vals.argsort()[::-1]          # 4. largest variance first
vals, vecs = vals[order], vecs[:, order]

print("explained ratio:", vals / vals.sum())
print("PC1:", vecs[:, 0])
Xrec = Xc @ vecs[:, :1] @ vecs[:, :1].T            # keep PC1, reconstruct
print("recon MSE:", np.mean(np.sum((Xc - Xrec)**2, 1)))

from scratch
explained ratio: [0.885 0.115]
PC1: [-0.806 -0.592]      # sign is arbitrary — same axis as sklearn
recon MSE: 0.58           # exactly the dropped eigenvalue

scikit-learn wraps those four steps and lands in the same place:

pca = PCA().fit(X)
print(pca.explained_variance_ratio_)   # [0.885 0.115]
print(pca.components_[0])               # [0.806 0.592]

Using PCA well

Because PCA chases variance, standardise features first when they're on different scales, or one large-range feature will masquerade as the main component. Remember it's purely linear and unsupervised: it finds directions of spread, which aren't always the directions that matter for a downstream label, and it can't unfold curved structure (for that, reach for kernel PCA, t-SNE, or UMAP). Used well, PCA is the workhorse for compression, denoising, visualising high-dimensional data in 2-D, and decorrelating features before another model.

You are like your neighbours

To classify a new point, k-NN looks at the k training examples nearest to it and lets them vote. That's the whole algorithm. It's lazy and instance-based: nothing is learned up front — the model simply stores every training point and does all its work at prediction time, measuring distances and counting votes.

Two overlapping classes

Fifty points in two classes that overlap in the middle — so the boundary between them is genuinely fuzzy, and the choice of k will matter.

What "nearest" means

Nearness needs a metric. The usual choice is straight-line (Euclidean) distance:

Other metrics — Manhattan (sum of absolute differences), cosine (angle between vectors) — are swappable, and the choice changes which points count as neighbours. Whatever the metric, every prediction is just a distance computation away.

Count the k closest

Sort the training points by distance to the query, take the closest k, and predict the class that appears most among them. With k = 1 you simply copy the single nearest point's label; with larger k you average over a neighbourhood. k is usually chosen odd for two classes so a tie can't happen, and votes can be distance-weighted so closer neighbours count more.

One query, three values of k

Classify the query point q = (4, 3.5) against six labelled points — and watch the answer change as k grows.

Regions of influence

Colour every spot on the plane by how k-NN would classify it and you get the decision regions below. Click anywhere to drop a query point — its k nearest neighbours light up and cast their votes. Slide k and watch the boundary change shape.

From jagged to smooth

k is the dial between overfitting and underfitting. At k = 1 the boundary is jagged and every noisy point carves out its own little island — high variance, it memorises. As k grows the boundary smooths and steadies, until a very large k washes out real structure and just predicts the overall majority — high bias. On this data, leave-one-out accuracy bears that out: k = 1 → 78%, k = 5 → 82% (the sweet spot), k = 15 → 76%. The right k is found by cross-validation, not by eye.

Distance is unforgiving

Because everything rests on distance, k-NN is acutely sensitive to feature scale: a feature measured in thousands will dominate one measured in fractions, no matter which is informative. Always standardise features first. Two refinements help further — weighting each neighbour's vote by 1/distance so closer points count more, and choosing a metric that suits the data. None of these change the core idea; they only change how distance is measured and tallied.

See it for yourself. Class here is decided entirely by feature x₁ (small scale, horizontal axis); x₂ is irrelevant noise on a 30× larger scale (vertical). Toggle the distance and watch which neighbours get chosen — and the accuracy:

Where "nearest" stops meaning much

k-NN's Achilles heel is high dimensions. As you add features, points spread out and all pairwise distances drift toward the same value — the nearest neighbour is barely closer than the farthest, so "nearest" carries little information. The plot tracks the contrast (farthest − nearest)/nearest as dimensions grow: it collapses fast.

This is why k-NN shines in 2–20 well-chosen features and struggles in raw high-dimensional spaces — often you reduce dimensions first (say with PCA) before reaching for it.

From scratch, then scikit-learn

There's nothing to fit — prediction is distance, sort, vote. Here it is in five lines of NumPy, with scikit-learn predicting identically below:

# --- k-NN from scratch ---
def knn_predict(Xtr, ytr, Q, k):
    out = []
    for q in Q:
        d   = ((Xtr - q)**2).sum(1)         # squared distance to every point
        idx = np.argsort(d)[:k]             # indices of the k nearest
        out.append(np.bincount(ytr[idx]).argmax())   # majority vote
    return np.array(out)

# leave-one-out accuracy (train accuracy would read 100% at k=1 — a lie)
def loo(k):
    ok = 0
    for i in range(len(X)):
        m = np.arange(len(X)) != i
        ok += knn_predict(X[m], y[m], X[i:i+1], k)[0] == y[i]
    return ok / len(X)

for k in (1, 5, 15): print(f"k={k}: LOO acc {loo(k):.2f}")

from scratch
k=1:  LOO acc 0.78        # jagged — overfits the noise
k=5:  LOO acc 0.82        # the sweet spot
k=15: LOO acc 0.76        # too smooth — underfits

scikit-learn does the same distance-and-vote, and predicts identically:

from sklearn.neighbors import KNeighborsClassifier
clf = KNeighborsClassifier(n_neighbors=5).fit(X, y)
print(clf.predict(Q))                       # [0 0 1 1 0]  — same as knn_predict(X, y, Q, 5)

Using k-NN well

k-NN is a superb baseline — it makes no assumptions about the shape of the boundary and a handful of lines implements it. The trade-offs: it must store all the data and search it at every prediction, so it's slow at scale (real systems use KD-trees, ball-trees, or approximate-nearest-neighbour indexes); it demands scaled features; and it degrades in high dimensions. Choose k by cross-validation, standardise first, and consider distance weighting.

Pick the most probable class

Given a new point, ask: which class most probably produced it? Naive Bayes answers with Bayes' rule, scoring every class by prior × likelihood and choosing the largest. The word "naive" is the trick that makes it tractable — it pretends every feature is independent of the others, so their likelihoods simply multiply.

Two classes in the plane

Sixty points, two classes, each a blob with its own centre and spread. Naive Bayes will model each class as a pair of bell curves — one per axis.

Flip the conditional

We want P(class | features) but it's easy to model P(features | class) — how a class generates data. Bayes' rule flips one into the other:

The denominator P(x) is the same for every class, so for choosing a class we can drop it and just compare the numerators:

Pretend the features don't talk

The likelihood P(x | c) of a whole feature vector is hard — features interact. Naive Bayes assumes they're conditionally independent given the class, so the joint likelihood factorises into a product of one-feature likelihoods:

This is almost never literally true — that's why it's "naive." But it turns an impossible joint density into a handful of easy 1-D densities, and the resulting classifier is often excellent anyway (chapter 8 explains why).

A bell curve per feature, per class

For continuous features the usual choice is a Gaussian: for each class and each feature, fit a mean and variance from the training data, then read the likelihood off the bell curve.

Fitting is trivial — just the per-class mean and variance of each feature. To classify, plug the point into every class's bell curves, multiply by the prior, and compare. In practice we add the logs instead of multiplying the tiny numbers, which turns the product into a sum and avoids underflow.

Drag the query across the two bell curves below — class A (μ = 2) and class B (μ = 5) — to read each class’s likelihood. Slide the prior and watch it tip the posterior and shift the decision point (the amber line, where the two classes are equally probable):

One feature, one query

The cleanest case: one feature, two classes fitted as bell curves — class A with mean 2, class B with mean 5, both with σ = 1. Classify the query x = 3, and watch the prior change the answer.

Bell curves carve the plane

Each class becomes a pair of Gaussians, drawn below as contour ellipses (the naive assumption makes them axis-aligned — no diagonal tilt). The shaded regions show which class wins the prior × likelihood contest at every point. Click anywhere to drop a query and read its posterior probabilities. Because the two classes have different variances, the boundary bends — it's a curve, not a line.

Wrong probabilities, right decisions

The independence assumption is usually false, so the probabilities Naive Bayes reports are often poorly calibrated — too confident, pushed toward 0 or 1. Yet its classifications are frequently spot-on. The reason: to pick the right class you don't need accurate probabilities, only the correct winner. Correlated features double-count evidence and inflate the magnitudes, but they usually inflate them in the same direction the true class points — so the argmax survives even when the numbers don't. Trust the label; distrust the exact percentage.

One template, three flavours

Swap the per-feature likelihood to fit the data type. Gaussian NB (above) for continuous features. Multinomial NB for counts — the classic spam filter, where features are word frequencies. Bernoulli NB for binary present/absent features. The recipe — prior × ∏ likelihoods, take the argmax — never changes.

One catch with counts: if a word never appeared in the training emails of a class, its likelihood is 0, which annihilates the entire product. Laplace smoothing fixes this by adding a pseudo-count of 1 to every tally, so nothing is ever impossible — just unlikely.

From scratch, then scikit-learn

Fitting is one pass for the per-class means, variances and priors; predicting is a sum of log-likelihoods. Here it is in NumPy, with scikit-learn's GaussianNB agreeing on every point:

# --- Gaussian Naive Bayes from scratch ---
def fit_gnb(X, y):
    stats = {}
    for c in np.unique(y):
        Xc = X[y == c]
        stats[c] = (Xc.mean(0), Xc.var(0), len(Xc)/len(X))   # mean, var, prior
    return stats

def predict_gnb(stats, Q):
    out = []
    for q in Q:
        score = {}
        for c, (mu, var, prior) in stats.items():
            # log posterior = log prior + sum of log-Gaussian likelihoods
            score[c] = np.log(prior) - 0.5*np.sum(np.log(2*np.pi*var) + (q-mu)**2/var)
        out.append(max(score, key=score.get))
    return np.array(out)

stats = fit_gnb(X, y)
print("train accuracy:", (predict_gnb(stats, X) == y).mean())

from scratch
class 0:  mean [-1.46 -0.31]  var [0.58 0.86]
class 1:  mean [ 1.14  0.83]  var [1.29 0.65]
train accuracy: 0.917

scikit-learn computes the very same statistics and agrees on every point of a test grid:

from sklearn.naive_bayes import GaussianNB
gnb = GaussianNB().fit(X, y)
print((predict_gnb(stats, grid) == gnb.predict(grid)).mean())   # 1.0 — identical everywhere
print(gnb.score(X, y))                                          # 0.917

Using Naive Bayes well

Naive Bayes is the speed champion: training is a single pass to count means and variances, prediction is a few additions, and it copes effortlessly with thousands of features — which is why it remained the default text and spam classifier for decades. It needs very little data to get going and never overfits in the dramatic way a deep tree can. The price is the independence assumption and the over-confident probabilities it produces; reach for it as a fast, strong baseline, especially on high-dimensional or text data, and calibrate the scores if you need to trust them.

Many wrong trees, one right answer

Ask one expert and you get one opinion, biases and all. Ask a thousand independent experts and average them, and the individual errors — pointing in random directions — largely cancel, leaving the shared signal. A random forest applies this to decision trees: build a crowd of them, each trained a little differently, and aggregate their votes. The maths is the bias–variance decomposition: averaging many high-variance, low-bias predictors keeps the low bias while slashing the variance.

One tree is twitchy

A fully-grown tree drives its training error to zero by carving the plane into ever-smaller boxes — including boxes drawn purely around noise. The result is a jagged boundary that nails the training set and stumbles on new data: classic high variance. Resample the data slightly and a different set of splits wins, producing a visibly different tree. That instability is precisely what averaging is built to fix.

Trick one: bootstrap the data

Bagging = bootstrap aggregating. For each tree, draw a training set the same size as the original but sampled with replacement — some points appear several times, others not at all. Each tree therefore sees a slightly different world and grows differently. A nice side effect: each bootstrap sample leaves out about 1/e ≈ 37% of the points (the out-of-bag set), which we'll reuse as free validation data. Click resample to see one bootstrap draw — filled points were selected, hollow ones left out.

Trick two: hide most of the features

Bagging alone isn't enough: if one feature is very predictive, every tree splits on it first and they all end up similar — correlated trees don't cancel. The forest's signature trick is to show each split only a random subset of features (typically √d of them). Forced to use different features, the trees decorrelate, and decorrelated errors cancel far better when averaged.

That single change — choosing splits from a random handful of features — is what separates a random forest from plain bagged trees.

Compare the two on the same data — toggle between bagging only (every split sees both features) and a real forest (each split sees one random feature). Watch the accuracy, and how often the trees agree with one another:

Aggregate the crowd

To predict, run the point down every tree and combine: majority vote for classification, average for regression. No single tree needs to be good — it only needs to be right slightly more often than chance and wrong in its own idiosyncratic way. The aggregate is what's accurate.

From staircase to smooth

The true boundary here is the diagonal x + y = 0 (dashed), with 12% of labels flipped as noise. A single deep tree can only approximate that diagonal with a jagged axis-aligned staircase, and it wraps boxes around the noise. Add trees and watch the forest's averaged boundary settle toward the true diagonal — and its accuracy climb past the lone tree's.

Why more trees help

Plot test accuracy against the number of trees and the shape tells the whole story: a steep early climb as the first few trees cancel each other's mistakes, then a plateau — beyond a point, extra trees barely move the needle (and never cause overfitting, they only stabilise the estimate). The dashed line is the single deep tree; the forest clears it quickly and stays there.

Free validation, no holdout

Because every tree skips the ~37% of points outside its bootstrap sample, each point has a committee of trees that never saw it during training. Predicting every point using only those trees gives the out-of-bag (OOB) error — an honest estimate of test performance computed for free, no separate validation split required. It's one of the quiet conveniences that makes forests so practical.

Which features mattered

Every split reduces impurity by some amount. Add up the impurity each feature removes, averaged across all trees, and you get a feature importance ranking — a cheap, useful read on which inputs the forest leaned on. It comes out of training at no extra cost (be aware it can inflate the importance of high-cardinality features; permutation importance is a more robust alternative). Interpretability is the main thing you trade away versus a single tree: a forest of hundreds of trees is no longer something you can read off as a flowchart.

From scratch, then scikit-learn

A forest is a loop over a tree: bootstrap, restrict the features, grow, repeat — then vote. Reusing the build/predict tree from the trees dive (now passing a feature subset to each split), the whole ensemble is a few lines, and it matches scikit-learn:

# --- random forest from scratch (build/predict come from the trees dive) ---
def forest(X, y, n_trees=100, max_features=1):
    trees = []
    for _ in range(n_trees):
        i = np.random.randint(0, len(X), len(X))     # bootstrap sample (with replacement)
        trees.append(build(X[i], y[i], max_features)) # each split sees a random feature subset
    return trees

def forest_predict(trees, Q):
    votes = np.array([[predict(t, x) for t in trees] for x in Q])
    return (votes.mean(1) >= 0.5).astype(int)          # majority vote

f = forest(Xtr, ytr, n_trees=100)
print("single deep tree:", tree_acc, " | forest:", (forest_predict(f, Xte) == yte).mean())

from scratch
single deep tree: 0.76   |   forest: 0.84      # the crowd beats the expert

scikit-learn's RandomForestClassifier bundles bootstrap + random features + voting, and lands in the same place:

from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=100, max_features=1).fit(Xtr, ytr)
print(rf.score(Xte, yte))      # 0.84 — matches the from-scratch forest

Using forests well

Random forests are the workhorse default for tabular data: accurate out of the box, hard to overfit, indifferent to feature scaling, comfortable with mixed numeric and categorical inputs, and equipped with free OOB error and importances. The cost is interpretability and prediction speed — hundreds of trees take memory and time, and you lose the single tree's readable logic. More trees is always safe (it never overfits, just plateaus); the knob that matters most is the feature-subset size. When forests aren't enough, boosting — growing trees in sequence, each fixing the last one's mistakes — is the usual next step.

Each tree fixes the last mistake

Start with a dumb prediction — say, the average. It's wrong everywhere; record by how much (the residuals). Now train a small tree not on the original target but on those residuals, so it learns the pattern in what's left over. Add a fraction of it to your prediction. The errors shrink. Repeat hundreds of times, each tree chipping away at what remains, and the running total converges on the truth.

Parallel crowd vs relay team

Both are tree ensembles, but they attack opposite problems. A random forest grows deep, independent trees on bootstrap samples and averages them — the trees are interchangeable and the goal is to cut variance. Boosting grows shallow, dependent trees in sequence — each one's job is defined by its predecessors' mistakes — and the goal is to cut bias, turning underfit stumps into an accurate whole.

The core loop

For regression with squared error, one round is dead simple. Hold the current prediction F. Compute the residual for every point. Fit a small tree to those residuals. Add a shrunk copy of it to F. That's it.

The first tree explains the biggest, coarsest errors; later trees pick at finer and finer leftovers. Because each tree only ever sees residuals, it's always working on the part of the signal nobody has captured yet.

Shrink every step

That ν — the learning rate (or shrinkage) — scales down each tree's contribution, typically to 0.01–0.3. Why deliberately cripple every step? Because taking small steps and using many of them generalises far better than a few greedy leaps: no single tree can dominate, and the ensemble approaches the target gently instead of lurching onto noise. Learning rate and number of trees trade off directly — halve ν and you'll want roughly twice as many trees.

The three curves are the same data boosted at three learning rates. A big rate (ν = 1.0) plunges fastest but bottoms out higher and overfits within a few trees; a small rate (ν = 0.05) descends slowly but reaches the lowest, flattest minimum — the safest fit. Each dot marks where that rate's test error bottoms out.

A sum of small corrections

Stack the rounds and the model is just a weighted sum of trees — an additive model built greedily, one stage at a time (this is "forward stagewise" fitting):

F₀ is the initial guess (the mean for regression). Every h is a shallow tree fit to the running residual. Nothing is ever re-fit — each stage is frozen once added.

Watch the curve assemble

Forty-two noisy points trace a sine wave. The forest's prediction starts as a flat line at the mean (stage 0). Drag the stage slider and watch each shallow tree bend the fit a little more toward the data — the grey residual lines (how far each point sits from the current fit) shrink as the corrections accumulate.

When more trees hurt

Here is boosting's crucial difference from a forest. Track the error as stages add: training error falls relentlessly toward zero — given enough trees, boosting will memorise every point. But test error traces a U: it drops fast, bottoms out (here around stage 7), then climbs as the later trees start fitting noise instead of signal. A forest would simply plateau; boosting overfits. The fix is to stop at the bottom of the U — early stopping on a validation set.

Residuals are a gradient

The name isn't decoration. For squared-error loss L = ½(y − F)², the derivative with respect to the prediction is ∂L/∂F = −(y − F) — exactly the negative residual. So fitting trees to residuals is gradient descent, performed in function space: each tree is a step in the steepest-downhill direction of the loss.

That reframing is the power of the method: swap in a different loss and the "residuals" become its negative gradient. Use the logistic loss and gradient boosting classifies; use the absolute or Huber loss and it shrugs off outliers. The loop never changes — only what counts as an error.

Classification, the same way

For classification the trees predict in log-odds space rather than directly in probabilities. F₀ starts at the data's base log-odds, the loss is logistic, and its negative gradient — the "residual" each tree fits — is simply (actual − predicted probability). Run the sum through a sigmoid at the end to recover a probability. Same additive loop, same shrinkage, same overfitting caution; only the loss and the final squashing differ.

From scratch, then scikit-learn

The whole method is a loop around a regression tree: predict, take the residual, fit a shallow tree to it, add a shrunk copy. Reusing a regression tree, the core is a dozen lines — and it matches scikit-learn's GradientBoostingRegressor to the digit:

# --- gradient boosting from scratch (squared-error regression) ---
def gradient_boost(X, y, n_stages=60, lr=0.3, depth=3):
    F0 = y.mean()                                  # start at the mean
    F  = np.full(len(y), F0)
    trees = []
    for _ in range(n_stages):
        residual = y - F                           # = negative gradient of squared loss
        h = fit_tree(X, residual, max_depth=depth) # shallow tree on the residual
        F = F + lr * predict(h, X)                 # add a shrunk copy
        trees.append(h)
    return F0, trees

def gb_predict(F0, trees, X, lr):
    return F0 + lr * sum(predict(h, X) for h in trees)

F0, trees = gradient_boost(Xtr, ytr, n_stages=60, lr=0.3)
print("test MSE:", ((gb_predict(F0, trees, Xte, 0.3) - yte)**2).mean())

from scratch
test MSE: 0.094

scikit-learn wraps the identical stagewise loop and lands on the same number:

from sklearn.ensemble import GradientBoostingRegressor
gbr = GradientBoostingRegressor(n_estimators=60, learning_rate=0.3, max_depth=3).fit(Xtr, ytr)
print(((gbr.predict(Xte) - yte)**2).mean())     # 0.094 — same as from scratch

Using boosting well

Gradient boosting is, for most tabular problems, the most accurate model you can reach for — the libraries XGBoost, LightGBM, and CatBoost are tuned, regularised, parallelised versions that win competition after competition. The cost is care: unlike a forest, it will overfit, so you tune the learning rate, tree depth (usually 3–8), and number of trees together, leaning on early stopping against a validation set. Start with a small learning rate and many trees, keep the trees shallow, and stop when validation error turns up.

Clusters with soft edges

Imagine the data was generated by a few Gaussian "blobs": pick a blob at random, then draw a point from its bell curve. A Gaussian mixture model reverses that story — given the points, recover the blobs. Crucially the assignment is soft: a point near the border between two clusters belongs, say, 70% to one and 30% to the other. That probability is called a responsibility, and it's the whole difference from k-means.

Loosening every assumption

k-Means is actually a special, stripped-down GMM. Generalise it on three fronts and you get the full model:

In fact, k-means is the limit of a GMM with tiny equal spherical covariances and hard assignments — the same loop, with the probabilities collapsed to 0/1.

A weighted sum of bells

The density is literally a weighted sum of K Gaussians. Each has a mixing weight π_k (how big a slice of the data it owns, with the π's summing to 1), a mean μ_k, and a covariance Σ_k (its shape and tilt):

Fitting means finding the π's, μ's, and Σ's that make the observed data most likely. The snag: if we knew which cluster each point came from, the fit would be trivial — but we don't. EM breaks the deadlock.

Guess the memberships

Hold the current Gaussians fixed and ask, for each point, the posterior probability it came from cluster k — its responsibility γ. It's just Bayes' rule over the clusters:

A point sitting deep inside one ellipse gets responsibility ≈ 1 there; a point on the fence splits its responsibility. This is the Expectation step — soft assignment given the current guess.

Refit each Gaussian

Now hold the responsibilities fixed and refit every Gaussian — but each point contributes to a cluster only in proportion to its responsibility there. Each cluster's effective size is N_k = Σ_i γ_ik, and:

These are exactly the ordinary mean, covariance, and proportion — only weighted by responsibility. This is the Maximization step: the parameters that maximise the likelihood given the soft assignments.

Alternate until it settles

Neither step can run without the other's output, so EM simply alternates: E (responsibilities from current Gaussians) → M (Gaussians from current responsibilities) → repeat. Each full round is guaranteed never to decrease the data's log-likelihood, so the process climbs steadily to a (local) maximum and stops. It's the same assign-then-update rhythm as k-means — only soft.

Here is that climb for a full run on the data above — the log-likelihood after each E–M round. It rises at every step (it can never dip) and flattens as the fit settles into a local maximum:

Watch the ellipses snap on

Two Gaussians start as round guesses in the wrong places. Press E-step to colour each point by its current responsibility (a blend of the two cluster colours), then M-step to refit the ellipses to those soft assignments — or just hit Run. Watch the ellipses rotate and stretch onto the data while the log-likelihood climbs and levels off.

Where circles can't cope

Here are two long, parallel, diagonal clusters. Because k-means measures plain distance to a centre, it carves the plane with a straight line and slices across both cigars — getting barely 60% right. A Gaussian mixture, free to stretch and tilt its ellipses, locks onto the two diagonals and separates them almost perfectly. Toggle between them:

The fine print

How many clusters? Because a GMM is a probabilistic model, you can compare different k by their likelihood, penalised for complexity — the BIC or AIC — and pick the elbow. A few hazards to respect: EM finds only a local optimum, so it's run from several random starts (usually k-means is used to initialise it); a Gaussian can collapse onto a single point, sending its variance to zero and the likelihood to infinity, which a little regularisation prevents; and you can constrain the covariance (spherical, diagonal, or full) to trade flexibility for fewer parameters. More clusters always fit the training data better, so the penalty term is what stops you from using one Gaussian per point.

From scratch, then scikit-learn

EM is two short steps in a loop: responsibilities, then weighted refit. Here it is in NumPy, with scikit-learn's GaussianMixture recovering the same clusters and log-likelihood:

# --- Gaussian mixture via EM, from scratch ---
def em(X, K, iters=60):
    n = len(X)
    mu = X[np.random.choice(n, K, replace=False)]      # init means from points
    S  = [np.cov(X.T) for _ in range(K)]
    pi = np.ones(K) / K
    for _ in range(iters):
        # E-step: responsibilities (Bayes over clusters)
        R = np.array([pi[k]*mvn_pdf(X, mu[k], S[k]) for k in range(K)]).T
        R /= R.sum(1, keepdims=True)
        # M-step: weighted mean, covariance, weight
        Nk = R.sum(0)
        for k in range(K):
            mu[k] = (R[:, k:k+1] * X).sum(0) / Nk[k]
            d = X - mu[k]
            S[k] = (R[:, k:k+1] * d).T @ d / Nk[k]
        pi = Nk / n
    return mu, S, pi

mu, S, pi = em(X, 2)
print("means:", mu.round(2), " weights:", pi.round(2))

from scratch
means: [[-1.35  0.3 ] [ 1.34  0.31]]   weights: [0.38 0.62]
log-likelihood: -324.5   (rose every iteration)

scikit-learn runs the same EM (from several starts) and converges to the same fit:

from sklearn.mixture import GaussianMixture
gm = GaussianMixture(n_components=2, covariance_type="full").fit(X)
print(gm.means_.round(2))                 # same two centres
print(gm.score(X) * len(X))               # total log-likelihood ≈ -324.5

Using mixtures well

Gaussian mixtures shine when clusters genuinely overlap or aren't round, when you need soft memberships (a point can be 60/40), or when you want a generative density model you can sample from and score new points against — anomaly detection falls out naturally as "low probability under the mixture." The costs are the usual EM caveats: local optima, sensitivity to k and to initialisation, and the assumption that clusters are roughly Gaussian. Initialise with k-means, run several restarts, pick k by BIC, and regularise the covariances.

Two ways to be wrong

Suppose there's a true pattern and you only ever see noisy samples of it. Train a model and its expected error on fresh data breaks into three pieces: bias (how far your model's average prediction sits from the truth — error from wrong assumptions), variance (how much your prediction jumps around as the training set changes — error from over-sensitivity), and noise (the irreducible randomness no model can remove).

Stubborn vs jittery

Bias is the error baked in by a model too simple to represent the truth — a straight line trying to fit a curve will be wrong in the same direction no matter which sample you give it. It's stubborn and systematic. Variance is the error from a model so flexible it memorises the quirks of its training set — give it a different sample and it draws a wildly different function. It's jittery and sample-specific.

Where it comes from

For squared-error loss, take the expected error at a point x over all possible training sets. Add and subtract the average prediction f̄(x) = E[f̂(x)], and the cross-term vanishes, leaving three non-negative pieces:

Bias² measures how the average model misses the truth; variance measures how individual models scatter around their own average. Crucially, the same lever that lowers one tends to raise the other.

Concretely: fix one input x₀ (here the peak of the sine, where the truth is 1.0) and refit on many noisy samples. Each fit makes a slightly different prediction there — a cloud of predictions. How far the cloud's centre sits from the truth is the bias; how spread out the cloud is, is the variance. Slide the degree and watch the two trade off at this single point:

One knob, opposite effects

Model complexity is the master knob. Turn it up — more polynomial terms, deeper trees, smaller k in k-NN — and bias falls (the model can represent more) while variance rises (it also fits more noise). Turn it down and the reverse happens. Because total error is their sum, it traces a U: high on the left from bias, high on the right from variance, with a sweet spot in between. Finding that spot is what model selection is.

Watch the fit go from stiff to frantic

Eighteen noisy points (grey) are sampled from a true sine wave (dashed). Fit a polynomial and slide its degree. At degree 1 the line is too stiff — it misses the curve everywhere (high bias). Around degree 3 it tracks the truth. Push to degree 9–10 and it starts threading through individual points, wiggling between them (high variance). The right panel plots the bias²/variance/total decomposition with a marker at your current degree:

The fan of fits

Variance is abstract until you watch it. Below, the same model is fit to many different noisy samples of the same sine, and every resulting fit is drawn together. At low degree the fits bunch tightly — but all of them miss the curve the same way (that consistent miss is bias). Crank the degree and the fits fan out wildly; each one chases its own sample's noise (that spread is variance). The bold line is their average.

Diagnosing from data size

The two ailments leave different fingerprints as you feed in more data. Plot training and validation error against training-set size. A high-bias model: both curves climb/fall to meet at a high plateau with almost no gap — adding data won't save a model that's too simple. A high-variance model: a wide gap, training error near zero while validation error sits well above it — more data steadily narrows the gap. Toggle the two:

The whole lab, one lens

Almost every technique you've met is a move in this game:

Measuring it from scratch

You can estimate bias and variance empirically: refit the model on many resampled datasets, then at each point compare the average prediction to the truth (bias) and measure the spread of predictions (variance).

# --- empirical bias-variance decomposition ---
def bias_variance(degree, n_sets=300):
    preds = []
    for _ in range(n_sets):
        ys = f(X) + np.random.randn(len(X)) * SIGMA   # a fresh noisy sample
        c  = np.polyfit(X, ys, degree)                # refit the model
        preds.append(np.polyval(c, grid))
    preds = np.array(preds)
    mean  = preds.mean(0)                              # the average model
    bias2 = ((mean - f(grid))**2).mean()              # avg miss from truth
    var   = preds.var(0).mean()                        # spread of the fits
    return bias2, var, bias2 + var + SIGMA**2

for d in (1, 3, 9):
    print(d, [round(v, 3) for v in bias_variance(d)])

bias², variance, total
1 [0.206, 0.004, 0.242]      # underfit: bias dominates
3 [0.005, 0.008, 0.045]      # sweet spot: both small
9 [0.0,   0.023, 0.055]      # overfit: variance dominates

The model here is np.polyfit; in scikit-learn the same fit is make_pipeline(PolynomialFeatures(d), LinearRegression()). The decomposition itself is a measurement, not a model — it's how you'd quantify which ailment you have.

Living with the tradeoff

You never see bias and variance directly — you see their sum on a validation set. So diagnose by symptom: high training error means high bias (the model can't even fit what it's seen — go more complex, add features, reduce regularization); a large train-to-validation gap means high variance (go simpler, regularise, bag, or gather more data). Use cross-validation to find the bottom of the U rather than trusting training error, which only ever points toward more complexity.

Accuracy isn't enough

How you score a model decides what it becomes — optimise the wrong metric and you'll get a model that's excellent at the wrong thing. Accuracy (fraction correct) hides two very different mistakes: flagging something that isn't there, and missing something that is. In spam filtering, a false alarm buries a real email; in disease screening, a miss costs a life. These errors rarely cost the same, so we need metrics that pull them apart.

Four outcomes

Every binary prediction lands in one of four cells, comparing what's true against what you predicted. This 2×2 table is the source of every classification metric:

Accuracy is just (TP+TN)/total. Everything more useful comes from looking at the off-diagonal — the two ways of being wrong — separately.

Two questions, two metrics

Precision asks: of everything I flagged positive, what fraction really was? It punishes false alarms. Recall (sensitivity) asks: of all the real positives, what fraction did I catch? It punishes misses.

They pull in opposite directions. Flag everything and recall is perfect but precision craters; flag only your single most confident case and precision is perfect but recall is near zero. Which matters more is a domain decision: a spam filter guards precision (don't trash real mail); a cancer screen guards recall (don't miss a tumour).

One number, when you need it

To collapse the pair into a single score, use the F1 — the harmonic mean of precision and recall. The harmonic mean (not the plain average) punishes imbalance: score 1.0 precision and 0.0 recall and the plain average is 0.5, but F1 is 0, correctly calling it a failure.

F1 is high only when both are high. When the costs of the two errors differ, a weighted version (F_β) tilts toward precision or recall.

The surface below is F1 across every precision–recall combination — bright only near the top-right corner where both are high, and collapsing to zero along either edge, exactly where the plain average would still read up to 0.5. Drag the sliders and compare F1 against the plain average:

Slide the decision line

A classifier doesn't output a label — it outputs a score (often a probability). You get a label by thresholding it. Below, 60 positives (teal) and 60 negatives (grey) are placed by their score; everything to the right of the line is predicted positive. Drag the threshold and watch the confusion matrix and metrics move: raise it and precision climbs while recall falls; lower it and the reverse. There is no single "right" threshold — only the one matching your costs.

Every threshold at once

Rather than pick one threshold, sweep them all and plot the result. The ROC curve traces the true-positive rate (recall) against the false-positive rate as the threshold moves from high to low. A perfect classifier hugs the top-left corner; random guessing is the diagonal. The single-number summary is the AUC — area under the curve — the probability the model scores a random positive above a random negative. The dot marks the threshold you set above.

The view that survives imbalance

Plot precision against recall across all thresholds and you get the PR curve. Its summary is average precision (AP), the area beneath it. Why have both? ROC counts true negatives, which flood the picture when negatives vastly outnumber positives and can make a weak model look strong. The PR curve ignores true negatives entirely, so it tells the truth on imbalanced problems — exactly the case we turn to next.

When accuracy lies

Make the positives rare — here just 10% of the data — and accuracy becomes a trap: a model that predicts "negative" for everything scores 90% while catching exactly zero positives. ROC barely notices (its AUC stays high because true negatives are abundant), but the PR curve collapses, and average precision drops sharply. Toggle the two regimes and compare the PR curves:

An honest estimate

One train/test split gives one noisy number — get unlucky with the split and you'll misjudge the model. k-fold cross-validation fixes this: partition the data into k equal folds, then train on k−1 and validate on the held-out fold, rotating so every fold serves as validation exactly once. Average the k scores for a stabler estimate, and read their spread as your uncertainty.

This is also how you find the bottom of the bias–variance U: sweep a hyperparameter, cross-validate at each value, and pick the one with the best mean validation score — never the best training score, which only ever rewards more complexity.

From scratch, then scikit-learn

Every metric is simple arithmetic on the four confusion counts; AUC is a sweep over thresholds. Here it is by hand, matching scikit-learn:

# --- metrics from the four counts ---
def metrics(y, pred):
    TP = ((pred==1)&(y==1)).sum(); FP = ((pred==1)&(y==0)).sum()
    FN = ((pred==0)&(y==1)).sum(); TN = ((pred==0)&(y==0)).sum()
    precision = TP / (TP + FP)
    recall    = TP / (TP + FN)
    f1        = 2*precision*recall / (precision + recall)
    return precision, recall, f1

pred = (scores >= 0.5).astype(int)
print(metrics(y, pred))                       # (0.77, 0.78, 0.78)

# --- AUC: sweep every threshold, integrate TPR vs FPR ---
def auc(y, scores):
    order = np.argsort(-scores); y = y[order]
    tpr = np.cumsum(y) / y.sum()
    fpr = np.cumsum(1-y) / (1-y).sum()
    return np.trapz(tpr, fpr)
print(round(auc(y, scores), 3))               # 0.820

from scratch
(0.77, 0.78, 0.78)
0.820

scikit-learn bundles all of it:

from sklearn.metrics import precision_score, recall_score, f1_score, roc_auc_score
print(precision_score(y, pred), recall_score(y, pred), f1_score(y, pred))   # 0.77 0.78 0.78
print(roc_auc_score(y, scores))               # 0.820 — identical
from sklearn.model_selection import cross_val_score
print(cross_val_score(model, X, y, cv=5).mean())   # honest k-fold estimate

Choosing well

Start from the cost of each error. Balanced classes and symmetric costs? Accuracy is fine. Care about catching positives? Watch recall. Care about not crying wolf? Watch precision. Need one number? F1, or AUC for threshold-independent ranking quality. Working with rare positives? Trust the PR curve and average precision over accuracy and ROC. And always estimate with cross-validation, tuning to the mean validation score. Set the operating threshold last — deliberately, from the curve, to match what a miss and a false alarm actually cost you.

Keep neighbours together

The goal is a 2-D map where points that were close in the original high-dimensional space stay close, and points that were far end up apart. We don't care about preserving exact distances or directions — only the neighbourhood structure, who-is-near-whom. That single, humble objective is enough to pull tangled high-dimensional clusters apart into a picture you can actually look at.

When a linear shadow fails

Here are three clusters living in six dimensions. Their real separation sits in a few low-variance directions, while a couple of high-variance "nuisance" directions dominate the spread. PCA, chasing variance, projects onto the nuisance — and the clusters smear into one blob. t-SNE ignores that global variance and reads the local neighbourhoods instead, recovering the three groups cleanly. Toggle the two projections of the same data:

Distances become probabilities

t-SNE never works with raw distances. For each point i it lays a Gaussian bell over its neighbours and converts distances into a probability that i would "pick" j as a neighbour — near points get high probability, far points almost none:

The width σ_i is set per point so that each ends up with roughly the same effective number of neighbours — a quantity called the perplexity. Dense regions get a narrow bell, sparse regions a wide one. Symmetrising these gives a single similarity p_ij for every pair.

The heavy-tailed trick

In the 2-D map, similarities q_ij use not a Gaussian but a Student-t (Cauchy) distribution — one with a heavy tail:

This is the masterstroke. Two dimensions have far less room than fifty, so moderately-distant points would pile up — the "crowding problem." The heavy tail lets dissimilar points sit comfortably far apart in the map without paying a big penalty, so clusters push away from each other and gaps open up.

The two kernels are plotted below as similarity versus distance. Near zero they agree — close points are similar either way. But drag the distance outward and the Gaussian collapses toward nothing while the Student-t keeps a stubborn heavy tail, so a dissimilar pair can sit far apart in the map without the similarity (and the cost) blowing up:

Pull neighbours, push strangers

Now place the points randomly in 2-D and slide them until the map's similarities Q match the high-D similarities P. The mismatch is measured by the KL divergence, minimised by gradient descent:

The gradient resolves into two forces on every point: an attraction toward points it was close to in high-D, and a repulsion from everything else. Neighbours clump; strangers drift apart; clusters emerge.

From noise to clusters

The same six-dimensional, three-cluster data, started as a random 2-D blob. Press Run and watch gradient descent on the KL divergence pull the neighbourhoods together — the three clusters separate out of the noise over a few hundred iterations, then settle. Each colour is a true cluster (the algorithm never sees the colours).

The neighbourhood dial

Perplexity is t-SNE's main knob — loosely, how many neighbours each point tries to keep. Too small and the map fragments into tiny specks, fixating on the very nearest points; too large and distinct groups blur together as the bells overlap everything. There's no universally right value (5–50 is the usual range); it's worth trying a few. Compare three on the same data:

What the picture does not mean

t-SNE maps are seductive and easy to over-read. Three rules keep you honest:

The one thing you can trust: points that land together really were neighbours in high-D. Local structure is the signal; global layout is decoration.

The faster cousin

UMAP pursues the same goal — preserve neighbourhoods — but builds a weighted k-nearest-neighbour graph and then optimises a low-dimensional layout of that graph (with attractive and repulsive forces, much like a force-directed graph). In practice it's substantially faster on large datasets, tends to preserve a bit more global structure (so inter-cluster distances mean a little more than in t-SNE), and exposes its own neighbourhood knob, n_neighbors, plus a min_dist controlling how tightly points may pack. The same reading cautions apply.

From scratch, then the libraries

The heart of t-SNE is the two distributions and a gradient loop — the same one driving the animation above:

# --- t-SNE core (sketch) ---
P = high_dim_affinities(X, perplexity=12)   # Gaussian p_ij, perplexity-calibrated σ
Y = 1e-4 * np.random.randn(n, 2)            # random 2-D start
for it in range(1000):
    num = 1 / (1 + sqdist(Y))               # Student-t kernel
    np.fill_diagonal(num, 0)
    Q = num / num.sum()                      # low-dim affinities q_ij
    grad = 4 * ((P - Q) * num)[:, :, None] * (Y[:, None] - Y[None])
    Y -= lr * grad.sum(1)                    # pull neighbours, push strangers

in practice
from sklearn.manifold import TSNE
Y = TSNE(n_components=2, perplexity=12, init="pca").fit_transform(X)

import umap            # the UMAP library
Y = umap.UMAP(n_neighbors=15, min_dist=0.1).fit_transform(X)

Always scale features first (both methods use distances), and run PCA down to ~50 dimensions before t-SNE on very wide data — it denoises and speeds things up.

Using embeddings well

Reach for t-SNE or UMAP when you want to see the structure of high-dimensional data — to check whether classes separate, spot outliers, or explore clusters. Keep them out of your modelling pipeline: the coordinates aren't stable features, distances aren't trustworthy, and there's no clean way to project new points (UMAP does this better than t-SNE). Use PCA when you need a fast, linear, reversible reduction that preserves global variance; use t-SNE/UMAP when you need a faithful local picture to look at.

Let the data say how many, and what shape

k-Means partitions space into straight-edged cells around k centres — wonderful for round, well-separated blobs, hopeless for anything that curves or interlocks, and it always returns exactly the k you asked for, noise and all. The methods here ask different questions. DBSCAN asks "where is the data dense, and which dense points connect?" Hierarchical clustering asks "if I repeatedly merge the closest groups, what tree do I get?" Neither needs k up front, and DBSCAN doesn't even force every point into a cluster.

Two moons

Here are two interleaving crescents with a scatter of noise. k-Means, forced to carve space around two centres, slices straight through both moons — its boundary is a line, and the moons aren't linearly separable. DBSCAN walks along each crescent's density and recovers them whole, setting the stray points aside as noise. Toggle the two:

Core, border, noise

DBSCAN needs two numbers: a radius ε (eps) and a minimum count minPts. Around each point, look within radius ε; that's its neighbourhood. Then every point is one of three kinds:

Only core points can grow a cluster; border points get pulled in at the edges; noise is left unlabelled. That third category — admitting some points belong to nothing — is what makes DBSCAN robust to outliers.

Density-reachable

The algorithm is a flood-fill through dense regions. Pick an unvisited point; if it's core, start a new cluster and absorb its entire ε-neighbourhood. For each new core point absorbed, absorb its neighbourhood too, rippling outward until no more core points are reachable. Points reachable through a chain of core neighbourhoods are density-reachable and land in the same cluster — which is exactly how a cluster can snake along a curved crescent. When the ripple stops, find the next unvisited point and start again. Non-core points that never get reached become noise.

Turn the two dials

The same moons-and-noise data. ε sets how far the flood-fill reaches; minPts sets how dense a region must be to count as core. Too small an ε and clusters shatter into fragments with everything flagged as noise; too large and separate clusters fuse into one. Find the sweet spot where the two moons emerge and the scatter stays grey (noise). Core points are drawn solid, noise hollow:

What density buys, and costs

See the squeeze directly: two tight clusters on the left, one sparse cluster on the right, all sharing one ε (minPts fixed at 4). At a small ε the dense pair resolves cleanly but the sparse cluster dissolves into noise; raise ε to rescue the sparse cluster and the two dense clusters fuse into one. No single value serves all three densities:

A common recipe for ε: plot each point's distance to its k-th nearest neighbour, sort, and look for the "elbow" — the knee of that curve is a good ε. Variants like HDBSCAN remove the single-ε limitation by working across density scales.

Build the whole tree

Agglomerative clustering takes the opposite tack to "pick k." Start with every point as its own cluster, then repeatedly merge the two closest clusters — closest by some rule — recording each merge and the distance at which it happened. After n−1 merges everything is one cluster, and the record of merges is a tree: the dendrogram. You don't choose k while building it; you build the entire hierarchy once, then cut it wherever you like.

What "closest" means

"Distance between two clusters" can be defined several ways, and the choice — the linkage — shapes the result:

Merge, then cut

Sixteen points in four little groups. The tree shows the merge order bottom-up — short joins for nearby points, tall joins for distant groups. Slide the cut line: every branch it crosses becomes one cluster, so a low cut gives many small clusters and a high cut gives few big ones. The big vertical gaps between merges are the natural places to cut. Switch the linkage and watch the tree rebuild:

From scratch, then scikit-learn

DBSCAN is a neighbourhood query plus a flood-fill:

# --- DBSCAN from scratch ---
def dbscan(X, eps, min_pts):
    labels = [0]*len(X)           # 0 = unvisited, -1 = noise, >0 = cluster
    def region(p): return [q for q in range(len(X)) if dist(X[p],X[q]) <= eps]
    C = 0
    for p in range(len(X)):
        if labels[p] != 0: continue
        nbrs = region(p)
        if len(nbrs) < min_pts: labels[p] = -1; continue   # not core -> noise (for now)
        C += 1; labels[p] = C; seeds = list(nbrs)
        i = 0
        while i < len(seeds):                    # flood-fill the dense region
            q = seeds[i]; i += 1
            if labels[q] == -1: labels[q] = C     # border point
            elif labels[q] == 0:
                labels[q] = C
                qn = region(q)
                if len(qn) >= min_pts: seeds += qn  # q is core -> keep expanding
    return labels

Agglomerative clustering is repeated nearest-pair merging:

# --- agglomerative (average linkage) ---
clusters = [[i] for i in range(len(X))]
while len(clusters) > 1:
    a, b = closest_pair(clusters)          # min average cross-distance
    clusters.append(clusters[a] + clusters[b])
    del clusters[max(a,b)]; del clusters[min(a,b)]   # record height = that distance

in practice
from sklearn.cluster import DBSCAN, AgglomerativeClustering
DBSCAN(eps=0.15, min_samples=4).fit_predict(X)        # -1 marks noise
AgglomerativeClustering(n_clusters=3, linkage="ward").fit_predict(X)
from scipy.cluster.hierarchy import linkage, dendrogram   # for the tree + plot

Choosing a clusterer

Reach for DBSCAN when clusters are irregularly shaped, the count is unknown, and outliers need flagging — spatial data, anomaly detection — but standardise features and expect to tune ε, and avoid it in high dimensions. Reach for hierarchical clustering when you want to see structure at every scale or there's a real nested taxonomy (gene families, document topics); the dendrogram is its own diagnostic, though it costs more than linear time. Keep k-Means for large, roughly-round, well-separated blobs where speed matters, and Gaussian mixtures when clusters overlap and you want soft, probabilistic membership.

Garbage in, garbage out

A learning algorithm only ever sees the numbers you feed it. If those numbers are on wild scales, miscoded, skewed, or riddled with holes, no amount of model tuning rescues them — and conversely, thoughtful features can make a simple linear model shine. Feature engineering is where domain knowledge enters the pipeline, and in practice it's where most of the time (and most of the gains) live. It's also where the subtlest bugs hide, because a preprocessing mistake doesn't crash — it quietly inflates your score.

Put features on equal footing

Any method that measures distance (k-NN, k-Means, SVM, DBSCAN, t-SNE) or follows gradients sees a feature measured in the thousands as overwhelmingly more important than one measured in fractions — purely because of its units. Below, the class is decided entirely by the small-scale feature on the x-axis, while the y-axis feature is large-scale noise. In raw units the noise dominates every distance, so k-NN picks neighbours in a useless horizontal band and misclassifies; standardised (each feature centred and divided by its spread), the signal feature gets an equal vote and the neighbours snap to the right class:

Standardisation (z-score) is the default; min–max scaling squeezes to [0, 1]; robust scaling uses the median and interquartile range to shrug off outliers. Tree-based models (decision trees, random forests, gradient boosting) are the exception — they split on thresholds, so monotonic rescaling doesn't change them at all.

Turning labels into numbers

Models eat numbers, so categorical variables must be encoded — and the choice matters:

The classic error is ordinal-encoding a nominal variable: labelling {red, green, blue} as {0, 1, 2} tells the model green is "between" red and blue and that blue is "three times" red — relationships that don't exist. One-hot avoids inventing that false order, at the cost of extra columns.

Let a linear model bend

A linear model can only draw straight boundaries — but it draws them in whatever feature space you give it. Add x² and a model linear in the parameters fits a parabola; add the product x₁·x₂ and it captures an interaction (the effect of one feature depending on another). This is exactly the kernel trick made explicit, and the lever behind the regularization and SVM dives. Other workhorse transforms: logarithms for multiplicative quantities (prices, counts), ratios and differences that encode domain meaning, and binning a continuous variable into ranges when the relationship is steppy rather than smooth.

Here is the simplest case made visible. One feature x, two classes: the inner class sits near zero, the outer class on both flanks — no single threshold on x can split them. But add the engineered feature x² and the outer points (large |x|) lift high while the inner ones stay low, so a single straight line now cuts cleanly between them. Drag the lift to fold in the x² feature:

Tame a long tail

Many real quantities — income, population, response times — are heavily right-skewed: most values are small, a few are enormous. That long tail distorts means, dominates distances, and violates the roughly-symmetric assumptions many models prefer. A log transform compresses the tail and pulls the distribution toward symmetry. Toggle it on the skewed feature below and watch the lump on the left spread into a bell:

The Box–Cox and Yeo–Johnson transforms generalise this, searching for the power that best symmetrises a feature automatically.

Filling the holes

Real datasets have gaps, and most models can't take a blank. The fixes, roughly in order of sophistication:

Crucially, imputation statistics (that mean, that median) must be computed on the training set alone — which brings us to the trap that ruins more analyses than any other.

The curse of dimensionality

More features is not more better. As dimensions grow, points spread apart until everything is roughly equidistant — distances lose meaning, and distance-based methods degrade. Worse, with enough columns some will correlate with the target by chance alone, so the model fits noise. Defences: drop near-constant and highly-correlated columns, select features by a real statistical or model-based criterion, or reduce dimensionality with PCA. But feature selection itself is dangerous if done naively — as the next chapter shows in the most vivid way possible.

The score that lies

Leakage is when information from outside the training set sneaks into it — and it produces validation scores that look spectacular and mean nothing. The data below is pure noise: random features, random labels, zero real signal, so honest accuracy must be 50%. The leaky pipeline picks the features most correlated with the label using the whole dataset, then cross-validates — and reports soaring accuracy, climbing as you add more noise features (more lottery tickets to find a spurious winner). The proper pipeline selects features inside each fold and stays pinned at chance. Slide the feature count:

The fixes are discipline, not cleverness: split first, then fit every preprocessing step — scaling, imputation, feature selection, encoding — on the training fold only, and apply the learned parameters to validation. Watch too for target leakage (a feature that secretly encodes the answer, like "date account closed" predicting churn) and temporal leakage (training on future data to predict the past).

Make leakage impossible

The reliable cure is to bind every preprocessing step and the model into a single pipeline object, then cross-validate the whole thing. Because the pipeline re-fits its preprocessing on each training fold automatically, the validation fold stays untouched — leakage becomes structurally impossible rather than something you must remember to avoid. A pipeline also guarantees that the exact transformations learned in training are reproduced, in order, at prediction time, so training and serving can't silently drift apart.

From scratch, then scikit-learn

The golden rule of scaling in two lines: learn the statistics on train, apply them to test — never the reverse.

# --- standardise: fit on TRAIN only ---
mu, sigma = X_train.mean(0), X_train.std(0)     # learned from training data alone
X_train_s = (X_train - mu) / sigma
X_test_s  = (X_test  - mu) / sigma              # SAME mu, sigma — no peeking at test

In practice, let a pipeline enforce it so you can't get it wrong:

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer

pre = ColumnTransformer([
    ("num", Pipeline([("impute", SimpleImputer(strategy="median")),
                       ("scale",  StandardScaler())]),        num_cols),
    ("cat", OneHotEncoder(handle_unknown="ignore"),           cat_cols),
])
model = Pipeline([("pre", pre), ("clf", LogisticRegression())])
cross_val_score(model, X, y, cv=5)      # preprocessing re-fit per fold -> no leakage

why it's safe
each fold: fit impute+scale+encode on train rows only, then transform val rows
-> the validation data never influences the statistics used to transform it

A working order

A sane sequence: split off a test set first; explore and clean only the training data; impute, encode, transform skew, and scale — all fit on train; engineer domain features and interactions; select or reduce if dimensionality bites; wrap it all in a pipeline and cross-validate. Resist the urge to peek at the test set until the very end, and be suspicious of any feature or score that looks too good — it usually is. The unglamorous truth of applied ML is that careful preprocessing, not exotic models, is what most reliably moves the needle.

The knobs you set, not learn

If you picked hyperparameters by training error, you'd always choose maximum flexibility — depth-infinity trees, k=1 neighbours, zero regularization — and overfit completely. The training set can't judge its own complexity. So tuning needs a separate yardstick: data the model didn't train on. The whole craft is (1) get an honest estimate of how each setting generalises, (2) search the space of settings efficiently, and (3) never contaminate your final estimate with the choices you made along the way.

Train, validate, test

The classic discipline splits data three ways, each with one job:

The test set is sacred: every time you peek at it to make a decision, it stops being an honest estimate of the future and becomes just another thing you've overfit. Lock it away until the model is final.

A better validation signal

A single validation split is noisy, and it spends data you'd rather train on. k-fold cross-validation (from the evaluation dive) fixes both: for each hyperparameter setting, train on k−1 folds and validate on the held-out fold, rotating through all k, and average. That mean is a lower-variance estimate of how the setting generalises, and every row gets used for both training and validation across the folds. Tuning then means: compute the cross-validated score for each candidate setting, and keep the best.

Try every combination

The simplest search lays a grid over the hyperparameter ranges and evaluates every cell by cross-validation. It's exhaustive and trivially parallel — but it suffers the curse of dimensionality. Five hyperparameters at six values each is 6⁵ = 7,776 fits, times k folds. And it spends that budget rigidly: every value of every hyperparameter gets tested the same number of times, whether it matters or not. Which leads to a surprising failure mode.

Why random beats grid

Here's the response surface for a model with two hyperparameters — but only the horizontal one actually matters (the bright ridge runs vertically; the vertical knob barely moves the score). A 3×3 grid spends 9 evaluations but tests only 3 distinct values of the important knob. Random search spends the same 9 evaluations on 9 distinct values of it — far more likely to land on the ridge. Toggle, and resample the random draw:

The lesson (Bergstra & Bengio, 2012): when only a few hyperparameters matter — which is usual — random search finds good settings faster, because it never wastes its budget re-testing the same value of an important knob. It's the sensible default before reaching for anything cleverer.

Sweep one knob

To tune a single hyperparameter, plot the cross-validated score against it. Here is k for k-NN on a noisy two-class problem: small k is flexible (training score near perfect, but it's memorising noise), large k is over-smoothed (both scores sag). The training score and the validation score tell different stories — and tuning means picking the peak of the validation curve, never the training one. Slide k and watch:

This is the bias–variance tradeoff seen through the tuning lens: the gap between the two curves is overfitting, and the validation peak is the sweet spot grid or random search is hunting for.

Beyond brute force

Here's Bayesian optimization on a single hyperparameter. The true score curve is hidden; the optimizer fits a surrogate (the mean line, with a shaded uncertainty band that's wide where it hasn't looked and tight near samples), then picks its next evaluation where mean-plus-uncertainty peaks — balancing exploiting good regions against exploring unknown ones. Sample a few times and watch it home in on the peak:

(the dashed line marks where it will sample next)

Don't grade your own tuning

Here's the subtle trap. If you tune over a grid with cross-validation and then report that best cross-validated score as your performance, you've cheated — you picked the maximum of many noisy estimates, so it's optimistically biased upward. The data below has no real signal (true accuracy 50%), yet selecting the best of twenty hyperparameter settings reports well above chance. Nested cross-validation fixes it: an inner loop tunes, an outer loop evaluates the entire tuning procedure on data it never saw — restoring an honest estimate.

Putting it in order

A dependable sequence: hold out a test set and lock it away; build a pipeline so every preprocessing step re-fits inside each fold; pick a search strategy (random search as default, grid for one or two knobs, Bayesian when fits are costly); search hyperparameters with cross-validation on the training data; refit the best configuration on all the training data; and only then, once, score it on the test set. Spend your budget on the hyperparameters that actually move the metric — a few usually dominate — and use early stopping to cut off hopeless runs. Report the test score, not the tuning score.

From the loop to scikit-learn

Tuning is conceptually a loop over settings, scoring each by cross-validation:

# --- the idea, by hand ---
best_score, best_k = -1, None
for k in range(1, 26):                       # candidate hyperparameters
    score = cross_val_score(KNN(k), X_train, y_train, cv=5).mean()
    if score > best_score:
        best_score, best_k = score, k        # keep the best on VALIDATION folds

scikit-learn wraps this — and crucially lets you tune a whole pipeline, so preprocessing re-fits per fold and can't leak:

from sklearn.model_selection import GridSearchCV, RandomizedSearchCV, cross_val_score
search = RandomizedSearchCV(pipeline, param_distributions, n_iter=40, cv=5)
search.fit(X_train, y_train)
print(search.best_params_)
final_score = search.score(X_test, y_test)   # test set, opened once

# nested CV for an unbiased estimate of the whole procedure:
honest = cross_val_score(search, X, y, cv=5).mean()   # outer loop wraps the search

key idea
the test set (or the outer fold) is touched only to score a fully-decided model
- never to choose anything. choosing happens on training/validation data alone.

Tuning well

Start simple: sensible defaults and a single validation split will get you most of the way, and many hyperparameters barely matter. Reach for cross-validated random search when you want to do it properly, grid search when there are only one or two knobs, and Bayesian or halving methods when each fit is expensive. Above all, keep the test set untouched and be honest about what your reported number measures — a score you tuned against is not a score you can trust. The model is rarely the bottleneck; disciplined evaluation is.

The wisdom of crowds

Ask one expert to guess the weight of an ox and they'll be off; average a few hundred fair guesses and the crowd lands almost exactly right. The individual errors point in different directions and cancel. Machine-learning ensembles bottle this effect: train many models, combine their outputs, and the combination is steadier and more accurate than its members — provided the members are decent (better than chance) and diverse (wrong in different ways). Every method here is a different answer to "how do we build diverse-but-decent members, and how do we combine them?"

Averaging shrinks variance

Here's the engine of the whole field. Each thin curve is a depth-3 tree fit to a bootstrap of the same noisy data — high-variance on its own, lurching wildly from sample to sample. But average them (the bold curve) and the wiggles cancel toward the true signal. Slide the number of trees and watch the average settle:

The math is exact: average N models whose errors each have variance σ² and correlation ρ, and the ensemble's variance is ρσ² + (1−ρ)σ²/N. The second term melts away as N grows — that's the averaging win. But the first term, ρσ², is a floor set by how correlated the models are. Independent models (ρ→0) average down to nothing; identical models (ρ=1) average to no improvement at all. Everything hinges on ρ.

See that law directly. The curve is ensemble variance (as a fraction of one model's) against the number of members; the dashed line is the ρσ² floor. Drag ρ: at zero, independent errors drive the variance down as 1/N toward nothing; raise it and the curve slams into the floor, where adding members stops helping no matter how many you pile on:

Parallel: bootstrap then aggregate

Bagging (bootstrap aggregating) is the direct application: build each member on a different bootstrap resample of the training data, then average (regression) or vote (classification). The resampling is what makes the members differ, so their variance averages down — while the bias stays put. It shines on high-variance, low-bias learners like deep trees, which is exactly why a random forest is bagged deep trees plus an extra trick (random feature subsets at each split) to push ρ even lower. Members train independently, so bagging is trivially parallel.

Sequential: fix the last mistake

Boosting takes the opposite tack. Instead of averaging independent models, it builds them in sequence, each one focused on the examples its predecessors got wrong — reweighting errors (AdaBoost) or fitting the residual gradient (gradient boosting). Members are deliberately weak (shallow trees, "stumps") and the ensemble grows by small steps, so boosting chips away at bias, turning a pile of underfit learners into a sharp predictor. The price: because each step depends on the last, boosting is sequential, and its relentless error-chasing can overfit if you let it run too long.

Two cures for two diseases

Neither dominates: bagging is the safer default and parallelises; boosting often wins accuracy contests but needs more careful tuning. Both have their own dives — this one is about what they share, and what comes next.

Let them vote

The simplest combiner doesn't even need resampling: take several decent classifiers and let them vote. Below, each voter is a weak depth-2 tree (each gets about 80% on its own); as you add voters, the majority-vote boundary sharpens around the true wavy boundary and the ensemble climbs past every individual member. Hard voting counts labels; soft voting averages predicted probabilities and usually does a touch better. Slide in the voters:

Voting is the cheapest way to combine genuinely different model types — say a tree, a k-NN, and a logistic regression — whose errors are decorrelated precisely because their inductive biases differ.

Learn how to combine

Why combine with a fixed rule like a vote when you could learn the best combination? Stacking trains a second-level meta-model whose inputs are the predictions of the base models. The meta-model discovers that, say, the k-NN is trustworthy here and the tree there, and weights them accordingly. The one subtlety — and it's the whole game — is that the meta-model must train on out-of-fold predictions: each base model predicts on data it didn't train on (via cross-validation), or the meta-model just learns to trust whichever base model overfit hardest. That's leakage, and the out-of-fold scheme is what prevents it.

The crowd must disagree

The Condorcet jury theorem makes the diversity condition precise. Suppose N independent voters each get the answer right with probability p, and you take the majority. If p > ½, the ensemble's accuracy races toward 100% as you add voters; if p < ½, it races toward zero — a committee of below-chance voters is worse than one. Drag p across the ½ line and watch the crowd's fate flip:

The catch is the word independent. Real models share data and assumptions, so their errors correlate — and correlation is exactly the ρ floor from chapter 2. This is why every good ensemble method is, at heart, a diversity-manufacturing machine: bootstrap samples, random feature subsets, different algorithms, different error-chasing. Maximise disagreement among decent members, and the crowd grows wise.

The price of a committee

From scratch, then scikit-learn

Bagging and voting are a few lines — resample, fit, combine:

# --- bagging from scratch ---
models = []
for _ in range(n_estimators):
    idx = rng.choice(len(X), len(X), replace=True)   # bootstrap sample
    models.append(clone(base).fit(X[idx], y[idx]))

def predict(X):
    P = np.array([m.predict(X) for m in models])      # every model votes
    return (P.mean(0) > 0.5).astype(int)              # majority vote

scikit-learn provides the combiners directly, including stacking with built-in out-of-fold prediction:

from sklearn.ensemble import VotingClassifier, StackingClassifier
from sklearn.linear_model import LogisticRegression

vote = VotingClassifier([("tree", tree), ("knn", knn), ("svc", svc)],
                        voting="soft")               # average probabilities

stack = StackingClassifier(
    estimators=[("tree", tree), ("knn", knn), ("svc", svc)],
    final_estimator=LogisticRegression(),            # meta-model
    cv=5)                                            # out-of-fold preds -> no leakage

the rule of thumb
diverse base models + a simple meta-model (or a plain vote) > any single model
- and the diversity matters far more than the cleverness of the combiner

Combining well

Start with a strong single model; reach for an ensemble when you need the last few points of accuracy and can pay for them. Bag high-variance learners, boost high-bias ones, and when you have several good-but-different models already, vote or stack them rather than picking one. Above all, engineer diversity — different data, different features, different algorithms — because a committee of clones is just an expensive single model. The leaderboards of tabular machine learning are owned by ensembles for a reason.

One pipeline, end to end

The whole process is a pipeline with a fixed order, and the discipline is mostly about not skipping ahead — especially not touching the test set until the very end. Here's the path this dive walks, and the dive each step draws on:

Understand before you model

Before any algorithm, look at the data. How many rows (569) and features (30)? What's the class balance — here 357 benign to 212 malignant, roughly 63/37, imbalanced enough that accuracy alone will mislead (a model that screams "benign" every time already scores 63%). Are features on wildly different scales? They are — radii in the tens, areas in the hundreds, smoothness in the hundredths — so any distance- or gradient-based model will need standardising. Are there missing values or obvious leaks? Catching all this now, on the training data, is the work the feature engineering dive is about.

Lock the test set away

The very first action — before exploring deeply, before fitting a scaler — is to split off a test set (here 20%, stratified to preserve the 63/37 balance) and not look at it again until the end. Every decision from here lives on the training data, judged by cross-validation. This is the single rule that separates an honest performance estimate from a fantasy, and it's the spine of the evaluation and hyperparameter-tuning dives. Touch the test set early and every number after it is contaminated.

Inside a pipeline, fit on train

These features are all numeric, so preprocessing here is mainly standardisation — but the rule is universal: bind every step (scale, and for other datasets impute and one-hot-encode) into a pipeline so it re-fits on each training fold and never sees validation data. Had there been categories, they'd be one-hot encoded; had there been gaps, imputed with a training-set statistic. The pipeline is what makes leakage structurally impossible, the hard-won lesson of the feature engineering dive.

Earn the right to be complex

Always establish a dumb baseline first, so you know what "good" even means. The majority-class baseline here — always guess benign — scores 63.2%. That's the number every real model must beat to justify itself. A baseline also catches catastrophic bugs early: if your fancy model can't beat "always guess the common class," something is broken in the data or the pipeline, not the algorithm.

Let the candidates audition

Now run several models through identical 5-fold cross-validation on the training set and compare honestly. Each is a full pipeline (scaling + estimator). The dashed line is the baseline; every real model towers over it, and a handful cluster near the top — on this clean dataset the differences are small, but the SVM edges ahead:

When scores are this close, other factors break the tie: training cost, interpretability, calibration. Here we carry the SVM forward, but logistic regression — nearly as accurate and far more interpretable — would be a defensible choice too.

Search the knobs

With a model chosen, tune its hyperparameters by cross-validated search — for the RBF SVM, that's the regularization C and the kernel width γ, swept on a log scale. Here the search confirmed the defaults were already near-optimal (C=1, γ=scale), nudging cross-validated accuracy to 97.8%. That's a real outcome worth internalising: tuning doesn't always produce fireworks, and a search that confirms sensible defaults has still done its job. This is the hyperparameter-tuning dive applied for real — and note we're still only touching training data.

Open the test set — once

Only now, with the model fully decided, do we unlock the test set (114 tumors the model has never seen) and score it once. But "score" means more than accuracy on an imbalanced medical problem — a missed malignancy is far costlier than a false alarm. Below is the final model on the test set; slide the decision threshold and watch the confusion matrix and the operating point trade missed cancers against false alarms:

At the default threshold the model catches 40 of 42 malignancies and lands around 96% accuracy, at an ROC-AUC of 0.993. Lower the threshold and you catch the last malignancies — at the price of more false alarms. Which threshold to ship is not a machine-learning question but a clinical one; the model's job is to make the tradeoff explicit, exactly as the evaluation dive argued.

What did it learn, and can we trust it?

A number isn't a deliverable; understanding is. Inspect which features drive the predictions (for this dataset, the "worst" nucleus measurements — largest radius, most concave points — dominate, which is biologically sensible and a good sign the model isn't exploiting an artefact). The interpretable logistic model from the leaderboard — nearly as accurate as the SVM — gives a clean window in: its standardised coefficients show which measurements push a prediction toward malignant, and they line up with intuition (bigger, rougher, more-concave nuclei):

Check where it errs: the two missed malignancies are the cases to study. State the caveats plainly — the estimate holds only for data like the training distribution, and a 0.993 AUC on a curated benchmark will not survive contact with messy real-world imaging unchanged. Honesty about limits is part of the result.

The whole thing, end to end

The entire workflow is a few dozen lines — and notice the test set appears exactly once, on the last line:

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2,
                                          stratify=y, random_state=0)   # split FIRST

pipe = make_pipeline(StandardScaler(), SVC(probability=True))           # preprocess + model, bound together

grid = {"svc__C": [0.1, 1, 10, 100], "svc__gamma": ["scale", 1e-2, 1e-3]}
search = GridSearchCV(pipe, grid, cv=5, scoring="accuracy")             # tune on TRAIN via CV
search.fit(X_tr, y_tr)

print(cross_val_score(search.best_estimator_, X_tr, y_tr, cv=5).mean()) # honest CV estimate
final = search.best_estimator_.fit(X_tr, y_tr)                          # refit on all training data
print(classification_report(y_te, final.predict(X_te)))                 # TEST SET — opened once, at the end

final test performance
accuracy 0.965   ROC-AUC 0.993
malignant: precision 0.95  recall 0.95   (caught 40 of 42; 2 missed)
benign:    precision 0.97  recall 0.97

The map of the whole field

That's the arc, and every step rests on a dive: look and preprocess (feature engineering), split and evaluate (model evaluation), compare across the model zoo (every algorithm dive), tune (hyperparameter tuning), and reason about why it generalises (bias–variance, ensembles). The model was almost the easy part; the discipline around it — splitting first, fitting on train, baselining, cross-validating, choosing the right metric, opening the test set once — is what turns a number into a result you can stand behind.