Random Forests — Averaging Away the Variance

Two Sources of Randomness

A single deep tree is low-bias and high-variance — it pinned a client's default probability only to an SD of 0.26 across resamples. A random forest averages many deep trees to cancel that variance, but averaging only helps if the trees are de-correlated, so Breiman adds two independent sources of randomness: bagging (each tree on a bootstrap resample) and the random subspace (only $m$ of the $p$ features considered at each split). Going from one tree to a 100-tree forest lifts test AUC from 0.652 to 0.764.

Out-of-bag error — a validation set for free

Because each bootstrap omits about 37% of rows, every tree carries a free test set: averaging each row's prediction over only the trees that did not see it gives an honest generalisation estimate with no validation split. It works — OOB AUC 0.7612 against test AUC 0.7636 from scratch, 0.7725 against 0.7737 in R's ranger. One wrinkle worth knowing: a small forest cannot give every row an OOB estimate. At 5 trees about 10% of rows are in-bag everywhere and have none at all, so they must be excluded rather than scored — including them drags the small-forest end of the curve down by ~0.03 AUC and overstates how much OOB improves with tree count.

What de-correlation actually buys

The $\texttt{max\_features}$ knob is the interesting one, and the usual explanation is that restricting features de-correlates the trees. For an average of $B$ trees with variance $\sigma^2$ and pairwise correlation $\rho$, the ensemble variance is $\rho\sigma^2 + (1-\rho)\sigma^2/B$, so no number of trees can average away the $\rho\sigma^2$ floor — only de-correlation lowers it. Measured directly, ensemble variance falls 44% as features are restricted, from 0.0429 using all 23 down to 0.0241 at one, tracking the AUC gain exactly.

Two things that measurement makes visible and the usual telling does not. First, individual trees move the other way: per-tree variance rises from 0.061 to 0.111 as more features are allowed, so restricting features makes each tree noisier and the ensemble better — that is the trade the random subspace makes. Second, the obvious statistic fails to show any of this: the raw pairwise correlation between trees' predicted probabilities is nearly flat (0.375–0.392) across the whole sweep, because it is dominated by the signal every tree captures. The variance of the ensemble is the quantity the decomposition is written about, and the one to look at.

Feature importance, and why the default misleads

A practical warning the example makes concrete. Adding a pure-noise continuous feature to the data, the default impurity importance ranks it 3rd of 25 — a continuous variable offers many thresholds, so it accumulates spurious impurity decrease. Permutation importance puts it at −0.0012 AUC, i.e. zero, which is the honest answer. R reproduces this exactly, ranking the same noise feature 2nd by impurity and near zero by permutation. Prefer permutation importance, especially with continuous or high-cardinality features.

Extra-Trees, and regression

Extremely randomised trees push the idea further by drawing split thresholds at random rather than optimising them — more de-correlation, and a faster fit. Here it is a wash on accuracy (0.7753 against the forest's 0.7752 in Python; 0.7762 against 0.7737 in R) at about two-thirds the fit time. On the California housing regression the forest reaches an out-of-sample RMSE of 0.508 against a single tree's 0.666 and a linear regression's 0.737 — smooth from averaging, yet still bent, which is the whole point.

Notebooks

Downloads

Random-Forest Module — Source Code

"""
rf.py -- Random Forests from scratch (Breiman 2001).

Backs the notebooks in  Z-ML Trees-Random Forests.  Self-contained (a lean deep-tree grower with per-node
feature subsampling, so it does not depend on the CART notebook's module).

A single deep tree is a low-bias, HIGH-variance learner (see the CART notebook). A random forest averages
many deep trees to cancel that variance -- but averaging only helps if the trees are DE-CORRELATED, so
Breiman adds two independent sources of randomness:

    1. BAGGING  -- each tree is grown on a bootstrap resample of the rows.
    2. RANDOM SUBSPACE -- at every split only a random subset of  max_features  columns is considered.

For an average of B trees each with variance sigma^2 and pairwise correlation rho, the ensemble variance is
    rho sigma^2 + (1-rho)/B sigma^2,
so raising B kills the second term while LOWERING rho (via max_features) attacks the first -- the whole
point of the random subspace trick. Because each bootstrap leaves out ~37% of rows, every tree has an
"out-of-bag" (OOB) test set for free: averaging each row's predictions over only the trees that did NOT
see it gives an unbiased generalisation estimate with no separate validation split.
"""

import numpy as np


class _Node:
    __slots__ = ("feat", "thr", "left", "right", "value", "n", "imp")
    def __init__(self, value, n, imp):
        self.feat = None; self.thr = None; self.left = None; self.right = None
        self.value = value; self.n = n; self.imp = imp


class _Tree:
    """deep CART tree that considers a random subset of max_features features at each split."""
    def __init__(self, mf, min_samples_leaf, max_depth, rng, regression, n_classes):
        self.mf = mf; self.msl = min_samples_leaf; self.max_depth = max_depth
        self.rng = rng; self.regression = regression; self.K = n_classes
        self.importances = None

    def _imp(self, y):
        if self.regression: return np.mean((y - y.mean()) ** 2) if len(y) else 0.0
        c = np.bincount(y, minlength=self.K).astype(float); p = c / c.sum(); return 1.0 - (p * p).sum()

    def _leaf(self, y):
        if self.regression: return float(y.mean())
        return np.bincount(y, minlength=self.K).astype(float) / len(y)

    def _best_split(self, X, y):
        n, p = X.shape; parent = self._imp(y); best = (0.0, None, None)
        feats = self.rng.choice(p, min(self.mf, p), replace=False)
        for j in feats:
            xs = X[:, j]; order = np.argsort(xs, kind="quicksort"); xj = xs[order]; yj = y[order]
            cut = np.where(np.diff(xj) > 0)[0]
            if len(cut) == 0: continue
            if self.regression:
                cs = np.cumsum(yj); cs2 = np.cumsum(yj * yj); tot = cs[-1]; tot2 = cs2[-1]
                nL = cut + 1.0; nR = n - nL
                child = ((cs2[cut] - cs[cut] ** 2 / nL) + ((tot2 - cs2[cut]) - (tot - cs[cut]) ** 2 / nR)) / n
            else:
                oh = np.zeros((n, self.K)); oh[np.arange(n), yj] = 1.0; cc = np.cumsum(oh, axis=0)
                nL = (cut + 1.0)[:, None]; nR = n - nL; cL = cc[cut]; cR = cc[-1] - cL
                giniL = 1.0 - ((cL / nL) ** 2).sum(1); giniR = 1.0 - ((cR / nR) ** 2).sum(1)
                child = (nL[:, 0] * giniL + nR[:, 0] * giniR) / n
            gain = parent - child
            ok = (np.asarray(nL).reshape(-1) >= self.msl) & (np.asarray(nR).reshape(-1) >= self.msl)
            if not ok.any(): continue
            gain = np.where(ok, gain, -np.inf); k = int(np.argmax(gain))
            if gain[k] > best[0]:
                best = (float(gain[k]), int(j), float((xj[cut[k]] + xj[cut[k] + 1]) / 2.0))
        return best

    def _build(self, X, y, depth):
        node = _Node(self._leaf(y), len(y), self._imp(y))
        if len(y) < 2 * self.msl or node.imp == 0.0 or (self.max_depth is not None and depth >= self.max_depth):
            return node
        gain, feat, thr = self._best_split(X, y)
        if feat is None or gain <= 0: return node
        m = X[:, feat] <= thr; node.feat = feat; node.thr = thr
        self.importances[feat] += node.imp * node.n - self._imp(y[m]) * m.sum() - self._imp(y[~m]) * (~m).sum()
        node.left = self._build(X[m], y[m], depth + 1); node.right = self._build(X[~m], y[~m], depth + 1)
        return node

    def fit(self, X, y):
        self.importances = np.zeros(X.shape[1]); self.root = self._build(X, y, 0); return self

    def predict_proba(self, X):
        out = []
        for x in X:
            nd = self.root
            while nd.feat is not None: nd = nd.left if x[nd.feat] <= nd.thr else nd.right
            out.append(nd.value)
        return np.array(out)


class RandomForest:
    def __init__(self, n_estimators=100, max_features="sqrt", min_samples_leaf=1,
                 max_depth=None, bootstrap=True, regression=False, random_state=0):
        self.n_estimators = n_estimators; self.max_features = max_features
        self.min_samples_leaf = min_samples_leaf; self.max_depth = max_depth
        self.bootstrap = bootstrap; self.regression = regression; self.random_state = random_state

    def _mf(self, p):
        if self.max_features == "sqrt": return max(1, int(np.sqrt(p)))
        if self.max_features == "log2": return max(1, int(np.log2(p)))
        if isinstance(self.max_features, float): return max(1, int(self.max_features * p))
        if self.max_features is None: return p
        return int(self.max_features)

    def fit(self, X, y):
        X = np.asarray(X, float); y = np.asarray(y); n, p = X.shape
        if not self.regression:
            self.classes_, y = np.unique(y, return_inverse=True); self.K = len(self.classes_)
        else:
            self.K = 1
        mf = self._mf(p); rng = np.random.default_rng(self.random_state)
        self.trees = []; imp = np.zeros(p)
        oob_sum = np.zeros((n, self.K)) if not self.regression else np.zeros(n)
        oob_cnt = np.zeros(n)
        for b in range(self.n_estimators):
            idx = rng.integers(0, n, n) if self.bootstrap else np.arange(n)
            tr = _Tree(mf, self.min_samples_leaf, self.max_depth,
                       np.random.default_rng(rng.integers(1 << 31)), self.regression, self.K).fit(X[idx], y[idx])
            self.trees.append(tr); imp += tr.importances
            if self.bootstrap:
                oob = np.ones(n, bool); oob[np.unique(idx)] = False
                if oob.any():
                    pr = tr.predict_proba(X[oob])
                    if self.regression: oob_sum[oob] += pr.reshape(-1)
                    else: oob_sum[oob] += pr
                    oob_cnt[oob] += 1
        s = imp.sum(); self.feature_importances_ = imp / s if s > 0 else imp
        self.n_features_ = p
        if self.bootstrap:
            seen = oob_cnt > 0
            if self.regression:
                self.oob_prediction_ = np.where(seen, oob_sum / np.maximum(oob_cnt, 1), np.nan)
            else:
                self.oob_decision_function_ = oob_sum / np.maximum(oob_cnt, 1)[:, None]
                self.oob_seen_ = seen
        return self

    def predict_proba(self, X):
        X = np.asarray(X, float)
        return np.mean([t.predict_proba(X) for t in self.trees], axis=0)

    def predict(self, X):
        if self.regression:
            X = np.asarray(X, float)
            return np.mean([t.predict_proba(X).reshape(-1) for t in self.trees], axis=0)
        return self.classes_[np.argmax(self.predict_proba(X), axis=1)]

References

• Breiman, L. (2001). Random forests. Machine Learning 45(1), 5–32. — bagging plus the random subspace, and the variance decomposition
• Breiman, L. (1996). Bagging predictors. Machine Learning 24(2), 123–140. — and out-of-bag estimation
• Geurts, P., Ernst, D. & Wehenkel, L. (2006). Extremely randomized trees. Machine Learning 63(1), 3–42. — random thresholds in place of optimised ones
• Strobl, C., Boulesteix, A.-L., Zeileis, A. & Hothorn, T. (2007). Bias in random forest variable importance measures. BMC Bioinformatics 8, 25. — why impurity importance favours continuous and high-cardinality features
• Wright, M. N. & Ziegler, A. (2017). ranger: a fast implementation of random forests. Journal of Statistical Software 77(1). — the R implementation used here