A model can fail in two opposite ways: by being too simple to capture the pattern, or by being complex enough to memorise the noise. Telling them apart is the single most useful diagnostic skill in machine learning, and it needs two numbers rather than one.
The two failures
import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
rng = np.random.default_rng(3)
n = 60
# The truth is a gentle curve. The data has noise on top.
x = np.sort(rng.uniform(0, 1, n))
truth = np.sin(2 * np.pi * x)
y = truth + rng.normal(0, 0.28, n)
x_tr, x_te, y_tr, y_te = train_test_split(x, y, test_size=0.4, random_state=0)
rows = []
for degree in [0, 1, 3, 9, 15, 25]:
coefficients = np.polyfit(x_tr, y_tr, degree)
rows.append({
'degree': degree,
'train_rmse': round(float(np.sqrt(mean_squared_error(
y_tr, np.polyval(coefficients, x_tr)))), 3),
'test_rmse': round(float(np.sqrt(mean_squared_error(
y_te, np.polyval(coefficients, x_te)))), 3),
})
table = pd.DataFrame(rows)
table['gap'] = (table['test_rmse'] - table['train_rmse']).round(3)
print(table.to_string(index=False))
print()
print(f'the noise floor is {0.28:.2f} - no model can beat it')
print()
print('degree 0-1: both errors high -> UNDERFITTING')
print('degree 3: both near the floor -> about right')
print('degree 15+: train tiny, test huge -> OVERFITTING')
Open this example in the editor and change it. Nothing you type is saved.
The diagnosis comes from the pair of numbers, never from one:
- Both high - underfitting. The model is not flexible enough, or the features do not contain the signal.
- Train low, test high - overfitting. The model has memorised the training set.
- Both low - you are done.
- Train high, test low - almost always a bug in the split, or a tiny test set.
Reporting only the test score cannot distinguish "this problem is hard" from "this model memorised".
Cross-validation
A single train/test split gives one number, and that number is noisy. Cross-validation gives several.
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import (cross_val_score, StratifiedKFold,
train_test_split)
from sklearn.metrics import accuracy_score
rng = np.random.default_rng(7)
n = 600
X = rng.normal(size=(n, 6))
y = (X[:, 0] + 0.8 * X[:, 1] + rng.normal(0, 1.1, n) > 0).astype(int)
# Five different single splits of the same data.
singles = []
for seed in range(5):
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.25, random_state=seed, stratify=y)
model = RandomForestClassifier(n_estimators=150, random_state=0).fit(X_tr, y_tr)
singles.append(accuracy_score(y_te, model.predict(X_te)))
print('five single splits of the SAME data:')
print(' ', [round(s, 3) for s in singles])
print(f' range {min(singles):.3f} to {max(singles):.3f} - a spread of '
f'{max(singles) - min(singles):.3f}')
print()
scores = cross_val_score(
RandomForestClassifier(n_estimators=150, random_state=0), X, y,
cv=StratifiedKFold(5, shuffle=True, random_state=0))
print(f'5-fold cross-validation: {scores.round(3)}')
print(f' mean {scores.mean():.3f} std {scores.std():.3f}')
print(f' report as {scores.mean():.3f} +/- {2 * scores.std():.3f}')
Open this example in the editor and change it. Nothing you type is saved.
The spread between single splits is the thing to notice. A model reported as 0.84 when a different random seed would have said 0.79 is not reporting a property of the model.
Report the mean and the spread, and when comparing two models, compare them on the same folds:
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(11)
n = 2000
X = rng.normal(size=(n, 8))
y = ((X[:, 0] + 0.9 * (X[:, 1] > 0) * (X[:, 2] > 0) + rng.normal(0, 1, n)) > 0.5).astype(int)
folds = StratifiedKFold(5, shuffle=True, random_state=0)
models = {
'logistic': Pipeline([('s', StandardScaler()), ('m', LogisticRegression(max_iter=1000))]),
'forest': RandomForestClassifier(n_estimators=200, random_state=0),
'boosting': HistGradientBoostingClassifier(random_state=0),
}
rows = []
for name, model in models.items():
scores = cross_val_score(model, X, y, cv=folds, scoring='roc_auc')
rows.append({'model': name, 'mean_auc': round(scores.mean(), 4),
'std': round(scores.std(), 4),
'folds': ' '.join(f'{s:.3f}' for s in scores)})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('The same folds for every model, so the comparison is paired')
print('and the fold-to-fold variation cancels.')
Open this example in the editor and change it. Nothing you type is saved.
Which splitter to use is the same question as in the split lesson: StratifiedKFold for a classification label, GroupKFold when rows share an entity, TimeSeriesSplit when the data is ordered in time.
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
X = np.arange(20).reshape(-1, 1)
for fold, (train, test) in enumerate(TimeSeriesSplit(n_splits=4).split(X), 1):
print(f'fold {fold}: train {train.min():>2}-{train.max():>2} '
f'test {test.min():>2}-{test.max():>2}')
print()
print('Each fold trains only on the past. A shuffled KFold would let')
print('the model see the future, and the score would be a fiction.')
Open this example in the editor and change it. Nothing you type is saved.
The learning curve
A learning curve answers a different question: would more data help?
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import learning_curve, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(13)
n = 4000
X = rng.normal(size=(n, 10))
y = ((1.2 * X[:, 0] + 0.9 * (X[:, 1] > 0) * (X[:, 2] > 0)
+ rng.normal(0, 1.0, n)) > 0.6).astype(int)
for name, model in [
('underfitting (1 feature)',
Pipeline([('s', StandardScaler()), ('m', LogisticRegression(max_iter=1000))])),
('overfitting (deep forest)',
RandomForestClassifier(n_estimators=100, random_state=0)),
]:
data = X[:, :1] if 'underfitting' in name else X
sizes, train_scores, test_scores = learning_curve(
model, data, y, train_sizes=np.linspace(0.05, 1.0, 6),
cv=StratifiedKFold(4, shuffle=True, random_state=0), scoring='roc_auc')
table = pd.DataFrame({
'n': sizes,
'train': train_scores.mean(axis=1).round(3),
'cv': test_scores.mean(axis=1).round(3),
})
table['gap'] = (table['train'] - table['cv']).round(3)
print(name)
print(table.to_string(index=False))
print()
Open this example in the editor and change it. Nothing you type is saved.
How to read the two shapes:
- Both curves flat and close together, at a low level - underfitting. More data will not help; a more flexible model or better features will.
- A wide and slowly closing gap - overfitting. More data will help, and so will regularisation.
- Both curves high and converged - you are at the limit of what this data supports.
That first case is the valuable one, because "we need more data" is an expensive conclusion and the learning curve is how you find out whether it is true.
Regularisation, and the bias-variance trade
import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score, KFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
rng = np.random.default_rng(17)
n = 120
x = rng.uniform(0, 1, n)
y = np.sin(2 * np.pi * x) + rng.normal(0, 0.3, n)
X = x.reshape(-1, 1)
rows = []
for alpha in [1e-6, 1e-3, 1e-1, 1.0, 10.0, 100.0]:
model = Pipeline([
('poly', PolynomialFeatures(degree=15)),
('scale', StandardScaler()),
('model', Ridge(alpha=alpha)),
])
scores = -cross_val_score(model, X, y, cv=KFold(5, shuffle=True, random_state=0),
scoring='neg_root_mean_squared_error')
fitted = model.fit(X, y)
train_rmse = float(np.sqrt(np.mean((y - fitted.predict(X)) ** 2)))
rows.append({'alpha': alpha, 'train_rmse': round(train_rmse, 3),
'cv_rmse': round(scores.mean(), 3)})
table = pd.DataFrame(rows)
table['gap'] = (table['cv_rmse'] - table['train_rmse']).round(3)
print(table.to_string(index=False))
print()
print(f'noise floor {0.30:.2f}')
print('Low alpha: a degree-15 polynomial fits the noise.')
print('High alpha: the model is pulled flat and underfits.')
print('The best alpha sits between, near the noise floor.')
Open this example in the editor and change it. Nothing you type is saved.
Regularisation is a dial between the two failures. Every model family has one:
- Linear: the
alphaof ridge or lasso. - Trees:
max_depth,min_samples_leaf,ccp_alpha. - Boosting:
learning_ratewithmax_iter, plus early stopping. - Neural networks: weight decay, dropout, early stopping.
More data moves the whole curve down; regularisation moves you along it.
A worked example
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import (StratifiedKFold, cross_validate,
learning_curve, train_test_split)
from sklearn.metrics import roc_auc_score
rng = np.random.default_rng(seed=71)
n = 3000
# A moderately hard problem: real signal, real noise, an interaction.
X = rng.normal(size=(n, 12))
signal = (1.1 * X[:, 0] - 0.7 * X[:, 1] + 1.2 * (X[:, 2] > 0) * (X[:, 3] > 0)
+ 0.6 * np.abs(X[:, 4]))
y = (signal + rng.normal(0, 1.3, n) > 1.0).astype(int)
X_rest, X_test, y_rest, y_test = train_test_split(
X, y, test_size=0.2, random_state=0, stratify=y)
folds = StratifiedKFold(5, shuffle=True, random_state=0)
print(f'{len(X_rest):,} for development, {len(X_test):,} held out')
print(f'positive rate {y.mean():.1%}')
print()
# ---- 1. Sweep the regularisation dial, measuring BOTH scores. ------
rows = []
for leaves, min_leaf in [(3, 50), (7, 25), (15, 10), (31, 5), (63, 1), (127, 1)]:
model = HistGradientBoostingClassifier(
max_leaf_nodes=leaves, min_samples_leaf=min_leaf,
learning_rate=0.1, max_iter=200, early_stopping=False, random_state=0)
result = cross_validate(model, X_rest, y_rest, cv=folds, scoring='roc_auc',
return_train_score=True)
rows.append({
'max_leaf_nodes': leaves, 'min_samples_leaf': min_leaf,
'train_auc': round(result['train_score'].mean(), 3),
'cv_auc': round(result['test_score'].mean(), 3),
'cv_std': round(result['test_score'].std(), 3),
})
table = pd.DataFrame(rows)
table['gap'] = (table['train_auc'] - table['cv_auc']).round(3)
print(table.to_string(index=False))
print()
best = table.loc[table['cv_auc'].idxmax()]
print(f'best cross-validated AUC {best["cv_auc"]:.3f} at '
f'max_leaf_nodes={int(best["max_leaf_nodes"])}')
print(f'the largest model trains to {table["train_auc"].max():.3f} and '
f'cross-validates at {table.loc[table["train_auc"].idxmax(), "cv_auc"]:.3f} '
f'- a gap of {table["gap"].max():.3f}')
print()
# ---- 2. The learning curve: would more data help? ------------------
model = HistGradientBoostingClassifier(
max_leaf_nodes=int(best['max_leaf_nodes']),
min_samples_leaf=int(best['min_samples_leaf']),
learning_rate=0.1, max_iter=200, early_stopping=False, random_state=0)
sizes, train_scores, cv_scores = learning_curve(
model, X_rest, y_rest, train_sizes=np.linspace(0.1, 1.0, 6),
cv=folds, scoring='roc_auc')
curve = pd.DataFrame({
'train_rows': sizes.astype(int),
'train_auc': train_scores.mean(axis=1).round(3),
'cv_auc': cv_scores.mean(axis=1).round(3),
})
curve['gap'] = (curve['train_auc'] - curve['cv_auc']).round(3)
print(curve.to_string(index=False))
print()
improvement = curve['cv_auc'].iloc[-1] - curve['cv_auc'].iloc[-2]
print(f'the last doubling of data bought {improvement:+.3f} AUC')
print(f'the gap is {"still closing" if curve["gap"].iloc[-1] < curve["gap"].iloc[0] else "not closing"}, '
f'{curve["gap"].iloc[0]:.3f} -> {curve["gap"].iloc[-1]:.3f}')
print()
if improvement > 0.005:
print('More data would still help: the curve has not flattened.')
else:
print('More data would help little: the curve has flattened, and the')
print('remaining error is noise or a missing feature.')
print()
# ---- 3. Early stopping, instead of guessing the iteration count. ---
with_stopping = HistGradientBoostingClassifier(
max_leaf_nodes=int(best['max_leaf_nodes']),
min_samples_leaf=int(best['min_samples_leaf']),
learning_rate=0.05, max_iter=2000, early_stopping=True,
validation_fraction=0.15, n_iter_no_change=30, random_state=0).fit(X_rest, y_rest)
print(f'early stopping used {with_stopping.n_iter_} of 2000 iterations')
print()
# ---- 4. The test set, once. ---------------------------------------
final = HistGradientBoostingClassifier(
max_leaf_nodes=int(best['max_leaf_nodes']),
min_samples_leaf=int(best['min_samples_leaf']),
learning_rate=0.05, max_iter=with_stopping.n_iter_,
early_stopping=False, random_state=0).fit(X_rest, y_rest)
train_auc = roc_auc_score(y_rest, final.predict_proba(X_rest)[:, 1])
test_auc = roc_auc_score(y_test, final.predict_proba(X_test)[:, 1])
print(f'train AUC {train_auc:.3f}')
print(f'TEST AUC {test_auc:.3f} (cross-validated estimate was {best["cv_auc"]:.3f})')
print(f'gap {train_auc - test_auc:.3f}')
print()
print('The cross-validated estimate and the held-out test score agree')
print('to within the fold-to-fold spread, which is what you want: it')
print('means the selection process did not overfit the folds.')
Open this example in the editor and change it. Nothing you type is saved.
Three outputs there answer three different questions. The sweep shows where the model stops generalising, and the gap column shows it directly. The learning curve says whether collecting more data is worth doing. And the agreement between the cross-validated estimate and the held-out score is evidence that the tuning itself did not overfit - which is the thing a single number cannot tell you.
When it goes wrong
| Symptom | Cause |
|---|---|
| Only the test score is reported | The gap is the diagnosis; report both |
| Both scores are poor | Underfitting, or the features lack the signal |
| Training score is near perfect | Memorisation; regularise or simplify |
| Scores swing with the random seed | A single split; use cross-validation |
| Cross-validation beats the test set badly | The folds were used for selection too |
| More data did not help | The curve had flattened; it was the model or the features |
| Cross-validation on a time series looks great | Use TimeSeriesSplit |
| The gap is small but the score is low | The problem may simply be hard |
A check you can run
Take your current model and print its training score alongside its test score.
If the training score is near perfect, the model is memorising and the test score will fall further in production. If both are poor, no amount of tuning will help and the next step is features, not hyperparameters. Those two numbers take one extra line and they tell you which of two completely different problems you have.