4.4.1Random Forests

Asking many people for an estimate and averaging their answers usually gives a better result than trusting any one of them, because their individual errors partly cancel out. A random forest applies the same idea to decision trees: it grows many different trees and combines their predictions.

In practice, random forests are among the most dependable models for tabular data. They need little preparation, work with features on any scale, capture interactions between features, and perform reasonably well with their default settings.

That makes a random forest a good first strong model after a simple baseline, especially when contest time is short. Later, it is also a useful point of comparison for gradient boosting, which can be more accurate but needs more care.

A fully grown decision tree fits its training data closely, but it has high variance: a small change in the training data can change its splits and many of its predictions. A random forest reduces this variance by averaging many trees, each grown on a random resample of the training data with a random choice of candidate features at every split. This module explains why averaging works, how the trees are made to differ, and the out-of-bag estimates and feature importances that forests provide.

Ensembles and averaging

Averaging helps because the errors of different members partly cancel. Let an ensemble have $B$ members, and let $e_b$ be the error of member $b$ on a random sample: its prediction minus the true label. Assume that every error has mean 0 and variance $\sigma^2 = \mathbb{E}[e_b^2]$, where $\mathbb{E}$ denotes the average over samples, and that every pair of errors has correlation $\rho = \mathbb{E}[e_b e_{b'}] / \sigma^2$ for $b \ne b'$. The correlation lies between −1 and 1; it is 1 when two errors are always equal and 0 when they have no linear relationship.

The ensemble's error is the average $\bar{e} = \frac{1}{B} \sum_{b=1}^{B} e_b$, and its variance is

$$\mathbb{E}[\bar{e}^2] = \rho \sigma^2 + \frac{1 - \rho}{B} \sigma^2.$$

The table evaluates this formula for members with $\sigma^2 = 1$.

Independent errors ($\rho = 0$) give a variance proportional to $1/B$. For correlated errors, the second term vanishes as $B$ grows but $\rho \sigma^2$ remains: with $\rho = 0.3$, ten members reduce the variance from 1 to 0.37, and no number of members reduces it below 0.3. The correlation therefore limits what averaging can achieve, and a random forest is designed to lower $\rho$ without making each tree much less accurate.

Voting behaves in the same way. If three classifiers are each correct with probability 0.7 and err independently, the majority is correct with probability $0.7^3 + 3 \times 0.7^2 \times 0.3 = 0.784$, and 25 such classifiers raise it to 0.983. Classifiers trained on the same data make correlated errors, so the real gain is smaller.

Bootstrap samples and bagging

Trees grown by CART on the same data differ at most where splits tie, so averaging them requires different training data for each tree. The bootstrap creates such data from one training set.

The code draws a bootstrap sample from 10 samples, numbered 0 to 9, and counts how often each one was drawn.

import numpy as np

rng = np.random.default_rng(0)
n = 10
idx = rng.integers(0, n, size=n)        # n draws with replacement
counts = np.bincount(idx, minlength=n)  # times each sample was drawn

Sample 0 was drawn three times and sample 8 twice, so the 10 draws contain only 7 distinct samples. Samples 4, 7 and 9 were never drawn: they are out-of-bag.

Each draw misses a given sample with probability $1 - 1/n$. The $n$ draws are independent, so all of them miss it with probability $(1 - 1/n)^n$, and the expected fraction of distinct samples is $1 - (1 - 1/n)^n$.

The limit follows from $(1 + x/n)^n \to e^x$ as $n \to \infty$, with $x = -1$. For a training set of a few hundred samples or more, a bootstrap sample therefore contains about 63.2% of the distinct samples, and 36.8% are out-of-bag; a simulation agrees:

n = 1000
fractions = [len(np.unique(rng.integers(0, n, size=n))) / n
             for _ in range(10_000)]
print(round(np.mean(fractions), 4))     # 0.6324

The code compares one unrestricted tree with 100 bagged trees on the breast cancer training split of Decision Trees, scoring each by 5-fold cross-validation: the mean accuracy on five validation folds of the training set.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import BaggingClassifier
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.tree import DecisionTreeClassifier

cancer = load_breast_cancer()
X_train, X_val, y_train, y_val = train_test_split(
    cancer.data, cancer.target, test_size=0.3,
    stratify=cancer.target, random_state=0,
)
tree = DecisionTreeClassifier(random_state=0)
bagging = BaggingClassifier(DecisionTreeClassifier(),
                            n_estimators=100, random_state=0)
for name, model in [("tree", tree), ("bagging", bagging)]:
    scores = cross_val_score(model, X_train, y_train, cv=5)
    print(name, scores.round(3), round(scores.mean(), 3))

tree [0.95  0.888 0.912 0.899 0.975] 0.925
bagging [0.962 0.95  0.962 0.924 0.975] 0.955

Bagging raises the mean accuracy from 0.925 to 0.955, and the worst fold from 0.888 to 0.924.

Random feature subsets

Bagged trees remain strongly correlated, because the strongest features win the upper splits in most bootstrap samples. When the bagging model above is fitted to the whole training split, worst perimeter is the root feature of 50 of its 100 trees, and only five different features appear at a root.

The number $m$ is the hyperparameter max_features, introduced in Decision Trees, and $m = d$ gives bagging. By default the scikit-learn classifier uses $m = \sqrt{d}$ rounded down, 5 of the 30 breast cancer features. Its 100 trees then have 14 different root features, and worst perimeter is the root of only 19.

function random_forest(X, y, n_trees, m):
    n = number of samples in X
    trees = []
    repeat n_trees times:
        rows = n random row indices, drawn with replacement
        tree = full CART tree grown on X[rows], y[rows],
               trying only m random features at each split
        trees.append(tree)
    return trees

function forest_predict(trees, x):
    if regression:
        return mean of tree.predict(x) over all trees
    else:  # classification
        return the class with the highest mean probability

Measuring the correlation

The formula can be tested on the validation split. For a sample with label $y \in \{0, 1\}$, tree $b$ predicts a probability $p_b$ of class 1 (benign) and has error $e_b = p_b - y$. The code estimates $\sigma^2$ and $\rho$ from these errors for three values of max_features and compares the formula with the forest's actual squared error.

from sklearn.ensemble import RandomForestClassifier

B = 100
for max_features in [None, "sqrt", 1]:
    forest = RandomForestClassifier(n_estimators=B,
                                    max_features=max_features,
                                    random_state=0)
    forest.fit(X_train, y_train)
    # Row b holds the errors of tree b on the validation samples
    E = np.array([t.predict_proba(X_val)[:, 1] - y_val
                  for t in forest.estimators_])
    sigma2 = (E ** 2).mean()
    # Mean correlation over all pairs of trees (upper triangle)
    rho = np.corrcoef(E)[np.triu_indices(B, k=1)].mean()
    formula = rho * sigma2 + (1 - rho) * sigma2 / B
    p = forest.predict_proba(X_val)[:, 1]
    print(max_features, round(sigma2, 4), round(rho, 3),
          round(formula, 4), round(((p - y_val) ** 2).mean(), 4),
          round(forest.score(X_val, y_val), 3))

The formula nearly equals the forest's squared error in every row. From bagging to "sqrt", the correlation falls from 0.527 to 0.460 while each tree becomes only slightly worse, and the forest's accuracy rises from 0.924 to 0.953. With one feature per split, the correlation falls further, but the trees become so much worse that the forest loses accuracy. max_features thus trades the accuracy of single trees against their correlation. The scikit-learn User Guide recommends "sqrt" for classification and all features (1.0) for regression as defaults, to be checked by cross-validation.

The number of trees

The code records the 5-fold cross-validated accuracy of forests of increasing size, their fitting time and the time to predict the 171 validation samples.

import time

for B in [1, 10, 50, 100, 200, 500, 1000]:
    forest = RandomForestClassifier(n_estimators=B, random_state=0)
    cv_acc = cross_val_score(forest, X_train, y_train, cv=5).mean()
    start = time.perf_counter()
    forest.fit(X_train, y_train)
    fit_s = time.perf_counter() - start
    start = time.perf_counter()
    forest.predict(X_val)
    pred_ms = 1000 * (time.perf_counter() - start)
    print(B, round(cv_acc, 3), f"{fit_s:.3f} s", f"{pred_ms:.1f} ms")

The accuracy rises steeply up to about 50 trees and then stays between 0.960 and 0.965, a difference of two training samples. It does not decrease as trees are added: Breiman (2001) proved that the error of a forest converges to a limit as the number of trees grows, so more trees do not cause overfitting. The cost, however, grows in proportion to the number of trees, so n_estimators should be just large enough for the score to settle.

The trees are independent, so n_jobs=-1 fits them and predicts with them in parallel on all processor cores, without changing the result. On a synthetic dataset of 20,000 samples and 40 features, n_jobs=-1 makes fitting 200 trees more than five times faster; on the small breast cancer data, the overhead of parallelism cancels most of the gain.

Out-of-bag evaluation

Each OOB prediction uses only trees trained without the sample, so the OOB score estimates performance on new data without a validation split or refitting.

forest = RandomForestClassifier(oob_score=True, random_state=0)
forest.fit(X_train, y_train)
print(round(forest.oob_score_, 3))                # 0.97
print(forest.oob_decision_function_[:2].round(3))
# [[0.033 0.967]
#  [0.03  0.97 ]]

in_bag = np.zeros((100, len(X_train)), dtype=bool)
for b, idx in enumerate(forest.estimators_samples_):
    in_bag[b, idx] = True
oob_share = (~in_bag).mean(axis=0)  # share of trees, per sample
print(round(oob_share.mean(), 3), oob_share.min(), oob_share.max())
# 0.369 0.21 0.54

cv = cross_val_score(RandomForestClassifier(random_state=0),
                     X_train, y_train, cv=5)
print(round(cv.mean(), 3))                        # 0.96
print(round(forest.score(X_val, y_val), 3))       # 0.953

oob_decision_function_ holds the OOB class probabilities and estimators_samples_ the indices drawn for each tree. A sample is out-of-bag for 36.9% of the trees on average, as the bootstrap calculation predicts. The OOB accuracy of 0.970 is close to the cross-validated accuracy of 0.960 and the validation accuracy of 0.953, but needs one fit instead of five.

The OOB score is accuracy for a classifier and $R^2$ for a regressor; oob_score also accepts a function metric(y_true, y_pred) for the task's metric. OOB estimates require bootstrap=True, and with few trees they are noisy, because each prediction averages only about a third of the trees.

Feature importances

The impurity-based importance (MDI) of a tree is defined in Decision Trees, and a forest's feature_importances_ is the mean over its trees. Random feature subsets spread importance across correlated features: the single tree in that module gives worst perimeter 0.8 and worst radius and worst area 0, whereas the forest gives them 0.171, 0.128 and 0.066. MDI is still computed from training data and still favours features with many distinct values.

from sklearn.inspection import permutation_importance

names = list(cancer.feature_names)
mdi = forest.feature_importances_
per_tree = [t.feature_importances_ for t in forest.estimators_]
print(np.allclose(mdi, np.mean(per_tree, axis=0)))  # True
for name in ["worst perimeter", "worst radius", "worst area"]:
    print(name, round(mdi[names.index(name)], 3))
# worst perimeter 0.171
# worst radius 0.128
# worst area 0.066

result = permutation_importance(forest, X_val, y_val,
                                n_repeats=10, random_state=0)
for j in result.importances_mean.argsort()[::-1][:5]:
    print(names[j], round(result.importances_mean[j], 3),
          round(result.importances_std[j], 3))

worst radius 0.018 0.005
worst concave points 0.017 0.01
worst area 0.016 0.006
worst texture 0.013 0.005
mean texture 0.009 0.004

Shuffling worst radius lowers the validation accuracy by 0.018 on average, about 3 of the 171 samples, with a standard deviation of 0.005 over the 10 shuffles; the scoring argument selects another metric. Because the drops are measured on validation data, a feature that only helped to fit noise in the training set receives an importance near 0.

All five drops are small, although the forest's validation accuracy is 0.953. Many breast cancer features are strongly correlated, and when one is shuffled the forest obtains similar information from the others, so correlated features can all appear unimportant; the scikit-learn example on multicollinear features in the resources shows this on the same dataset.

In scikit-learn

RandomForestClassifier and RandomForestRegressor in sklearn.ensemble accept the hyperparameters of a single tree, such as criterion, max_depth and min_samples_leaf, with the same defaults, which grow full trees; the API documentation warns that such trees can be very large on some datasets. The table lists max_features, whose default differs from that of a single tree, and the hyperparameters specific to forests.

predict_proba returns the mean of the trees' class probabilities and predict the most probable class; in Breiman's original method, as the User Guide notes, each tree instead votes for one class.

proba = np.mean([t.predict_proba(X_val) for t in forest.estimators_],
                axis=0)
print(np.allclose(forest.predict_proba(X_val), proba))  # True

RandomForestRegressor averages its trees' predictions, each a mean of training labels, so like a regression tree it never predicts outside the range of its training labels: a forest fitted to $y = 2x$ for $x$ between 0 and 10 predicts 19.88 at both $x = 20$ and $x = 100$. Forests need no feature scaling, and in scikit-learn 1.9.1 they accept NaN values in X.

Extremely randomised trees

ExtraTreesClassifier and ExtraTreesRegressor add further randomness. Each split considers a random subset of features, but for each candidate feature it draws a threshold at random instead of searching, and it uses the best of these random splits. By default each tree sees the whole training set (bootstrap=False). According to the User Guide, this usually reduces variance a little more than a random forest, at the cost of slightly more bias.

from sklearn.ensemble import ExtraTreesClassifier

for model in [RandomForestClassifier(random_state=0),
              ExtraTreesClassifier(random_state=0)]:
    scores = cross_val_score(model, X_train, y_train, cv=5)
    print(type(model).__name__, round(scores.mean(), 3))
# RandomForestClassifier 0.96
# ExtraTreesClassifier 0.975

Random forests in olympiad tasks

A random forest is a strong first model for tabular data and for fixed-length feature vectors, such as embeddings from a pretrained network. It needs no scaling, captures interactions and performs reasonably with its defaults, so it gives a reference score soon after the baseline.

• Compare it with other models. In Approaching a Typical Problem, a default forest reached a validation balanced accuracy of 0.944, against 0.991 for logistic regression. Gradient boosting, covered in XGBoost, is the main alternative among tree ensembles.
• Tune little. The User Guide names n_estimators and max_features as the main parameters to adjust. The OOB score evaluates each setting with one fit; Hyperparameter Tuning covers systematic searches.
• Respect the time limit. Measure fitting and prediction time early, as advised in Reproducibility and time limits. Fewer trees, a smaller max_depth or max_samples, a larger min_samples_leaf and n_jobs=-1 all reduce it.
• Fix the seed. With random_state set, the graded notebook rebuilds the same forest, as shown in Scikit-learn Basics.