Random Forests — averaging away the variance¶
Breiman (2001), Random Forests¶
The CART notebook ended on a problem: a single deep tree is low-bias but high-variance — refit on a bootstrap resample and a client's predicted default probability swung by an SD of ~0.26, and averaging 25 such trees lifted test AUC from 0.65 to 0.76. The random forest is the disciplined version of that idea, and it rests on one insight about averaging.
Averaging $B$ trees each of variance $\sigma^2$ with pairwise correlation $\rho$ gives an ensemble variance $$\rho\,\sigma^2 + \frac{1-\rho}{B}\,\sigma^2.$$ More trees ($B\uparrow$) kills the second term; but the first term floors out at $\rho\sigma^2$. So the trees must be de-correlated. Breiman injects two independent sources of randomness:
- Bagging — grow each tree on a bootstrap resample of the rows.
- Random subspace — at every split, consider only a random subset of
max_featurescolumns.
The random subspace is what lowers $\rho$: it stops every tree from locking onto the same dominant feature (here PAY_1). A free bonus falls out — each bootstrap omits ~37% of rows, so every row has an out-of-bag test set built in. We build all of this from scratch (rf.py), verify it against scikit-learn, expose the well-known feature-importance bias, and finish with Extremely Randomized Trees. Data: Taiwan credit-card default, California housing for regression.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, time
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier, RandomForestRegressor
from sklearn.tree import DecisionTreeClassifier
from sklearn.inspection import permutation_importance
from sklearn.metrics import roc_auc_score, mean_squared_error
import rf
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"
d = pd.read_csv("credit_default.csv"); feat=[c for c in d.columns if c!="default"]
X=d[feat].to_numpy(float); y=d["default"].to_numpy(int)
Xtr,Xte,ytr,yte = train_test_split(X,y,test_size=0.3,random_state=0,stratify=y)
# a subsample keeps the from-scratch forest fast; scikit-learn runs on the full training set
Xs,ys = Xtr[:7000], ytr[:7000]
print(f"credit default: train {Xtr.shape[0]:,} (from-scratch demo uses {Xs.shape[0]:,}), test {Xte.shape[0]:,}; {X.shape[1]} features")
credit default: train 21,000 (from-scratch demo uses 7,000), test 9,000; 23 features
1. From one tree to a forest¶
Our from-scratch RandomForest bags the rows and subsamples $\sqrt{p}$ features per split. Against a single deep tree on the same data it should show the ensemble's signature: a large jump in test AUC from averaging de-correlated trees. (ROC-AUC = the probability the model ranks a random defaulter above a random non-defaulter; 0.5 = chance, 1 = perfect; threshold-free and robust to class imbalance — defined in the CART notebook.)
one = DecisionTreeClassifier(min_samples_leaf=5, random_state=0).fit(Xs,ys)
t=time.time(); forest = rf.RandomForest(n_estimators=100, max_features="sqrt", min_samples_leaf=5, random_state=0).fit(Xs,ys); el=time.time()-t
auc1 = roc_auc_score(yte, one.predict_proba(Xte)[:,1]); aucF = roc_auc_score(yte, forest.predict_proba(Xte)[:,1])
print(f"single deep tree : test AUC {auc1:.4f}")
print(f"random forest (100 trees, from scratch, {el:.1f}s): test AUC {aucF:.4f} (+{aucF-auc1:.3f})")
print("Same trees, but de-correlated and averaged -- the variance cancels and AUC jumps.")
single deep tree : test AUC 0.6520 random forest (100 trees, from scratch, 9.0s): test AUC 0.7636 (+0.112) Same trees, but de-correlated and averaged -- the variance cancels and AUC jumps.
2. Out-of-bag error — a validation set for free¶
Because each tree omits ~37% of the rows, we can score every row using only the trees that never saw it. That OOB estimate needs no held-out split, and it tracks the true test error closely. We check OOB-AUC against test-AUC, and watch both converge as trees are added (using scikit-learn's warm_start to grow the forest incrementally).
oob_auc = roc_auc_score(ys[forest.oob_seen_], forest.oob_decision_function_[forest.oob_seen_,1])
print(f"from-scratch forest: OOB AUC {oob_auc:.4f} vs test AUC {aucF:.4f} -> OOB is an honest free estimate")
# convergence of OOB and test AUC with number of trees (scikit-learn, warm_start)
import warnings
ns=[5,10,20,40,80,150,300]; oob=[]; te=[]; cov=[]
rfw=RandomForestClassifier(warm_start=True, oob_score=True, max_features="sqrt", min_samples_leaf=5, random_state=0, n_jobs=-1)
for nt in ns:
with warnings.catch_warnings(): # the coverage this warns about is reported below
warnings.simplefilter("ignore", UserWarning)
rfw.set_params(n_estimators=nt).fit(Xtr,ytr)
D = rfw.oob_decision_function_
# a small forest leaves some rows out-of-bag for NO tree; they have no OOB estimate and
# must be excluded rather than scored as if they did.
ok = ~np.isnan(D).any(1) & (D.sum(1) > 0); cov.append(ok.mean())
oob.append(roc_auc_score(ytr[ok], D[ok,1])); te.append(roc_auc_score(yte, rfw.predict_proba(Xte)[:,1]))
fig,ax=plt.subplots(figsize=(7.5,4))
ax.plot(ns,oob,"o-",color=ORANGE,label="OOB AUC"); ax.plot(ns,te,"s-",color=BLUE,label="test AUC")
ax.set_xlabel("number of trees"); ax.set_ylabel("ROC-AUC"); ax.set_title("OOB tracks test error, and both converge with more trees")
ax.legend(frameon=False); plt.tight_layout(); plt.show()
print("OOB AUC sits right on the test AUC curve: adding trees never hurts (unlike depth), it only reduces variance until it plateaus.")
print(f"OOB coverage: " + ", ".join(f"{n} trees {c:.1%}" for n, c in zip(ns, cov)))
print("A small forest cannot give every row an OOB estimate -- with 5 trees about a tenth of rows are")
print("in-bag everywhere -- so those rows are excluded rather than scored. Including them would drag the")
print("left end of the curve down by ~0.03 AUC and overstate how much OOB improves with tree count.")
from-scratch forest: OOB AUC 0.7612 vs test AUC 0.7636 -> OOB is an honest free estimate
OOB AUC sits right on the test AUC curve: adding trees never hurts (unlike depth), it only reduces variance until it plateaus. OOB coverage: 5 trees 89.9%, 10 trees 98.9%, 20 trees 100.0%, 40 trees 100.0%, 80 trees 100.0%, 150 trees 100.0%, 300 trees 100.0% A small forest cannot give every row an OOB estimate -- with 5 trees about a tenth of rows are in-bag everywhere -- so those rows are excluded rather than scored. Including them would drag the left end of the curve down by ~0.03 AUC and overstate how much OOB improves with tree count.
3. The de-correlation knob — max_features¶
max_features is the dial between strength and de-correlation. Set it to $p$ (all features) and every tree greedily grabs PAY_1 first — strong trees, but highly correlated, so averaging helps little (this is plain bagging). Shrink it and trees are forced to explore other features — weaker individually, but far less correlated. The test-AUC curve typically peaks at an intermediate value; $\sqrt{p}$ is Breiman's default.
p=X.shape[1]; mfs=[1,2,3,4,int(np.sqrt(p)),8,12,16,p]; mfs=sorted(set(mfs)); aucm=[]
for mf in mfs:
m=RandomForestClassifier(n_estimators=200,max_features=mf,min_samples_leaf=5,random_state=0,n_jobs=-1).fit(Xtr,ytr)
aucm.append(roc_auc_score(yte,m.predict_proba(Xte)[:,1]))
fig,ax=plt.subplots(figsize=(7.5,4))
ax.plot(mfs,aucm,"o-",color=GREEN); ax.axvline(int(np.sqrt(p)),color=RED,ls="--",label=f"sqrt(p)={int(np.sqrt(p))} (default)")
ax.axvline(p,color="grey",ls=":",label=f"p={p} (= bagging)")
ax.set_xlabel("max_features per split"); ax.set_ylabel("test AUC"); ax.set_title("Too many features -> correlated trees; the sweet spot is intermediate")
ax.legend(frameon=False); plt.tight_layout(); plt.show()
print(f"Best max_features = {mfs[int(np.argmax(aucm))]} (AUC {max(aucm):.4f}); using all {p} features (pure bagging) is worse.")
# "de-correlation matters" is so far an inference from AUC. Measure what the variance formula
# is actually about -- the variance of the ENSEMBLE prediction -- as max_features varies.
def spread(mf, B=60):
f = RandomForestClassifier(n_estimators=B, max_features=mf, min_samples_leaf=5,
random_state=0, n_jobs=-1).fit(Xtr, ytr)
Pm = np.column_stack([t.predict_proba(Xte)[:, 1] for t in f.estimators_])
R = np.corrcoef(Pm, rowvar=False)
return (float(R[np.triu_indices_from(R, 1)].mean()), # pairwise correlation between trees
float(Pm.var(0).mean()), # mean variance of a single tree
float(Pm.mean(1).var())) # variance of the ensemble
stats = [spread(mf) for mf in mfs]
print(f"{'max_features':>13s}{'tree corr':>11s}{'per-tree var':>14s}{'ensemble var':>14s}{'test AUC':>10s}")
for mf, (r, v1, ve), a in zip(mfs, stats, aucm):
print(f"{mf:13d}{r:11.3f}{v1:14.4f}{ve:14.4f}{a:10.4f}")
ve_lo, ve_hi = stats[0][2], stats[-1][2]
print()
print(f"The ensemble variance is what the formula predicts, and it behaves: {ve_hi:.4f} using all {p} features")
print(f"(pure bagging) down to {ve_lo:.4f} at max_features={mfs[0]} -- a {100*(1-ve_lo/ve_hi):.0f}% reduction, tracking the AUC gain.")
print("Note the individual trees move the OTHER way: restricting features makes each tree noisier, and the")
print("ensemble still wins. That is the trade the random subspace makes.")
print()
print("One caution, since it is the obvious thing to reach for: the raw pairwise CORRELATION between trees'")
print("predicted probabilities is nearly flat across the sweep. It is dominated by the signal every tree")
print("captures -- they all track PAY_1 -- so it is a poor probe of the randomisation. The variance of the")
print("ensemble is the quantity the rho*sigma^2 + (1-rho)*sigma^2/B decomposition is written about, and the")
print("one to look at.")
Best max_features = 3 (AUC 0.7761); using all 23 features (pure bagging) is worse.
max_features tree corr per-tree var ensemble var test AUC
1 0.383 0.0613 0.0241 0.7761
2 0.386 0.0795 0.0316 0.7753
3 0.392 0.0884 0.0355 0.7761
4 0.387 0.0933 0.0371 0.7747
8 0.390 0.1026 0.0411 0.7733
12 0.385 0.1062 0.0420 0.7701
16 0.381 0.1088 0.0426 0.7687
23 0.375 0.1111 0.0429 0.7657
The ensemble variance is what the formula predicts, and it behaves: 0.0429 using all 23 features
(pure bagging) down to 0.0241 at max_features=1 -- a 44% reduction, tracking the AUC gain.
Note the individual trees move the OTHER way: restricting features makes each tree noisier, and the
ensemble still wins. That is the trade the random subspace makes.
One caution, since it is the obvious thing to reach for: the raw pairwise CORRELATION between trees'
predicted probabilities is nearly flat across the sweep. It is dominated by the signal every tree
captures -- they all track PAY_1 -- so it is a poor probe of the randomisation. The variance of the
ensemble is the quantity the rho*sigma^2 + (1-rho)*sigma^2/B decomposition is written about, and the
one to look at.
4. From-scratch vs scikit-learn, and full-scale¶
The from-scratch forest should match RandomForestClassifier (same bagging + subspace rule). We confirm it on the subsample, then let scikit-learn run on the full training set for the production number.
sk = RandomForestClassifier(n_estimators=100,max_features="sqrt",min_samples_leaf=5,random_state=0,n_jobs=-1).fit(Xs,ys)
print(f"on {Xs.shape[0]:,} rows: from-scratch AUC {aucF:.4f} scikit-learn AUC {roc_auc_score(yte,sk.predict_proba(Xte)[:,1]):.4f}")
skf = RandomForestClassifier(n_estimators=400,max_features="sqrt",min_samples_leaf=5,random_state=0,n_jobs=-1).fit(Xtr,ytr)
print(f"full training set (400 trees): scikit-learn test AUC {roc_auc_score(yte,skf.predict_proba(Xte)[:,1]):.4f}")
print(f"reminder -- a single pruned tree (CART notebook) reached ~0.74; the forest clears it comfortably.")
on 7,000 rows: from-scratch AUC 0.7636 scikit-learn AUC 0.7668
full training set (400 trees): scikit-learn test AUC 0.7752 reminder -- a single pruned tree (CART notebook) reached ~0.74; the forest clears it comfortably.
5. Feature importance — and why the default is biased¶
Random forests give two importance measures, and the difference matters. Impurity (Gini) importance sums each feature's split gains — but it is biased toward high-cardinality / continuous features, which offer more thresholds to overfit, even when they carry no signal. Permutation importance instead shuffles a feature on held-out data and measures the drop in performance — model-agnostic and far more trustworthy. We prove the bias by injecting two pure-noise features (one continuous, one binary) and watching only the impurity measure light up the continuous one.
rng=np.random.default_rng(0)
Xn_tr=np.column_stack([Xtr, rng.normal(size=len(Xtr)), rng.integers(0,2,len(Xtr))])
Xn_te=np.column_stack([Xte, rng.normal(size=len(Xte)), rng.integers(0,2,len(Xte))])
fn=feat+["NOISE_cont","NOISE_binary"]
m=RandomForestClassifier(n_estimators=300,max_features="sqrt",min_samples_leaf=5,random_state=0,n_jobs=-1).fit(Xn_tr,ytr)
imp_gini=pd.Series(m.feature_importances_,index=fn)
perm=permutation_importance(m,Xn_te,yte,n_repeats=10,random_state=0,scoring="roc_auc",n_jobs=-1)
imp_perm=pd.Series(perm.importances_mean,index=fn)
fig,ax=plt.subplots(1,2,figsize=(14,5))
imp_gini.sort_values().tail(10).plot.barh(ax=ax[0],color=[RED if "NOISE" in k else BLUE for k in imp_gini.sort_values().tail(10).index])
ax[0].set_title("Impurity (Gini) importance — inflates NOISE_cont"); ax[0].set_xlabel("importance")
imp_perm.sort_values().tail(10).plot.barh(ax=ax[1],color=[RED if "NOISE" in k else GREEN for k in imp_perm.sort_values().tail(10).index])
ax[1].set_title("Permutation importance — noise correctly ~0"); ax[1].set_xlabel("AUC drop when shuffled")
plt.tight_layout(); plt.show()
print(f"Impurity importance ranks the pure-noise continuous feature #{list(imp_gini.sort_values(ascending=False).index).index('NOISE_cont')+1} of {len(fn)} -- a spurious signal from its many thresholds.")
print(f"Permutation importance puts it at {imp_perm['NOISE_cont']:+.4f} AUC (i.e. ~0): the honest measure. Prefer permutation importance, especially with continuous features.")
Impurity importance ranks the pure-noise continuous feature #3 of 25 -- a spurious signal from its many thresholds. Permutation importance puts it at -0.0012 AUC (i.e. ~0): the honest measure. Prefer permutation importance, especially with continuous features.
6. Extremely Randomized Trees (ExtraTrees)¶
Geurts, Ernst & Wehenkel (2006) push the randomness further: instead of searching for the best threshold on each candidate feature, pick the threshold at random. This de-correlates trees even more (lowering $\rho$ further) and is faster — sometimes trading a little bias for a useful variance cut. We benchmark it against the random forest.
rf400=RandomForestClassifier(n_estimators=400,max_features="sqrt",min_samples_leaf=5,random_state=0,n_jobs=-1)
et400=ExtraTreesClassifier(n_estimators=400,max_features="sqrt",min_samples_leaf=5,random_state=0,n_jobs=-1)
for name,mdl in [("Random Forest",rf400),("Extra-Trees",et400)]:
t=time.time(); mdl.fit(Xtr,ytr); ft=time.time()-t
print(f"{name:15s}: test AUC {roc_auc_score(yte,mdl.predict_proba(Xte)[:,1]):.4f} (fit {ft:.1f}s)")
print("Extra-Trees uses random thresholds -> more de-correlation and a faster fit; here it is competitive with the RF.")
Random Forest : test AUC 0.7752 (fit 1.2s)
Extra-Trees : test AUC 0.7753 (fit 0.7s) Extra-Trees uses random thresholds -> more de-correlation and a faster fit; here it is competitive with the RF.
7. How the forest improves — graphically¶
Numbers aside, it helps to see what the forest buys over a single tree and over the parametric baseline the CART notebook used. Two one-feature views over continuous predictors (a tree can only step, never curve, on a discrete feature like PAY_1, so we use continuous ones here): classification — $P(\text{default})$ against the credit limit, with the empirical proportions (quantile bins), logistic regression, one tree, and the forest; and regression — California median value against median income, with the linear fit, one tree, and the forest. In both, averaging hundreds of bootstrapped trees turns the single tree's coarse steps into a smooth curve, while still bending where the straight-line logit/linear baseline cannot.
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.tree import DecisionTreeRegressor
GREY="#a0aec0"
fig,ax=plt.subplots(1,2,figsize=(14,4.8))
# classification: P(default) vs CREDIT LIMIT (continuous) -- quantile-binned empirical rate
lim=Xtr[:,feat.index("LIMIT_BAL")]/1000.0 # NT$ thousands
qb=np.quantile(lim,np.linspace(0,1,16)); bidx=np.clip(np.digitize(lim,qb)-1,0,14)
ctr=np.array([lim[bidx==k].mean() for k in range(15)]); emp=np.array([ytr[bidx==k].mean() for k in range(15)])
g=np.linspace(lim.min(),np.percentile(lim,99),300)
lo_=LogisticRegression(max_iter=1000).fit(lim[:,None],ytr).predict_proba(g[:,None])[:,1]
t1=DecisionTreeClassifier(max_depth=4,random_state=0).fit(lim[:,None],ytr).predict_proba(g[:,None])[:,1]
f1=RandomForestClassifier(n_estimators=400,min_samples_leaf=200,random_state=0,n_jobs=-1).fit(lim[:,None],ytr).predict_proba(g[:,None])[:,1]
ax[0].scatter(ctr,emp,s=45,color="k",zorder=3,label="empirical (quantile bins)")
ax[0].plot(g,lo_,color=ORANGE,lw=2,label="logistic"); ax[0].plot(g,t1,color=GREY,lw=2,label="single tree"); ax[0].plot(g,f1,color=GREEN,lw=2,label="random forest")
ax[0].set_xlabel("credit limit (NT$ thousands)"); ax[0].set_ylabel("P(default)"); ax[0].set_title("Classification: P(default) vs credit limit"); ax[0].legend(frameon=False,fontsize=8)
# regression: California value vs median income
h=pd.read_csv("cali_housing.csv"); hf=[c for c in h.columns if c!="MedHouseVal"]
Xh=h[hf].to_numpy(float); yh=h["MedHouseVal"].to_numpy(float)
Xhtr,Xhte,yhtr,yhte=train_test_split(Xh,yh,test_size=0.3,random_state=0)
xi=Xhtr[:,hf.index("MedInc")]; gh=np.linspace(xi.min(),np.percentile(xi,99),300)
ols=LinearRegression().fit(xi[:,None],yhtr).predict(gh[:,None])
trg=DecisionTreeRegressor(max_depth=4,random_state=0).fit(xi[:,None],yhtr).predict(gh[:,None])
frg=RandomForestRegressor(n_estimators=400,min_samples_leaf=40,random_state=0,n_jobs=-1).fit(xi[:,None],yhtr).predict(gh[:,None])
ax[1].scatter(xi,yhtr,s=4,alpha=0.07,color="grey")
ax[1].plot(gh,ols,color=ORANGE,lw=2,label="linear regression"); ax[1].plot(gh,trg,color=GREY,lw=2,label="single tree"); ax[1].plot(gh,frg,color=GREEN,lw=2,label="random forest")
ax[1].set_xlabel("median income"); ax[1].set_ylabel("median house value ($100k)"); ax[1].set_title("Regression: value vs income"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
# full multivariate out-of-sample error for the record
ols_f=LinearRegression().fit(Xhtr,yhtr); tr_f=DecisionTreeRegressor(max_depth=8,random_state=0).fit(Xhtr,yhtr)
rf_f=RandomForestRegressor(n_estimators=300,max_features=1/3.,min_samples_leaf=3,random_state=0,n_jobs=-1).fit(Xhtr,yhtr)
print(f"California OOS RMSE ($100k): linear regression {mean_squared_error(yhte,ols_f.predict(Xhte))**.5:.3f} | single tree {mean_squared_error(yhte,tr_f.predict(Xhte))**.5:.3f} | random forest {mean_squared_error(yhte,rf_f.predict(Xhte))**.5:.3f}")
print("The forest curve is smooth (variance averaged) AND bent (nonlinearity kept) -- and it has the lowest out-of-sample error,")
print("beating both the single tree and the straight-line/logit baseline on classification and regression alike.")
California OOS RMSE ($100k): linear regression 0.737 | single tree 0.666 | random forest 0.508 The forest curve is smooth (variance averaged) AND bent (nonlinearity kept) -- and it has the lowest out-of-sample error, beating both the single tree and the straight-line/logit baseline on classification and regression alike.
8. Summary¶
The random forest turns the single tree's fatal flaw — variance — into a strength by averaging many de-correlated deep trees. We built it from scratch (bagging + random subspace + OOB) and matched scikit-learn: on credit default it cleared the single pruned tree's AUC comfortably, the out-of-bag estimate tracked the test error for free, and max_features traced the strength-vs-de-correlation trade-off with an intermediate optimum. We also showed the trap in the default impurity importance (it inflates high-cardinality noise) and the fix (permutation importance), and that Extra-Trees buys extra de-correlation with random thresholds.
The forest reduces error by cutting variance — it averages independent trees, each already low-bias. The next notebook takes the opposite route: boosting adds small, dependent trees in sequence, each correcting the last, to cut bias (Friedman 2001) — and then the production libraries XGBoost / LightGBM / CatBoost. The R companion fits this same forest with ranger (fast RF with permutation importance built in). The Bayesian cousin is BART (its own notebook in this subsection), which regularises many small trees with priors and returns a full posterior — credible intervals on every prediction, which the forest cannot give.