1. The AutoML promise¶
Around 2018–2022, a wave of AutoML tools promised to automate the tedious middle of machine learning: pick the model family, tune the hyperparameters, report the winner. auto-sklearn, TPOT, H2O AutoML, and pycaret all offered a one-line compare_models() experience, and earlier editions of this chapter taught pycaret on the same temperature dataset used in lesson 3.7.
That wave receded. Most academic AutoML tools are no longer maintained, and pycaret is retired from this course. Two things survived and are worth teaching:
- Hyperparameter-optimization libraries. Automated search over model settings is still standard practice; Optuna is the current default library for it.
- Strong tabular defaults. Gradient-boosted trees (scikit-learn’s
HistGradientBoostingRegressor, LightGBM, XGBoost) win on feature tables so consistently that model search is rarely the bottleneck anymore. Pick a boosted tree, tune it a little, and spend the saved time on data quality and evaluation.
This lesson covers both survivors, then looks at what replaced AutoML in 2026: code-writing agents, and the verification skills they demand from you.
2. Hyperparameter search that survived: grid search vs Optuna¶
We reuse the exact dataset from lesson 3.7 — the same generator, the same features, the same split — so numbers are comparable across the two lessons. The model is HistGradientBoostingRegressor, and the quantity we optimize is the 5-fold cross-validated MAE on the training set.
import numpy as np
import pandas as pd
def make_daily_temps(start="2012-01-01", end="2019-12-31", seed=42):
"""Synthetic Seattle-like daily maximum temperature record (degrees F).
Seasonal climatology + a weak warming trend + AR(1) weather noise.
Generated in-notebook so the lesson does not depend on a remote file.
"""
rng = np.random.default_rng(seed)
dates = pd.date_range(start, end, freq="D")
day_of_year = dates.dayofyear.to_numpy()
climatology = 62.0 - 15.0 * np.cos(2 * np.pi * (day_of_year - 203) / 365.25)
trend = 0.05 * np.arange(len(dates)) / 365.25
noise = np.zeros(len(dates))
for i in range(1, len(dates)):
noise[i] = 0.65 * noise[i - 1] + rng.normal(0.0, 3.0)
df = pd.DataFrame(
{
"date": dates,
"average": np.round(climatology, 1), # historical average for that calendar day
"actual": np.round(climatology + trend + noise, 1),
}
)
df["temp_1"] = df["actual"].shift(1) # yesterday's max
df["temp_2"] = df["actual"].shift(2) # two days ago
df["month"] = df["date"].dt.month
df["day"] = df["date"].dt.day
df["doy_sin"] = np.sin(2 * np.pi * day_of_year / 365.25)
df["doy_cos"] = np.cos(2 * np.pi * day_of_year / 365.25)
return df.dropna().reset_index(drop=True)
df = make_daily_temps()
features = ["temp_1", "temp_2", "average", "month", "day", "doy_sin", "doy_cos"]
X = df[features]
y = df["actual"]
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
print("Training features:", X_train.shape, " Testing features:", X_test.shape)Training features: (2190, 7) Testing features: (730, 7)
Grid search¶
GridSearchCV tries every combination in a fixed lattice of values. With 3 depths, 3 learning rates, and 2 leaf-node limits, that is 18 candidates, each cross-validated 5 times: 90 fits.
import time
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error
from sklearn.model_selection import GridSearchCV, KFold, cross_val_score
cv = KFold(n_splits=5, shuffle=True, random_state=42)
param_grid = {
"max_depth": [2, 4, 8],
"learning_rate": [0.03, 0.1, 0.3],
"max_leaf_nodes": [15, 31],
}
grid = GridSearchCV(
HistGradientBoostingRegressor(random_state=42),
param_grid,
cv=cv,
scoring="neg_mean_absolute_error",
)
t0 = time.perf_counter()
grid.fit(X_train, y_train)
grid_time = time.perf_counter() - t0
n_candidates = len(grid.cv_results_["params"])
grid_fits = n_candidates * cv.get_n_splits()
grid_cv_mae = -grid.best_score_
grid_test_mae = mean_absolute_error(y_test, grid.best_estimator_.predict(X_test))
print("Best params:", grid.best_params_)
print(f"Best CV MAE: {grid_cv_mae:.3f} F")
print(f"Fits: {grid_fits} ({n_candidates} candidates x {cv.get_n_splits()} folds), wall time {grid_time:.1f} s")Best params: {'learning_rate': 0.1, 'max_depth': 2, 'max_leaf_nodes': 15}
Best CV MAE: 2.458 F
Fits: 90 (18 candidates x 5 folds), wall time 6.8 s
Optuna¶
Optuna samples the hyperparameter space instead of walking a lattice. Its TPE sampler builds a model of which regions score well and concentrates later trials there. The objective function is anything you can compute: here, the same 5-fold CV MAE.
import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)
def objective(trial):
params = {
"max_depth": trial.suggest_int("max_depth", 2, 10),
"learning_rate": trial.suggest_float("learning_rate", 0.01, 0.4, log=True),
"max_leaf_nodes": trial.suggest_int("max_leaf_nodes", 10, 60),
}
model = HistGradientBoostingRegressor(random_state=42, **params)
scores = -cross_val_score(model, X_train, y_train, cv=cv, scoring="neg_mean_absolute_error")
return scores.mean()
study = optuna.create_study(direction="minimize", sampler=optuna.samplers.TPESampler(seed=42))
t0 = time.perf_counter()
study.optimize(objective, n_trials=30)
optuna_time = time.perf_counter() - t0
optuna_fits = 30 * cv.get_n_splits()
optuna_cv_mae = study.best_value
print("Best params:", study.best_params)
print(f"Best CV MAE: {optuna_cv_mae:.3f} F")
print(f"Fits: {optuna_fits} (30 trials x {cv.get_n_splits()} folds), wall time {optuna_time:.1f} s")Best params: {'max_depth': 3, 'learning_rate': 0.05958491008634781, 'max_leaf_nodes': 35}
Best CV MAE: 2.438 F
Fits: 150 (30 trials x 5 folds), wall time 12.6 s
# Refit each winner on the full training set, score once on the test set
optuna_model = HistGradientBoostingRegressor(random_state=42, **study.best_params)
optuna_model.fit(X_train, y_train)
optuna_test_mae = mean_absolute_error(y_test, optuna_model.predict(X_test))
comparison = pd.DataFrame(
{
"best CV MAE (F)": [grid_cv_mae, optuna_cv_mae],
"test MAE (F)": [grid_test_mae, optuna_test_mae],
"wall time (s)": [grid_time, optuna_time],
"fits": [grid_fits, optuna_fits],
},
index=["GridSearchCV", "Optuna (TPE)"],
)
comparison.round(3)Grid search spends its budget on a fixed lattice: every point gets tried whether or not its neighborhood already looked bad. Optuna samples the continuous space and adapts to past trials, so it can land between lattice points and skip dead regions. On a problem this small the difference is minutes; on a deep-learning search with hour-long fits, it is days. Either way, both methods report a cross-validated score on training data and touch the test set exactly once, at the end.
3. 2026: agents write the pipeline, you verify it¶
Model exploration today often happens conversationally: you describe the dataset to an LLM agent, and it writes the pipeline code, runs it, and reports a score. The search problem AutoML tried to solve became a verification problem. Machine-written code fails in the same ways as hurried human code — leaked preprocessing, dropped classes, scores computed on the wrong split — only more fluently, wrapped in tidy comments and confident output.
Below is a modeling script of the kind an AI assistant will happily produce. It runs, it prints a strong score, and it is wrong in three distinct ways. Find them before opening the solution.
# --- AI-generated modeling script: do NOT trust it yet ---
import mlgeo_synth
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
# Whole-rock geochemistry: oxide wt%, density, magnetic susceptibility -> rock type
df = mlgeo_synth.geochem_table(n=4000, seed=7)
df = df[df["label"].isin(["granite", "basalt"])] # remove sparse label noise
feature_cols = ["SIO2", "AL2O3", "FEO", "MGO", "CAO", "NA2O", "K2O", "density_g_cm3", "mag_susc_si"]
X = df[feature_cols].to_numpy()
y = df["label"].to_numpy()
scaler = StandardScaler() # normalize features
X = scaler.fit_transform(X)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=7)
models = {
"logistic_regression": LogisticRegression(max_iter=2000),
"random_forest": RandomForestClassifier(n_estimators=200, random_state=7),
"hist_gradient_boosting": HistGradientBoostingClassifier(random_state=7),
}
best_name, best_score = None, -1.0
for name, clf in models.items():
clf.fit(X_train, y_train)
score = clf.score(X_train, y_train) # evaluate each candidate
print(f"{name}: accuracy = {score:.3f}")
if score > best_score:
best_name, best_score = name, score
print(f"\nbest model: {best_name}, accuracy = {best_score:.3f}")logistic_regression: accuracy = 1.000
random_forest: accuracy = 1.000
hist_gradient_boosting: accuracy = 1.000
best model: logistic_regression, accuracy = 1.000
Solution
- The minority class is silently dropped.
df[df["label"].isin(["granite", "basalt"])]discards every andesite sample under a comment about “label noise”. A 3-class problem becomes an easier 2-class one, the reported accuracy applies to a different task than the one posed, and in production this would be a silent scientific error: the model can never predict andesite. - The scaler is fit before the split.
StandardScaler().fit_transform(X)on the full matrix computes means and variances using the test rows, then the split happens. Test-set statistics leak into the training features. The effect is small here, but the pattern is exactly the leakage discussed in earlier lessons, and with other preprocessors (imputation, target encoding) it can be large. - Model selection uses training-set accuracy.
clf.score(X_train, y_train)rewards memorization, so this comparison systematically favors the most overfit candidate. Here every model scores ~1.0 on data it has already seen, so the comparison cannot distinguish them at all. The printed “best model accuracy” says nothing about generalization.
The tell-tale smell: a near-perfect score appearing with no baseline and no held-out evaluation. Real results come with a baseline to beat and a test set touched once.
Corrected pipeline¶
The fixes: keep all three classes; split first, with stratification so the andesite minority appears in both subsets; put the scaler inside a Pipeline so it is fit only on training folds; report a majority-class baseline before any model; select among candidates by cross-validated macro-F1 on the training set (macro-F1 weights the minority class equally); and touch the test set once, at the end, with per-class scores.
from sklearn.dummy import DummyClassifier
from sklearn.metrics import accuracy_score, classification_report, f1_score
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
df = mlgeo_synth.geochem_table(n=4000, seed=7) # all three classes kept
print(df["label"].value_counts(), "\n")
X = df[feature_cols]
y = df["label"]
# Split FIRST, stratified so class proportions match in train and test
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.25, random_state=7, stratify=y
)
# Baseline before any model
dummy = DummyClassifier(strategy="most_frequent").fit(X_train, y_train)
print(f"Majority-class baseline: accuracy = {dummy.score(X_test, y_test):.3f}, "
f"macro-F1 = {f1_score(y_test, dummy.predict(X_test), average='macro'):.3f}\n")
candidates = {
"logistic_regression": make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)),
"random_forest": make_pipeline(StandardScaler(), RandomForestClassifier(n_estimators=200, random_state=7)),
"hist_gradient_boosting": make_pipeline(StandardScaler(), HistGradientBoostingClassifier(random_state=7)),
}
best_name, best_cv = None, -1.0
for name, pipe in candidates.items():
scores = cross_val_score(pipe, X_train, y_train, cv=5, scoring="f1_macro")
print(f"{name}: CV macro-F1 = {scores.mean():.3f} +/- {scores.std():.3f}")
if scores.mean() > best_cv:
best_name, best_cv = name, scores.mean()
# One look at the test set, for the selected model only
final = candidates[best_name].fit(X_train, y_train)
y_pred = final.predict(X_test)
print(f"\nSelected model: {best_name}")
print(f"Test accuracy = {accuracy_score(y_test, y_pred):.3f}, "
f"test macro-F1 = {f1_score(y_test, y_pred, average='macro'):.3f}\n")
print(classification_report(y_test, y_pred))label
granite 2205
basalt 1397
andesite 398
Name: count, dtype: int64
Majority-class baseline: accuracy = 0.551, macro-F1 = 0.237
logistic_regression: CV macro-F1 = 0.999 +/- 0.002
random_forest: CV macro-F1 = 0.999 +/- 0.003
hist_gradient_boosting: CV macro-F1 = 0.999 +/- 0.001
Selected model: hist_gradient_boosting
Test accuracy = 0.999, test macro-F1 = 0.998
precision recall f1-score support
andesite 1.00 0.99 0.99 100
basalt 1.00 1.00 1.00 349
granite 1.00 1.00 1.00 551
accuracy 1.00 1000
macro avg 1.00 1.00 1.00 1000
weighted avg 1.00 1.00 1.00 1000
The corrected pipeline also reports a high score — these rock types are genuinely well separated in oxide space — but the number now makes a different claim. It covers all three classes, including the andesite the AI script could never predict; it is measured on held-out data instead of memorized data; and it stands against a 0.55 majority baseline. Same digits, different meaning. On a harder dataset the two workflows diverge: the AI script’s training score stays near 1.0 no matter what, while the honest number drops to tell you so.
4. What to keep¶
Automation moved. In 2020 it lived in search algorithms — AutoML looping over models and hyperparameters. In 2026 it lives in code-writing agents that produce the whole pipeline on request. Verification did not move. A trivial baseline, a leak-free split, cross-validation inside the training set, and one look at a held-out test set catch machine-written mistakes exactly as they catch human ones. Automate the search; never the checking.