Lab 4 Exercise - Cross Validation

Author

Josh Davis/Lucas Schwengber

Learning goals

By the end of this lab you should be able to:

  1. Explain the fundamental reason why cross validation is needed, and why it almost always works better than using training error as a model selection technique.
  2. Explain why training error is not unbiased for the risk. It is too optimistic, because the same data is used to fit and to evaluate the model, and this gets worse as model complexity grows.
  3. State what CV error estimates. If implemented correctly, CV error is unbiased for the expected risk of a model trained on \(N - N/K\) points (under the condition of no distribution shift). Because the final model is trained on all \(N\) points, CV error slightly overestimates its risk.
  4. State the principles required for CV to actually work: the validation data must not be used in training, and it must be distributed like the data the model will be used on. For IID data, a random train-val split achieves this.
  5. Describe how the sizes of the training and validation sets affect the bias and variance of the error estimate.
  6. Describe how the number of folds \(K\) affects the bias and variance of the risk estimate, and the computational burden. Larger \(K\) lowers bias and increases computational burden. The effect on variance is not monotone, which is why \(K = 5\) or \(10\) is common.
Part Learning goals
Setup
1. Training error vs. true risk 1, 2
2. What does CV estimate? 3
3. When CV breaks 4
4. Size of the validation set 5
5. Number of folds 6
6. Challenge 1, 3
ImportantHow to read this lab

Everything you have to do is in a red TODO box like this one. Text outside the boxes is background you need to understand the TODOs.

Each TODO says in words what to do, followed by functions that may be useful. Each TODO has a number (TODO 1.1, TODO 1.2, …). In your notebook, put each TODO under a markdown heading with its number, for example ### TODO 1.2, and write your answers to the questions in a markdown cell under their TODO.

Setup

In Parts 1-5 you will work in a simulated “practice world”: a binary classification problem in which we know the true probability \(p(\boldsymbol{x}) = P(y = 1 \mid \boldsymbol{x})\). This means that for any fitted model \(\hat f\) that predicts probabilities \(q(\boldsymbol{x})\), we can compute its true risk (expected log loss on a new observation) exactly, \[ \text{Err}_\mathcal{T} = \mathbb{E}_{\boldsymbol{x}}\left[-p(\boldsymbol{x}) \log q(\boldsymbol{x}) - (1 - p(\boldsymbol{x})) \log (1 - q(\boldsymbol{x}))\right], \] where \(\mathcal{T}\) denotes the training set \(\hat f\) was fit on. We can then check how well training error, validation error and cross validation error estimate it. In real life you never get to do this, which is why we use a simulation.

You are given four files:

ImportantTODO (Setup)

TODO 0.1. Create a new notebook. Copy the code below into its first cell and run it. It downloads the four files into the folder of your notebook and imports everything you need for the whole lab.

import urllib.request
import numpy as np
import matplotlib.pyplot as plt
import lightgbm as lgb
from scipy.stats import sem
from sklearn.dummy import DummyClassifier
from sklearn.metrics import log_loss, make_scorer
from sklearn.model_selection import (KFold, GroupKFold, train_test_split,
                                     cross_val_score, cross_validate, GridSearchCV)

base_url = "https://stat154.berkeley.edu/fall-2026/datasets/cv_lab"
for f in ["simulator.py", "X_train.npy", "y_train.npy", "X_test.npy"]:
    urllib.request.urlretrieve(f"{base_url}/{f}", f)

from simulator import simulate, true_risk, bayes_risk

# Log loss as a scikit-learn scorer. labels=[0, 1] makes it work even when a validation
# fold contains only one class (the built-in "neg_log_loss" returns nan there; see Part 3).
log_loss_scorer = make_scorer(log_loss, greater_is_better=False,
                              response_method="predict_proba", labels=[0, 1])

The simulator gives you:

  • simulate(n, seed) - draws a training set X, y of n IID observations with 10 features.
  • true_risk(model) - the true risk \(\text{Err}_\mathcal{T}\) of a fitted model. It averages the formula above over 20,000 new feature vectors from the same distribution, which the model never saw. This is the log loss on a huge test set, with each label replaced by its true probability to remove label noise.
  • bayes_risk() - the smallest possible risk, achieved by predicting \(p(\boldsymbol{x})\) itself.
NoteGradient boosting and LightGBM in a nutshell

The model in this lab is a gradient boosted ensemble of decision trees. You do not need the details. Think of it as a flexible black box whose complexity is controlled by hyperparameters, which is exactly when cross validation is needed.

  • Decision tree. Splits the feature space into rectangles with a sequence of yes/no questions such as “is \(x_3 > 0.7\)?”, and predicts a constant in each rectangle (a leaf).
  • Boosting. Instead of one big tree, add many small trees one at a time. Each new tree is fit to correct the errors the current ensemble still makes on the training data (formally, it is fit to the negative gradient of the loss, hence gradient boosting). The sum of the trees is a log odds \(F(\boldsymbol{x})\), turned into a probability by \(q(\boldsymbol{x}) = 1 / (1 + e^{-F(\boldsymbol{x})})\).
  • LightGBM is a fast and widely used implementation of gradient boosting.
Hyperparameter Meaning Larger value gives a model that is
min_child_samples minimum number of training observations in a leaf simpler
n_estimators number of trees more complex
learning_rate how much each new tree’s correction counts more complex

n_estimators and learning_rate trade off against each other: a smaller learning rate needs more trees.

ImportantTODO (Setup)

TODO 0.2. Write a function make_model(min_child_samples=50) that returns a new, unfitted LightGBM classifier with:

  • min_child_samples set to the function’s argument,
  • n_jobs=1, which is about 10 times faster than the default on small datasets like ours,
  • verbosity=-1, which silences LightGBM’s warnings,
  • LightGBM’s defaults for every other hyperparameter.

Check: make_model() should display as LGBMClassifier(min_child_samples=50, n_jobs=1, verbosity=-1). You will use make_model in Parts 1-5.

May be useful: lgb.LGBMClassifier.

Predicted probabilities of a fitted model come from model.predict_proba(X)[:, 1].

Simulation noise. Whenever you average a quantity over \(R\) simulated datasets, also report its standard error \(\text{sd} / \sqrt{R}\), which sem (imported earlier from scipy.stats) computes. A difference smaller than about two standard errors could be simulation noise.

1 Training error vs. true risk

ImportantTODO (Part 1)

TODO 1.1. Simulate one training set of 500 observations with seed 0.

May be useful: simulate.

TODO 1.2. For each value of min_child_samples in [1, 2, 5, 10, 20, 40, 80, 160, 320], compute these three numbers and store each in its own list:

  • the training error: fit a model with that min_child_samples on all 500 observations, and compute its log loss on the same 500 observations,
  • the 5-fold CV error, using shuffled folds with random_state=0,
  • the true risk of the model fit on all 500 observations.

May be useful: make_model, .fit, .predict_proba, log_loss, KFold, cross_val_score with scoring=log_loss_scorer, true_risk. Note that cross_val_score returns one score per fold, and the scores are negative log losses.

TODO 1.3. Plot the three curves against min_child_samples, with a log scale on the x axis, a horizontal dashed line at the Bayes risk, and a legend.

May be useful: plt.plot, plt.xscale, plt.axhline, bayes_risk, plt.legend.

TODO 1.4. Answer in a markdown cell:

  1. Which value of min_child_samples minimizes the training error? The CV error? The true risk?
  2. How does the gap between training error and true risk change as the model becomes more complex? Why?
  3. What would happen if you selected min_child_samples by minimizing training error?

2 What does CV estimate?

Here you check what CV error is an unbiased estimator of. We use \(K = 2\) so that each training fold has only \(N/2 = 250\) observations, which makes the effect easy to see. Use make_model(), i.e. min_child_samples=50, throughout this part.

ImportantTODO (Part 2)

TODO 2.1. Repeat the following for \(R = 50\) datasets, r = 0, ..., 49:

  • simulate a dataset of 500 observations with seed 1000 + r,
  • compute \(\text{Err}_\mathcal{T}\): the true risk of the model fit on all 500 observations,
  • compute the 2-fold CV error, using shuffled folds with random_state=r,
  • compute the true risk of each of the 2 models that CV fit on the training folds (250 observations each), and average the two.

At the end you should have three arrays of length 50.

May be useful: cross_validate with return_estimator=True. It returns a dictionary with the fold scores under "test_score" and the fitted fold models under "estimator".

TODO 2.2. Print the average of each of the three arrays over the 50 datasets, each with its standard error.

May be useful: np.mean, sem.

TODO 2.3. Make a scatter plot with \(\text{Err}_\mathcal{T}\) on the x axis and the CV error on the y axis (one point per dataset), and add the line \(y = x\). Print the correlation between the two.

May be useful: plt.scatter, plt.plot, np.corrcoef.

TODO 2.4. Answer in a markdown cell:

  1. Which of the two true risks does the CV error match on average?
  2. Is CV error an optimistic or a pessimistic estimate of the risk of the model fit on all 500 observations? Why?
  3. If one dataset has a lower CV error than another, does its model tend to have a lower true risk? Use your scatter plot and correlation.
  4. CV averages over models fit on different training sets. Based on your scatter plot, does CV estimate \(\text{Err}_\mathcal{T}\), the risk of the model fit on your dataset, or the expected risk \(\text{Err} = \mathbb{E}_\mathcal{T}[\text{Err}_\mathcal{T}]\) over random datasets of the same size? (Compare with ESL Section 7.12.) What does this mean for hyperparameter tuning: is the value CV chooses the best for your dataset, or for a typical dataset like yours?

3 When CV breaks

CV works when the validation data is distributed like the data the model will be used on and was not used in training. Each of the two experiments below breaks one of these principles.

  • Experiment B: duplicated observations. Sometimes the same unit appears more than once in a dataset, for example two measurements of the same patient.
  • Experiment A: data that arrives sorted. Real data files are often sorted, for example by label.

Use make_model(), i.e. min_child_samples=50, and log_loss_scorer throughout this part.

ImportantTODO (Part 3)

TODO 3.1 (Experiment A). Repeat the following for \(R = 20\) datasets, r = 0, ..., 19:

  • simulate a dataset of 500 observations with seed 1000 + r,
  • sort its rows by label, so that all observations with \(y = 0\) come first,
  • compute three numbers:
    • the 5-fold CV error without shuffling the folds,
    • the 5-fold CV error with shuffled folds (random_state=r),
    • the true risk of the model fit on all 500 sorted observations.

May be useful: np.argsort(y, kind="stable") and indexing X[order]. KFold(5) does not shuffle by default. Some unshuffled folds contain only one class, which log_loss_scorer handles thanks to labels=[0, 1].

TODO 3.2 (Experiment B). Repeat the following for r = 0, ..., 19:

  • simulate a dataset of 250 observations with seed 1000 + r,
  • repeat every row twice, so that the dataset has 500 rows and the two copies of each row are next to each other,
  • make an array groups of length 500 recording which of the 250 original observations each row came from: 0, 0, 1, 1, 2, 2, ...,
  • compute three numbers:
    • the 5-fold CV error with shuffled folds (random_state=r),
    • the 5-fold CV error with folds that keep all copies of an observation in the same fold,
    • the true risk of the model fit on all 500 rows.

May be useful: np.repeat (with axis=0 for X), np.arange, GroupKFold, the groups argument of cross_val_score.

TODO 3.3. For each experiment, print the average of each of its three numbers over the 20 datasets, each with its standard error.

TODO 3.4. Answer in a markdown cell:

  1. In each experiment, which CV estimate is far from the true risk, and in which direction?
  2. Which of the two principles does each broken version violate? Explain how.
  3. Why do shuffling and GroupKFold fix the problems?

4 Size of the validation set

Now consider a single train-validation split instead of CV. What we want to estimate is \(\text{Err}_\mathcal{T}\), the true risk of the model fit on all 500 observations. Use make_model(), i.e. min_child_samples=50, throughout this part.

ImportantTODO (Part 4)

TODO 4.1. Repeat the following for \(R = 100\) datasets, r = 0, ..., 99:

  • simulate a dataset of 500 observations with seed 1000 + r,
  • compute \(\text{Err}_\mathcal{T}\): the true risk of the model fit on all 500 observations,
  • for each validation fraction in [0.1, 0.2, 0.3, 0.5, 0.7]:
    • randomly split the data into a training part and a validation part containing that fraction of the observations, with random_state=r,
    • fit a model on the training part only,
    • compute its log loss on the validation part: the validation error.

At the end you should have 100 values of \(\text{Err}_\mathcal{T}\), and 100 validation errors for each fraction.

May be useful: train_test_split with test_size, a dictionary {fraction: list of validation errors}.

TODO 4.2. For each fraction, compute and print:

  • the bias: the average over the 100 datasets of (validation error \(-\) \(\text{Err}_\mathcal{T}\)), with its standard error,
  • the standard deviation of the validation error across the 100 datasets.

May be useful: sem, np.std with ddof=1.

TODO 4.3. Make two plots side by side, both with the validation fraction on the x axis:

  • left: the bias, with error bars of \(\pm 2\) standard errors, and a horizontal line at 0,
  • right: the standard deviation.

May be useful: plt.subplots(1, 2), ax.errorbar, ax.axhline.

TODO 4.4. Answer in a markdown cell:

  1. How do the bias and the standard deviation change as the validation set gets larger? Explain each trend.
  2. Which validation fraction would you choose here, and why? Why can’t a single split make both small?

5 Number of folds

Use make_model(), i.e. min_child_samples=50, throughout this part.

ImportantTODO (Part 5)

TODO 5.1. Repeat the following for \(R = 100\) datasets, r = 0, ..., 99:

  • simulate a dataset of 500 observations with seed 1000 + r,
  • compute \(\text{Err}_\mathcal{T}\): the true risk of the model fit on all 500 observations,
  • for each \(K\) in [2, 5, 10, 20], compute:
    • the \(K\)-fold CV error with shuffled folds (random_state=r),
    • the computation time of the CV: the total time spent fitting and scoring over the \(K\) folds.

May be useful: cross_validate. Its dictionary also contains the time of each fold under "fit_time" and "score_time".

TODO 5.2. For each \(K\), compute and print:

  • the bias: the average over the 100 datasets of (CV error \(-\) \(\text{Err}_\mathcal{T}\)), with its standard error,
  • the standard deviation of the CV error across the 100 datasets,
  • the average computation time per dataset, in seconds.

TODO 5.3. Make three plots side by side, all with \(K\) on the x axis: the bias with error bars of \(\pm 2\) standard errors and a horizontal line at 0, the standard deviation, and the time.

May be useful: plt.subplots(1, 3).

TODO 5.4. Answer in a markdown cell:

  1. How does the bias change with \(K\)? Relate this to what you found in Part 2.
  2. How does the standard deviation change with \(K\)? Does it keep decreasing?
  3. How does the computation time grow with \(K\)? Given all three, why are \(K = 5\) or \(K = 10\) common choices?
TipOptional

ESL (Section 7.10.1) says that leave-one-out CV (\(K = N\)) “can have high variance because the N ‘training sets’ are so similar to one another”. Check this by repeating TODO 5.1 and 5.2 with datasets of 100 observations, min_child_samples=10, \(R = 50\), and \(K\) in [5, 10, 100] (\(K = 100\) is leave-one-out). Do your standard deviations agree with ESL?

6 Challenge

The challenge data X_train.npy, y_train.npy (2,000 observations) and X_test.npy (10,000 observations) are IID draws from a different, unknown distribution: not the simulator. As in real life, there is no true_risk. I have hidden y_test from you. Your task is to predict the probabilities \(P(y = 1 \mid \boldsymbol{x})\) for X_test as well as possible, using only X_train and y_train. Read Allowed methods and How to submit predictions below before you start.

For every cross validation in this part, use 5 shuffled folds with random_state=0, and log_loss_scorer.

ImportantTODO (Part 6)

TODO 6.1. Load the data. Load the three challenge files into X_train, y_train and X_test.

May be useful: np.load.

TODO 6.2. Baselines. Compute and print the 5-fold CV error of two models:

  • DummyClassifier(strategy="prior"), which predicts the training-fold mean of y for every observation,
  • lgb.LGBMClassifier(n_jobs=1, verbosity=-1), which is LightGBM with its default hyperparameters.

TODO 6.3. Tune one hyperparameter. Starting from lgb.LGBMClassifier(n_jobs=1, verbosity=-1), use a grid search with 5-fold CV to choose min_child_samples from a grid of your choice, for example [5, 10, 20, 40, 80, 160, 320, 640]. Plot the CV error against min_child_samples (log scale on the x axis), and print the best value and its CV error.

May be useful: GridSearchCV with param_grid, scoring and cv. After .fit, it has the attributes cv_results_["mean_test_score"] (the scores of every grid value), best_params_ and best_score_. As with cross_val_score, the scores are negative log losses.

TODO 6.4. Tune several hyperparameters jointly. Repeat the grid search, now also tuning n_estimators and learning_rate together with min_child_samples. Print the best combination and its CV error. A grid search tries every combination of values and fits 5 models for each, so keep the total number of combinations under about 50.

May be useful: param_grid takes a dictionary with one list of values per hyperparameter.

TODO 6.5. Final predictions. Predict the probabilities for X_test with your best model, check that there are 10,000 of them, and save them to pred_{codename}.npy, with your group’s codename in place of {codename} (see How to submit predictions).

May be useful: the best_estimator_ attribute of your grid search, which is the best model already refit on all of X_train, y_train (because refit=True by default). np.save.

TODO 6.6. Answer in a markdown cell:

  1. Which hyperparameters did you choose, and what is their CV error? How much better is it than the baselines?
  2. Two things make your CV error differ from the test error of your final model. First, your final model is trained on 2,000 observations, while each CV model was trained on 1,600. Second, you chose the hyperparameters with the lowest CV error out of many. In which direction does each push? Write down a guess of your test error, and compare it with the leaderboard once it is published. (We will come back to the second effect in the next lab.)

6.1 Allowed methods

You may only use lgb.LGBMClassifier, with any hyperparameters you like, trained only on X_train and y_train. The simulator comes from a different distribution, so drawing extra data from it will not help you.

6.2 How to submit predictions

Submit a binary .npy file containing a one-dimensional numpy array of predicted probabilities (numbers between 0 and 1), whose length is the same as the number of rows in X_test:

np.save("pred_{codename}.npy", pred)

The name of the file should be "pred_{codename}.npy" where {codename} is a secret codename of your group’s choosing, that you will be able to use later to identify your submission in the leaderboard.

6.3 How your predictions will be scored

Your predictions will be compared against the real y_test using average log loss: \[ \text{your loss} = -\frac{1}{N}\sum_{n=1}^N \left( y_n \log \text{pred}_n + (1 - y_n) \log (1 - \text{pred}_n) \right). \]

A leaderboard will be published with every group’s codename alongside its loss on the test set, ranked from low to high. No other identifying information will be shared. Your rank on the leaderboard will not factor into your grade in any way.

Submission

ImportantTODO (Submission)

Upload these two files on Gradescope by Tuesday at 9 PM:

  1. your notebook, with every TODO in Parts 1-6 done and every question answered,
  2. pred_{codename}.npy from TODO 6.5.

The grading will be based on completion of Parts 1-6.

LLM usage policy

You may use auto-completion for filling out small code snippets (affect around 2-3 lines of code) and documenting your code. Small code snippets are things like “add a label for the axis of the plot”, “change the color of the points”, “add an argument for the number of interactions in the fit”. These are small code operations which do not involve significant planning nor automated changes to the logic of the code without your supervision. All edits must be verified and accepted manually, and you should make sure you know what the code is doing.

Things that are explicitly not allowed are major edits to the code (affect > 10 lines) and ones which involve automating huge parts of the development of the logic structure of your analysis. For example, “write a function to extract features from the data”, or “run a test set validation loop and report the best model”.

Delegate to language models only things which you know exactly how you would have done it yourself or would be comfortable googling/asking for help from someone who already took the course/knows how to code, but who cannot do the learning for you.