B-forest 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 of each day; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them 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 with python3 randomforest.py, from any directory — it writes its output files into the current directory. Requires NumPy and Matplotlib only, and uses 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
"""Three bootstrap regression trees and their average at Krems: the
forest fits no worse than its trees on average.
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 of each day; this
script downloads the records for 2024 from the GeoSphere data hub
(dataset klima-v2-1d) and writes them 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 bootstrap
sample of the 366 days, and the random forest predicts by averaging
the three tree predictions. With a single feature there is no room
for feature subsampling, so the forest's randomness comes from the
bootstrap alone, as in bagging.
The demo checks the entry's central claims: each tree is a piecewise
constant map with few pieces; the bootstrap makes the trees differ;
and the squared-error risk of the averaged prediction on the full
dataset never exceeds the average of the trees' squared-error risks
(the algebraic identity behind variance reduction by averaging).
Deterministic: the bootstrap samples use fixed seeds. Self-contained:
numpy + matplotlib only (stdlib urllib for the download).
Blocks
------
[B-fetch] Download the 366 daily temperature pairs at Krems for 2024
and write them 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 bootstrap
sample (fixed seeds); check that each tree is piecewise
constant with at most four pieces, fits its bootstrap sample
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 squared-error risk on all
366 days is smaller than the average of the trees' risks.
Outputs
-------
randomforest_weather.csv : date, tmin, tmax for the 366 days of 2024
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.png : preview (checking only) -- scatterplot of
the days, the three tree curves, and the
forest curve
"""
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 temperature pairs at Krems for 2024 and write them 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&station_ids=3805"
"&start=2024-01-01&end=2024-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 2024",
len(tmin) == 366)
feb1 = stamps.index("2024-02-01")
check("[B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4)",
tmin[feb1] == -3.8 and tmax[feb1] == 10.4)
[ok] [B-fetch] 366 daily temperature pairs downloaded for 2024 [ok] [B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4)
Train three depth-2 regression trees, each on a bootstrap sample (fixed seeds); check that each tree is piecewise constant with at most four pieces, fits its bootstrap sample 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)) # bootstrap sample
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 bootstrap sample 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 bootstrap 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 bootstrap sample better than a constant [ok] [B-trees] the bootstrap 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 squared-error risk on all 366 days is smaller than the average of the trees' risks.
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 smaller squared-error risk "
"than the trees on average",
forest_risk < np.mean(tree_risks) - 1e-9)
print(f" tree risks {[round(r, 2) for r in tree_risks]}, "
f"forest risk {forest_risk:.2f}")
fig, ax = plt.subplots(figsize=(6.4, 4.4))
ax.plot(tmin, tmax, "o", ms=2, color="0.6", label="days of 2024 (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 bootstrap trees and their average (random forest)")
ax.legend(frameon=False, fontsize=8)
fig.tight_layout()
fig.savefig(OUT_DIR / "randomforest.png", dpi=150)
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 ""))
[ok] [B-forest] the forest curve is the pointwise mean of the tree curves [ok] [B-forest] the averaged prediction has smaller squared-error risk than the trees on average tree risks [np.float64(22.32), np.float64(22.36), np.float64(22.52)], forest risk 18.84 7/7 checks passed
B-forest writes when the script runs