Examples
In this example we generate data, write \(K\)-fold cross-validation from scratch, and use it to select a polynomial degree.
import numpy as np
import matplotlib.pyplot as plt
RNG = np.random.default_rng(322)
SIGMA = 0.15
def f(x):
return np.sin(2 * x) + 2 * (x - 0.5) ** 3 - 0.5 * x
def sample(n, rng=RNG):
x = rng.uniform(0, 1, n)
return x, f(x) + rng.normal(0, SIGMA, n)
x, y = sample(80)
print(f"n = {len(x)}, x in [{x.min():.2f}, {x.max():.2f}], y in [{y.min():.2f}, {y.max():.2f}]")n = 80, x in [0.00, 0.97], y in [-0.30, 1.00]Polynomial regression is linear regression on a transformed design matrix. The model is non-linear in \(x\), but it is linear in the parameters, so least squares applies without any change.
def design(x, degree):
return np.vander(x, degree + 1, increasing=True)
def fit(x, y, degree):
coef, *_ = np.linalg.lstsq(design(x, degree), y, rcond=None)
return coef
def mse(coef, x, y):
pred = design(x, len(coef) - 1) @ coef
return float(np.mean((pred - y) ** 2))
coef = fit(x, y, 3)
print("coefficients:", np.round(coef, 3))
print("training MSE:", round(mse(coef, x, y), 4))lstsq solves the normal equations \(X^\top X\beta = X^\top y\) in a stable way. Do not form \((X^\top X)^{-1}\) yourself.
coefficients: [-0.279 3.238 -4.031 1.718]
training MSE: 0.0239The training MSE always falls when we increase the degree. So we cannot use it to pick the degree. This is the whole problem.
degrees = np.arange(1, 11)
train_err = [mse(fit(x, y, d), x, y) for d in degrees]
fig, ax = plt.subplots()
ax.plot(degrees, train_err, "o-")
ax.axhline(SIGMA**2, ls="--", lw=1.2, color="#a8452f")
ax.annotate("irreducible error $\\sigma^2$", (degrees[-1], SIGMA**2),
xytext=(-4, 6), textcoords="offset points",
ha="right", color="#a8452f", fontsize=8.5)
ax.set(xlabel="polynomial degree", ylabel="training MSE", yscale="log")
ax.set_title("Training error")
plt.show()
def kfold_indices(n, K, rng=RNG):
"""Shuffle, then cut into K blocks of nearly equal size."""
return np.array_split(rng.permutation(n), K)
def cross_val_mse(x, y, degree, K=10, rng=RNG):
folds = kfold_indices(len(x), K, rng)
total, count = 0.0, 0
for i in range(K):
val = folds[i]
train = np.concatenate([folds[j] for j in range(K) if j != i])
coef = fit(x[train], y[train], degree)
total += mse(coef, x[val], y[val]) * len(val)
count += len(val)
return total / count
print("10-fold CV MSE, degree 3:", round(cross_val_mse(x, y, 3), 4))i is held out. Everything else is used for training.
10-fold CV MSE, degree 3: 0.0261Now we run this for every degree and take the smallest.
cv_err = [cross_val_mse(x, y, d, K=10) for d in degrees]
best = degrees[int(np.argmin(cv_err))]
print(f"selected degree: {best} CV MSE: {min(cv_err):.4f}")selected degree: 4 CV MSE: 0.0248fig, ax = plt.subplots()
ax.plot(degrees, train_err, "o-", label="training MSE", alpha=0.55)
ax.plot(degrees, cv_err, "o-", label="10-fold CV MSE")
ax.axvline(best, color="#1f5f7a", lw=1, ls=":", alpha=0.8)
ax.axhline(SIGMA**2, ls="--", lw=1.2, color="#a8452f", label="$\\sigma^2$")
ax.set(xlabel="polynomial degree", ylabel="MSE", yscale="log")
ax.set_title(f"Selected degree: {best}")
ax.legend()
plt.show()
In practice we do not write the loop ourselves. Note that the model and the polynomial transformation go into one Pipeline, so that the transformation is also fitted inside each fold.
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
X = x.reshape(-1, 1)
cv = KFold(10, shuffle=True, random_state=0)
scores = []
for d in degrees:
model = make_pipeline(PolynomialFeatures(d), LinearRegression())
s = cross_val_score(model, X, y, cv=cv, scoring="neg_mean_squared_error")
scores.append(-s.mean())
print("selected degree:", degrees[int(np.argmin(scores))])selected degree: 4For a model that is linear in its parameters the leave-one-out error has a closed form. We fit once, read the diagonal of the hat matrix, and divide.
def loocv_mse(x, y, degree):
X = design(x, degree)
Q, _ = np.linalg.qr(X)
h = np.sum(Q**2, axis=1)
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
resid = y - X @ coef
return float(np.mean((resid / (1 - h)) ** 2))
loo = [loocv_mse(x, y, d) for d in degrees]
print("LOOCV selects degree:", degrees[int(np.argmin(loo))])
print("brute force check (degree 3):",
round(cross_val_mse(x, y, 3, K=len(x)), 6), "vs", round(loocv_mse(x, y, 3), 6))LOOCV selects degree: 3
brute force check (degree 3): 0.025934 vs 0.025934The two numbers agree to six digits. That is the point of the check. The identity is exact, not an approximation.
The value min(cv_err) above is not an estimate of the test error of the selected model. It is the smallest of ten correlated noisy numbers, and we chose it because it was small. See nested cross-validation.
The panel below makes the same point in a different way. Nothing changes except the shuffle.