B-splits writes when the script runsDictionary of Applied Machine Learning · random forest
Numerical companion to the entry random forest: it recomputes what the entry states and prints one line per check
Numerical companion to the glossary entry 'randomforest' (random forest). The GeoSphere Austria weather station Krems (station id 3805) records the minimum and maximum air temperature and the sunshine hours of each day; this script downloads the records for 2020 from the GeoSphere data hub (dataset klima-v2-1d) and writes the temperatures to randomforest_weather.csv. Each day is a data point whose feature is the minimum temperature tmin and whose label is the maximum temperature tmax.
Run it without installing anything:uv run https://dictionaryofml.org/terms/randomforest.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download randomforest.py into a folder and run uv run randomforest.py there. With NumPy and Matplotlib already installed, python3 randomforest.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download randomforest.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""Random forest at Krems: three trees and their average, which features
the splits use, and what withholding features costs.
Purpose
-------
Numerical companion to the glossary entry 'randomforest' (random
forest). The GeoSphere Austria weather station Krems (station id 3805)
records the minimum and maximum air temperature and the sunshine hours
of each day; this script downloads the records for 2020 from the
GeoSphere data hub (dataset klima-v2-1d) and writes the temperatures to
randomforest_weather.csv. Each day is a data point whose feature is the
minimum temperature tmin and whose label is the maximum temperature tmax.
Three depth-2 regression trees are trained, each on its own random
selection (with replacement) of the 366 days, and the random forest
predicts by averaging the three tree predictions. A second experiment
adds three more features of the day and grows forests whose base
learners are also offered only a random pair of the four features, as
the entry's random forest does.
The demo checks the entry's central claims: each tree is a piecewise
constant map with few pieces; the random selection makes the trees
differ; the mean squared prediction error of the averaged prediction on
the full dataset never exceeds the average of the trees' errors (the
algebraic identity behind variance reduction by averaging); the feature
importance scores are shares of the splits; and withholding features
lowers the correlation between the trees while raising the variance of
each of them, which is the trade the entry describes.
Deterministic: every random selection uses a fixed seed. Self-contained:
numpy + matplotlib only (stdlib urllib for the download).
Blocks
------
[B-fetch] Download the 366 daily records at Krems for 2020 and write the
temperature pairs to randomforest_weather.csv; check the count
and one pinned value against the archive.
[B-trees] Train three depth-2 regression trees, each on a random
selection of the days (fixed seeds); check that each tree is
piecewise constant with at most four pieces, fits its
selection better than a constant, and that the trees differ
pairwise.
[B-forest] Average the three trees into the random forest prediction;
check that the forest curve is the pointwise mean of the tree
curves and that the forest's mean squared prediction error on all
366 days is smaller than the average of the trees' errors.
[B-splits] Grow a 60-tree forest on four features of the day, record the
feature tested at every split, and report the importance of a
feature as the fraction of all splits that test it. Each tree
sees randomly selected days and a randomly selected pair of
the four features. Written to randomforest_importance.csv.
[B-trade] On 120 training sets of 200 days each, grow two trees per
training set, once with all four features offered and once
with a random pair; estimate at every day the variance of a
single tree's prediction and the correlation between the two
trees, as bagging.py does. Withholding features lowers the
correlation and raises the variance -- the trade. Written to
randomforest_trade.csv.
Outputs
-------
randomforest_weather.csv : date, tmin, tmax for the 366 days of 2020
randomforest_tree1.csv : tmin, pred -- first tree on a tmin grid
randomforest_tree2.csv : tmin, pred -- second tree on the grid
randomforest_tree3.csv : tmin, pred -- third tree on the grid
randomforest_forest.csv : tmin, pred -- the averaged (forest) curve
randomforest_importance.csv : feature, number of splits testing it, that
count as a fraction of all splits, and the
bar position in the entry's figure
randomforest_trade.csv : features offered, mean correlation between two
trees, mean variance of one tree
randomforest.png : preview (checking only) -- the three trees and
the forest, and the importance scores
"""
import json
import urllib.request
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
OUT_DIR = Path(__file__).parent
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
Download the 366 daily records at Krems for 2020 and write the temperature pairs to randomforest_weather.csv; check the count and one pinned value against the archive.
URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/"
"klima-v2-1d?parameters=tlmin,tlmax,so_h&station_ids=3805"
"&start=2020-01-01&end=2020-12-31")
with urllib.request.urlopen(URL, timeout=120) as resp:
payload = json.load(resp)
params = payload["features"][0]["properties"]["parameters"]
stamps = [t[:10] for t in payload["timestamps"]]
tmin = np.array(params["tlmin"]["data"], dtype=float)
tmax = np.array(params["tlmax"]["data"], dtype=float)
with open(OUT_DIR / "randomforest_weather.csv", "w") as f:
f.write("date,tmin,tmax\n")
for d, a, b in zip(stamps, tmin, tmax):
f.write(f"{d},{a:.1f},{b:.1f}\n")
check("[B-fetch] 366 daily temperature pairs downloaded for 2020",
len(tmin) == 366)
feb1 = stamps.index("2020-02-01")
check("[B-fetch] the record matches the archive (Feb 1: 4.2 to 17.0)",
tmin[feb1] == 4.2 and tmax[feb1] == 17.0)
[ok] [B-fetch] 366 daily temperature pairs downloaded for 2020 [ok] [B-fetch] the record matches the archive (Feb 1: 4.2 to 17.0)
Train three depth-2 regression trees, each on a random selection of the days (fixed seeds); check that each tree is piecewise constant with at most four pieces, fits its selection better than a constant, and that the trees differ pairwise.
def best_split(x, y):
"""Threshold minimizing the summed squared error of the two parts,
or None if no split with at least 5 data points per part exists."""
order = np.argsort(x)
xs, ys = x[order], y[order]
best, best_sse = None, np.inf
for i in range(5, len(xs) - 5 + 1):
if xs[i - 1] == xs[i]:
continue
left, right = ys[:i], ys[i:]
sse = ((left - left.mean()) ** 2).sum() \
+ ((right - right.mean()) ** 2).sum()
if sse < best_sse:
best_sse, best = sse, (xs[i - 1] + xs[i]) / 2
return best
def depth2_tree(x, y):
"""Thresholds and leaf means of a depth-2 regression tree."""
t0 = best_split(x, y)
cuts = []
for side in (x <= t0, x > t0):
t = best_split(x[side], y[side])
if t is not None:
cuts.append(t)
cuts = sorted(cuts + [t0])
edges = [-np.inf] + cuts + [np.inf]
means = [y[(x > lo) & (x <= hi)].mean()
for lo, hi in zip(edges[:-1], edges[1:])]
return np.array(cuts), np.array(means)
def predict(cuts, means, xq):
return means[np.searchsorted(cuts, xq)]
grid = np.linspace(tmin.min() - 1.0, tmin.max() + 1.0, 400)
trees, curves = [], []
for seed in (1, 2, 3):
rng = np.random.default_rng(seed)
idx = rng.integers(0, len(tmin), len(tmin)) # random selection, with replacement
cuts, means = depth2_tree(tmin[idx], tmax[idx])
trees.append((idx, cuts, means))
curves.append(predict(cuts, means, grid))
np.savetxt(OUT_DIR / f"randomforest_tree{seed}.csv",
np.stack([grid, curves[-1]], 1),
delimiter=",", header="tmin,pred", comments="", fmt="%.2f")
check("[B-trees] each tree is piecewise constant with at most 4 pieces",
all(len(means) <= 4 for _, _, means in trees))
check("[B-trees] each tree fits its own selection of the days better than a constant",
all(((tmax[idx] - predict(cuts, means, tmin[idx])) ** 2).mean()
< tmax[idx].var() for idx, cuts, means in trees))
check("[B-trees] the random selection makes the trees differ pairwise",
all(np.max(np.abs(curves[i] - curves[j])) > 0.5
for i in range(3) for j in range(i + 1, 3)))
[ok] [B-trees] each tree is piecewise constant with at most 4 pieces [ok] [B-trees] each tree fits its own selection of the days better than a constant [ok] [B-trees] the random selection makes the trees differ pairwise
Average the three trees into the random forest prediction; check that the forest curve is the pointwise mean of the tree curves and that the forest's mean squared prediction error on all 366 days is smaller than the average of the trees' errors.
forest = np.mean(curves, axis=0)
np.savetxt(OUT_DIR / "randomforest_forest.csv",
np.stack([grid, forest], 1),
delimiter=",", header="tmin,pred", comments="", fmt="%.2f")
tree_risks = [((tmax - predict(cuts, means, tmin)) ** 2).mean()
for _, cuts, means in trees]
forest_pred = np.mean([predict(cuts, means, tmin)
for _, cuts, means in trees], axis=0)
forest_risk = ((tmax - forest_pred) ** 2).mean()
check("[B-forest] the forest curve is the pointwise mean of the tree curves",
np.allclose(forest, np.mean(curves, axis=0)))
check("[B-forest] the averaged prediction has a smaller mean squared prediction error "
"than the trees on average",
forest_risk < np.mean(tree_risks) - 1e-9)
print(f" tree errors {[round(r, 2) for r in tree_risks]}, "
f"forest error {forest_risk:.2f}")
fig, (ax, axb) = plt.subplots(1, 2, figsize=(10.4, 4.2),
gridspec_kw={"width_ratios": [1.5, 1]})
ax.plot(tmin, tmax, "o", ms=2, color="0.6", label="days of 2020 (Krems)")
for c, (style, name) in zip(curves, [(":", "tree 1"), ("--", "tree 2"),
("-.", "tree 3")]):
ax.plot(grid, c, style, lw=1.2, label=name)
ax.plot(grid, forest, "-", lw=2.2, color="black", label="random forest")
ax.set_xlabel("minimum temperature of the day (°C)")
ax.set_ylabel("maximum temperature of the day (°C)")
ax.set_title("three trees and their average (random forest)")
ax.legend(frameon=False, fontsize=8)
[ok] [B-forest] the forest curve is the pointwise mean of the tree curves [ok] [B-forest] the averaged prediction has a smaller mean squared prediction error than the trees on average tree errors [np.float64(25.2), np.float64(26.51), np.float64(25.29)], forest error 24.10
Grow a 60-tree forest on four features of the day, record the feature tested at every split, and report the importance of a feature as the fraction of all splits that test it. Each tree sees randomly selected days and a randomly selected pair of the four features. Written to randomforest_importance.csv.
# A forest is grown on four features of the day, and every split it makes is
# recorded. The importance of a feature is the fraction of all splits in the
# forest that test it -- read straight off the trees, with no impurity
# bookkeeping. Each base learner sees randomly selected days AND a randomly
# selected pair of the four features, so a feature cannot be used by a tree
# that was not offered it.
sun = np.array(params["so_h"]["data"], dtype=float)
doy = np.arange(1.0, len(tmax) + 1.0)
FEATURE_NAMES = ["minimum temperature", "sunshine hours",
"day of the year", "maximum temperature of the previous day"]
X = np.column_stack([tmin[1:], sun[1:], doy[1:], tmax[:-1]])
Y = tmax[1:]
check("[B-splits] no feature has a missing value over the 365 usable days",
X.shape == (365, 4) and not np.isnan(X).any() and not np.isnan(Y).any())
def fit_tree(xs, ys, cols, depth, splits=None):
"""Depth-limited regression tree over the feature columns in `cols`.
Returns a nested tuple, ("leaf", mean) or ("node", j, cut, left,
right), and appends the column index of every split to `splits`."""
if depth == 0 or len(ys) < 8:
return ("leaf", float(ys.mean()))
best = None
for j in cols:
v = xs[:, j]
o = np.argsort(v)
vs, ws = v[o], ys[o]
for i in range(4, len(vs) - 4):
if vs[i] == vs[i - 1]:
continue
lo, hi = ws[:i], ws[i:]
sse = ((lo - lo.mean()) ** 2).sum() + ((hi - hi.mean()) ** 2).sum()
if best is None or sse < best[0]:
best = (sse, j, 0.5 * (vs[i] + vs[i - 1]))
if best is None:
return ("leaf", float(ys.mean()))
_, j, cut = best
if splits is not None:
splits.append(j)
left = xs[:, j] <= cut
return ("node", j, cut,
fit_tree(xs[left], ys[left], cols, depth - 1, splits),
fit_tree(xs[~left], ys[~left], cols, depth - 1, splits))
def tree_predict(node, xs):
"""Evaluate a fitted tree on every row of xs."""
out = np.empty(len(xs))
for r in range(len(xs)):
n = node
while n[0] == "node":
n = n[3] if xs[r, n[1]] <= n[2] else n[4]
out[r] = n[1]
return out
srng = np.random.default_rng(11)
split_cols = []
N_TREES, N_FEAT = 60, 2
for _ in range(N_TREES):
rows = srng.integers(0, len(Y), len(Y)) # randomly selected days
cols = srng.choice(4, size=N_FEAT, replace=False) # and features
fit_tree(X[rows], Y[rows], cols, 3, splits=split_cols)
counts = np.bincount(split_cols, minlength=4)
frac = counts / counts.sum()
order = np.argsort(-frac)
print(f" {len(split_cols)} splits over {N_TREES} trees; "
f"importance = fraction of splits testing the feature")
for j in order:
print(f" {FEATURE_NAMES[j]:42s} {counts[j]:4d} {frac[j]*100:5.1f}%")
check("[B-splits] the importance scores are fractions that sum to one",
abs(frac.sum() - 1.0) < 1e-12 and bool(np.all(frac >= 0)))
check("[B-splits] every feature was offered to some tree, so a zero score "
"would mean unused rather than unavailable", bool(np.all(counts > 0)))
check(f"[B-splits] the two temperature features together take "
f"{(frac[0] + frac[3]) * 100:.0f}% of the splits, more than the "
f"calendar day at {frac[2] * 100:.0f}%", frac[0] + frac[3] > frac[2])
# rank = bar position in the entry's figure, 1 for the smallest share
with open(OUT_DIR / "randomforest_importance.csv", "w") as f:
f.write("feature,splits,fraction,rank\n")
for pos, j in enumerate(order):
f.write(f"{FEATURE_NAMES[j]},{counts[j]},{frac[j]:.4f},{len(order) - pos}\n")
# ---- preview, second panel: the importance shares of [B-splits]
axb.barh([FEATURE_NAMES[j].replace(" of the previous day", "\nof the previous day")
for j in order][::-1], frac[order][::-1], color="0.8", edgecolor="k")
axb.set_xlabel("share of all splits that test the feature")
axb.set_ylabel("feature of the day")
axb.set_title("share of splits per feature")
axb.tick_params(axis="y", labelsize=7)
fig.tight_layout()
fig.savefig(OUT_DIR / "randomforest.png", dpi=150)
[ok] [B-splits] no feature has a missing value over the 365 usable days
420 splits over 60 trees; importance = fraction of splits testing the feature
maximum temperature of the previous day 157 37.4%
day of the year 120 28.6%
minimum temperature 72 17.1%
sunshine hours 71 16.9%
[ok] [B-splits] the importance scores are fractions that sum to one
[ok] [B-splits] every feature was offered to some tree, so a zero score would mean unused rather than unavailable
[ok] [B-splits] the two temperature features together take 55% of the splits, more than the calendar day at 29%
features offered: all four rho = 0.102 sigma^2 = 4.82
features offered: a random pair rho = 0.022 sigma^2 = 8.81
B-splits writes when the script runsOn 120 training sets of 200 days each, grow two trees per training set, once with all four features offered and once with a random pair; estimate at every day the variance of a single tree's prediction and the correlation between the two trees, as bagging.py does. Withholding features lowers the correlation and raises the variance -- the trade. Written to randomforest_trade.csv.
# The entry says restricting the features is a trade: it lowers the
# correlation rho between the base learners and raises the variance sigma^2
# of each of them. Measured here the way bagging.py measures the variance
# identity: 120 training sets of 200 days, two trees on each, the same two
# random selections of days for both settings, so the only difference is
# whether a tree is offered all four features or a random pair. At every
# day, sigma^2 is the variance of one tree's prediction over (training set,
# tree) and rho the correlation between the two trees of a training set.
trng = np.random.default_rng(13)
N_TS, TS_SIZE = 120, 200
ALL, PAIR = np.arange(4), None
tr_pred = {"all four": np.empty((N_TS, 2, len(Y))),
"a random pair": np.empty((N_TS, 2, len(Y)))}
for t in range(N_TS):
ts = trng.choice(len(Y), size=TS_SIZE, replace=False)
for b in range(2):
sel = ts[trng.integers(0, TS_SIZE, TS_SIZE)] # one selection ...
pair = trng.choice(4, size=2, replace=False) # ... one feature pair
tr_pred["all four"][t, b] = tree_predict(fit_tree(X[sel], Y[sel], ALL, 3), X)
tr_pred["a random pair"][t, b] = tree_predict(fit_tree(X[sel], Y[sel], pair, 3), X)
trade = {}
for name, pr in tr_pred.items():
c = pr - pr.mean(axis=(0, 1)) # centre per day
sigma2 = (c ** 2).mean(axis=(0, 1))
rho = (c[:, 0] * c[:, 1]).mean(axis=0) / sigma2
trade[name] = (float(rho.mean()), float(sigma2.mean()))
print(f" features offered: {name:14s} rho = {rho.mean():.3f} "
f"sigma^2 = {sigma2.mean():.2f}")
r4, s4 = trade["all four"]; r2, s2 = trade["a random pair"]
check(f"[B-trade] offering a random pair instead of all four features lowers "
f"the correlation between two trees, {r4:.3f} -> {r2:.3f}", r2 < r4)
check(f"[B-trade] and raises the variance of a single tree's prediction, "
f"{s4:.2f} -> {s2:.2f}", s2 > s4)
with open(OUT_DIR / "randomforest_trade.csv", "w") as f:
f.write("features_offered,rho,sigma2\n")
for name, (rho_m, s2_m) in trade.items():
f.write(f"{name},{rho_m:.4f},{s2_m:.4f}\n")
failed = [name for name, ok in report if not ok]
print(f"{len(report) - len(failed)}/{len(report)} checks passed"
+ (f", FAILED: {failed}" if failed else ""))
if failed:
raise SystemExit(1)
[ok] [B-trade] offering a random pair instead of all four features lowers the correlation between two trees, 0.102 -> 0.022 [ok] [B-trade] and raises the variance of a single tree's prediction, 4.82 -> 8.81 13/13 checks passed