ml-003 · Random forest vs gradient boosting — Python

Problem statement

Open this practice in GitHub Codespaces runs in your browser; open this folder under practices/ after it starts

Simulate a nonlinear tabular regression task in-script (no files): n = 2000 rows, 8 features drawn Uniform(-2, 2), outcome y = 4*sin(1.5*f1) + 3*f2*f3 + 2*f4^2 - 1.5*f5 + noise, with f6-f8 pure noise features, seed 42. On one fixed train/validation split (80/20), compare RandomForestRegressor (all defaults) against XGBRegressor at learning rates 0.30 and 0.05 across ensemble sizes 25-400, printing a validation-RMSE table. Then fit XGBoost once more with early_stopping_rounds and report the best iteration, and finally compare the two models’ top-3 feature importances. You know it is correct when: both models beat the predict-the- mean baseline, the RF column barely moves across sizes while the XGBoost columns move a lot, early stopping halts before the round budget, and the noise features f6-f8 sit at the bottom of both importance rankings.

Learning objectives

After completing this practice you should be able to:

• explain bagging vs boosting and predict how each responds to more trees
• fit RandomForestRegressor and XGBRegressor on the same split and compare them fairly with a validation-RMSE curve
• explain why RF is robust to its defaults while boosting requires the learning rate and the number of rounds to be tuned together
• use XGBoost early stopping (early_stopping_rounds + eval_set) and read best_iteration
• compare impurity-based and gain-based feature importances and explain why they can disagree on interaction terms

Concepts and definitions

See concepts for the shared background. Python-specific notes: since xgboost 2.0, early_stopping_rounds is a constructor argument of XGBRegressor (passing it to .fit() raises an error), and after early stopping .predict() automatically uses the best iteration. sklearn’s RF feature_importances_ is impurity-based; XGBoost’s is gain-based by default.

Input / output description

Item	Description
Input	Nothing external — data simulated in-script with numpy.random.default_rng(42): X (2000 x 8 DataFrame, columns f1-f8), y (Series)
Output	Printed report: baseline RMSE; a 5-row table of validation RMSE (rf_default, xgb_lr_030, xgb_lr_005) by n_estimators; early-stopping best iteration + RMSE; top-3 importances per model and their overlap count

Before-coding checklist

• I can restate the problem without looking at it
• I know what the output should look like for a tiny example
• I have identified the tricky part
• I know how I will check my answer

My attempt

Work here. Run cells interactively; nothing below overwrites this section.

# Milestone 1: simulate the data (n=2000, 8 uniform features, nonlinear
# signal on f1-f5 only, seed 42) and make an 80/20 train/validation split.

# Milestone 2: baseline RMSE (predict the training mean), then fit
# RandomForestRegressor with defaults at n_estimators in [25, 50, 100, 200,
# 400] and record validation RMSE for each size.

# Milestone 3: same sizes for XGBRegressor (max_depth=3) at learning_rate
# 0.30 and 0.05; assemble the three curves into one comparison table.

# Milestone 4: refit XGBoost with n_estimators=500 and
# early_stopping_rounds=30 monitored on the validation set; report
# best_iteration and the resulting RMSE.

# Milestone 5: compare top-3 feature importances (RF impurity vs XGB gain)
# and check how much importance lands on the noise features f6-f8.

Hint 1

TipHint 1 — a nudge

Write small functions, not one script: make_data, split_data, rmse, then one function per model that takes (X_tr, y_tr, X_va, y_va, n_estimators, ...) and returns a validation RMSE. The comparison table is then just a loop over [25, 50, 100, 200, 400] building a DataFrame. Fix random_state/seed everywhere so reruns match.

Hint 2

TipHint 2 — almost a spoiler

For early stopping the modern xgboost API is: XGBRegressor(n_estimators=500, learning_rate=0.15, max_depth=3, early_stopping_rounds=30, eval_metric="rmse") and then model.fit(X_tr, y_tr, eval_set=[(X_va, y_va)], verbose=False); read model.best_iteration afterwards. For the importance comparison, wrap both feature_importances_ arrays in pd.Series(..., index=X.columns), sort descending, and intersect the top-3 index sets.

Complete solution

NoteOfficial solution — open only after a real attempt

The official, tested script lives at solutions/solution.py. Run it with:

python practices/ml-003-forest-vs-xgboost/solutions/solution.py

Core code (the script adds a main() that prints the full report):

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor

SEED = 42
FEATURE_NAMES = [f"f{i}" for i in range(1, 9)]
ENSEMBLE_SIZES = [25, 50, 100, 200, 400]

def make_data(n_rows=2000, seed=SEED):
    rng = np.random.default_rng(seed)
    X = pd.DataFrame(rng.uniform(-2.0, 2.0, size=(n_rows, 8)),
                     columns=FEATURE_NAMES)
    signal = (4.0 * np.sin(1.5 * X["f1"]) + 3.0 * X["f2"] * X["f3"]
              + 2.0 * X["f4"] ** 2 - 1.5 * X["f5"])
    y = signal + rng.normal(0.0, 1.5, size=n_rows)
    return X, pd.Series(np.asarray(y), name="y")

def split_data(X, y, valid_frac=0.2, seed=SEED):
    rng = np.random.default_rng(seed)
    idx = rng.permutation(len(X))
    n_valid = int(len(X) * valid_frac)
    va, tr = idx[:n_valid], idx[n_valid:]
    return X.iloc[tr], X.iloc[va], y.iloc[tr], y.iloc[va]

def rmse(y_true, y_pred):
    return float(np.sqrt(np.mean((np.asarray(y_true, dtype=float)
                                  - np.asarray(y_pred, dtype=float)) ** 2)))

def rf_validation_rmse(X_tr, y_tr, X_va, y_va, n_estimators):
    model = RandomForestRegressor(n_estimators=n_estimators,
                                  random_state=SEED, n_jobs=-1)
    model.fit(X_tr, y_tr)
    return rmse(y_va, model.predict(X_va))

def xgb_validation_rmse(X_tr, y_tr, X_va, y_va, n_estimators, learning_rate):
    model = XGBRegressor(n_estimators=n_estimators,
                         learning_rate=learning_rate, max_depth=3,
                         random_state=SEED, n_jobs=-1, verbosity=0)
    model.fit(X_tr, y_tr)
    return rmse(y_va, model.predict(X_va))

def validation_curve_table(X_tr, y_tr, X_va, y_va, sizes=ENSEMBLE_SIZES):
    return pd.DataFrame([
        {"n_estimators": s,
         "rf_default": rf_validation_rmse(X_tr, y_tr, X_va, y_va, s),
         "xgb_lr_030": xgb_validation_rmse(X_tr, y_tr, X_va, y_va, s, 0.30),
         "xgb_lr_005": xgb_validation_rmse(X_tr, y_tr, X_va, y_va, s, 0.05)}
        for s in sizes])

def fit_xgb_early_stopping(X_tr, y_tr, X_va, y_va, learning_rate=0.15,
                           max_rounds=500, patience=30):
    model = XGBRegressor(n_estimators=max_rounds,
                         learning_rate=learning_rate, max_depth=3,
                         random_state=SEED, n_jobs=-1, verbosity=0,
                         early_stopping_rounds=patience, eval_metric="rmse")
    model.fit(X_tr, y_tr, eval_set=[(X_va, y_va)], verbose=False)
    return model, int(model.best_iteration), rmse(y_va, model.predict(X_va))

def top_features(importances, k=3):
    return list(importances.sort_values(ascending=False).index[:k])

def importance_overlap(imp_a, imp_b, k=3):
    return len(set(top_features(imp_a, k)) & set(top_features(imp_b, k)))

Line-by-line explanation

NoteWhat each part does

• make_data — one default_rng(seed) generator produces both the feature matrix and the noise, so a single seed reproduces everything. The signal mixes a sine (smooth nonlinearity), a product f2*f3 (pure interaction — invisible to any single-feature look), a square, and one linear term; f6-f8 never enter, giving the importance comparison a ground truth.
• split_data — permutes indices with an independent generator and slices the first 20% as validation. Returning views via .iloc keeps the original row indices, which the tests use to verify the split is disjoint.
• rmse — coerces both inputs to float arrays so it accepts Series or arrays interchangeably; returns a plain float for clean printing.
• rf_validation_rmse — deliberately touches nothing but n_estimators: the point of the exercise is what defaults buy you. n_jobs=-1 because bagging is embarrassingly parallel.
• xgb_validation_rmse — max_depth=3 gives classic shallow weak learners; the learning rate is a parameter because the whole lesson is comparing two of them under a shared round budget.
• validation_curve_table — one row per ensemble size; refitting from scratch per size keeps the code obvious (see Design decisions for the faster alternative).
• fit_xgb_early_stopping — early_stopping_rounds goes in the constructor (post-2.0 API); fit(..., eval_set=[(X_va, y_va)]) supplies the monitored set. best_iteration is the round with the lowest validation RMSE, and subsequent predict calls use it automatically.
• top_features / importance_overlap — sort a named Series descending, take the first k names, and count the intersection of the two top-k sets: a crude but honest agreement measure between impurity and gain rankings.

Design decisions

• One fixed split, not CV. The lesson is about shapes of curves, not precise error bars; a single 80/20 split keeps every fit comparable and the runtime in seconds. Extension 1 adds CV.
• RF at pure defaults vs XGBoost at two learning rates. That asymmetry is the point: the forest gets no tuning on purpose, so the table shows its worst case against boosting’s sensitivity.
• max_depth=3 for XGBoost. Shallow trees make the learning-rate/rounds trade-off vivid; the depth-6 default would blur it by adding a second moving part.
• Refit per ensemble size instead of one fit plus predict(iteration_range=...) (XGBoost) or warm_start (RF): O(sizes) redundant work, but each table cell is exactly “this model at this size”, with no shared state to reason about.
• Early stopping at lr=0.15, budget 500, patience 30. Chosen so the demo genuinely stops early (~iteration 400) rather than exhausting the budget, which is the behavior the section is teaching.

Common mistakes

• Passing early_stopping_rounds to .fit() — that was the pre-2.0 API and now raises a TypeError; it belongs in the XGBRegressor constructor.
• Comparing the models at a single n_estimators value. At 25 trees XGBoost lr=0.05 looks hopeless and at 400 it looks fine; only the curve tells the real story.
• Judging boosting rounds on training RMSE, which decreases monotonically by construction — always monitor a held-out set.
• Forgetting to fix random_state on the models as well as the data, then chasing phantom differences between reruns.
• Expecting the interaction features f2 and f3 to top both importance lists: the interaction’s credit is split between them (and shared with the strong main effects), so rankings differ between models — that is a property of importances, not a bug in your code.

Alternative solution

NoteAnother way to do it

sklearn alone can do both sides: HistGradientBoostingRegressor(max_iter=500, learning_rate=0.15, early_stopping=True, validation_fraction=0.2, n_iter_no_change=30) is a fast LightGBM-style booster with built-in early stopping — no explicit eval_set plumbing — and permutation_importance gives a model-agnostic importance you can apply identically to both models, which is a cleaner comparison than impurity vs gain. Prefer this stack when you want to stay inside sklearn pipelines or avoid the xgboost dependency; prefer native XGBoost when you need its regularization knobs (min_child_weight, subsample, lambda) or exact control over the monitored metric.

Extensions and harder variants

• Replace the single split with 5-fold CV and add standard errors to the table; check whether the RF-vs-XGB gap survives the noise.
• Sweep learning rate over {0.01, 0.05, 0.1, 0.3} with early stopping and plot best-iteration vs learning rate — verify the roughly inverse relationship.
• Add a strongly correlated duplicate of f1 and watch what happens to both importance rankings (importance splitting across correlated features).
• Raise the noise sd from 1.5 to 4 and re-run: does boosting’s edge over RF shrink when the signal-to-noise ratio drops?
• Compare permutation_importance on the validation set against the built-in importances for both models; explain any rank changes.

Review questions

Answer out loud before expanding.

CautionQ1. Why does the RF column barely change from 25 to 400 trees while the XGBoost columns change a lot?

RF trees are i.i.d. draws averaged together: after a few dozen trees the average has essentially converged, so extra trees only shave a little variance. In boosting, each round adds capacity — the ensemble is still moving toward (or past) the optimum — so validation RMSE depends strongly on the round count until early stopping or overfitting kicks in.

CautionQ2. Why must early_stopping_rounds be paired with an eval_set, and what does best_iteration mean?

Early stopping needs an out-of-sample signal to know when improvement has stalled; without a monitored set there is nothing to stop on (training loss never stalls). best_iteration is the 0-based round index with the best monitored metric; predict uses exactly that many rounds afterwards, so the final model is the best validation-set model, not the last one trained.

CautionQ3. XGBoost at lr=0.05 with 25 rounds has RMSE near 4.7 — is that overfitting or underfitting, and how do you tell?

Underfitting: 25 steps of size 0.05 have moved the model only a fraction of the way down the loss, so it is still close to predicting the mean (baseline RMSE about 6). The tell is that both training and validation error are high, and validation error keeps improving as rounds are added — an overfit model would show low training error and worsening validation error.

CautionQ4. The validation RMSE reported at the early-stopping best iteration is slightly optimistic. Why?

The validation set was used to choose the number of rounds, so the reported minimum is selected over ~500 correlated tries — a winner’s-curse effect. It is a fine model-selection signal but a biased estimate of generalization error; honest reporting needs a third untouched test set or nested CV.

CautionQ5. Both models put some importance on the pure-noise features f6-f8. Why is it not exactly zero, and which model gives them less?

Trees will occasionally split on noise features by chance — any split that happens to reduce impurity in-sample earns importance credit, and deep RF trees (grown to purity) take many such spurious splits. XGBoost’s shallow, regularized trees have only 7 splits each to spend and gain-weighting concentrates credit on real signal, so the noise share is typically smaller (here roughly 7% vs 11% for RF).

Personal reflection

What was hard? What would you do differently? Date your entries.

References and related practices

• concepts — shared background for this practice
• r.qmd — the same exercise with ranger + xgboost’s xgb.train
• ml-001 (train/validation split) — prerequisite; the split logic here is the same idea
• ml-002 (as-of price join) — sibling demo practice in this lab
• Hastie, Tibshirani, Friedman, The Elements of Statistical Learning, ch. 10 (boosting) and ch. 15 (random forests)
• XGBoost docs: “Notes on Parameter Tuning” and the early-stopping section of the Python API reference