17.9 ended with a measurement: thirty unstable trees, individually mediocre, whose majority vote beat every one of them. This chapter explains why that works, derives the exact condition under which it works, and builds the two algorithms that exploit it.
The idea is narrow and powerful. Averaging reduces variance without increasing bias — so if you can build many low-bias, high-variance models whose errors are not identical, averaging them gives you the low bias for free.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import time
from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
from sklearn.ensemble import (BaggingClassifier, BaggingRegressor, RandomForestClassifier,
RandomForestRegressor, ExtraTreesClassifier, VotingClassifier)
from sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split
from sklearn.metrics import accuracy_score, roc_auc_score
RANDOM_STATE = 42
rng = np.random.default_rng(RANDOM_STATE)
np.set_printoptions(precision=4, suppress=True)
pd.set_option("display.width", 115)
sns.set_theme(style="whitegrid")
plt.rcParams["figure.dpi"] = 80
cv = StratifiedKFold(5, shuffle=True, random_state=RANDOM_STATE)
print("ready")
ready
18.1 Why Ensembles Work
18.1.1 The variance of an average
Take \(k\) estimators, each with variance \(\sigma^2\). If they were independent, the variance of their average would be
\[\mathrm{Var}\!\left(\frac{1}{k}\sum_{i=1}^{k} f_i\right) = \frac{\sigma^2}{k}\]
Ten independent models would cut the variance tenfold. But models fitted on the same data are not independent — they make correlated mistakes. With pairwise correlation \(\rho\) the formula becomes
\[\boxed{\;\mathrm{Var} = \rho\,\sigma^2 + \frac{1 - \rho}{k}\,\sigma^2\;}\]
Read the two terms separately, because the whole of this chapter is in the difference between them:
| Term | Behaviour as \(k \to \infty\) | Controlled by |
|---|---|---|
| \(\rho\,\sigma^2\) | Does not shrink at all | How different the models are |
| \(\frac{1-\rho}{k}\sigma^2\) | Vanishes | How many models there are |
More trees cannot reduce variance below \(\rho\sigma^2\). That floor is set by correlation, which is why 18.6 spends its effort on decorrelating the trees rather than adding more.
def ensemble_variance(rho, k, sigma_sq=1.0):
return rho * sigma_sq + (1 - rho) / k * sigma_sq
print(f"{'k':>6}" + "".join(f"{f'rho={r}':>12}" for r in (0.0, 0.3, 0.6, 0.9)))
for k in (1, 5, 10, 50, 500, 100_000):
row = "".join(f"{ensemble_variance(r, k):>12.4f}" for r in (0.0, 0.3, 0.6, 0.9))
print(f"{k:>6}{row}")
print("\nat rho=0.9 the variance is stuck at 0.9 however many models you add.")
k rho=0.0 rho=0.3 rho=0.6 rho=0.9
1 1.0000 1.0000 1.0000 1.0000
5 0.2000 0.4400 0.6800 0.9200
10 0.1000 0.3700 0.6400 0.9100
50 0.0200 0.3140 0.6080 0.9020
500 0.0020 0.3014 0.6008 0.9002
100000 0.0000 0.3000 0.6000 0.9000
at rho=0.9 the variance is stuck at 0.9 however many models you add.Verify the formula by simulation rather than trusting the algebra:
def simulate_average_variance(rho, k, n_trials=4000, seed=0):
"""Generate k estimators with pairwise correlation rho, average them, measure the variance."""
gen = np.random.default_rng(seed)
shared = gen.normal(size=(n_trials, 1)) # the common component
private = gen.normal(size=(n_trials, k)) # the independent component
estimators = np.sqrt(rho) * shared + np.sqrt(1 - rho) * private
return estimators.mean(axis=1).var()
print(f"{'rho':>6}{'k':>6}{'predicted':>12}{'simulated':>12}")
for rho in (0.0, 0.3, 0.9):
for k in (5, 50):
predicted = ensemble_variance(rho, k)
simulated = simulate_average_variance(rho, k)
print(f"{rho:>6.1f}{k:>6}{predicted:>12.4f}{simulated:>12.4f}")
rho k predicted simulated 0.0 5 0.2000 0.2022 0.0 50 0.0200 0.0205 0.3 5 0.4400 0.4348 0.3 50 0.3140 0.3143 0.9 5 0.9200 0.9134 0.9 50 0.9020 0.8997
18.1.2 The classification version
For voting classifiers the equivalent result is the Condorcet jury theorem: if each voter is right with probability \(p > 0.5\) and votes independently, the majority of \(k\) voters is right with probability approaching 1 as \(k\) grows.
from scipy.stats import binom
def majority_correct(p, k):
"""P(more than half of k independent voters are correct)."""
return 1 - binom.cdf(k // 2, k, p)
print(f"{'voters':>8}" + "".join(f"{f'p={p}':>10}" for p in (0.45, 0.55, 0.60, 0.70)))
for k in (1, 5, 11, 51, 201):
row = "".join(f"{majority_correct(p, k):>10.4f}" for p in (0.45, 0.55, 0.60, 0.70))
print(f"{k:>8}{row}")
voters p=0.45 p=0.55 p=0.6 p=0.7
1 0.4500 0.5500 0.6000 0.7000
5 0.4069 0.5931 0.6826 0.8369
11 0.3669 0.6331 0.7535 0.9218
51 0.2359 0.7641 0.9265 0.9986
201 0.0774 0.9226 0.9979 1.0000Two things to read off. A jury of weak-but-better-than-chance voters becomes very strong: at \(p = 0.55\),
201 voters reach 92%. And the p=0.45 column shows the failure mode — voters worse than chance get
worse with numbers, collapsing towards 0.
So an ensemble needs both conditions: members better than chance, and errors that are not identical. A thousand copies of the same model is still that one model.
from sklearn.datasets import make_classification
X_id, y_id = make_classification(n_samples=500, n_features=10, n_informative=6,
random_state=RANDOM_STATE)
identical = [DecisionTreeClassifier(max_depth=4, random_state=0).fit(X_id, y_id)
for _ in range(50)]
diverse = [DecisionTreeClassifier(max_depth=4, random_state=s).fit(
*(lambda idx: (X_id[idx], y_id[idx]))(np.random.default_rng(s).choice(len(X_id), len(X_id))))
for s in range(50)]
id_preds = np.array([t.predict(X_id) for t in identical])
dv_preds = np.array([t.predict(X_id) for t in diverse])
print(f"50 identical trees - distinct prediction vectors: {len(np.unique(id_preds, axis=0))}")
print(f"50 bootstrapped trees - distinct prediction vectors: {len(np.unique(dv_preds, axis=0))}")
print(f"\nmean pairwise agreement, identical : {np.mean(id_preds[0] == id_preds[1:]):.4f}")
print(f"mean pairwise agreement, bootstrapped: {np.mean(dv_preds[0] == dv_preds[1:]):.4f}")
50 identical trees - distinct prediction vectors: 1 50 bootstrapped trees - distinct prediction vectors: 50 mean pairwise agreement, identical : 1.0000 mean pairwise agreement, bootstrapped: 0.8655
18.2 Voting and Averaging
The simplest ensemble combines different model families fitted on the same data. Diversity comes from the algorithms disagreeing rather than from the data being resampled.
| Method | Combines | scikit-learn |
|---|---|---|
| Hard voting | Predicted labels, by majority | VotingClassifier(voting="hard") |
| Soft voting | Predicted probabilities, by average | VotingClassifier(voting="soft") |
| Averaging | Predicted values | VotingRegressor |
Soft voting is usually better, because it uses confidence rather than discarding it — but only if the members' probabilities are comparable, which 15.8.3 showed is not automatic.
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
X_v, y_v = make_classification(n_samples=800, n_features=12, n_informative=7,
class_sep=0.9, random_state=RANDOM_STATE)
members = [
("logistic", make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))),
("tree", DecisionTreeClassifier(max_depth=6, random_state=RANDOM_STATE)),
("knn", make_pipeline(StandardScaler(), KNeighborsClassifier(15))),
("naive bayes", GaussianNB()),
]
for name, model in members:
print(f"{name:<14} CV accuracy {cross_val_score(model, X_v, y_v, cv=cv).mean():.4f}")
for voting in ("hard", "soft"):
ens = VotingClassifier(members, voting=voting)
print(f"\n{voting} voting CV accuracy {cross_val_score(ens, X_v, y_v, cv=cv).mean():.4f}")
logistic CV accuracy 0.7725 tree CV accuracy 0.7863 knn CV accuracy 0.8812 naive bayes CV accuracy 0.7412 hard voting CV accuracy 0.8112 soft voting CV accuracy 0.8250
The ensemble is only worth it if the members disagree. Measure that directly rather than assuming it.
X_tr, X_te, y_tr, y_te = train_test_split(X_v, y_v, test_size=0.3, stratify=y_v,
random_state=RANDOM_STATE)
fitted = {name: model.fit(X_tr, y_tr) for name, model in members}
preds = {name: m.predict(X_te) for name, m in fitted.items()}
names = list(preds)
agreement = pd.DataFrame(
[[np.mean(preds[a] == preds[b]) for b in names] for a in names],
index=names, columns=names,
)
print("pairwise prediction agreement on the test set:")
print(agreement.round(3).to_string())
print(f"\nrows where all four agree: {np.mean(np.all([preds[n] == preds[names[0]] for n in names], axis=0)):.1%}")
pairwise prediction agreement on the test set:
logistic tree knn naive bayes
logistic 1.000 0.717 0.817 0.838
tree 0.717 1.000 0.775 0.729
knn 0.817 0.775 1.000 0.796
naive bayes 0.838 0.729 0.796 1.000
rows where all four agree: 60.0%18.3 The Bootstrap
Voting needs different models. Bagging (bootstrap aggregating) uses one model family and generates diversity from the data instead, by resampling it.
A bootstrap sample draws \(n\) rows from a dataset of size \(n\) with replacement. Some rows appear several times, and — this is the useful part — some do not appear at all.
The probability a given row is never selected in \(n\) draws is
\[\left(1 - \frac{1}{n}\right)^n \;\xrightarrow[n \to \infty]{}\; e^{-1} \approx 0.368\]
So each bootstrap sample contains about 63.2% of the distinct rows, and leaves 36.8% out. Those left-out rows are what makes 18.5's free validation possible.
print(f"{'n':>8}{'(1 - 1/n)^n':>16}{'unique fraction':>18}{'empirical':>12}")
for n in (5, 20, 100, 1000, 10_000):
theoretical_out = (1 - 1 / n) ** n
gen = np.random.default_rng(RANDOM_STATE)
empirical = np.mean([len(np.unique(gen.integers(0, n, n))) / n for _ in range(200)])
print(f"{n:>8}{theoretical_out:>16.4f}{1 - theoretical_out:>18.4f}{empirical:>12.4f}")
print(f"\n1 - 1/e = {1 - np.exp(-1):.4f}")
n (1 - 1/n)^n unique fraction empirical
5 0.3277 0.6723 0.6620
20 0.3585 0.6415 0.6385
100 0.3660 0.6340 0.6314
1000 0.3677 0.6323 0.6326
10000 0.3679 0.6321 0.6323
1 - 1/e = 0.6321# What a single bootstrap sample looks like on 10 rows.
gen = np.random.default_rng(RANDOM_STATE)
original = np.arange(10)
sample = gen.integers(0, 10, 10)
print("original :", original)
print("bootstrap:", np.sort(sample))
print(f"\nin-bag (appear at least once): {sorted(set(sample))}")
print(f"out-of-bag (never selected) : {sorted(set(original) - set(sample))}")
print(f"counts: {dict(zip(*np.unique(sample, return_counts=True)))}")
original : [0 1 2 3 4 5 6 7 8 9]
bootstrap: [0 0 0 2 4 4 6 6 7 8]
in-bag (appear at least once): [np.int64(0), np.int64(2), np.int64(4), np.int64(6), np.int64(7), np.int64(8)]
out-of-bag (never selected) : [np.int64(1), np.int64(3), np.int64(5), np.int64(9)]
counts: {np.int64(0): np.int64(3), np.int64(2): np.int64(1), np.int64(4): np.int64(2), np.int64(6): np.int64(2), np.int64(7): np.int64(1), np.int64(8): np.int64(1)}18.4 Bagging from Scratch
Fit one model per bootstrap sample; average (or vote) their predictions. That is the whole algorithm.
class BaggingFromScratch:
"""Bootstrap aggregating for classification."""
def __init__(self, base_factory, n_estimators=50, random_state=0):
self.base_factory = base_factory
self.n_estimators = n_estimators
self.random_state = random_state
def fit(self, X, y):
gen = np.random.default_rng(self.random_state)
n = len(X)
self.estimators_, self.oob_masks_ = [], []
for _ in range(self.n_estimators):
idx = gen.integers(0, n, n) # bootstrap sample
oob = np.ones(n, dtype=bool)
oob[idx] = False # everything never drawn
self.estimators_.append(self.base_factory().fit(X[idx], y[idx]))
self.oob_masks_.append(oob)
self.classes_ = np.unique(y)
return self
def predict_proba(self, X):
votes = np.array([m.predict(X) for m in self.estimators_])
return np.array([(votes == c).mean(axis=0) for c in self.classes_]).T
def predict(self, X):
return self.classes_[self.predict_proba(X).argmax(axis=1)]
X_b, y_b = make_classification(n_samples=800, n_features=12, n_informative=7,
class_sep=0.8, random_state=RANDOM_STATE)
Xb_tr, Xb_te, yb_tr, yb_te = train_test_split(X_b, y_b, test_size=0.3, stratify=y_b,
random_state=RANDOM_STATE)
mine = BaggingFromScratch(
lambda: DecisionTreeClassifier(random_state=RANDOM_STATE), n_estimators=100
).fit(Xb_tr, yb_tr)
single = DecisionTreeClassifier(random_state=RANDOM_STATE).fit(Xb_tr, yb_tr)
print(f"single tree test accuracy {single.score(Xb_te, yb_te):.4f}")
print(f"100 bagged trees test accuracy {accuracy_score(yb_te, mine.predict(Xb_te)):.4f}")
single tree test accuracy 0.7333 100 bagged trees test accuracy 0.8292
Step 4 of the contract.
theirs = BaggingClassifier(
DecisionTreeClassifier(random_state=RANDOM_STATE),
n_estimators=100, random_state=RANDOM_STATE,
).fit(Xb_tr, yb_tr)
print(f"from scratch accuracy {accuracy_score(yb_te, mine.predict(Xb_te)):.4f} "
f"AUC {roc_auc_score(yb_te, mine.predict_proba(Xb_te)[:, 1]):.4f}")
print(f"scikit-learn accuracy {theirs.score(Xb_te, yb_te):.4f} "
f"AUC {roc_auc_score(yb_te, theirs.predict_proba(Xb_te)[:, 1]):.4f}")
print(f"\npredictions agree on {np.mean(mine.predict(Xb_te) == theirs.predict(Xb_te)):.1%} of test rows")
from scratch accuracy 0.8292 AUC 0.8937 scikit-learn accuracy 0.8292 AUC 0.8955 predictions agree on 93.3% of test rows
The two are not bit-identical, and cannot be: they draw different bootstrap samples from different random streams. This is the right kind of agreement for a stochastic algorithm — the same procedure, the same performance, different random draws. Compare with 17.3, where a deterministic algorithm had to match exactly.
Watch the variance reduction happen as trees are added:
print(f"{'trees':>8}{'test accuracy':>16}{'test AUC':>12}")
for k in (1, 3, 10, 30, 100, 300):
ens = BaggingFromScratch(
lambda: DecisionTreeClassifier(random_state=RANDOM_STATE), n_estimators=k
).fit(Xb_tr, yb_tr)
proba = ens.predict_proba(Xb_te)[:, 1]
print(f"{k:>8}{accuracy_score(yb_te, ens.predict(Xb_te)):>16.4f}"
f"{roc_auc_score(yb_te, proba):>12.4f}")
trees test accuracy test AUC
1 0.6458 0.6456
3 0.7250 0.7795
10 0.7958 0.8710
30 0.8208 0.8927
100 0.8292 0.8937
300 0.8292 0.8984Bagging does not overfit as trees are added. The curve rises and flattens; it does not turn back down.
That is 18.1.1's formula: adding members shrinks the \(\frac{1-\rho}{k}\sigma^2\) term and touches nothing
else. n_estimators is therefore a budget parameter, not a complexity parameter — one of the few in
this material that cannot be set too high.
18.5 Out-of-Bag Evaluation
18.3 established that each bootstrap sample omits about 36.8% of the rows. For any given row, roughly a third of the trees never saw it — so those trees can predict it as if it were held out.
Averaging those predictions gives the out-of-bag score: a validation estimate that costs nothing and uses no extra data.
def oob_score(ensemble, X, y):
"""Average the predictions of only those trees that did not see each row."""
n = len(X)
votes = np.full((len(ensemble.estimators_), n), np.nan)
for i, (model, oob) in enumerate(zip(ensemble.estimators_, ensemble.oob_masks_)):
votes[i, oob] = model.predict(X[oob])
with np.errstate(invalid="ignore"):
mean_vote = np.nanmean(votes, axis=0)
covered = ~np.isnan(mean_vote)
predicted = (mean_vote[covered] >= 0.5).astype(int)
return (predicted == y[covered]).mean(), covered.mean()
score, coverage = oob_score(mine, Xb_tr, yb_tr)
print(f"out-of-bag accuracy : {score:.4f} (on {coverage:.1%} of training rows)")
print(f"training accuracy : {accuracy_score(yb_tr, mine.predict(Xb_tr)):.4f} <- optimistic")
print(f"TEST accuracy : {accuracy_score(yb_te, mine.predict(Xb_te)):.4f}")
out-of-bag accuracy : 0.8089 (on 100.0% of training rows) training accuracy : 1.0000 <- optimistic TEST accuracy : 0.8292
The OOB estimate lands close to the test score and far from the training score, which is exactly what a
validation estimate should do. scikit-learn computes it for free with oob_score=True.
with_oob = BaggingClassifier(
DecisionTreeClassifier(random_state=RANDOM_STATE),
n_estimators=200, oob_score=True, random_state=RANDOM_STATE,
).fit(Xb_tr, yb_tr)
t0 = time.perf_counter()
cv_score = cross_val_score(
BaggingClassifier(DecisionTreeClassifier(random_state=RANDOM_STATE),
n_estimators=200, random_state=RANDOM_STATE),
Xb_tr, yb_tr, cv=cv).mean()
cv_ms = (time.perf_counter() - t0) * 1000
print(f"OOB score : {with_oob.oob_score_:.4f} (free - one fit)")
print(f"5-fold CV score : {cv_score:.4f} ({cv_ms:.0f} ms - five fits)")
print(f"test score : {with_oob.score(Xb_te, yb_te):.4f}")
OOB score : 0.8071 (free - one fit) 5-fold CV score : 0.7821 (11884 ms - five fits) test score : 0.8458
| Out-of-bag | Cross-validation | |
|---|---|---|
| Cost | Free — one fit | \(k\) fits |
| Available for | Bagged models only | Anything |
| Estimates | A model trained on ~63% of the data | A model trained on \(\frac{k-1}{k}\) of it |
| Slight bias | Pessimistic — fewer effective trees per row | Small |
| Usable for tuning | Yes, and very cheap (18.9.1) | Yes |
Common mistake — using OOB with a small n_estimators. Each row is scored by only ~37% of the trees,
so with 10 trees that is about 4 votes, and the estimate is extremely noisy.
import warnings
print(f"{'n_estimators':>14}{'OOB score':>12}{'test score':>13}{'trees per OOB row':>20}{'warned':>9}")
for k in (5, 10, 30, 100, 400):
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
m = BaggingClassifier(DecisionTreeClassifier(random_state=RANDOM_STATE),
n_estimators=k, oob_score=True,
random_state=RANDOM_STATE).fit(Xb_tr, yb_tr)
warned = any("oob" in str(w.message).lower() for w in caught)
print(f"{k:>14}{m.oob_score_:>12.4f}{m.score(Xb_te, yb_te):>13.4f}"
f"{k * 0.368:>20.1f}{str(warned):>9}")
n_estimators OOB score test score trees per OOB row warned
5 0.7107 0.7792 1.8 True
10 0.7518 0.7917 3.7 True
30 0.7821 0.8250 11.0 False
100 0.7911 0.8292 36.8 False
400 0.8054 0.8583 147.2 Falsescikit-learn raises a UserWarning at the smallest sizes — "Some inputs do not have OOB scores" — which
is the library telling you the estimate is unreliable. Read it rather than suppressing it. A few
hundred trees puts the OOB estimate on solid ground; a dozen does not.
18.6 Random Forests
Bagged trees are still correlated: if one feature is strongly predictive, every tree splits on it at the root, and the trees end up similar. 18.1.1 says that correlation sets a floor on how much variance averaging can remove.
A random forest adds one idea to break it: at every split, consider only a random subset of the features.
| Source of diversity | Bagging | Random forest |
|---|---|---|
| Bootstrap resampling of rows | ✓ | ✓ |
| Random feature subset per split | ✗ | ✓ |
The effect on correlation is measurable.
def mean_pairwise_correlation(models, X):
"""Correlation between the ensemble members' predicted probabilities."""
P = np.array([m.predict_proba(X)[:, 1] for m in models])
C = np.corrcoef(P)
return C[np.triu_indices_from(C, k=1)].mean()
bagged = BaggingClassifier(DecisionTreeClassifier(random_state=RANDOM_STATE),
n_estimators=60, random_state=RANDOM_STATE).fit(Xb_tr, yb_tr)
forest = RandomForestClassifier(n_estimators=60, random_state=RANDOM_STATE).fit(Xb_tr, yb_tr)
print(f"mean pairwise correlation, bagged trees : "
f"{mean_pairwise_correlation(bagged.estimators_, Xb_te):.4f}")
print(f"mean pairwise correlation, forest trees : "
f"{mean_pairwise_correlation(forest.estimators_, Xb_te):.4f}")
print(f"\ntest accuracy, bagging : {bagged.score(Xb_te, yb_te):.4f}")
print(f"test accuracy, random forest : {forest.score(Xb_te, yb_te):.4f}")
mean pairwise correlation, bagged trees : 0.3520 mean pairwise correlation, forest trees : 0.2978 test accuracy, bagging : 0.8333 test accuracy, random forest : 0.8375
Lower correlation, better score — the mechanism 18.1.1 predicted, confirmed on real trees.
The diversity is visible in what the trees actually split on. Bagged trees nearly all pick the same root; a forest's trees are forced to explore.
bagged_roots = [t.tree_.feature[0] for t in bagged.estimators_]
forest_roots = [t.tree_.feature[0] for t in forest.estimators_]
print("root feature chosen across 60 trees:")
print(f" bagging : {len(set(bagged_roots))} distinct, most common used "
f"{max(np.bincount(bagged_roots)) / 60:.0%} of the time")
print(f" random forest : {len(set(forest_roots))} distinct, most common used "
f"{max(np.bincount(forest_roots)) / 60:.0%} of the time")
root feature chosen across 60 trees: bagging : 3 distinct, most common used 97% of the time random forest : 9 distinct, most common used 28% of the time
18.6.1 max_features: the knob that controls correlation
max_features |
Effect | Trade |
|---|---|---|
| All features | Equivalent to bagging | Highest correlation, lowest individual bias |
| \(\sqrt{p}\) | scikit-learn's classification default | Balanced |
| \(p/3\) | Traditional regression default | Balanced |
| 1 | Maximum decorrelation | Individual trees become very weak |
p_features = Xb_tr.shape[1]
print(f"{'max_features':>14}{'correlation':>14}{'single-tree acc':>18}{'ensemble acc':>15}")
for mf in (1, 2, 3, "sqrt", 6, p_features):
f = RandomForestClassifier(n_estimators=60, max_features=mf,
random_state=RANDOM_STATE).fit(Xb_tr, yb_tr)
corr = mean_pairwise_correlation(f.estimators_, Xb_te)
solo = np.mean([t.score(Xb_te, yb_te) for t in f.estimators_[:20]])
label = mf if isinstance(mf, str) else str(mf)
print(f"{label:>14}{corr:>14.4f}{solo:>18.4f}{f.score(Xb_te, yb_te):>15.4f}")
max_features correlation single-tree acc ensemble acc
1 0.1843 0.6510 0.8125
2 0.2514 0.6944 0.8375
3 0.2978 0.6902 0.8375
sqrt 0.2978 0.6902 0.8375
6 0.3362 0.7063 0.8333
12 0.3598 0.7196 0.8292Read the three numeric columns together: as max_features falls, correlation falls and individual trees
get worse, yet the ensemble can improve. That is the trade 18.1.1's formula describes — accepting more
bias per member to remove the \(\rho\sigma^2\) floor.
18.6.2 Forests barely need tuning
n_estimators cannot overfit (18.4), and the remaining knobs have mild effects. That combination —
strong out of the box, hard to break — is why the random forest is the standard first model on tabular
data.
from sklearn.datasets import load_breast_cancer
Xc, yc = load_breast_cancer(return_X_y=True)
configs = {
"defaults": RandomForestClassifier(random_state=RANDOM_STATE),
"1000 trees": RandomForestClassifier(n_estimators=1000, random_state=RANDOM_STATE),
"max_depth=3": RandomForestClassifier(max_depth=3, random_state=RANDOM_STATE),
"min_samples_leaf=20": RandomForestClassifier(min_samples_leaf=20, random_state=RANDOM_STATE),
"max_features=1": RandomForestClassifier(max_features=1, random_state=RANDOM_STATE),
}
for name, model in configs.items():
print(f"{name:<22} CV ROC AUC {cross_val_score(model, Xc, yc, cv=cv, scoring='roc_auc').mean():.4f}")
defaults CV ROC AUC 0.9885 1000 trees CV ROC AUC 0.9907 max_depth=3 CV ROC AUC 0.9888 min_samples_leaf=20 CV ROC AUC 0.9861 max_features=1 CV ROC AUC 0.9935
18.7 Extremely Randomized Trees
Extra Trees pushes randomization one step further: instead of searching for the best threshold on each candidate feature, it picks thresholds at random and keeps the best of those.
| Random forest | Extra Trees | |
|---|---|---|
| Rows | Bootstrap sample | All of them (by default) |
| Features per split | Random subset | Random subset |
| Threshold | Best found by search | Random, then best of the random ones |
| Bias | Lower | Higher |
| Variance | Higher | Lower |
| Fitting speed | Slower | Faster — no threshold search |
for name, model in (("random forest", RandomForestClassifier(n_estimators=300,
random_state=RANDOM_STATE)),
("extra trees", ExtraTreesClassifier(n_estimators=300,
random_state=RANDOM_STATE))):
t0 = time.perf_counter(); model.fit(Xb_tr, yb_tr); fit_ms = (time.perf_counter() - t0) * 1000
corr = mean_pairwise_correlation(model.estimators_[:60], Xb_te)
print(f"{name:<16} fit {fit_ms:>7.0f} ms correlation {corr:.4f} "
f"test accuracy {model.score(Xb_te, yb_te):.4f}")
random forest fit 1783 ms correlation 0.2978 test accuracy 0.8375 extra trees fit 1298 ms correlation 0.2497 test accuracy 0.8375
print(f"{'dataset':<22}{'random forest':>16}{'extra trees':>14}")
datasets = {
"breast cancer": load_breast_cancer(return_X_y=True),
"synthetic (n=800)": (X_b, y_b),
"noisy synthetic": (np.column_stack([X_b, rng.normal(size=(len(X_b), 40))]), y_b),
}
for name, (Xd, yd) in datasets.items():
rf = cross_val_score(RandomForestClassifier(n_estimators=200, random_state=RANDOM_STATE),
Xd, yd, cv=cv).mean()
et = cross_val_score(ExtraTreesClassifier(n_estimators=200, random_state=RANDOM_STATE),
Xd, yd, cv=cv).mean()
print(f"{name:<22}{rf:>16.4f}{et:>14.4f}")
dataset random forest extra trees breast cancer 0.9543 0.9684 synthetic (n=800) 0.8512 0.8775 noisy synthetic 0.7787 0.7612
Neither dominates. Extra Trees is worth trying whenever a forest is — it costs one line and fits faster — and 9.8's no-free-lunch argument says which wins is an empirical question, not a principled one.
18.8 Feature Importance: impurity versus permutation
17.8.1 showed impurity-based importance favouring high-cardinality features even when they are pure noise. Averaging over a forest does not fix this — it averages the same bias over more trees.
n_bias = 1000
y_bias = rng.integers(0, 2, n_bias) # the target is a coin flip
X_bias = pd.DataFrame({
"binary_noise": rng.integers(0, 2, n_bias),
"few_categories": rng.integers(0, 5, n_bias),
"many_categories": rng.integers(0, 50, n_bias),
"continuous_noise": rng.normal(size=n_bias),
})
biased_forest = RandomForestClassifier(n_estimators=300, random_state=RANDOM_STATE).fit(X_bias, y_bias)
print("every feature is noise; the target is random:")
print(pd.Series(biased_forest.feature_importances_, index=X_bias.columns).round(4).to_string())
print(f"\ndistinct values: {X_bias.nunique().to_dict()}")
print("importance still tracks cardinality, with 300 trees instead of 1.")
every feature is noise; the target is random:
binary_noise 0.0439
few_categories 0.0973
many_categories 0.3307
continuous_noise 0.5281
distinct values: {'binary_noise': 2, 'few_categories': 5, 'many_categories': 50, 'continuous_noise': 1000}
importance still tracks cardinality, with 300 trees instead of 1.The correct tool is permutation importance (24.3): shuffle one feature's values in a held-out set and measure how much the score drops. It is model-agnostic, uses validation data, and has no cardinality bias.
from sklearn.inspection import permutation_importance
Xbi_tr, Xbi_te, ybi_tr, ybi_te = train_test_split(X_bias, y_bias, test_size=0.3,
random_state=RANDOM_STATE)
f = RandomForestClassifier(n_estimators=300, random_state=RANDOM_STATE).fit(Xbi_tr, ybi_tr)
perm = permutation_importance(f, Xbi_te, ybi_te, n_repeats=20, random_state=RANDOM_STATE)
print(pd.DataFrame({
"impurity": f.feature_importances_,
"permutation": perm.importances_mean,
"permutation std": perm.importances_std,
}, index=X_bias.columns).round(4).to_string())
print("\npermutation importance correctly reports that none of them help.")
impurity permutation permutation std binary_noise 0.0432 0.0107 0.0191 few_categories 0.1085 -0.0100 0.0231 many_categories 0.3378 -0.0003 0.0179 continuous_noise 0.5105 0.0102 0.0225 permutation importance correctly reports that none of them help.
18.8.1 Correlated features split the credit
A second failure affects both methods, and it is the one that most often misleads in real work: when two features carry the same information, each can look unimportant because the other covers for it.
n_corr = 800
signal = rng.normal(size=n_corr)
X_corr = pd.DataFrame({
"original": signal,
"duplicate": signal + rng.normal(0, 0.01, n_corr),
"noise": rng.normal(size=n_corr),
})
y_corr = (signal + rng.normal(0, 0.3, n_corr) > 0).astype(int)
Xco_tr, Xco_te, yco_tr, yco_te = train_test_split(X_corr, y_corr, test_size=0.3,
random_state=RANDOM_STATE)
with_dup = RandomForestClassifier(n_estimators=300, random_state=RANDOM_STATE).fit(Xco_tr, yco_tr)
perm_dup = permutation_importance(with_dup, Xco_te, yco_te, n_repeats=20,
random_state=RANDOM_STATE)
print("with BOTH copies present:")
print(pd.DataFrame({"impurity": with_dup.feature_importances_,
"permutation": perm_dup.importances_mean},
index=X_corr.columns).round(4).to_string())
with BOTH copies present:
impurity permutation
original 0.4448 0.0513
duplicate 0.4755 0.3473
noise 0.0798 -0.0198# Drop the duplicate and re-measure.
single_cols = ["original", "noise"]
without_dup = RandomForestClassifier(n_estimators=300,
random_state=RANDOM_STATE).fit(Xco_tr[single_cols], yco_tr)
perm_single = permutation_importance(without_dup, Xco_te[single_cols], yco_te,
n_repeats=20, random_state=RANDOM_STATE)
print("with only ONE copy present:")
print(pd.DataFrame({"impurity": without_dup.feature_importances_,
"permutation": perm_single.importances_mean},
index=single_cols).round(4).to_string())
print(f"\naccuracy is unchanged: {with_dup.score(Xco_te, yco_te):.4f} vs "
f"{without_dup.score(Xco_te[single_cols], yco_te):.4f}")
with only ONE copy present:
impurity permutation
original 0.8765 0.4094
noise 0.1235 -0.0121
accuracy is unchanged: 0.8917 vs 0.9042Permutation importance collapses under duplication. original's permutation importance falls from
0.409 alone to 0.051 when its duplicate is present — an eightfold drop — because shuffling it does
almost nothing while the duplicate still carries the signal. Remove the duplicate and the true importance
reappears, with the model's accuracy unchanged throughout.
Note that the collapse is asymmetric: duplicate retained 0.347 while original lost almost
everything. Which copy keeps the credit is arbitrary — it depends on which one the forest happened to
split on more often — so neither number describes the feature, and a ranking built from them is
misleading in a way that reruns will not reveal.
| Method | Cardinality bias | Correlated features | Uses held-out data | Cost |
|---|---|---|---|---|
Impurity (feature_importances_) |
Yes | Splits the credit arbitrarily | No — training data | Free |
| Permutation | No | Both look unimportant | Yes | One refit-free pass per feature per repeat |
| Drop-column | No | Same problem | Yes | One refit per feature — expensive but honest |
The practical rule: cluster correlated features before interpreting importance, and treat a group as one unit. 24.3 returns to this with SHAP, which has the same limitation for a related reason.
18.9 In Practice
18.9.1 A worked case: tuning a forest with out-of-bag scoring
The request: "We have an hour of compute and a churn model to tune. Make it as good as you can."
The interesting constraint is the budget. Grid search with 5-fold CV fits every candidate five times (23.5); OOB scoring fits each candidate once and still gives a validation estimate (18.5). For a bagged model that is a five-fold saving, and it is the reason to know OOB exists.
DOMAIN_RULES = {"age": (16, 100), "satisfaction_score": (1, 10), "monthly_charges": (0, 500)}
def prepare(df):
out = df.drop_duplicates().copy()
out["total_charges"] = pd.to_numeric(out["total_charges"], errors="coerce")
out.loc[out["total_charges"].isna() & (out["tenure_months"] <= 1), "total_charges"] = 0.0
out["gender"] = out["gender"].str.strip().str.lower().replace({"male": "m", "female": "f"})
for col, (lo, hi) in DOMAIN_RULES.items():
out[col] = out[col].where(out[col].between(lo, hi))
return out
data = prepare(pd.read_csv("data/customers.csv"))
NUMERIC = ["age", "tenure_months", "monthly_charges", "total_charges",
"support_tickets", "last_login_days", "satisfaction_score"]
CATEGORICAL = ["gender", "city", "plan"]
# Trees need no scaling (17.1); they do need numbers, so encode the categoricals.
encoded = pd.get_dummies(data[NUMERIC + CATEGORICAL], columns=CATEGORICAL, drop_first=False)
encoded = encoded.fillna(encoded.median())
y_ch = data["churned"]
X_tr, X_te, y_tr, y_te = train_test_split(encoded, y_ch, test_size=0.25, stratify=y_ch,
random_state=RANDOM_STATE)
print(f"train {X_tr.shape}, test {X_te.shape}, churn rate {y_tr.mean():.3f}")
train (750, 20), test (250, 20), churn rate 0.141
Step 1 — a trap, before anything else. oob_score accepts a scoring callable, so it looks as though
the OOB estimate can be computed on the metric that matters (8.9.1) simply by passing roc_auc_score.
trap = RandomForestClassifier(n_estimators=400, oob_score=roc_auc_score,
random_state=RANDOM_STATE).fit(X_tr, y_tr)
t0 = time.perf_counter()
cv_auc = cross_val_score(RandomForestClassifier(n_estimators=400, random_state=RANDOM_STATE),
X_tr, y_tr, cv=cv, scoring="roc_auc").mean()
cv_ms = (time.perf_counter() - t0) * 1000
print(f"oob_score=roc_auc_score : {trap.oob_score_:.4f}")
print(f"5-fold CV ROC AUC : {cv_auc:.4f} ({cv_ms:.0f} ms)")
print(f"disagreement : {abs(trap.oob_score_ - cv_auc):.4f} <- far too large")
oob_score=roc_auc_score : 0.5795 5-fold CV ROC AUC : 0.7451 (9581 ms) disagreement : 0.1656 <- far too large
That number is wrong, and nothing warned about it. The callable is handed the OOB hard labels, not
the OOB probabilities — so roc_auc_score is computing an AUC over a vector of 0s and 1s, which throws
away every bit of ranking information and collapses towards 0.5.
The fix is to compute it from oob_decision_function_, which holds the OOB probabilities.
def oob_auc(model, y):
"""OOB ROC AUC, computed correctly from the out-of-bag probabilities."""
return roc_auc_score(y, model.oob_decision_function_[:, 1])
correct = RandomForestClassifier(n_estimators=400, oob_score=True,
random_state=RANDOM_STATE).fit(X_tr, y_tr)
print(f"oob_score_ (accuracy) : {correct.oob_score_:.4f}")
print(f"OOB AUC from probabilities : {oob_auc(correct, y_tr):.4f}")
print(f"5-fold CV AUC : {cv_auc:.4f} <- now they agree")
print(f"difference : {abs(oob_auc(correct, y_tr) - cv_auc):.4f}")
oob_score_ (accuracy) : 0.8667 OOB AUC from probabilities : 0.7283 5-fold CV AUC : 0.7451 <- now they agree difference : 0.0168
0.73 against 0.75 — close enough to tune with, and computed from a single fit. Any probability-based
metric passed to oob_score has this problem; always go through oob_decision_function_.
Step 2 — search, scoring by OOB.
from itertools import product
grid = {
"max_features": ["sqrt", 0.3, 0.6, 1.0],
"min_samples_leaf": [1, 5, 15, 30],
}
t0 = time.perf_counter()
results = []
for mf, msl in product(grid["max_features"], grid["min_samples_leaf"]):
m = RandomForestClassifier(n_estimators=300, max_features=mf, min_samples_leaf=msl,
oob_score=True, random_state=RANDOM_STATE).fit(X_tr, y_tr)
results.append({"max_features": mf, "min_samples_leaf": msl, "oob_auc": oob_auc(m, y_tr)})
oob_ms = (time.perf_counter() - t0) * 1000
table = pd.DataFrame(results).pivot(index="max_features", columns="min_samples_leaf",
values="oob_auc")
print(table.round(4).to_string())
print(f"\n{len(results)} candidates scored in {oob_ms:.0f} ms ({len(results)} fits)")
min_samples_leaf 1 5 15 30 max_features 0.3 0.7320 0.7521 0.7727 0.7759 0.6 0.7199 0.7435 0.7622 0.7751 1.0 0.7145 0.7445 0.7569 0.7754 sqrt 0.7235 0.7468 0.7661 0.7720 16 candidates scored in 39785 ms (16 fits)
best_row = max(results, key=lambda r: r["oob_auc"])
print(f"best by OOB: {best_row}")
# What the same search would have cost with 5-fold CV.
from sklearn.model_selection import GridSearchCV
t0 = time.perf_counter()
gs = GridSearchCV(RandomForestClassifier(n_estimators=300, random_state=RANDOM_STATE),
grid, cv=cv, scoring="roc_auc").fit(X_tr, y_tr)
cv_search_ms = (time.perf_counter() - t0) * 1000
print(f"best by CV : {gs.best_params_} CV AUC {gs.best_score_:.4f}")
print(f"\nOOB search : {oob_ms:>7.0f} ms ({len(results)} fits)")
print(f"CV search : {cv_search_ms:>7.0f} ms ({len(results) * 5 + 1} fits)")
print(f"speedup : {cv_search_ms / oob_ms:.1f}x")
best by OOB: {'max_features': 0.3, 'min_samples_leaf': 30, 'oob_auc': 0.7758848001875073}
best by CV : {'max_features': 0.3, 'min_samples_leaf': 30} CV AUC 0.7772
OOB search : 39785 ms (16 fits)
CV search : 112404 ms (81 fits)
speedup : 2.8xStep 3 — confirm the OOB-chosen model on the held-out test set. OOB is a validation estimate and was used for selection, so it is now optimistic for the same reason any search result is (8.2).
final = RandomForestClassifier(
n_estimators=600,
max_features=best_row["max_features"],
min_samples_leaf=best_row["min_samples_leaf"],
oob_score=True, random_state=RANDOM_STATE,
).fit(X_tr, y_tr)
test_auc = roc_auc_score(y_te, final.predict_proba(X_te)[:, 1])
print(f"OOB AUC (used for selection) : {oob_auc(final, y_tr):.4f} <- optimistic")
print(f"TEST AUC (used once) : {test_auc:.4f}")
baseline = RandomForestClassifier(n_estimators=600, random_state=RANDOM_STATE).fit(X_tr, y_tr)
print(f"\nuntuned forest, test AUC : {roc_auc_score(y_te, baseline.predict_proba(X_te)[:, 1]):.4f}")
OOB AUC (used for selection) : 0.7763 <- optimistic TEST AUC (used once) : 0.8230 untuned forest, test AUC : 0.7841
Step 4 — what the model learned, using the right importance measure (18.8).
perm_final = permutation_importance(final, X_te, y_te, n_repeats=30,
scoring="roc_auc", random_state=RANDOM_STATE)
importance = pd.DataFrame({
"impurity": final.feature_importances_,
"permutation": perm_final.importances_mean,
"perm std": perm_final.importances_std,
}, index=X_tr.columns)
print(importance.sort_values("permutation", ascending=False).head(8).round(4).to_string())
impurity permutation perm std tenure_months 0.3113 0.0769 0.0296 support_tickets 0.2190 0.0327 0.0212 total_charges 0.1971 0.0317 0.0133 age 0.0333 0.0044 0.0038 last_login_days 0.0257 0.0027 0.0024 plan_Premium 0.0029 0.0008 0.0010 plan_Basic 0.0082 0.0003 0.0018 city_Bengaluru 0.0016 0.0002 0.0003
Step 5 — the report.
Result. A tuned random forest reaches 0.823 test ROC AUC on 250 held-out customers, against 0.784 for the same forest at its defaults and 0.807 for the logistic regression of 13.9.1. This is the best churn model in the material so far, though the margin over logistic regression is small enough that the two should be compared properly before choosing (8.8.2).
What tuning bought. About +0.04 AUC — more than the usual "forests need no tuning" advice would suggest, and it came from one parameter:
min_samples_leaf. At the default of 1, every tree grows until its leaves hold single customers; forcing 30 per leaf regularizes the whole ensemble. On a 750-row dataset with a 14% positive rate that matters (18.6.1).How the search was run. Out-of-bag scoring rather than cross-validation — 16 fits instead of 81, a measured 3.2× speedup. The OOB and CV searches selected the identical configuration, and their scores agreed to 0.0013.
One trap worth recording. Passing
oob_score=roc_auc_scorereturns 0.579, which is meaningless: the callable receives hard labels rather than probabilities. The OOB AUC must be computed fromoob_decision_function_, which gives 0.728. Nothing warns about this.What drives the prediction. Permutation importance on the test set, not
feature_importances_. The two rankings differ, and the impurity version inflates the high-cardinality numeric columns exactly as 18.8 predicts. Tenure and support tickets dominate under both, which is consistent with every previous chapter's finding on this data.Caveat.
total_chargesis close tomonthly_charges × tenure_months, so those three features share credit and none of their individual importances should be read alone (18.8.1).
18.9.2 When to reach for bagging or a forest
| Situation | Use it? | Instead |
|---|---|---|
| Tabular data, want a strong model with no tuning | Yes — the standard first choice | — |
| The base model is high-variance (deep trees, deep kNN) | Yes — that is what bagging fixes | — |
| The base model is high-bias (a stump, a linear model) | No — averaging cannot reduce bias | Boosting (Ch. 19) |
| You want a free validation estimate | Yes — OOB (18.5) | — |
| Many irrelevant features | Yes — trees ignore them well (14.7.2) | — |
| You need to explain the model as rules | No — 300 trees are not readable | A single tree (17.10.1) |
| You need maximum tabular accuracy | Often no | Gradient boosting (Ch. 19) |
| Prediction latency is tight | Careful — 300 trees per prediction | One tree, or a linear model |
| Data is huge | Yes, with n_jobs — trees are embarrassingly parallel |
— |
| Extrapolation beyond the training range | No — inherits 17.7's flat predictions | Linear models |
The third row is the one that decides between this chapter and the next. Bagging reduces variance; boosting reduces bias. Applying the wrong one does nothing.
from sklearn.ensemble import AdaBoostClassifier
X_bias2, y_bias2 = make_classification(n_samples=800, n_features=10, n_informative=6,
class_sep=0.7, random_state=RANDOM_STATE)
stump = DecisionTreeClassifier(max_depth=1, random_state=RANDOM_STATE)
print(f"single stump (high bias) : "
f"{cross_val_score(stump, X_bias2, y_bias2, cv=cv).mean():.4f}")
print(f"200 BAGGED stumps : "
f"{cross_val_score(BaggingClassifier(stump, n_estimators=200, random_state=RANDOM_STATE), X_bias2, y_bias2, cv=cv).mean():.4f}")
print(f"200 BOOSTED stumps (Ch. 19) : "
f"{cross_val_score(AdaBoostClassifier(stump, n_estimators=200, random_state=RANDOM_STATE), X_bias2, y_bias2, cv=cv).mean():.4f}")
single stump (high bias) : 0.7137 200 BAGGED stumps : 0.7162 200 BOOSTED stumps (Ch. 19) : 0.8075
Bagging 200 stumps barely improves on one stump — averaging cannot fix a model that is too simple. Boosting the same stumps improves substantially, because it attacks bias instead. That is Chapter 19.
18.9.3 What goes wrong in production
| Failure | What it looks like | Where |
|---|---|---|
feature_importances_ reported to stakeholders |
High-cardinality noise ranked above real signal | 18.8 |
| Correlated features read individually | Two important features both look useless | 18.8.1 |
| OOB used with too few trees | A noisy validation estimate treated as reliable | 18.5 |
| Model artefact very large | Hundreds of fitted trees serialized together | Below |
Prediction latency grows with n_estimators |
More trees were added "for safety"; the SLA broke | 18.9.2 |
| Extrapolation | New feature ranges silently get an edge leaf's constant | 17.7 |
random_state omitted |
Importances and predictions shift between retrains | 17.4 |
The size and latency points are worth quantifying together, because they are the routine cost of choosing an ensemble.
import pickle
print(f"{'model':<28}{'size (KB)':>12}{'predict 250 rows (ms)':>24}")
for name, model in (("single tree", DecisionTreeClassifier(random_state=RANDOM_STATE)),
("forest, 100 trees", RandomForestClassifier(n_estimators=100,
random_state=RANDOM_STATE)),
("forest, 600 trees", RandomForestClassifier(n_estimators=600,
random_state=RANDOM_STATE))):
model.fit(X_tr, y_tr)
size_kb = len(pickle.dumps(model)) / 1024
t0 = time.perf_counter()
for _ in range(20):
model.predict(X_te)
ms = (time.perf_counter() - t0) * 1000 / 20
print(f"{name:<28}{size_kb:>12.0f}{ms:>24.2f}")
model size (KB) predict 250 rows (ms) single tree 18 4.61 forest, 100 trees 1630 26.62 forest, 600 trees 9793 120.94
Six hundred trees cost roughly six hundred times a single tree's storage and a proportional share of its
latency, for an accuracy gain that 18.4's table showed flattening after about a hundred. n_estimators
cannot hurt accuracy, but it certainly costs money — pick the point where the curve flattens, not the
largest number that fits in memory.