Cross-validation

What this chapter covers

A single train/test split gives you one accuracy number, and that number is part luck. Reshuffle the split and it moves. On a small dataset it can move a lot — enough to make a worse model look better, which is the exact mistake you're trying to avoid when you measure at all. Cross-validation is the fix: instead of holding out one slice and trusting it, hold out every slice in turn, score each one, and average. You use all your data for testing and all of it for training, just never at the same time.

We build it from scratch in NumPy — the fold bookkeeping, the score loop, the mean and standard deviation — then hand the exact same splits to scikit-learn and check that our fold scores match its cross_val_score to the digit. Along the way you watch the folds rotate, the per-fold scores land, and the running mean tighten with an error bar that tells you how much to trust it. The payoff number at the end is the whole argument: one split's estimate wanders by four points of accuracy across seeds; the five-fold estimate barely moves.

The data is the Pima Indians Diabetes set again — 768 patients, eight measurements each, a binary diagnosis — with a shallow decision tree as the classifier. The tree is beside the point here. Cross-validation is the thing we build, and it works the same way around any model you drop into it.

A bit of history

The idea of holding data back to check yourself is old and obvious, but the version we use has two specific parents, both from the mid-seventies. Mervyn Stone's 1974 paper in the Journal of the Royal Statistical Society — "Cross-validatory choice and assessment of statistical predictions" — is the one that gave it the name and the theory. Stone laid out the rotation we still run: leave a point (or a group) out, predict it from the rest, do that for every point, and use the aggregate to both choose a model and estimate how well it will do. The next year Seymour Geisser published the predictive sample reuse method in the Journal of the American Statistical Association, coming at the same idea from a predictive angle and generalizing it. Between the two of them cross-validation went from folklore to a named procedure with a justification.

What makes the timing interesting is that they were arguing against the grain. The statistical culture of the day leaned on closed-form criteria — adjusted measures, information criteria — that estimate out-of-sample error from a formula on the training fit. Stone and Geisser said: don't estimate it, measure it, by actually predicting data the model didn't see. That's a computer's answer, and it only got more attractive as computers got cheaper. Today cross-validation is the default way anyone reports a model's accuracy on data that isn't huge, and every cross_val_score call is running Stone's rotation.

The intuition

Here's the move. Cut the data into k equal parts — call them folds. Now run k rounds. In each round, one fold sits out as the validation set and the other k-1 folds are the training set; you fit a fresh model on the training folds and score it on the held-out one. Rotate so that every fold is the held-out one exactly once. That gives you k scores, one per round, each measured on data the model in that round never touched. Average them and you have your estimate.

The picture is the whole idea. Each row below is one round; the highlighted cell is the fold held out that round; the grey cells are what the model trained on. Read down the diagonal and you watch the held-out fold march across the data until every part has been the test set once.

Every patient is a test case exactly once and a training case k-1 times. Nobody is wasted, and no score is ever measured on data the model trained on. That's the trick that makes the estimate honest and stable at the same time.

The math

Split the N rows into k disjoint folds. Let $F_j$ be the set of row indices in fold j, and let $\hat{y}_i$ be the prediction for row i from the model trained on every fold except the one containing i. The score for fold j is just the accuracy on that fold:

The cross-validation estimate is the average of the k fold scores:

That average is the number you report. But an average hides disagreement, and the disagreement is information — it tells you how much the estimate would wobble if the folds had fallen differently. So report the spread too, the standard deviation across the folds:

Two numbers, not one. The mean is your best guess at out-of-sample accuracy; the standard deviation is how loudly the folds argued about it. A model that scores $0.75 \pm 0.01$ and one that scores $0.75 \pm 0.06$ are not the same result, even though their headline is identical, and the chapter's whole point is that you can't see the difference without the second number.

What it's good at, what it isn't

The good side is that you stop throwing data away. On a small set — and 768 rows is small — a 20% test split means you're estimating your accuracy from 150 patients and training on the other 600, and both of those numbers being small hurts. Cross-validation trains on 80% and tests on 100% across the rounds, so every row does double duty. You get a lower-variance estimate and a full-data model, and you get the spread for free, which is the part people skip and then wonder why their leaderboard is noise. For model selection on limited data it's not really optional; it's how you tell a real improvement from a lucky split.

The cost is compute, and it's exactly k times. Five-fold cross-validation fits your model five times; leave-one-out fits it N times, once per row. On a shallow tree over 768 rows that's nothing, but on anything that takes minutes to train the multiplier is real, and it's why you see three- and five-fold in practice far more than ten. The other trap is subtler: if you tune anything by looking at the cross-validation score — the tree's depth, a decision cutoff, which columns you keep — then that score has seen the validation data and is optimistic again. The fix is nested cross-validation, and it's where this chapter points next.

The data

Same patients as chapter one: 768 people, eight measurements each — glucose, BMI, age, blood pressure and so on — and a label that says whether each was diagnosed with diabetes. What matters here isn't the features, it's the balance of the label. Only about 35% of the patients are positive, so the classes are lopsided, and that lopsidedness is the thing that will bite a naive fold split in a minute.

The classifier is a decision tree capped at depth three — shallow on purpose. It needs no feature scaling, trains instantly, and with a fixed random seed it's deterministic, so the only thing that can move a fold's score is the data in that fold. That's what lets us prove our folds match scikit-learn's exactly: same data, same model, same number.

Build it, one function at a time

Cross-validation is bookkeeping around a scoring loop, and it's easiest to build from the inside out. Start with the thing every fold needs — accuracy:

def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth — unitless, 0 to 1."""
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return float((y_true == y_pred).mean())

Now the folds. Plain k-fold just cuts the row indices into k contiguous blocks. The only wrinkle is the remainder when N doesn't divide by k: the first few folds get one extra row each. That's the same rule scikit-learn follows, and matching it is what makes the scores line up later.

def kfold_indices(n, k):
    """Split range(n) into k contiguous folds; return (train, val) index pairs.

    No shuffling: fold 0 is the first block of rows, fold 1 the next, and so
    on. When n isn't divisible by k the first n % k folds get one extra row —
    the same rule scikit-learn's KFold(shuffle=False) uses, so the splits line
    up exactly.
    """
    sizes = np.full(k, n // k, dtype=int)
    sizes[: n % k] += 1                      # spread the remainder over the first folds
    fold_id = np.repeat(np.arange(k), sizes)  # which fold each row belongs to
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]

Contiguous blocks have a problem on imbalanced data, though. If the rows aren't shuffled with respect to the label — and you can never assume they are — a block can end up with too many positives or too few. Stratified k-fold fixes that by splitting each class separately, so every fold inherits the class balance of the whole set:

def stratified_kfold_indices(y, k):
    """Like kfold_indices, but each fold keeps the class balance of the whole.

    Assign folds within each class separately so a fold can't end up with too
    few of a rare label. This reproduces scikit-learn's StratifiedKFold
    (shuffle=False): for every class the counts are spread across folds as
    evenly as possible, then handed out in the order the rows appear.
    """
    y = np.asarray(y)
    classes, y_enc = np.unique(y, return_inverse=True)
    # how many of each class land in each fold, balanced by the k-way stride
    ordered = np.sort(y_enc)
    per_fold = np.array([np.bincount(ordered[i::k], minlength=len(classes))
                         for i in range(k)])
    fold_id = np.empty(len(y), dtype=int)
    for c in range(len(classes)):
        fold_id[y_enc == c] = np.repeat(np.arange(k), per_fold[:, c])
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]

The allocation line is the whole trick: for each class, spread its members across the folds as evenly as the counts allow, then hand them out in order. It reproduces scikit-learn's StratifiedKFold exactly, which we'll lean on when we check the numbers.

Leave-one-out is k-fold taken to its limit — one row per fold, N folds. It needs no new code; it's just k = N:

def loo_indices(n):
    """Leave-one-out: the k = n extreme. Each row is its own validation fold."""
    return kfold_indices(n, n)

With the folds in hand, the driver is a short loop: for each split, fit a fresh model on the training rows, predict the validation rows, record the accuracy. The estimator comes in as a callable so this loop doesn't care whether it's wrapping a tree, a logistic regression, or a neural net.

def cross_val_score(fit_predict, X, y, splits):
    """Train and score once per fold; return the k raw scores.

    `fit_predict(X_train, y_train, X_val)` fits a fresh model on the training
    rows and returns predictions for the validation rows — so this loop is
    estimator-agnostic. A new model is fit for every fold; nothing leaks across
    them.
    """
    X, y = np.asarray(X), np.asarray(y)
    scores = np.empty(len(splits))
    for j, (train, val) in enumerate(splits):
        y_pred = fit_predict(X[train], y[train], X[val])
        scores[j] = accuracy(y[val], y_pred)
    return scores

Last, collapse the k scores into the two numbers you actually report — the mean and its spread:

def cv_summary(scores):
    """The headline number and its uncertainty: mean and std of the fold scores.

    The std here is the population std over the k folds (ddof=0) — the same
    thing `scores.std()` gives. Report both: the mean is your estimate, the std
    is how far the folds disagreed about it.
    """
    scores = np.asarray(scores)
    return float(scores.mean()), float(scores.std())

Watch it work

This is five-fold cross-validation running for real on the 768 patients, one frame per round. The top strip is the split for that round — four folds in grey training the model, one in orange held out. The lower panel is where the fold scores land as they come in: each point is one fold's accuracy, the dashed line is the running mean, and the shaded band is one standard deviation on either side of it. The metric is accuracy — the fraction of held-out patients the tree classifies correctly, unitless, from 0 to 1.

Press play. Watch the held-out fold march left to right while the running mean settles and the band forms. Reset and run it as often as you like — this is plain k-fold with no shuffling, so it's fully deterministic; the same folds fall the same way every time.

Now look at how much those five points disagree. Fold four scores 0.830 and fold two scores 0.669 — sixteen points of accuracy apart, on the same data, the same model, just a different slice held out. The mean lands at 0.750 but the band is wide, $\pm 0.057$, and that width is the honest signal that a single slice of this data is a shaky thing to trust. If you'd run one train/test split and happened to draw fold four as your test set, you'd have reported 0.83 and believed it.

Here are the five stratified fold scores on their own, with the mean line and the $\pm 1$ standard deviation band drawn in. This is the version you'd actually report:

Notice the band is far tighter than the plain run's — $\pm 0.012$ instead of $\pm 0.057$. That's stratification earning its keep. Remember the label is only 35% positive; when you cut the rows into plain contiguous blocks, the positive rate per fold swings from 0.255 to 0.416, so some folds are much easier or harder than others and the scores scatter. Stratified folds hold the positive rate near 0.35 in every fold, so the folds are comparable and the estimate steadies. Here's the same five folds both ways:

The stratified points cluster; the plain points spray across the axis. Same data, same model, same k — the only difference is whether the folds kept the class balance. On imbalanced data, stratify. It costs nothing and it's the difference between an estimate you can quote and one that's mostly noise.

The full implementation

The whole thing, top to bottom — six small functions and a loader, no scikit-learn anywhere in the file. This is what the animation above actually ran:

"""Cross-validation, built from scratch in pure NumPy.

The idea: instead of trusting one train/test split, cut the data into k folds,
train on k-1 of them and score on the one you left out, and rotate until every
fold has been the held-out set exactly once. Average the k scores for a stable
estimate, and keep their spread so you know how much to trust it.

This file is only the CV machinery — the fold bookkeeping, the score loop, the
mean/std. The estimator is handed in as a `fit_predict` callable, so nothing
here imports scikit-learn. Every function appears in the chapter one step at a
time (the `# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: accuracy
def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth — unitless, 0 to 1."""
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return float((y_true == y_pred).mean())
# endregion


# region: kfold_indices
def kfold_indices(n, k):
    """Split range(n) into k contiguous folds; return (train, val) index pairs.

    No shuffling: fold 0 is the first block of rows, fold 1 the next, and so
    on. When n isn't divisible by k the first n % k folds get one extra row —
    the same rule scikit-learn's KFold(shuffle=False) uses, so the splits line
    up exactly.
    """
    sizes = np.full(k, n // k, dtype=int)
    sizes[: n % k] += 1                      # spread the remainder over the first folds
    fold_id = np.repeat(np.arange(k), sizes)  # which fold each row belongs to
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]
# endregion


# region: stratified_kfold_indices
def stratified_kfold_indices(y, k):
    """Like kfold_indices, but each fold keeps the class balance of the whole.

    Assign folds within each class separately so a fold can't end up with too
    few of a rare label. This reproduces scikit-learn's StratifiedKFold
    (shuffle=False): for every class the counts are spread across folds as
    evenly as possible, then handed out in the order the rows appear.
    """
    y = np.asarray(y)
    classes, y_enc = np.unique(y, return_inverse=True)
    # how many of each class land in each fold, balanced by the k-way stride
    ordered = np.sort(y_enc)
    per_fold = np.array([np.bincount(ordered[i::k], minlength=len(classes))
                         for i in range(k)])
    fold_id = np.empty(len(y), dtype=int)
    for c in range(len(classes)):
        fold_id[y_enc == c] = np.repeat(np.arange(k), per_fold[:, c])
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]
# endregion


# region: loo_indices
def loo_indices(n):
    """Leave-one-out: the k = n extreme. Each row is its own validation fold."""
    return kfold_indices(n, n)
# endregion


# region: cross_val_score
def cross_val_score(fit_predict, X, y, splits):
    """Train and score once per fold; return the k raw scores.

    `fit_predict(X_train, y_train, X_val)` fits a fresh model on the training
    rows and returns predictions for the validation rows — so this loop is
    estimator-agnostic. A new model is fit for every fold; nothing leaks across
    them.
    """
    X, y = np.asarray(X), np.asarray(y)
    scores = np.empty(len(splits))
    for j, (train, val) in enumerate(splits):
        y_pred = fit_predict(X[train], y[train], X[val])
        scores[j] = accuracy(y[val], y_pred)
    return scores
# endregion


# region: cv_summary
def cv_summary(scores):
    """The headline number and its uncertainty: mean and std of the fold scores.

    The std here is the population std over the k folds (ddof=0) — the same
    thing `scores.std()` gives. Report both: the mean is your estimate, the std
    is how far the folds disagreed about it.
    """
    scores = np.asarray(scores)
    return float(scores.mean()), float(scores.std())
# endregion


FEATURES = [
    "Pregnancies", "Glucose", "BloodPressure", "SkinThickness",
    "Insulin", "BMI", "DiabetesPedigreeFunction", "Age",
]


def load_xy(path="../data/diabetes.csv"):
    """Pima Indians Diabetes: 768 patients, 8 features, binary Outcome."""
    df = pd.read_csv(path)
    return df[FEATURES].to_numpy(float), df["Outcome"].to_numpy(int)

The library version

Nobody writes their own fold loop in production, and once you've written it once you shouldn't either. scikit-learn's model_selection module has the splitters and the scorer, and they're the exact counterparts of what we built: KFold for kfold_indices, StratifiedKFold for stratified_kfold_indices, LeaveOneOut for the k = N case, and cross_val_score for the driver loop.

def make_model():
    """A shallow, deterministic classifier: depth-3 tree, fixed seed.

    Shallow so it needs no feature scaling and trains in milliseconds, seeded
    so the only thing that moves a fold's score is the data in that fold — not
    the model. That's what makes the scratch-vs-library scores match to the
    digit.
    """
    return DecisionTreeClassifier(max_depth=3, random_state=0)

def sklearn_kfold(X, y, k):
    """cross_val_score over a plain k-fold split (no shuffle)."""
    return cross_val_score(make_model(), X, y, cv=KFold(n_splits=k, shuffle=False))


def sklearn_stratified(X, y, k):
    """cross_val_score over a stratified k-fold split (no shuffle)."""
    return cross_val_score(make_model(), X, y,
                           cv=StratifiedKFold(n_splits=k, shuffle=False))


def sklearn_loo(X, y):
    """cross_val_score with leave-one-out (k = n folds)."""
    return cross_val_score(make_model(), X, y, cv=LeaveOneOut())

One thing worth flagging, because it trips people up: if you call cross_val_score(model, X, y, cv=5) with a bare integer on a classifier, scikit-learn quietly gives you stratified folds, not plain ones. That's a good default, but it means the plain-vs-stratified difference we just watched is invisible unless you pass the splitter object yourself. We pass it explicitly in both cases so the comparison is real.

Scratch versus library

The claim is that our from-scratch folds and scikit-learn's are the same split, so the fold scores should be identical, not merely close. They are. Here's the mean of the five fold scores, ours against the library's, for both the plain and the stratified split — the bars land exactly on top of each other because the underlying numbers match to the digit:

The generator asserts this equality on every run — plain, stratified, and leave-one-out all match scikit-learn fold for fold — so if the two ever drifted, the data wouldn't build. Matching the library isn't the point in itself; it's the proof that the from-scratch version is the real algorithm and not a lookalike.

Now the number that justifies the whole chapter. I ran two estimators of the tree's accuracy thirty times each, with thirty different seeds. The first is a single 80/20 train/test split — one held-out slice, one number. The second is the five-fold cross-validation mean. Same data, same model, same thirty seeds; the only difference is one estimate looks at one slice and the other averages five. Here's every run, both ways:

Both estimators center on the same place — 0.7385 for the single split, 0.7386 for the cross-validation mean, so cross-validation isn't biased, it's not cheating the accuracy up or down. What changes is the scatter. The single split ranges from 0.649 to 0.799 across seeds, a standard deviation of 0.040 — draw a different random test set and your reported accuracy swings by fourteen points between the unlucky seed and the lucky one. The cross-validation mean ranges from 0.718 to 0.763, a standard deviation of 0.010 — a quarter of the wobble. That's the trade in one picture: same answer on average, four times steadier. The single-split number isn't wrong, it's just imprecise, and on a small dataset imprecise is how you ship the worse model because it drew the kinder split.

Takeaways

Use cross-validation for any accuracy number you're going to make a decision on, and make the decision on the mean while quoting the standard deviation next to it. The mean is your estimate; the spread is your confidence, and reporting one without the other is how model comparisons turn into coin flips. If you take one habit from this chapter, take that: two numbers, always. I've watched teams pick model A over model B on a 0.02 gap that was smaller than either model's fold-to-fold spread, which means they picked at random and dressed it up as progress.

Reach for it hardest exactly when data is scarce, because that's when a single split is least trustworthy and cross-validation buys you the most. On imbalanced labels, stratify — it's free and it's the difference we watched between a tight estimate and a noisy one. Five folds is the sensible default; go to ten when the model is cheap and you want a touch more precision, and to leave-one-out only on genuinely tiny sets, where its N-times cost is affordable and you can't spare a single row for a test fold. The cost you're always paying is compute, k times over, so keep that in the back of your mind when the model gets expensive.

The one thing cross-validation can't do by itself is protect a score you've tuned against. The moment you pick a hyperparameter, a threshold, or a feature set by its cross-validation score, that score has peeked at the validation data and gone optimistic — the honest fix is to wrap the tuning in an outer cross-validation loop, which is nested cross-validation. That's the natural next step, and it's exactly where hyperparameter tuning picks up: now that you can measure a model reliably, you can start searching for the best one without fooling yourself.