# take the first five cars from the mtcars data frame
cars5 <- head(mtcars, 5)
# stack the two inputs as columns, one row per car
X_cars5 <- as.matrix(cars5[, c("hp", "wt")])
# stack the outcomes into one vector
y_cars5 <- cars5$mpg
# show the input matrix and the outcome vector side by side
cbind(X_cars5, mpg = y_cars5)
#> hp wt mpg
#> Mazda RX4 110 2.620 21.0
#> Mazda RX4 Wag 110 2.875 21.0
#> Datsun 710 93 2.320 22.8
#> Hornet 4 Drive 110 3.215 21.4
#> Hornet Sportabout 175 3.440 18.7Validating Predictions
Last modified: 2026-10-27 23:24:14 (PDT)
Example 26 (The expected prediction of a line on cylinders) Treat the 392 cars of the Auto data [6] as the whole world of cars. A training set is 30 cars drawn at random from it, and the learning algorithm fits the line of mpg on cylinders. Draw 2000 training sets, fit a line to each, and record each line’s intercept, slope and prediction at 4 cylinders. Both tabs draw with the SyncRNG package [10], a random number generator with the same algorithm in R and Python, from the same seed, so they pick the same cars:
[R code]
import numpy as np
import pandas as pd
from SyncRNG import SyncRNG
auto_af = pd.read_csv("data/auto.csv").dropna(subset=["horsepower"])
cyl_af = auto_af["cylinders"].to_numpy(float)
mpg_af = auto_af["mpg"].to_numpy(float)
rng_af = SyncRNG(seed=1)
coefs_af = np.empty((2000, 2))
for s_af in range(2000):
# One training set: the first 30 rows of a shuffle of all the rows (Python counts from 0), and the line fitted to it.
pick_af = rng_af.shuffle(list(range(len(mpg_af))))[:30]
slope_af, intercept_af = np.polyfit(cyl_af[pick_af], mpg_af[pick_af], 1)
coefs_af[s_af] = intercept_af, slope_af
preds_af = coefs_af[:, 0] + coefs_af[:, 1] * 4
# The expected prediction at 4 cylinders is the mean of those predictions.
# The Monte Carlo standard error of that mean is the spread of the predictions over sqrt(2000).
mc_se_af = preds_af.std(ddof=1) / np.sqrt(2000)
print(preds_af.mean().round(2), preds_af.min().round(2), preds_af.max().round(2), mc_se_af.round(3))
#> 28.72 23.89 32.74 0.031[R code]
auto_af <- read.csv("data/auto.csv")
auto_af <- na.omit(auto_af[, c("mpg", "cylinders", "horsepower")])
# SyncRNG is a random number generator with the same algorithm
# in R and Python, so the draws match the Python tab.
library(SyncRNG)
rng_af <- SyncRNG(seed = 1)
coefs_af <- t(replicate(2000, {
# One training set: the first 30 rows of a shuffle of all the rows,
# and the line fitted to it.
rows <- rng_af$shuffle(seq_len(nrow(auto_af)))[1:30]
coef(lm(mpg ~ cylinders, data = auto_af[rows, ]))
}))
preds_af <- coefs_af[, 1] + coefs_af[, 2] * 4
# The expected prediction at 4 cylinders is the mean of those predictions.
# The Monte Carlo standard error of that mean is the spread of the
# predictions over sqrt(2000).
mc_se_af <- sd(preds_af) / sqrt(2000)
summary_af <- c(
mean = mean(preds_af),
min = min(preds_af),
max = max(preds_af)
)
c(round(summary_af, 2), mc_se = round(mc_se_af, 3))
#> mean min max mc_se
#> 28.720 23.890 32.740 0.031No single training set produces the expected prediction. The 2000 lines predict anywhere from 23.9 to 32.7 miles per gallon for a four-cylinder car, and their mean, 28.7, is the expected prediction \({\mu_{h}\mathopen{}\left(4\right)\mathclose{}}\) of Equation 30, with a 95% Monte Carlo confidence interval of 28.66 to 28.78 (1.96 Monte Carlo standard errors either side). The two tabs print the same numbers up to rounding.
Figure 1 plots 100 of the 2000 lines, their pointwise average, and the true mean function:
[R code]
import numpy as np
import matplotlib.pyplot as plt
grid_af = np.arange(3, 9)
# The pointwise expected prediction is the mean of the fitted lines at each cylinder count.
avg_af = (coefs_af[:, [0]] + coefs_af[:, [1]] * grid_af).mean(axis=0)
# The true mean function is the mean mpg of the cars with each cylinder count.
truth_af = auto_af.groupby("cylinders")["mpg"].mean()
fig, ax = plt.subplots(figsize=(6, 3.5), layout="constrained")
for intercept_af, slope_af in coefs_af[:100]:
ax.plot(grid_af, intercept_af + slope_af * grid_af, color="grey", alpha=0.3, linewidth=0.8)
# One legend entry for the 100 grey lines, with a solid line to stand for them.
ax.plot([], [], color="grey", label="One fitted line (100 shown)")
ax.plot(grid_af, avg_af, color="steelblue", linewidth=2.5, label="Expected prediction")
ax.plot(truth_af.index, truth_af.values, color="darkorange", linewidth=2.5, marker="o", label="True mean function")
ax.set_xlabel("cylinders")
ax.set_ylabel("miles per gallon")
fig.legend(loc="outside lower center", ncol=3, frameon=False)
plt.show()mpg against cylinders.
[R code]
library(ggplot2)
grid_af <- 3:8
# One row per fitted line and cylinder count;
# the pointwise expected prediction is the mean over lines.
lines_af <- data.frame(
line = rep(seq_len(nrow(coefs_af)), each = length(grid_af)),
cylinders = grid_af,
fit = as.vector(t(coefs_af[, 1] + outer(coefs_af[, 2], grid_af)))
)
avg_af <- aggregate(fit ~ cylinders, data = lines_af, FUN = mean)
# The true mean function is the mean mpg of the cars with each cylinder count.
truth_af <- aggregate(mpg ~ cylinders, data = auto_af, FUN = mean)
line_name_af <- "One fitted line (100 shown)"
avg_name_af <- "Expected prediction"
truth_name_af <- "True mean function"
ggplot() +
geom_line(
data = subset(lines_af, line <= 100),
aes(cylinders, fit, group = line, colour = line_name_af),
alpha = 0.3
) +
geom_line(
data = avg_af,
aes(cylinders, fit, colour = avg_name_af),
linewidth = 1.5
) +
geom_line(
data = truth_af,
aes(cylinders, mpg, colour = truth_name_af),
linewidth = 1.5
) +
geom_point(
data = truth_af,
aes(cylinders, mpg, colour = truth_name_af),
size = 3
) +
scale_colour_manual(
values = c("grey60", "steelblue", "darkorange"),
breaks = c(line_name_af, avg_name_af, truth_name_af)
) +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
labs(x = "cylinders", y = "miles per gallon", colour = NULL) +
theme_minimal() +
theme(legend.position = "bottom")
Example 30 (Training error and generalization error by degree) Real data never reveal \(\operatorname{Err}_{\mathcal{D}_{\text{train}}}\), but a simulation can, because it can draw as many new observations as we like. Here the population is known: \(X\sim \operatorname{Uniform}(0, 1)\) and \(Y\mid X= x\sim {\operatorname{N}\mathopen{}\left(\sin(2 \pi x), \sigma^2\right)\mathclose{}}\), with \(\sigma= 0.3\). We draw one training set of \(n= 30\) observations, fit a polynomial in \(x\) of each degree from 1 to 10 by least squares, and score each fit on the training set and on \(10{,}000\) new observations. By Theorem 6, the second score has expectation \(\operatorname{Err}_{\mathcal{D}_{\text{train}}}\) given the training set, and with \(10{,}000\) observations it is a precise estimate of it.
sim_sigma <- 0.3
sim_truth <- function(x) sin(2 * pi * x)
sim_draw <- function(n) {
x <- stats::runif(n)
tibble::tibble(x = x, y = sim_truth(x) + stats::rnorm(n, sd = sim_sigma))
}
set.seed(1)
sim_train <- sim_draw(30)
sim_new <- sim_draw(10000)
sim_degrees <- 1:10
sim_errors <- tibble::tibble(degree = sim_degrees) |>
dplyr::mutate(
fit = purrr::map(degree, \(d) lm(y ~ poly(x, d), data = sim_train)),
training = purrr::map_dbl(fit, rmse, data = sim_train),
generalization = purrr::map_dbl(fit, rmse, data = sim_new)
) |>
dplyr::select(-fit)
sim_errors[R code]
sim_errors |>
tidyr::pivot_longer(
c(training, generalization),
names_to = "error",
values_to = "rmse"
) |>
ggplot2::ggplot(ggplot2::aes(degree, rmse, colour = error, shape = error)) +
ggplot2::geom_hline(yintercept = sim_sigma, linetype = "dashed") +
ggplot2::geom_line() +
ggplot2::geom_point(size = 2) +
ggplot2::scale_x_continuous(breaks = sim_degrees) +
ggplot2::scale_y_log10() +
ggplot2::labs(
x = "polynomial degree",
y = "root mean squared error (log scale)",
colour = NULL,
shape = NULL
) +
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "top")
[R code]
sim_best_generalization <-
sim_errors$degree[which.min(sim_errors$generalization)]
sim_best_training <- sim_errors$degree[which.min(sim_errors$training)]
sim_n_below <- c(
training = sum(sim_errors$training < sim_sigma),
generalization = sum(sim_errors$generalization < sim_sigma)
)
c(
best_by_training = sim_best_training,
best_by_generalization = sim_best_generalization,
n_below_sigma2 = sim_n_below
)
#> best_by_training best_by_generalization
#> 10 5
#> n_below_sigma2.training n_below_sigma2.generalization
#> 8 0The training error is smallest at degree 10 and never rises as the degree grows, because each polynomial family contains the one before it. The estimated generalization error is smallest at degree 5. No fitted model can have generalization pRMSE below \(\sigma= 0.3\). Yet the training sRMSE falls below \(\sigma\) at 8 of the ten degrees, while the estimated generalization error falls below it at 0 of them: the training error is optimistic.
Example 34 (Numerical example) This example (Example 34) uses R’s built-in mtcars dataset (\(n=32\) cars) to predict fuel efficiency (mpg) from vehicle weight (wt). It uses one random split to compare linear, quadratic, cubic, and quartic models (mpg ~ wt, mpg ~ wt + I(wt^2), mpg ~ wt + I(wt^2) + I(wt^3), and mpg ~ wt + I(wt^2) + I(wt^3) + I(wt^4)) on a validation set, then reports the chosen model’s test sRMSE on untouched test data [6, p. 213].
set.seed(108)
n <- nrow(mtcars)
n_train <- floor(0.6 * n)
n_valid <- floor(0.2 * n)
idx_train <- sample.int(n, size = n_train)
idx_remaining <- setdiff(seq_len(n), idx_train)
idx_valid <- sample(idx_remaining, size = n_valid)
idx_test <- setdiff(idx_remaining, idx_valid)
split_sizes <- tibble::tibble(
split = c("training", "validation", "test"),
n = c(length(idx_train), length(idx_valid), length(idx_test))
)
split_sizes[R code]
train_dat <- mtcars[idx_train, ]
valid_dat <- mtcars[idx_valid, ]
test_dat <- mtcars[idx_test, ]
model_linear <- lm(mpg ~ wt, data = train_dat)
model_quadratic <- lm(mpg ~ wt + I(wt^2), data = train_dat)
model_cubic <- lm(mpg ~ wt + I(wt^2) + I(wt^3), data = train_dat)
model_quartic <- lm(
mpg ~ wt + I(wt^2) + I(wt^3) + I(wt^4),
data = train_dat
)
candidate_models <- list(
linear = model_linear,
quadratic = model_quadratic,
cubic = model_cubic,
quartic = model_quartic
)
validation_results <- tibble::tibble(
model = names(candidate_models)
) |>
dplyr::mutate(
training_RMSE = purrr::map_dbl(
model,
\(model_name) rmse(candidate_models[[model_name]], train_dat)
),
validation_RMSE = purrr::map_dbl(
model,
\(model_name) rmse(candidate_models[[model_name]], valid_dat)
)
)
model_labels <- c(
linear = "linear model",
quadratic = "quadratic model",
cubic = "cubic model",
quartic = "quartic model"
)
chosen_model_name <-
validation_results |>
dplyr::arrange(validation_RMSE) |>
dplyr::slice(1) |>
dplyr::pull(model)
chosen_model_label <- model_labels[[chosen_model_name]]
chosen_model <- candidate_models[[chosen_model_name]]
chosen_model_validation_rmse <- rmse(chosen_model, valid_dat)
performance_comparison <- tibble::tibble(
split = c("validation", "test"),
model = chosen_model_name,
RMSE = c(
chosen_model_validation_rmse,
rmse(chosen_model, test_dat)
)
)
chosen_model_name
#> [1] "linear"[R code]
partition_colors <- c(
training = "#1b9e77",
validation = "#d95f02",
test = "#7570b3"
)
partition_shapes <- c(
training = 16,
validation = 17,
test = 15
)
wt_range <- range(c(train_dat$wt, valid_dat$wt, test_dat$wt))
mpg_range <- range(c(train_dat$mpg, valid_dat$mpg, test_dat$mpg))
wt_grid <- seq(wt_range[1], wt_range[2], length.out = 200)
partition_data <- dplyr::bind_rows(
training = train_dat,
validation = valid_dat,
test = test_dat,
.id = "series"
)
plot_partitions_with_fit <- function(model, fit_label) {
series_levels <- c(names(partition_colors), fit_label)
partition_data <- partition_data |>
dplyr::mutate(
series = factor(partition_data$series, levels = series_levels)
)
grid_data <- tibble::tibble(wt = wt_grid)
fitted_curve <- tibble::tibble(
wt = wt_grid,
mpg = predict(model, newdata = grid_data),
series = factor(fit_label, levels = series_levels)
)
# nolint start: object_usage_linter.
ggplot2::ggplot(
partition_data,
ggplot2::aes(
x = wt,
y = mpg,
color = series,
shape = series
)
) +
ggplot2::geom_point() +
ggplot2::geom_line(
ggplot2::aes(linetype = series),
data = fitted_curve,
linewidth = 1
) +
ggplot2::scale_color_manual(
limits = series_levels,
values = c(partition_colors, stats::setNames("blue", fit_label))
) +
ggplot2::scale_shape_manual(
limits = series_levels,
values = c(partition_shapes, stats::setNames(NA, fit_label))
) +
ggplot2::scale_linetype_manual(
limits = series_levels,
values = stats::setNames(c(0, 0, 0, 1), series_levels)
) +
ggplot2::coord_cartesian(xlim = wt_range, ylim = mpg_range) +
ggplot2::labs(color = NULL, shape = NULL, linetype = NULL) +
ggplot2::theme_minimal() +
ggplot2::theme(
legend.position = "inside",
legend.position.inside = c(0.85, 0.8)
)
# nolint end
}
rbind(wt = wt_range, mpg = mpg_range)
#> [,1] [,2]
#> wt 1.513 5.424
#> mpg 10.400 33.900
The chosen model is the linear model, because it has the lowest validation sRMSE of the four. This split uses different cars from Exercise 1, and its validation set has only 6 cars, so the model it picks can differ from the one chosen there.
Source
The validation-set lab of James et al. [6, p. 213].
Example 35 (Cross-validation versus the generalization error) This example applies \(k\)-fold cross-validation with \(k= 5\) to the simulated training set of Example 30, using only those 30 observations. The function assign_folds() assigns each observation to a fold at random, cv_squared_errors() fits the model once per fold without that fold and returns each observation’s squared prediction error, and cv_mse() averages them, which is Equation 52.
assign_folds <- function(n, k) {
sample(rep_len(seq_len(k), n))
}
cv_squared_errors <- function(formula, data, folds) {
y <- stats::model.response(stats::model.frame(formula, data))
sq_errors <- numeric(nrow(data))
for (j in unique(folds)) {
in_fold <- folds == j
fit <- lm(formula, data = data[!in_fold, ])
pred <- predict(fit, newdata = data[in_fold, ])
sq_errors[in_fold] <- (pred - y[in_fold])^2
}
sq_errors
}
cv_mse <- function(formula, data, folds) {
mean(cv_squared_errors(formula, data, folds))
}
set.seed(5)
sim_folds <- assign_folds(nrow(sim_train), k = 5)
table(fold = sim_folds)
#> fold
#> 1 2 3 4 5
#> 6 6 6 6 6The same fold assignment is used for every degree, so the degrees are compared on the same splits. The table and figure report the square root of each cross-validation estimate, so that it is in the same unit as the training sRMSE and generalization pRMSE of Example 30.
[R code]
[R code]
sim_errors_cv |>
tidyr::pivot_longer(
c(training, generalization, cv_5),
names_to = "error",
values_to = "rmse"
) |>
dplyr::mutate(
error = factor(
error,
levels = c("training", "cv_5", "generalization"),
labels = c("training", "5-fold CV", "generalization")
)
) |>
ggplot2::ggplot(ggplot2::aes(degree, rmse, colour = error, shape = error)) +
ggplot2::geom_hline(yintercept = sim_sigma, linetype = "dashed") +
ggplot2::geom_line() +
ggplot2::geom_point(size = 2) +
ggplot2::scale_x_continuous(breaks = sim_degrees) +
ggplot2::scale_y_log10() +
ggplot2::labs(
x = "polynomial degree",
y = "root mean squared error (log scale)",
colour = NULL,
shape = NULL
) +
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "top")
[R code]
sim_best_cv <- sim_errors_cv$degree[which.min(sim_errors_cv$cv_5)]
sim_n_fit <- nrow(sim_train) - max(table(sim_folds))
sim_top <- dplyr::filter(sim_errors_cv, degree == max(sim_degrees))
sim_cv_ratio_top <- sim_top$cv_5 / sim_top$generalization
c(
best_by_cv = sim_best_cv,
best_by_generalization = sim_best_generalization,
n_per_fit = sim_n_fit,
cv_over_generalization_at_top_degree = sim_cv_ratio_top
)
#> best_by_cv best_by_generalization
#> 3.00000 5.00000
#> n_per_fit cv_over_generalization_at_top_degree
#> 24.00000 4.07316The cross-validation estimate is smallest at degree 3, and the estimated generalization error at degree 5. Unlike the training error, the cross-validation estimate does not keep falling as the degree grows, so it does not simply reward higher degrees. It is not an estimate of \(\operatorname{Err}_{\mathcal{D}_{\text{train}}}\) for the fit to all 30 observations, though (Remark 29). Each of its fits uses only 24 observations. At degree 10, polynomials fitted to 24 points predict their held-out points much worse than the fit to all 30 predicts new data: the cross-validation estimate is 4.1 times the estimated generalization pRMSE.
Part of that gap is expected: by Remark 29, with \(g\) the least-squares fit of a degree-10 polynomial, the cross-validation estimate targets \({\operatorname{Err}^{g}\mathopen{}\left(24\right)\mathclose{}}\), not the expected generalization error of a fit to 30 points. A simulation of many training sets of each size shows how much the size alone matters at this degree.
[R code]
sim_generalization_draws <- function(
n_train, degree, n_reps, draw, new_data, score
) {
purrr::map_dbl(seq_len(n_reps), \(r) {
fit <- lm(y ~ poly(x, degree), data = draw(n_train))
score(fit, new_data)
})
}
sim_n_reps <- 500
set.seed(24)
sim_size_effect <- tibble::tibble(
n_train = c(sim_n_fit, nrow(sim_train))
) |>
dplyr::mutate(
draws = purrr::map(
n_train,
\(n) {
sim_generalization_draws(
n, max(sim_degrees), sim_n_reps,
draw = sim_draw, new_data = sim_new, score = prediction_mse
)
}
),
median = purrr::map_dbl(draws, stats::median),
quantile_90 = purrr::map_dbl(draws, \(v) stats::quantile(v, 0.9)),
mean = purrr::map_dbl(draws, mean),
largest = purrr::map_dbl(draws, max),
top_two_share = purrr::map_dbl(
draws,
\(v) sum(sort(v, decreasing = TRUE)[1:2]) / sum(v)
)
) |>
dplyr::select(-draws)
sim_size_effectOver 500 simulated training sets of each size, degree-10 fits to 24 points have median generalization pMSE 1.14, against 0.27 for fits to 30 points, and 90th percentiles of 207.31 and 40.65. The sample means are no guide to \({\operatorname{Err}^{g}\mathopen{}\left(24\right)\mathclose{}}\) here: the largest single value among the fits to 24 points is 18,000,000, from one fit that is wildly wrong away from the data, and the two largest values supply 97% and 51% of the two means. So fits of this degree to 24 points are much worse than fits to 30 points, and a large cross-validation estimate here is no surprise. How much of this particular gap comes from the training-set size, and how much from the chance of this one data set, a single run cannot say. The estimate also jumps around from one degree to the next, because it rests on only 30 observations rather than on \(10{,}000\) new ones.
Example 39 (The one-standard-error rule for the simulated polynomials) This example applies 10-fold cross-validation to the simulated training set of Example 30, records the sMSE in each fold, and applies Definition 67 to the ten degrees.
[R code]
set.seed(10)
sim_folds_ten <- assign_folds(nrow(sim_train), k = 10)
sim_cv_ten <- tibble::tibble(degree = sim_degrees) |>
dplyr::mutate(
fold_mse = purrr::map(
degree,
\(d) {
sq_errors <- cv_squared_errors(y ~ poly(x, d), sim_train, sim_folds_ten)
as.numeric(tapply(sq_errors, sim_folds_ten, mean))
}
),
cv = purrr::map_dbl(fold_mse, mean),
se = purrr::map_dbl(fold_mse, \(v) sd(v) / sqrt(length(v)))
) |>
dplyr::select(-fold_mse)
sim_min <- dplyr::slice_min(sim_cv_ten, cv, n = 1)
sim_one_se <-
sim_cv_ten |>
dplyr::filter(cv <= sim_min$cv + sim_min$se) |>
dplyr::slice_min(degree, n = 1)
sim_cv_tenWith 30 observations and 10 folds, each fold holds 3 observations, so all ten folds have equal size and cv is \(\widehat{\operatorname{Err}}_{\text{CV}(10)}\) (Remark 29).
[R code]
sim_cv_ten |>
ggplot2::ggplot(ggplot2::aes(degree, cv)) +
ggplot2::geom_hline(
yintercept = sim_min$cv + sim_min$se,
linetype = "dashed"
) +
ggplot2::geom_errorbar(
ggplot2::aes(ymin = cv - se, ymax = cv + se),
width = 0.2
) +
ggplot2::geom_point(size = 2) +
ggplot2::scale_x_continuous(breaks = sim_degrees) +
ggplot2::scale_y_log10() +
ggplot2::labs(
x = "polynomial degree",
y = "10-fold CV estimate (log scale)"
) +
ggplot2::theme_bw()
The smallest estimate is at degree 8. The lowest degree within one fold-based standard error of it is degree 3, which is the one-standard-error rule’s choice. The rule prefers the simpler model when the data cannot clearly tell it apart from the degree with the smallest estimate. Because the data are simulated, we can check both choices against the estimated generalization pRMSE of Example 30 for each: 0.335 for degree 8 and 0.315 for degree 3.