Regression, Trees, and the Metrics That Judge Them

21.1 What a model is, precisely

A model is a function with adjustable numbers in it. You choose the shape of the function, the data chooses the numbers, and a third thing — the loss function — decides what "chose well" means.

That third thing is the one people skip, and it is the one that determines almost everything about how your model behaves. Two models with the identical shape, fitted to the identical data, will disagree completely if they are minimising different losses. This chapter is organised around that fact.

21.2 The loss function is the whole argument

Suppose you are predicting tomorrow's demand for a product, and your model is off by 10 units. How bad is that? There is no answer in the data. It depends entirely on what being wrong costs you, and that is a business fact, not a statistical one.

If being off by 10 is twice as bad as being off by 5, you want absolute error. If being off by 10 is four times as bad as being off by 5 — because a large miss empties the shelf and sends the customer to a competitor — you want squared error. The loss function is where that judgement gets written down.

21.3 Squared, absolute, and Huber loss

Look at the shapes. The horizontal axis is the residual — how wrong you were, positive or negative. The vertical axis is how much you are punished for it.

Squared error is a parabola, so the punishment grows faster than the error does. At a residual of 1 it costs 1; at a residual of 3 it costs 9. This is why a single mis-keyed row can dominate an entire fit: that one row contributes more to the total loss than a hundred small errors put together, and the optimiser, which only sees the total, will contort the whole model to reduce it.

Absolute error is a V. Punishment grows in step with the error, so an outlier is just another data point with a large residual. It gets no special power.

import numpy as np

# One corrupted row. Watch what each loss does about it.
y     = np.array([10.0, 11.0, 9.0, 10.5, 10.0])
clean = np.array([10.0, 11.0, 9.0, 10.5, 10.0])
dirty = clean.copy()
dirty[4] = 60.0                      # a data-entry error: 10.0 -> 60.0

def mse_fit(v):   return v.mean()               # minimises squared error
def mae_fit(v):   return float(np.median(v))    # minimises absolute error

print(f"clean data:  mean {mse_fit(clean):5.2f}   median {mae_fit(clean):5.2f}")
print(f"one outlier: mean {mse_fit(dirty):5.2f}   median {mae_fit(dirty):5.2f}")
print(f"mean moved   {mse_fit(dirty) - mse_fit(clean):5.2f}")
print(f"median moved {mae_fit(dirty) - mae_fit(clean):5.2f}")

One corrupted value moves the squared-error fit by a full 10 units and the absolute-error fit by 0.5 — twenty times less. Note that the median does not stay perfectly still, which is the sort of detail that gets rounded off in textbooks. Replacing a value reshuffles the ordering, so the middle element shifts a little. The honest claim is "twenty times less sensitive", not "immune".

Huber loss is the obvious compromise, and it is worth knowing because it comes up in interviews and in robust forecasting. It is squared near zero and linear further out, switching at a point you choose called delta. Small residuals get the smooth gradient that optimisers like; large ones get the linear treatment that refuses to let outliers take over. In the figure it tracks the parabola near the origin and then straightens out to run parallel to the V.

21.4 Least squares, solved by hand

Linear regression assumes the relationship is a straight line (or a flat plane, in more dimensions) and then finds the line that minimises squared error. Because the loss is a parabola in the coefficients, and a parabola has exactly one bottom, there is a unique answer and it can be written down algebraically. No iteration, no learning rate, no random seed.

import numpy as np

# Ordinary least squares has a closed form. No iteration, no learning rate:
# beta = (X'X)^-1 X'y. Solving it directly shows what the model IS.
rng = np.random.default_rng(0)
n = 200
x = rng.normal(0, 1, n)
y = 3.0 + 2.5 * x + rng.normal(0, 0.8, n)   # true intercept 3, slope 2.5

X = np.column_stack([np.ones(n), x])        # column of 1s = the intercept
beta = np.linalg.solve(X.T @ X, X.T @ y)    # solve, never invert explicitly

print(f"intercept {beta[0]:.3f}   slope {beta[1]:.3f}")

# The residuals are what is left over. OLS makes their SQUARED sum smallest.
resid = y - X @ beta
print(f"sum of residuals   {resid.sum():.6f}")   # ~0 by construction
print(f"sum of squares     {(resid ** 2).sum():.2f}")

Two details in that output are worth pausing on. The residuals sum to zero — not approximately, but to machine precision. That is forced by the intercept term: the fitting procedure makes the residuals orthogonal to every column of the design matrix, and the intercept's column is all ones, so their sum must vanish. And notice we used np.linalg.solve rather than computing an inverse. Explicitly inverting a matrix is slower and numerically worse; solving the system directly is the right habit and interviewers notice.

21.5 The five assumptions, and which ones matter

Textbooks list assumptions behind linear regression and imply that violating any of them invalidates everything. That is not true, and knowing which ones actually bite is the difference between using the tool and reciting about it. The critical question is always: am I predicting, or am I explaining?

Normality is the most over-taught item on that list. It is not required for the coefficient estimates to be correct on average; it is required only for exact small-sample inference. With a few thousand rows the central limit theorem does the work and you can stop worrying about it. Interviewers who ask "what if the residuals aren't normal?" are usually checking whether you know that.

21.6 Reading a coefficient without lying

A fitted coefficient of 2.5 on price means: comparing two rows that differ by one unit of price and are identical in every other feature in the model, the prediction differs by 2.5. Every clause in that sentence is load-bearing.

• "in every other feature in the model" — not every other feature in the world. Omit a variable that drives both price and the outcome and the coefficient absorbs its influence. This is confounding, and no amount of data fixes it.
• "comparing two rows" — not "changing the price of one row". The model saw observational data; it knows about association. Reading a coefficient as the result of an intervention requires an experiment or a causal design, which is chapter 19's subject.
• "the prediction differs" — the prediction, not reality.

21.7 Multicollinearity and the VIF

When two features carry nearly the same information, the fitting procedure cannot tell how to divide the credit between them, because many different divisions produce almost the same predictions.

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(3)
n = 400
x1 = rng.normal(0, 1, n)
x2 = x1 + rng.normal(0, 0.01, n)      # x2 is essentially a copy of x1
y = 2.0 * x1 + rng.normal(0, 0.5, n)   # truth depends on x1 only

def vif(X, j):
    """Variance inflation factor: how much collinearity inflates a variance."""
    others = np.delete(X, j, axis=1)
    r2 = LinearRegression().fit(others, X[:, j]).score(others, X[:, j])
    return 1.0 / max(1e-12, 1.0 - r2)

X = np.column_stack([x1, x2])
fit = LinearRegression().fit(X, y)
print(f"coefficients: x1 {fit.coef_[0]:+.2f}   x2 {fit.coef_[1]:+.2f}")
print(f"VIF x1 {vif(X, 0):,.0f}   VIF x2 {vif(X, 1):,.0f}")
print(f"sum of the two coefficients: {fit.coef_.sum():+.2f}  (the truth is +2)")

# Note what survived and what did not: the SUM is right, the SPLIT is noise.
# Prediction is fine. Any sentence about "the effect of x1" is not.

Read that output carefully, because it contains the resolution of a common confusion. The individual coefficients are nonsense — wild in magnitude and unstable against tiny changes in the data. But their sum recovers the truth, and the model's predictions are perfectly good.

The variance inflation factor quantifies it: regress each feature on all the others and see how well it can be predicted. A VIF of 1 means no collinearity; above 5 is worth a look; above 10 is conventionally a problem. The numbers above are in the hundreds because the two columns are near-identical by construction.

21.8 R-squared, adjusted R-squared, and their limits

R-squared answers one narrow question: what fraction of the variance in the target did the model account for, compared with just predicting the mean every time? An R-squared of 0.7 means the model's errors are 30% as large, in squared terms, as the errors of that trivial baseline.

21.9 MAE, RMSE, MAPE, and which to report

These are the same three ideas as the loss functions in 21.3, now used for reporting rather than fitting. You do not have to use the same one for both, and often you should not: fit with squared error for the smooth optimisation, report MAE because stakeholders understand it.

21.10 Why you cannot fit a line to a yes/no

The target is now 0 or 1 — churned or not, fraud or not. The obvious move is to fit a straight line to it and call anything above 0.5 a yes. Three things go wrong, and they are worth stating because the fixes define logistic regression.

1. The output leaves the interval. A line is unbounded, so for extreme inputs it predicts 1.4 or −0.3. There is no such probability.
2. The effect cannot really be constant. A line says one more unit of tenure changes the churn probability by the same amount whether the customer was at 0.02 or 0.95. Near the ceiling there is almost nothing left to change.
3. Squared error is the wrong punishment. Predicting 0.9 when the truth is 0 should be far worse than predicting 0.6, and squaring does not express that strongly enough.

21.11 Odds, log-odds, and the sigmoid

The fix is to stop modelling the probability directly and model something unbounded instead. Take it in two steps.

Odds convert a probability into a ratio: p / (1 − p). A probability of 0.5 is odds of 1 ("evens"). A probability of 0.9 is odds of 9. As p approaches 1 the odds run off to infinity — so the ceiling is gone, but the floor at zero remains.

Log-odds take the logarithm of that, which pushes the floor down to minus infinity. Now the quantity is unbounded in both directions, which is exactly what a linear model can safely produce. Log-odds are also called the logit.

import numpy as np

# Logistic regression is linear in the LOG-ODDS, not in the probability.
# That single sentence explains every coefficient interpretation.
def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

intercept, beta_age = -3.0, 0.05

for age in (20, 40, 60, 80):
    z = intercept + beta_age * age          # the log-odds ("logit")
    p = sigmoid(z)
    odds = p / (1 - p)
    print(f"age {age}:  logit {z:+.2f}   odds {odds:.3f}   p {p:.3f}")

# A one-unit rise in age multiplies the ODDS by exp(beta), always.
# It does NOT add a fixed amount to the probability.
print(f"odds ratio per year of age: {np.exp(beta_age):.4f}")

The assertions in that snippet make the crucial point: a one-unit change in a feature has a constant multiplicative effect on the odds, and a non-constant additive effect on the probability. Near p = 0.5 the curve is steep and a small push moves the probability a lot; out in the tails it is nearly flat and the same push barely registers. That is precisely the behaviour objection 2 in 21.10 demanded.

21.12 Log loss, and what it punishes

Logistic regression is fitted by minimising log loss, also called binary cross-entropy. For a single row it is −log(p), where p is the probability the model assigned to the class that actually occurred.

Look again at the right-hand panel of Figure 24.1. Assign 0.99 to the true class and you pay almost nothing. Assign 0.5 — a shrug — and you pay 0.69. Assign 0.02 to something that then happens and you pay 3.9. The curve has no upper bound: as the assigned probability approaches zero the penalty approaches infinity.

21.13 Interpreting an odds ratio out loud

Exponentiating a coefficient gives an odds ratio. A coefficient of 0.05 gives e^0.05 = 1.051, so each additional unit multiplies the odds by 1.051 — about a 5% increase in the odds.

21.14 Separation, and other ways it breaks

Logistic regression has one dramatic failure mode that surprises people the first time. If some feature (or combination) splits the classes perfectly, the fitting procedure diverges: it can always reduce the loss further by making the coefficient larger, so the coefficient runs off toward infinity.

This is complete separation. You will see it as an absurd coefficient, a standard error even more absurd, and possibly a convergence warning. Counter-intuitively it happens most often on small datasets and after adding a feature that is a little too good — frequently one that leaked the answer.

21.15 Ridge: shrink everything

Here the objective stops being the loss. Ridge regression minimises the squared error plus a penalty proportional to the sum of the squared coefficients. There is a dial, alpha, controlling how much the penalty counts.

The effect is that a coefficient now has to earn its size. A feature that reduces error only slightly is not worth the penalty it incurs, so its coefficient is pulled toward zero. This trades a little bias for a large reduction in variance — chapter 18.6's trade-off, made adjustable.

21.16 Lasso: shrink some to nothing

Lasso is identical except the penalty uses absolute values instead of squares. That one change has a consequence out of all proportion to its size: lasso sets coefficients exactly to zero, while ridge only ever makes them small.

import numpy as np
from sklearn.linear_model import Ridge, Lasso
from sklearn.datasets import make_regression
from sklearn.preprocessing import StandardScaler

X, y = make_regression(n_samples=120, n_features=80, n_informative=6,
                       noise=40.0, random_state=1)
X = StandardScaler().fit_transform(X)   # penalties are scale-sensitive!

print(f"{'alpha':>8} {'ridge nonzero':>14} {'lasso nonzero':>14}")
for alpha in (0.1, 1.0, 10.0, 100.0):
    r = Ridge(alpha=alpha).fit(X, y)
    l = Lasso(alpha=alpha, max_iter=20000).fit(X, y)
    print(f"{alpha:>8} {int((r.coef_ != 0).sum()):>14} "
          f"{int((l.coef_ != 0).sum()):>14}")

# Ridge shrinks coefficients toward zero but never TO zero.
# Lasso sets them exactly to zero -- that is why it selects features.

The intuition for why is worth carrying around. The gradient of the squared penalty shrinks as the coefficient shrinks, so the pull toward zero fades away as you approach it — you converge on something tiny but nonzero. The gradient of the absolute-value penalty is constant right up to zero, so the pull never weakens and the coefficient is pushed all the way in.

21.17 Why ridge cannot reach zero, four ways

The previous section showed lasso zeroing coefficients and ridge not. That is a demonstration, and a demonstration is the weakest form of the argument: it shows that it happens without showing why, and it leaves the obvious objection standing — maybe ridge would zero too, with a large enough penalty. It would not, and the reason is worth four separate looks, because each one makes a different part of it obvious.

Argument 1 — the force that remains at the origin

Picture the coefficient sliding toward zero and ask what is still pushing it. The penalty’s derivative is that push.

# Ask one question: as a coefficient approaches zero, how hard is the penalty
# still pushing? Differentiate each penalty and watch the limit.

print(f"{'b':>12} {'d/db (b^2)':>14} {'d/db |b|':>10}")
for b in (1.0, 0.1, 0.01, 1e-4, 1e-8):
    print(f"{b:>12} {2 * b:>14.10f} {1.0:>10}")

# L2: the push is 2b, which VANISHES as b vanishes. The penalty stops caring
#     precisely where you need it to care most, so motion stops early.
# L1: the push is 1. Constant. It never relents, however close to zero you
#     get, so the coefficient is driven all the way in.
#
# Squaring makes small numbers negligible (0.01^2 = 0.0001), and a penalty
# that treats small coefficients as negligible has no reason to remove them.

The squared penalty’s push is 2b, which fades to nothing exactly where you need it most. Halve the coefficient and you halve the pressure to shrink it further, so the shrinking slows as it approaches zero and settles somewhere tiny but positive. The absolute-value penalty pushes with constant force 1, no matter how close to zero the coefficient gets, so nothing stops it arriving.

Argument 2 — the closed forms, checked rather than quoted

With one standardised feature both problems have exact solutions, and the whole phenomenon is visible in that single dimension. Rather than quote the formulas, the code below verifies them against brute-force search — four million candidate values, no optimiser, no calculus, no trust required.

import numpy as np

# One standardised feature, so X'X = 1 and both problems collapse to
#     minimise  0.5 * (b - z)^2  +  penalty(b)
# where z is the least-squares answer. The entire ridge/lasso difference is
# already visible in this single dimension.

def brute_force(z, lam, kind, lo=-4.0, hi=4.0, n=4_000_001):
    """Minimise by exhaustive search: no formula, no optimiser, no faith."""
    b = np.linspace(lo, hi, n)
    penalty = lam * np.abs(b) if kind == "lasso" else 0.5 * lam * b**2
    return b[np.argmin(0.5 * (b - z) ** 2 + penalty)]

lam = 1.0
print(f"{'z':>7} {'lasso':>9} {'soft-thr':>9} {'ridge':>9} {'z/(1+lam)':>10}")
for z in (0.2, 0.5, 0.999, 1.0, 1.5, 3.0):
    lasso_num = brute_force(z, lam, "lasso")
    lasso_thy = np.sign(z) * max(abs(z) - lam, 0.0)     # soft-thresholding
    ridge_num = brute_force(z, lam, "ridge")
    ridge_thy = z / (1.0 + lam)                         # proportional shrink
    print(f"{z:>7} {lasso_num:>9.5f} {lasso_thy:>9.5f} "
          f"{ridge_num:>9.5f} {ridge_thy:>10.5f}")

# The formulas reproduce the brute-force search, so they can be trusted.
# Now read what they say:
#     lasso   b = sign(z) * max(|z| - lam, 0)   SUBTRACTS -> can reach 0
#     ridge   b = z / (1 + lam)                 DIVIDES   -> cannot

The formulas reproduce the search exactly, so they can be believed. And once written down they settle the question by inspection:

Subtraction reaches zero and stops there. Division never reaches zero: halving a number repeatedly gets you as close as you like without ever arriving. That distinction — subtract versus divide — is the difference between a method that selects variables and a method that does not.

Argument 3 — an equation versus an interval

Optimality means no downhill direction is available. For a smooth penalty that condition is an equation; for a kinked one it is an interval, and everything follows from that.

import numpy as np

# Optimality means "no downhill direction available". For a smooth penalty
# that is an EQUATION; for a kinked one it is an INTERVAL. That is the
# difference, stated exactly.
#
# RIDGE. The objective is differentiable everywhere, so at the optimum
#     (b - z) + lam * b = 0     ->     b = z / (1 + lam)
# Setting b = 0 forces z = 0: ridge returns zero only when least squares
# already did. No value of lam changes this.

for lam in (1.0, 1e2, 1e6, 1e12):
    b = 0.4 / (1.0 + lam)
    print(f"lam={lam:>8.0e}   ridge b = {b:>10.3e}   exactly zero? {b == 0.0}")

# At lam = 1e12 the coefficient is 4e-13. Tiny, and still not zero -- the
# variable remains in the model, still collected, still served, still a
# dependency you have to maintain.

print()
# LASSO. |b| has no derivative at 0; it has a SET of valid slopes, the whole
# interval [-lam, lam]. Zero is optimal whenever the data's pull |z| fits
# inside that interval, so an entire RANGE of z collapses to zero rather
# than a single unlucky point.
lam = 1.0
for z in (0.4, 1.0, 1.6):
    b = np.sign(z) * max(abs(z) - lam, 0.0)
    print(f"z={z:>4}   |z| <= lam? {str(abs(z) <= lam):<5}   ->   b = {b:.2f}")

Ridge’s objective is differentiable everywhere, so the optimum satisfies (b − z) + λb = 0, giving b = z/(1 + λ). For that to equal zero, z must equal zero. Ridge returns zero only when least squares already had. No value of λ changes this — at λ = 10¹² the coefficient measures 4.0e-13, which is very small and is not zero. The variable is still in the model, still needs collecting at serving time, still breaks the pipeline when its source goes down.

Lasso has no derivative at zero. It has a set of valid slopes there — the whole interval [−λ, λ] — and zero is optimal whenever the data’s pull fits inside it. That is why an entire range of situations collapses to exactly zero, rather than one improbable point.

Argument 4 — the corner, made quantitative

The geometric telling is the famous one. Both methods can be written as “minimise squared error, subject to keeping the coefficients inside a ball”. The L1 ball is a diamond with corners on the axes; the L2 ball is a sphere, smooth everywhere. Expand the least-squares contours until they first touch the ball, and that contact point is the answer.

A corner sitting on an axis is a coefficient of zero, and corners stick out, so they get touched first. That is the standard argument and it is correct — but it is a picture, and pictures are not evidence. So measure it: project random points onto each ball and count how often the result lands on an axis.

import numpy as np

# The standard picture: constrain the coefficients inside a ball, expand the
# least-squares contours until they touch it. The L1 ball (a diamond) has
# CORNERS on the axes; the L2 ball (a circle) is smooth everywhere. A corner
# sitting on an axis IS a zero coefficient.
#
# "Corners stick out" is a picture, not evidence. So measure it: project many
# random points onto each ball and count how often the answer lands on an axis.

rng = np.random.default_rng(1)

def project_l1(v, r):                       # onto the diamond {||b||_1 <= r}
    if np.abs(v).sum() <= r:
        return v
    u = np.sort(np.abs(v))[::-1]
    c = np.cumsum(u)
    rho = np.nonzero(u * np.arange(1, len(v) + 1) > (c - r))[0][-1]
    return np.sign(v) * np.maximum(np.abs(v) - (c[rho] - r) / (rho + 1.0), 0)

def project_l2(v, r):                       # onto the sphere {||b||_2 <= r}
    n = np.linalg.norm(v)
    return v if n <= r else v * r / n

TRIALS = 20_000
print(f"{'p':>4} {'L1 -> a zero':>14} {'L2 -> a zero':>14}")
for p in (2, 5, 20):
    l1 = sum((project_l1(rng.normal(size=p) * 2, 1.0) == 0).any()
             for _ in range(TRIALS))
    l2 = sum((project_l2(rng.normal(size=p) * 2, 1.0) == 0).any()
             for _ in range(TRIALS))
    print(f"{p:>4} {100 * l1 / TRIALS:>13.1f}% {100 * l2 / TRIALS:>13.3f}%")

# The corner is not a rare special case. By 20 dimensions it is the rule,
# because the diamond is nearly all corner. The sphere produced a zero
# exactly zero times in 60,000 attempts: landing on an axis requires aiming
# at it perfectly, and nothing makes that happen.

The corner is not an edge case that textbook diagrams exaggerate. By twenty dimensions it is the rule, because a high-dimensional diamond is almost entirely corner — and real models have far more than twenty features. The sphere produced a zero not once in 60,000 attempts, because landing on an axis requires aiming at it perfectly and nothing in the geometry encourages that.

21.18 Elastic net, and choosing between the three

Elastic net uses both penalties at once, mixing them with a ratio you choose. It exists because lasso has a specific weakness: faced with a group of correlated features it tends to pick one arbitrarily and zero the rest, and which one it picks can flip with a small change in the data. Elastic net's ridge component keeps correlated features together, so the selection is more stable.

21.19 Standardise first, always

Every penalty above sums up coefficient sizes. A coefficient's size depends on the units of its feature. Therefore the penalty depends on your units — which is absurd, since nothing about the problem changes when you record revenue in pounds instead of pence.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=200, n_features=5, noise=10.0, random_state=0)

# Simulate a realistic mixed-unit table: one column in millions, one tiny.
X_mixed = X.copy()
X_mixed[:, 0] *= 1_000_000      # e.g. revenue in currency units
X_mixed[:, 1] /= 1_000          # e.g. a rate

a = Ridge(alpha=1.0).fit(X, y).coef_
b = Ridge(alpha=1.0).fit(X_mixed, y).coef_

print("same data, different units, SAME penalty:")
print(f"  col 0 (x1,000,000): {a[0]:10.3f}  ->  {b[0]:.8f}")
print(f"  col 1 (/1,000):     {a[1]:10.3f}  ->  {b[1]:10.3f}")

# Both columns are damaged, in opposite ways, and neither is what you meant.
# The big-unit column needs only a minuscule coefficient to do its work, so
# the penalty barely touches it -- it is effectively exempt. The small-unit
# column needs a LARGE coefficient to have the same influence, and a large
# coefficient is exactly what the penalty punishes, so the model gives up on
# the feature instead.

21.20 Stepwise selection, and why to distrust it

Stepwise selection adds the most helpful feature, then the next, and so on (forward), or starts with everything and drops the least useful (backward), stopping when some criterion stops improving. It is intuitive, it is widely taught, and it has a serious flaw.

There is a second, quieter problem: stepwise is greedy. Having taken the best single feature first, it can never find a pair that is powerful together but unimpressive individually. Lasso considers all coefficients simultaneously and does not have this blind spot, which is the main reason to prefer it.

21.21 Filter, wrapper, embedded

Every variable-selection method falls into one of three families. Knowing the family tells you the cost and the failure mode.

21.22 When to remove a variable

Selection methods tell you what is predictive. They cannot tell you what is allowed, or what is stable, and those are often the decisions that matter more.

21.23 One tree, and what a split optimises

A decision tree asks a sequence of yes/no questions about one feature at a time. Everything in the tree family — forests, boosting, the whole gradient-boosted industry — is built from this one object, so it pays to be exact about what it does.

To grow a tree, the algorithm tries every threshold on every feature, scores each candidate split by how much purer it makes the two resulting groups, and keeps the best one. Then it repeats on each group. That is the whole algorithm. The only remaining question is when to stop, which turns out to be the question that matters.

21.24 Gini, entropy, and variance reduction

"Purer" needs a definition. A node is pure if everything in it belongs to one class, and maximally impure if the classes are evenly mixed.

import numpy as np

def gini(labels):
    if len(labels) == 0:
        return 0.0
    p = np.mean(labels)
    return 2 * p * (1 - p)          # for 2 classes: 1 - p^2 - (1-p)^2

def entropy(labels):
    if len(labels) == 0:
        return 0.0
    p = np.mean(labels)
    if p in (0.0, 1.0):
        return 0.0
    return -(p * np.log2(p) + (1 - p) * np.log2(1 - p))

parent = np.array([1, 1, 1, 1, 0, 0, 0, 0])          # perfectly mixed
left   = np.array([1, 1, 1, 0])
right  = np.array([1, 0, 0, 0])

def weighted(fn, l, r):
    n = len(l) + len(r)
    return (len(l) / n) * fn(l) + (len(r) / n) * fn(r)

print(f"parent  gini {gini(parent):.3f}   entropy {entropy(parent):.3f}")
print(f"split   gini {weighted(gini, left, right):.3f}   "
      f"entropy {weighted(entropy, left, right):.3f}")
print(f"gain    gini {gini(parent) - weighted(gini, left, right):.3f}   "
      f"entropy {entropy(parent) - weighted(entropy, left, right):.3f}")

# A tree tries every split on every feature and keeps the biggest gain.
# That is the entire algorithm. Everything else is when to stop.

21.25 Pruning, depth, and the stopping rules

Left alone, a tree keeps splitting until every leaf is pure. It can always achieve this, if necessary by isolating single rows. At that point it has memorised the training set.

The left panel is the shape to remember: training error falls forever as depth grows, while held-out error falls, bottoms out, and then climbs again. Everything to the right of that minimum is the model learning noise. The measured numbers on the right make it concrete — an unrestricted tree reached depth 19 with 335 leaves, scored a perfect 1.0 AUC on its training data, and managed only 0.866 on data it had not seen.

21.26 Bagging and random forests

A single deep tree has low bias and terrible variance: it can express almost any shape, but shift a few training rows and you get a visibly different tree. Averaging is the classic remedy for variance, and that is exactly what a forest is.

Bagging trains many trees on bootstrap samples — each drawn with replacement, so each tree sees a slightly different dataset — and averages their predictions. The errors are partly independent, so averaging cancels some of them.

Random forests add one more idea, and it is the important one. At every split, each tree may only consider a random subset of the features. This sounds like sabotage. It works because bagged trees are still highly correlated: if one feature is strongly predictive, every tree splits on it first and they all end up similar, so averaging buys little. Forcing trees to sometimes ignore the best feature decorrelates them, and averaging decorrelated errors cancels far more.

21.27 Extra trees, and why more randomness helps

Extremely randomised trees push the idea one step further. A random forest picks the best threshold among a random subset of features; extra trees pick a random threshold too, and keep the best of those.

Each individual tree is therefore worse. The ensemble often is not, because the trees are even less correlated and there is more variance to cancel. They are also noticeably faster to train, since scanning for the optimal split point is the expensive part and they skip it. Above, extra trees scored 0.967 against the forest's 0.964 — a real if modest edge on this data.

21.28 Boosting: AdaBoost to gradient boosting

Bagging builds trees in parallel and averages them; every tree is an equal vote and none knows the others exist. Boosting builds trees in sequence, each one trained specifically to fix what the current ensemble is still getting wrong.

AdaBoost (1995) does this by reweighting rows: after each tree, examples that were misclassified get more weight, so the next tree pays them more attention. The final prediction is a weighted vote, with more accurate trees getting louder votes.

Gradient boosting reframes that idea in a way that generalises to any differentiable loss. Instead of reweighting rows, each new tree is fitted to the negative gradient of the loss with respect to the current predictions. For squared error the gradient is just the residual, so the picture is simple: each tree predicts the leftover error of everything before it, and you add a small fraction of it — the learning rate — to the running total.

21.29 XGBoost, LightGBM, CatBoost

All three are gradient boosting. They differ in engineering, and the differences do change which one you should reach for.

21.30 Feature importance, and how it misleads

Every tree library exposes feature_importances_, and it is one of the most-quoted and least-questioned numbers in applied ML. It is computed by adding up how much each feature reduced impurity across all splits. That has a bias built into it.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(7)
n = 3000
signal    = rng.integers(0, 2, n)              # genuinely predictive
noise_bin = rng.integers(0, 2, n)              # pure noise, 2 levels
noise_id  = rng.integers(0, 1000, n)           # pure noise, 1000 levels
y = signal ^ (rng.random(n) < 0.15).astype(int)

X = np.column_stack([signal, noise_bin, noise_id]).astype(float)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.4, random_state=0)
rf = RandomForestClassifier(n_estimators=300, random_state=0,
                            n_jobs=-1).fit(Xtr, ytr)

names = ["signal", "noise (2 levels)", "noise (1000 levels)"]
print("impurity importance (the .feature_importances_ everyone uses):")
for nm, v in zip(names, rf.feature_importances_):
    print(f"  {nm:<22} {v:.3f}")

pi = permutation_importance(rf, Xte, yte, n_repeats=20,
                            random_state=0, scoring="roc_auc")
print("permutation importance, measured on HELD-OUT data:")
for nm, v in zip(names, pi.importances_mean):
    print(f"  {nm:<22} {v:+.3f}")

21.31 The confusion matrix, slowly

Every classification metric in existence is built from four counts. Get these genuinely solid and the rest of this chapter is arithmetic; stay fuzzy on them and no amount of memorising formulas will help.

import numpy as np

# Everything in classification metrics is four numbers. Build them once,
# by hand, and every formula afterwards is obvious.
y_true = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0])
y_pred = np.array([1, 1, 0, 0, 1, 0, 0, 0, 0, 0])

TP = int(((y_true == 1) & (y_pred == 1)).sum())   # said yes, was yes
FP = int(((y_true == 0) & (y_pred == 1)).sum())   # said yes, was no  (false alarm)
FN = int(((y_true == 1) & (y_pred == 0)).sum())   # said no,  was yes (miss)
TN = int(((y_true == 0) & (y_pred == 0)).sum())   # said no,  was no

print(f"            predicted 1   predicted 0")
print(f"actual 1  {TP:>11}   {FN:>11}")
print(f"actual 0  {FP:>11}   {TN:>11}")

precision = TP / (TP + FP)     # of those I FLAGGED, how many were right?
recall    = TP / (TP + FN)     # of those that EXIST, how many did I catch?
print(f"precision {precision:.2f}   recall {recall:.2f}")

# Read the denominators. Precision divides by what the MODEL said.
# Recall divides by what REALITY contains. That is the whole difference.

Two more that appear constantly and are worth naming, because ROC curves are built from them. Sensitivity is another word for recall. Specificity is recall for the negative class: TN / (TN + FP), the fraction of genuine negatives correctly left alone. The false positive rate is 1 − specificity.

21.32 Precision, recall, and the threshold

Here is the thing that reframes everything. A classifier does not output a class. It outputs a score, and somebody — usually a library default nobody examined — chose to cut it at 0.5.

Change the cut and every one of those four counts changes with it. Lower the threshold and you flag more things: recall rises, precision falls. Raise it and the reverse. Precision and recall are not properties of a model. They are properties of a model at a chosen operating point, and quoting them without the threshold is meaningless.

import numpy as np
from sklearn.metrics import (accuracy_score, balanced_accuracy_score,
                             matthews_corrcoef, f1_score)

# 5% of rows are positive. A model that predicts "no" for everything.
y = np.array([0] * 95 + [1] * 5)
lazy = np.zeros(100, dtype=int)

print(f"accuracy          {accuracy_score(y, lazy):.3f}   <- looks great")
print(f"balanced accuracy {balanced_accuracy_score(y, lazy):.3f}   <- honest")
print(f"F1                {f1_score(y, lazy, zero_division=0):.3f}")
print(f"MCC               {matthews_corrcoef(y, lazy):.3f}   <- 0 = no skill")

# Accuracy rewards the majority class for existing. The other three do not.

21.33 F-beta, and choosing beta

You often need one number to compare models. F-beta combines precision and recall into one, weighting recall as beta times more important than precision.

def fbeta(precision, recall, beta):
    """F-beta is a weighted harmonic mean. beta = how many times more
    you care about recall than precision."""
    b2 = beta * beta
    return (1 + b2) * precision * recall / (b2 * precision + recall)

p, r = 0.80, 0.40      # precise but misses a lot

for beta in (0.5, 1.0, 2.0):
    print(f"beta {beta}: F = {fbeta(p, r, beta):.3f}")

# Sanity: beta=1 is the plain harmonic mean, sitting nearer the WORSE value.
print(f"arithmetic mean would flatter it: {(p + r) / 2:.3f}")
print(f"F1 refuses to:                    {fbeta(p, r, 1.0):.3f}")

It is a harmonic mean, not an arithmetic one, and that is deliberate: the harmonic mean is dragged toward the smaller of the two values. Precision 0.8 with recall 0.4 gives an arithmetic mean of 0.6 but an F1 of 0.53. You cannot rescue a bad recall with an excellent precision, which is exactly the behaviour you want from a summary score.

21.34 ROC-AUC versus PR-AUC

Both summarise a classifier across all thresholds at once, which is what makes them useful for comparing models before you have chosen an operating point. They differ in what they plot, and under imbalance they tell strikingly different stories.

The ROC curve plots recall against the false positive rate. ROC-AUC has a genuinely intuitive reading: take one random positive and one random negative; ROC-AUC is the probability the model scores the positive higher. 0.5 is a coin flip, 1.0 is perfect.

The PR curve plots precision against recall. Its baseline is not 0.5 — a random classifier scores the prevalence, so on 1% positives the floor is 0.01.

21.35 The metrics nobody teaches you

Beyond the standard set are several that solve specific problems well, and knowing one or two signals real experience.

21.36 Calibration, and why ranking is not enough

Every metric so far only cares about order. Push all your predicted probabilities through any monotonic squashing function and ROC-AUC, PR-AUC and precision@k are all completely unchanged.

That is fine if you only need a ranking. It is not fine the moment anyone multiplies your output by money.

Calibration asks the direct question: of all the cases you scored 0.10, did about 10% actually happen? Plot predicted probability against observed frequency in bins; a calibrated model sits on the diagonal. Measure it with the Brier score (mean squared error of the probabilities — 0.0379 for the logistic model used in 21.37) or with log loss.

21.37 Choosing a threshold with money

This is where the chapter has been heading. Metrics are proxies. The thing you actually want to minimise is cost, and if you can put numbers on the two kinds of mistake you can skip the proxy entirely.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=20000, n_features=8, n_informative=5,
                           weights=[0.95, 0.05], random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0)
model = LogisticRegression(max_iter=2000).fit(Xtr, ytr)
p = model.predict_proba(Xte)[:, 1]

# Put money on the outcomes. THIS is the conversation to have with the
# business, instead of arguing about which metric is "best".
COST_FN = 50.0     # a missed fraud: we refund the customer
COST_FP = 4.0      # a false alarm: manual review by an analyst

def total_cost(t):
    fn = int(((yte == 1) & (p < t)).sum())
    fp = int(((yte == 0) & (p >= t)).sum())
    return fn * COST_FN + fp * COST_FP

grid = np.linspace(0.01, 0.99, 99)
costs = [total_cost(t) for t in grid]
best = float(grid[int(np.argmin(costs))])

print(f"default threshold 0.50 costs {total_cost(0.5):>9,.0f}")
print(f"best threshold    {best:.2f} costs {total_cost(best):>9,.0f}")
print(f"saving:                      {total_cost(0.5) - total_cost(best):>9,.0f}")

# For well-calibrated probabilities the optimum is analytic:
print(f"theory says: {COST_FP / (COST_FP + COST_FN):.3f}")

21.38 A decision guide

Everything above, compressed into the questions actually worth asking.

21.39 Key takeaways

21.40 Interview drills

Fourteen questions in the order the chapter built them. Try to answer before opening each one — the answers include the traps interviewers are usually probing for.

Expand all