import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import permutation_importance
import xgboost as xgb
import warnings
warnings.filterwarnings('ignore')
from pathlib import Path
def data_path(name):
for base in (Path.cwd(), *Path.cwd().parents):
p = base / "data" / "processed" / name
if p.exists():
return p
raise FileNotFoundError(
f"Could not find data/processed/{name}. "
"Run `python scripts/prep_m5.py` from the repository root first."
)
RANDOM_STATE = 42
daily = pd.read_csv(data_path("m5_daily.csv"), parse_dates=["date"])
daily["dow"] = daily["date"].dt.dayofweek
daily["doy"] = daily["date"].dt.dayofyear
FEATURES = ["units_lag1", "units_lag7", "units_lag14", "units_lag28", "units_lag365",
"snap", "is_event", "dow", "doy", "avg_price"]
df = daily.dropna(subset=FEATURES).copy()
N_TEST = 365
train, test = df.iloc[:-N_TEST], df.iloc[-N_TEST:]
X_train, y_train = train[FEATURES].values, train["units"].values
X_test, y_test = test[FEATURES].values, test["units"].values
rmse = lambda a, b: float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) ** 2)))
scoreboard = {
"Benchmark (same day last year)": rmse(y_test, test.units_lag365.values),
"Single tree (Lab 4)": 324,
"20 trees averaged by hand (Lab 4)": 266,
}
for k, v in scoreboard.items():
print(f"{k:<36} {v:>6,.0f}")5 Lab 5: Tree Ensembles — Averaging, Then Boosting
Run this one live: the executable lab.
5.1 What we are doing
At the end of Lab 4 you built a bagged ensemble by accident. You drew twenty bootstrap samples, fitted a tree to each, averaged the predictions, and the RMSE fell from about 324 to 266 — with nothing tuned and no new information added. What you did not do is the one thing that turns bagging into a Random Forest.
Today we do that properly, then find out what else the real algorithm does, then replace averaging with something better.
Same series throughout: daily units for CA_1 FOODS. You have now fitted a GAM to it (Homework 2) and a single tree (Lab 4), so every number today lands on a scale you already know.
Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, scikit-learn, xgboost, and data/processed/m5_daily.csv.
5.2 Setup
Walkthrough
Benchmark (same day last year) 625
Single tree (Lab 4) 324
20 trees averaged by hand (Lab 4) 266The scoreboard carries forward. 625 is the benchmark we have used since Lecture 1. 324 was a single tree. 266 was your hand-built bagging, and those two numbers are hard-coded here deliberately — they are last week’s measured results, not refitted.
Everything today gets added to this list, and at the end you will read it as one table.
5.3 Step 1 — The real thing
rf = RandomForestRegressor(n_estimators=300, random_state=RANDOM_STATE, n_jobs=-1)
rf.fit(X_train, y_train)
scoreboard["Random Forest (300 trees)"] = rmse(y_test, rf.predict(X_test))
print(f"Random Forest: {scoreboard['Random Forest (300 trees)']:,.0f}")Better than your hand-built version, and much better than one tree. Two differences account for it: more trees, and a second source of randomness we are about to meet.
Walkthrough
Random Forest: 248Three hundred trees, no tuning, two lines of code.
248 against your hand-built 266, and against 324 for a single tree. So the library version buys about 7% over what you built yourself in three lines last week.
Sit with how small that is for a moment. Most of the gain was already yours. The idea — resample, refit, average — is where the improvement came from; the professional implementation adds a useful but modest amount on top. Two things account for that remainder: three hundred trees instead of twenty, and a second source of randomness we meet in Step 3.
That is worth saying because the reverse impression is common — that the value is in the library. It isn’t. It’s in the idea, and you already had it.
5.4 Step 2 — How many trees is enough?
The usual worry about a big model is overfitting. Adding trees to a forest is the rare knob where that worry does not apply — but it does stop paying.
counts = [5, 20, 50, 100, 300, 600]
curve = []
for n in counts:
r = RandomForestRegressor(n_estimators=n, random_state=RANDOM_STATE, n_jobs=-1).fit(X_train, y_train)
curve.append(rmse(y_test, r.predict(X_test)))
for n, v in zip(counts, curve):
print(f" n_estimators={n:<4} test RMSE {v:,.0f}")
fig, ax = plt.subplots(figsize=(6, 3.2))
ax.plot(counts, curve, marker="o")
ax.set_xscale("log"); ax.set_xlabel("n_estimators (log scale)"); ax.set_ylabel("test RMSE")
plt.tight_layout()Say what the curve does after about fifty trees, and what that means for how you would spend a compute budget. More trees never hurt accuracy — but every one costs time to fit and to predict.
Walkthrough
n_estimators=5 test RMSE 288
n_estimators=20 test RMSE 257
n_estimators=50 test RMSE 251
n_estimators=100 test RMSE 249
n_estimators=300 test RMSE 248
n_estimators=600 test RMSE 248
Try it
Pause and describe this curve in one sentence, then say what it implies for a compute budget.
It drops hard to about fifty trees, and then it is flat. Fifty trees gets you 251; six hundred gets you 248. You paid twelve times the compute for three units of RMSE — on a series where a typical day is around 2,900 units.
Now the property that makes this curve unusual, and please notice it, because Step 6 breaks it: the curve never turns up. More trees never made the forest worse. That is not the normal behavior of a model knob. Usually “more” eventually means overfitting, and here it simply does not, because each tree is fitted independently and averaging more independent things cannot make the average worse.
So n_estimators in a forest is not really a tuning parameter — it is a budget question. Pick the point where the curve goes flat, spend nothing beyond it, and never worry that you have chosen badly. Hold that thought for twenty minutes.
5.5 Step 3 — The “random” in Random Forest
Bagging alone gives every tree a different sample of rows. A Random Forest also gives each split a different random subset of columns — max_features. Lecture 5 argued this decorrelates the trees, so their errors cancel more completely when averaged.
Test it.
for mf in [None, 0.6, "sqrt", 0.2]:
r = RandomForestRegressor(n_estimators=300, max_features=mf,
random_state=RANDOM_STATE, n_jobs=-1).fit(X_train, y_train)
label = "all features" if mf is None else str(mf)
print(f" max_features={label:<13} test RMSE {rmse(y_test, r.predict(X_test)):,.0f}")
The textbook rule does not win here, and you should know why
max_features="sqrt" is the standard recommendation, and on this dataset it is worse than using every feature. That is not a bug — it is the rule meeting its own boundary condition.
You have 10 features, and one of them (dow) carries most of the signal. sqrt(10) is about 3, so at any given split most trees are not even allowed to look at day-of-week. Decorrelation is only a bargain when you have enough redundant, correlated features that losing a few at each split costs you little.
In Homework 3 you will run a forest on 27 features with no single dominant one. Check whether the ranking flips there. The lesson is the condition, not the default.
Walkthrough
Here is the second randomness. Bagging gives every tree a different sample of rows. A random forest also gives every split a different random subset of columns. That is the entire difference between the two algorithms — and it is why what you built last week was bagging and not a forest.
The argument from Lecture 5: if one feature is strongly predictive, every bagged tree splits on it first, so the trees are near-copies and their errors do not cancel. Deny each split most of the columns and the trees are forced to be different.
You have already seen that failure mode. Last week, twenty out of twenty bootstrap trees found dow at the root.
max_features=all features test RMSE 248
max_features=0.6 test RMSE 247
max_features=sqrt test RMSE 254
max_features=0.2 test RMSE 257max_features="sqrt" is the textbook default. It is the third-best of four here, six units worse than simply using every column. And 0.2 — aggressive decorrelation — is worse still.
This is not a bug and it is not a bad dataset. It is the rule meeting its own boundary condition, and the condition is stated in the callout: you have ten features, one of which carries most of the signal. The square root of ten is about three, so at a typical split most trees are not allowed to look at day-of-week — and day-of-week, as you measured last week, is roughly two thirds of the variance in this series.
Decorrelation is a trade. You give up some accuracy per tree to gain independence between trees. That trade is a bargain when you have many partly-redundant features and no single dominant one. It is a bad trade when you have ten features and one of them is the answer.
Notice the honest middle of the table: 0.6 is the best row, at 247. Mild subsampling helps by a hair. So the finding is not “decorrelation is wrong here,” it is “the standard dose is too strong for this feature set.”
In Homework 3 you run a forest on twenty-seven features with no single dominant one. Check whether this ranking flips. It should — and that difference between the lab and the homework is the actual lesson. What you are learning is the condition, not the default.
5.6 Step 4 — Validation for free
Each bootstrap sample leaves out about a third of the rows. Those out-of-bag rows never trained that tree, so they can score it — a validation estimate with no holdout set.
rf_oob = RandomForestRegressor(n_estimators=300, oob_score=True,
random_state=RANDOM_STATE, n_jobs=-1).fit(X_train, y_train)
oob_rmse = rmse(y_train, rf_oob.oob_prediction_)
print(f"OOB R² : {rf_oob.oob_score_:.3f}")
print(f"OOB RMSE : {oob_rmse:,.0f}")
print(f"test RMSE: {rmse(y_test, rf_oob.predict(X_test)):,.0f}")Compare the two numbers and account for the gap. The OOB estimate is computed on training-period rows scattered through 2012–2015; the test RMSE is a solid block of the final year. Which one would you quote to someone asking how the model will do next month, and why?
Walkthrough
A bootstrap sample of n rows drawn with replacement leaves out about a third of them. Those rows never trained that tree, so they can score it — for free, with no holdout.
OOB R² : 0.809
OOB RMSE : 297
test RMSE: 248
Try it
Pause. The OOB estimate is 297 and the test RMSE is 248, so OOB is pessimistic by about 20%. Work out why before I tell you, and then decide which of the two you would quote to someone asking how this model will do next month.
Two reasons, and they pull in the same direction.
First, an OOB prediction for a given row comes from only the third of trees that did not see it — about a hundred trees rather than three hundred. Step 2 told you exactly what a smaller forest costs. So OOB is scoring a systematically smaller forest than the one you are going to ship.
Second, the rows are different. OOB rows are scattered across 2012 to 2015, the whole training period, including the earlier, noisier stretch. The test set is one solid block of the final year.
So which do you quote? Neither, uncritically. OOB is the cheap sanity check you can compute without giving up any data — genuinely valuable when data is scarce. But it is scoring the wrong model on the wrong period for a forecasting question, and it is not a walk-forward estimate. If someone asks how this does next month, the only honest answer comes from evaluating on time you held out, in order — which is what the test column is, and what Homework 3 makes you do properly.
5.7 Step 5 — Two importances, two different answers
mdi = pd.Series(rf.feature_importances_, index=FEATURES)
perm = permutation_importance(rf, X_test, y_test, n_repeats=10,
random_state=RANDOM_STATE, n_jobs=-1,
scoring="neg_root_mean_squared_error")
perm_imp = pd.Series(perm.importances_mean, index=FEATURES)
comp = pd.DataFrame({"MDI": mdi.round(3), "permutation_ΔRMSE": perm_imp.round(1)})
comp["MDI_rank"] = comp.MDI.rank(ascending=False).astype(int)
comp["perm_rank"] = comp["permutation_ΔRMSE"].rank(ascending=False).astype(int)
print(comp.sort_values("MDI", ascending=False).to_string())Permutation importance is in units of RMSE: shuffle that column on the test set, and this is how much worse the forecast gets. That is directly interpretable in a way MDI is not.
Find the disagreements. Three are worth naming before you read on:
- one lag that MDI ranks above
units_lag1, while permutation ranks it far below snap, which MDI places near the bottom and permutation places near the topavg_priceagainstis_event: MDI ranks the price column above the event flag, while permutation scores price at exactly zero and the event flag at 22 units of RMSE
Your turn: plot the two rankings side by side as horizontal bar charts so the disagreements are visible at a glance.
# ---- your turn: uncomment and fill in the two blanks ----
# fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))
# mdi.sort_values().plot.barh(ax=axes[0], title="MDI (impurity)")
# ___.sort_values().plot.barh(ax=axes[1], title="Permutation (ΔRMSE)") # <- which Series?
# plt.tight_layout()
Walkthrough
MDI permutation_ΔRMSE MDI_rank perm_rank
dow 0.583 284.2 1 1
units_lag28 0.119 36.8 2 4
units_lag1 0.086 116.9 3 2
units_lag7 0.053 35.6 4 5
units_lag365 0.036 10.2 5 8
doy 0.036 18.1 5 7
units_lag14 0.031 7.7 7 9
snap 0.027 40.4 8 3
avg_price 0.017 0.0 9 10
is_event 0.012 22.3 10 6First, the units, because they change what the numbers mean. MDI is a share of impurity reduction — it sums to one and it is computed on the training data by the fitting process. Permutation importance is measured in RMSE: shuffle that column in the test set and see how much worse the forecast gets. So dow at 284 means “destroy day-of-week and you lose 284 units of accuracy.” You can say that sentence to a manager. You cannot say “0.583” to anyone.
Try it
Pause and find the three disagreements before I name them.
One. units_lag28 is second on MDI and fourth on permutation, at 36.8 — while units_lag1, which MDI ranks below it, is worth 116.9. Three times as much. MDI counts how often a feature was used for splitting; permutation measures what happens when you take it away. A feature used constantly for small refinements can lose to one used rarely for large ones.
Two, and this is the one that matters commercially. snap is eighth on MDI and third on permutation, worth 40 units of RMSE. SNAP benefits land on fixed days of the month and this is a food aisle: SNAP days average 3,201 units against 2,839 otherwise, about 13% higher, on 400 of 1,211 training days. MDI nearly buried a real, actionable effect.
Three. Look at avg_price against is_event. MDI ranks price above the event flag — 0.017 to 0.012. Permutation says price is worth exactly zero and the event flag is worth 22.
That inversion has a known mechanism, and Lecture 4’s footnote named it: MDI is biased toward high-cardinality features. avg_price takes 173 distinct values in the training set, so it offers 172 candidate split points, and with that many chances some will reduce impurity by luck alone. is_event is binary — one possible split, no opportunity to get lucky. MDI rewards opportunity; permutation rewards contribution.
The practical rule: MDI is free and permutation is honest. MDI comes out of the fit at no cost and is fine for a quick screen. Anything you are going to act on — or say out loud to a category manager — gets checked with permutation, on held-out data.
Try it
— fill in the one blank and plot the two rankings side by side. The blank is the permutation series. Come back when you can see the snap bar jump.
5.8 Step 6 — Stop averaging. Start correcting.
A forest builds many trees in parallel and averages them. Boosting builds them in sequence, each one fitted to what the previous ones got wrong.
The first four rows hold the learning rate fixed and add rounds — the same experiment you ran on the forest in Step 2. The last two change the rate instead.
for lr, n in [(0.1, 100), (0.1, 300), (0.1, 1000), (0.1, 3000),
(0.3, 100), (0.05, 1000)]:
m = xgb.XGBRegressor(n_estimators=n, learning_rate=lr, max_depth=4,
subsample=0.8, colsample_bytree=0.8,
random_state=RANDOM_STATE).fit(X_train, y_train)
print(f" lr={lr:<5} n={n:<5} train {rmse(y_train, m.predict(X_train)):5,.0f}"
f" test {rmse(y_test, m.predict(X_test)):5,.0f}")best = xgb.XGBRegressor(n_estimators=100, learning_rate=0.1, max_depth=4,
subsample=0.8, colsample_bytree=0.8,
random_state=RANDOM_STATE).fit(X_train, y_train)
scoreboard["XGBoost (lr=0.1, 100 trees)"] = rmse(y_test, best.predict(X_test))
print(pd.Series(scoreboard).round(0).astype(int).sort_values(ascending=False).to_string())Read the first four rows as a curve, then compare it with Step 2. Training error falls to zero. Test error bottoms out around 300 rounds and then climbs back up — 3,000 rounds is worse than 100. In Step 2 you added trees to a forest all the way to 600 and it never once hurt.
That is the practical difference between the two ensembles. A forest’s trees are fitted independently and averaged, so extra ones can only sharpen the average. Boosting’s trees are fitted to the previous ensemble’s residuals, so extra ones keep going — eventually correcting noise. n_estimators is a free choice in a forest and a tuning parameter in boosting.
The last two rows say the obvious complement: a big rate with few rounds and a small rate with many land in much the same place, because roughly it is the product that sets the budget. What you cannot do is set a rate and then let the rounds run.
Note which configuration we put on the scoreboard below — lr=0.1, n=100, the conventional default, not the row with the best test RMSE. Picking the winner out of a column computed on the test set is the mistake Lecture 1 warned about, and it does not stop being a mistake because the grid is small.
Walkthrough
Everything so far today has been the same idea: build many trees independently, average them. Boosting is a different idea. Trees are built in sequence, and each new one is fitted to what the ensemble so far got wrong. Nothing is averaged. Every tree is a correction.
lr=0.1 n=100 train 171 test 244
lr=0.1 n=300 train 95 test 243
lr=0.1 n=1000 train 19 test 250
lr=0.1 n=3000 train 0 test 253
lr=0.3 n=100 train 95 test 238
lr=0.05 n=1000 train 54 test 239Read the first four rows as a curve — the learning rate is fixed, only the number of rounds changes. That is the identical experiment you ran in Step 2 on the forest.
Training error goes 171, 95, 19, zero. Test error goes 244, 243, 250, 253. It bottoms out around three hundred rounds and then turns up.
Now put that next to Step 2, where you took a forest from five trees to six hundred and it never got worse once.
Same data, same features, same family of base learner. Opposite behavior on the same knob. The reason is structural, and it is worth stating carefully: a forest’s trees are fitted independently of one another, so an extra one only adds to an average and cannot corrupt what is already there. A boosted tree is fitted to the residuals of everything before it — so once the real signal is used up, the next tree fits noise, and it is added to the ensemble anyway, with nothing to average it away.
So: n_estimators in a forest is a budget. n_estimators in boosting is a tuning parameter, and getting it wrong costs you accuracy in a way no amount of forest-tuning ever will. That is also why boosting has a learning rate — it is the knob that decides how much damage each additional tree is allowed to do.
The last two rows are the complement. A big rate with few rounds (238) and a small rate with many (239) end up in much the same place: roughly, the product sets the budget. What you cannot do is fix a rate and let the rounds run.
Benchmark (same day last year) 625
Single tree (Lab 4) 324
20 trees averaged by hand (Lab 4) 266
Random Forest (300 trees) 248
XGBoost (lr=0.1, 100 trees) 244Now, look at which XGBoost row we put on that scoreboard. lr=0.1, n=100 — 244. Not the lr=0.3 row at 238, which was the best number in the grid.
That is deliberate, and it is not modesty. Every test-RMSE number in the grid above was computed on the year we held out. If I pick the winner out of that column and then report it as my result, I have used the test set to choose the model and I no longer have an honest estimate of anything. Lecture 1 said this about forecast evaluation and it does not stop being true because the grid has only six rows. So the scoreboard gets the conventional default, and the way to actually choose here is walk-forward validation inside the training period — Homework 4.
5.9 Before you leave
Your scoreboard should run from a benchmark near 625 down to roughly 245, across five models fitted to one series.
Discuss (the ranking): Boosting beat the forest here. Name one property of this dataset that favors boosting, and one situation where you would still reach for the forest first — thinking about tuning effort and about what happens when nobody retunes the model for a year.
Discuss (the importances): You have two importance tables that disagree about
snap. A category manager asks whether SNAP timing matters for stocking. Which table do you show them, what do you say, and what would you want to check before committing to an answer?Discuss (the whole arc): Across Labs 4 and 5 the RMSE fell from 324 to about 245 on the same series with the same features. What was actually purchased with that improvement, and what was given up along the way? Be specific about the thing you can no longer do that you could do in Step 1 of Lab 4.
Solutions posted to Canvas after class.
Walkthrough
Read the scoreboard as one object. 625 to 244 — a 61% reduction against the benchmark, on one series, with the same ten features throughout. Almost all of it arrived in two jumps: a single tree over the benchmark, and averaging over the single tree. The last two steps, the forest and boosting, bought 22 units between them.
Try it
Three questions on the discussion board.
The ranking. Boosting won, 244 to 248 — by about one and a half percent. Name one property of this dataset that favours it. Then name a situation where you would still reach for the forest first, thinking about two things specifically: how much tuning each one needed today, and what happens to each when nobody retunes it for a year. was fire-and-forget, and boosting’s best setting is a live decision that can go stale.*
The importances. You have two tables that disagree about snap — eighth on one, third on the other. A category manager asks whether SNAP timing matters for stocking. Which table do you show them, what exactly do you say, and what would you check before committing to an answer? Be concrete about that last part.
The arc. Across Labs 4 and 5 the error fell from 324 to 244 on the same series with the same features. What was purchased with that, and what was given up? Be specific about the thing you could do in Step 1 of Lab 4 that you cannot do now.
I’ll say the shape of that last answer, because it is the through-line of the whole course and I do not want it to be missed. In Lab 4, Step 1, you printed a tree and read a path out loud — weekday, four weeks ago was busy, yesterday was busy, predict 2,950. You cannot do that now. There is no path to read. Three hundred trees, or a hundred sequential corrections, and the best account you can give of the model is a permutation importance table.
You bought 25% accuracy with the entire explanation. Sometimes that is obviously right. Sometimes you are standing in front of a category manager and it is obviously wrong. Knowing which is which is the actual skill.
Solutions for the coded parts go up on Canvas after the deadline. Next week: regularization, and what to do when you have far more features than this.
5.10 What your run should produce
Your numbers will differ where a model is stochastic; these are the values this lab produces on the pinned environment.
Deterministic; random_state=42 throughout. Whole lab runs in about 10 seconds.
| Quantity | Value |
|---|---|
| Scoreboard carried in | benchmark 625 · single tree 324 · 20 trees averaged 266 |
| Random Forest, 300 trees | 248 |
| Tree-count curve | 5: 288 · 20: 257 · 50: 251 · 100: 249 · 300: 248 · 600: 248 (never rises) |
max_features |
all: 248 · 0.6: 247 · sqrt: 254 · 0.2: 257 |
| OOB | R² 0.809 · RMSE 297 vs test 248 |
| MDI top three | dow 0.583 · units_lag28 0.119 · units_lag1 0.086 |
| Permutation top three (ΔRMSE) | dow 284.2 · units_lag1 116.9 · snap 40.4 |
snap |
MDI rank 8 → permutation rank 3 (40.4 RMSE) |
avg_price vs is_event |
MDI 0.017 > 0.012; permutation 0.0 < 22.3 |
| Cardinality behind that inversion | avg_price 173 distinct values · is_event 2 |
| SNAP effect in the data | 3,201 units on SNAP days vs 2,839 otherwise (+13%), 400 of 1,211 days |
Boosting, fixed lr=0.1 |
n=100: 244 · n=300: 243 · n=1000: 250 · n=3000: 253 (turns up) |
| Boosting, rate changes | lr=0.3 n=100: 238 · lr=0.05 n=1000: 239 |
| Scoreboard entry (chosen as the default, not the winner) | XGBoost lr=0.1 n=100: 244 |