Sampling, Confidence Intervals and Hypothesis Testing

Context

A test score is computed from one held-out set of listings, so it carries the luck of which listings were held out. A confidence interval (a range built by a procedure that holds the true score in most samples) says how far that luck reaches, and a hypothesis test (a check of whether a result could be luck alone) asks whether the gap between two models scored on the same listings could be that luck.

A pricing team quotes a model's average miss to the business with the range it could stand for. A reviewer decides whether boosting's lower test MAE is a better model or a better draw of 600 listings, and a report that compares three models, seven neighborhoods or twenty settings has to say how many of its gaps luck alone would produce.

Without them, the failures read as findings. A gain of a few thousand is announced as a better model when it sits within the luck of 600 listings, two overlapping intervals are read as no difference while the paired comparison of the same listings shows one, a p-value (the chance of a gap at least that large if the two models were equally good) of 0.092 is read as proof that two models are equal, and the one comparison in twenty that fell below 0.05 is reported as the result. Lesson 3.6.3 kept the test rows for one final score, and this lesson measures how far that one score, and the gap between two of them, can be trusted.

Goal

• Compute a test MAE's standard error from its 600 per-listing errors, and show how far an MAE over 100 listings swings.
• Build a 95 percent bootstrap interval for a score with np.random.default_rng(42), set the normal approximation beside it, and state what the 95 percent belongs to.
• Compare two fixed models on the same test listings by paired differences and a paired interval, and explain why overlapping separate intervals are not the test.
• Run a sign-flip test by hand beside scipy.stats.ttest_rel, and explain precisely what a p-value says and what a large one does not.
• Measure how often tests of a difference that is not there fall below 0.05 when they are run 20 at a time.

Prerequisites

None of this is required, and the lesson stands on its own. Both are here if you want resampling and p-values drawn out on small examples before you start.

• YouTube - "Bootstrapping Main Ideas!!!" by StatQuest with Josh Starmer. A handful of measurements resampled with replacement over and over, the mean of each resample kept, and a standard error and a confidence interval read from how those means spread.
• YouTube - "p-values: What they are and how to interpret them" by StatQuest with Josh Starmer. Why a p-value is needed, what it is, the thresholds people set for calling one small, and why a p-value does not measure how large a difference is.

Table of Contents

• A Test Score Is a Sample
• The Bootstrap and a Confidence Interval
• Comparing Two Models on the Same Rows
• Hypothesis Tests and What a p-Value Says
• Many Comparisons

Theory

Nothing to install: SciPy came with scikit-learn in Lesson 3.1.1, and every library test this lesson runs comes from scipy.stats, which held Lesson 3.6.3's distributions and holds SciPy's statistical tests too. Every script is saved in ~/dev/bootcamp beside listings.csv and price_models.py, the module saved in Lesson 3.5.1, whose split, make_prep and fill every model here is built from. The first block is a module, held_out_errors.py, which every script imports. Every other Theory block goes into one script, is_real.py, each added under the one before it and run with python3 is_real.py, except one block in the second subtitle, which is a script of its own, coverage.py.

Nothing here chooses a model. The three models were each fitted once on the 2,400 training rows, with settings chosen before the test rows scored them, and this lesson measures how far their test scores, and the gaps between them, can be trusted. Every interval and every test is computed from the 600 test listings' errors, and no model is fitted again.

A Test Score Is a Sample

What. A test MAE is the mean of 600 absolute errors, one per held-out listing. Those 600 listings are one sample from the listings they were drawn from, and another 600 drawn the same way would give another mean, so the score carries the luck of the draw. Why. Lesson 3.1.2 moved the split's seed and watched a held-out MAE move by thousands. One split's errors already say how far a score can move, with no new split and no refit, and the standard error says it in one number. When. Whenever a test score is reported, and before any gap between two scores is read.

held_out_errors.py fits the three models of Sections 3.2 and 3.5 once each on the 2,400 training rows: Lesson 3.2.3's linear pipeline, Lesson 3.5.2's chosen forest with max_features=0.5 and its 100 trees, and Lesson 3.5.3's HistGradientBoosting at its defaults, the last two on the tree preparation, make_prep("passthrough"):

import numpy as np
from sklearn.ensemble import (
    HistGradientBoostingRegressor,
    RandomForestRegressor,
)
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from price_models import X_test, X_train, fill, make_prep, y_test, y_train

models = {
    "linear": make_pipeline(make_prep(fill()), LinearRegression()),
    "forest": make_pipeline(
        make_prep("passthrough"),
        RandomForestRegressor(max_features=0.5, random_state=42),
    ),
    "boosting": make_pipeline(
        make_prep("passthrough"),
        HistGradientBoostingRegressor(random_state=42),
    ),
}
errors = {}
for name, model in models.items():
    model.fit(X_train, y_train)  # once, on the 2,400 rows
    # one absolute error per test listing, in the test rows' order
    errors[name] = np.abs(y_test.to_numpy() - model.predict(X_test))

errors maps each model's name to a NumPy array of its 600 absolute errors, in the test rows' order, so position 0 is the same listing in all three. is_real.py opens with each model's MAE, the mean of its array:

import numpy as np
from held_out_errors import errors

linear = errors["linear"]
print(len(linear), round(linear.mean()))  # 600 56562
for name in ["forest", "boosting"]:
    print(name, round(errors[name].mean()))

forest 53802
boosting 52773

These are the 600 test listings on which Lessons 3.5.2 and 3.5.3 put the three models' R² at 0.766, 0.775 and 0.777, now as the average miss in currency, Lesson 3.1.4's MAE. The test rows came out of train_test_split shuffled, so every run of 100 consecutive errors is a random 100 of the 600. Lesson 2.2.1's reshape cuts the array into six such slices, no listing in two of them:

# six slices of 100 test listings, one MAE each
slices = linear.reshape(6, 100).mean(axis=1)
print([round(mae) for mae in slices])

[61392, 65649, 57997, 53257, 52821, 48252]

One model, one fit and six held-out sets of 100 listings put its MAE anywhere from 48,252 to 65,649. How far an MAE can move is what the standard error (the standard deviation a mean would show across many samples of the same size) measures, estimated from one sample as the values' standard deviation over the square root of their count:

spread = linear.std(ddof=1)
standard_error = spread / np.sqrt(len(linear))
print(round(spread), round(standard_error))  # 74222 3030
print(round(spread / np.sqrt(100)))          # 7422

ddof=1 divides by n minus 1, the sample version Lesson 2.2.5 named, because the 600 errors are a sample. They spread by 74,222 around their mean of 56,562, so the 600-listing MAE has a standard error of 3,030. An MAE over 100 listings has one of 7,422, about two and a half times as large, and the six slices, each an MAE over 100 listings, spread on that scale. The standard error shrinks with the square root of the count: six times the listings divides it by about 2.45, and halving it takes four times the listings. Lesson 3.1.2's ten seeds measured the same kind of movement by refitting on ten splits, and one split's errors measure it with no refit at all.

DON'T:

• DON'T report a test MAE as the model's error to the last unit. 56,562 is the mean of one sample of 600 listings, with a standard error of 3,030.
• DON'T read an MAE over a small held-out set as closely as one over a large set. Six sets of 100 listings put one fitted model anywhere from 48,252 to 65,649.

The Bootstrap and a Confidence Interval

What. A bootstrap resample of the errors is Lesson 3.5.2's bootstrap sample drawn from the 600 test errors: 600 errors drawn with replacement, so some listings count twice or more and some not at all. Each resample's MAE is a score a sample like this one could have given, and the 2.5th to the 97.5th percentile of 2,000 such MAEs is a 95 percent confidence interval (a range computed from a sample by a procedure built so that, over many samples, its ranges hold the true value 95 times in 100). Why. The standard error's formula serves a mean. The bootstrap needs no formula, so it gives an interval for a median as readily as for a mean, and because it resamples errors and not training rows, no model is fitted again. When. With every score that is reported, as the range beside the number, and for any score whose standard error has no formula in this module.

rng.integers, Lesson 2.2.1's, draws all 2,000 resamples at once as rows of positions, and Lesson 2.2.2's fancy indexing with a 2-D array of positions returns an array of the same shape, here 2,000 rows of 600 errors, one resample per row:

rng = np.random.default_rng(42)
# 2,000 resamples, each 600 positions drawn with replacement
picks = rng.integers(0, 600, size=(2000, 600))
boot_mae = linear[picks].mean(axis=1)
print(boot_mae.shape)           # (2000,)
low, high = np.percentile(boot_mae, [2.5, 97.5])
print(round(low), round(high))  # 50970 62951

95 percent of the 2,000 resampled MAEs fall from 50,970 to 62,951, the percentile interval of the linear pipeline's MAE. The standard error gives a second interval. A normal distribution, Lesson 2.2.1's bell-shaped spread, holds 95 percent of its values within 1.96 standard deviations of its center, so the normal approximation is the MAE plus and minus 1.96 standard errors:

margin = 1.96 * standard_error
print(round(linear.mean() - margin))  # 50623
print(round(linear.mean() + margin))  # 62501
print(round(boot_mae.std()))          # 3008

The two routes, one through a formula and one through resamples, land within 450 of each other at either end, and the resampled MAEs spread with a standard deviation of 3,008 against the formula's 3,030. The percentile interval reaches further above the MAE than below it, while the normal one is symmetric by its construction. The median absolute error, the miss half the listings stay under, has no standard error formula in this module, and the same 2,000 resamples give it an interval all the same:

boot_median = np.median(linear[picks], axis=1)
print(round(np.median(linear)))  # 39237
print(np.percentile(boot_median, [2.5, 97.5]).round())

[35743. 44424.]

Half the test listings are missed by less than 39,237, while the average miss is 56,562, and the median carries an interval of its own, 35,743 to 44,424, from resamples alone.

An interval computed from one sample either holds the true MAE (the model's MAE over all the listings the test rows were drawn from, the number a test score estimates) or does not, and nothing in the sample says which. The 95 belongs to the procedure, and the procedure can be checked where the true MAE is known. A market that is exactly these 600 listings has the linear errors' mean as its true MAE, and samples drawn from it with replacement are held-out sets from that market. coverage.py draws 1,000 samples of 600, builds each one's percentile interval from 2,000 resamples, and counts the intervals that hold the market's true MAE and those that hold the MAE of the next sample drawn after them. It is a script of its own, so the rest of is_real.py never waits on its two million resampled MAEs:

import numpy as np
from held_out_errors import errors


def coverage(market, n, samples=1000):
    # samples of n listings from a market whose true MAE is known:
    # how many intervals hold it, and how many hold the next MAE
    sim = np.random.default_rng(42)
    true_mae = market.mean()
    held = next_held = 0
    for _ in range(samples):
        sample = market[sim.integers(0, len(market), size=n)]
        boot = sample[sim.integers(0, n, size=(2000, n))].mean(axis=1)
        low, high = np.percentile(boot, [2.5, 97.5])
        held += low <= true_mae <= high
        next_mae = market[sim.integers(0, len(market), size=n)].mean()
        next_held += low <= next_mae <= high
    return int(held), int(next_held)


print(coverage(errors["linear"], 600))  # (927, 810)

927 of the 1,000 intervals held the true MAE. The procedure aims at 950, and on these errors it falls short: a percentile interval is an approximation, close and not exact. Each of the 1,000 either held the true MAE or missed it, and once an interval was computed nothing about it was left to chance. The second count says what an interval is not. It held the next sample's MAE 810 times in 1,000, because the next sample's MAE carries its own draw as well. The linear pipeline's 50,970 to 62,951 is a range of true MAEs these 600 listings are consistent with, built by a procedure that held a known true MAE about 93 times in 100 on errors like these, and it is never a range the next test score will land in 95 times in 100.

DON'T:

• DON'T say the true MAE lies between 50,970 and 62,951 with probability 0.95. The 95 belongs to the procedure: on a market whose true MAE was known, 927 of 1,000 intervals held it, and each one held it or did not.
• DON'T read the interval as the range the next test score will land in. It held the next sample's MAE 810 times in 1,000.
• DON'T report a score without its interval. As far as these 600 listings can tell, the linear pipeline's 56,562 is consistent with any true MAE from 50,970 to 62,951.

Comparing Two Models on the Same Rows

What. Two models scored on the same 600 listings give each listing a paired difference (one listing's error under the first model minus its error under the second). Their mean is the gap between the two MAEs, and resampling the listings, the same resamples for both models, gives that gap a paired interval. Why. A listing that is hard to price tends to be hard for both models, so its difficulty lands in both of its errors. The difference cancels it listing by listing, and the paired interval measures how the two models differ rather than how hard the drawn listings were. When. Whenever two fixed models are compared on the same held-out rows, as Lessons 3.5.2 and 3.5.3 left the forest and boosting set against the linear pipeline.

The rest of the lesson goes on in is_real.py, under the median's block. Boosting first:

boosting = errors["boosting"]
gain = linear - boosting        # above 0 where boosting missed less
print(round(gain.mean()))       # 3788
# 329 0.548
print((gain > 0).sum(), round((gain > 0).mean(), 3))
paired = gain[picks].mean(axis=1)
low, high = np.percentile(paired, [2.5, 97.5])
print(round(low), round(high))  # 365 7303

gain is positive where the linear pipeline missed by more, so its mean, 3,788, is the linear MAE less boosting's, up to rounding, and boosting priced 329 of the 600 listings closer, 54.8 percent. picks holds the bootstrap block's 2,000 resamples, so every resample takes the same listings for both models, and the gain's paired interval runs from 365 to 7,303, which excludes 0. Each model's own interval, from the same resamples:

for name in ["linear", "boosting"]:
    resampled = errors[name][picks].mean(axis=1)
    print(name, np.percentile(resampled, [2.5, 97.5]).round())

linear [50970. 62951.]
boosting [47340. 59193.]

The two separate intervals overlap from 50,970 to 59,193, most of each, and overlap is not the test. Each separate interval carries in full the luck of how hard the drawn listings were, and the two models share that luck, because they priced the same listings. The run that takes the pairing away draws a second set of resamples for boosting, so that in every resample its MAE comes from a draw of listings made apart from the linear pipeline's, the comparison Lesson 3.1.2's DON'T ruled out for two models scored on different splits:

other = rng.integers(0, 600, size=(2000, 600))
apart = linear[picks].mean(axis=1) - boosting[other].mean(axis=1)
print(np.percentile(apart, [2.5, 97.5]).round())
print(round(np.corrcoef(linear, boosting)[0, 1], 3))

[-4566. 12484.]
0.829

The mean gain is unchanged, and the interval runs from -4,566 to 12,484, about two and a half times as wide, with 0 inside it. np.corrcoef (NumPy's correlation, returned as a 2 by 2 table whose entry [0, 1] is Lesson 2.7.1's Pearson r between its two arguments) puts the two models' errors at 0.829: a listing one model misses by a lot, the other mostly misses by a lot too. Pairing takes that shared part out of every difference before the average, and the paired interval, narrower by that much, excludes 0. The forest, against the same linear pipeline:

forest = errors["forest"]
forest_gain = linear - forest
forest_paired = forest_gain[picks].mean(axis=1)
low, high = np.percentile(forest_paired, [2.5, 97.5])
print(round(forest_gain.mean()), round(low), round(high))
print(round((forest_gain > 0).mean(), 3))

2760 -393 6089
0.508

The forest priced 50.8 percent of the listings closer, and its gain of 2,760 has a paired interval from -393 to 6,089, which holds 0. Lesson 3.5.2 left the forest's 0.009 of R² over the linear pipeline as one draw of 600 listings. In currency it is 2,760, and these 600 listings cannot tell it from no gain at all, which answers on the test rows the question Lesson 3.6.2 raised when it set the same two models side by side on five folds.

Both comparisons measure models fixed before the test rows scored them. A model kept because its test interval looked best would be chosen by the test score, Lesson 3.1.2's second break, and its interval would have been spent on the choice, as the search winner's score was in Lesson 3.6.3.

DON'T:

• DON'T read overlapping intervals as no difference. The linear pipeline's 50,970 to 62,951 and boosting's 47,340 to 59,193 overlap from 50,970 to 59,193, and the paired interval of the gain, 365 to 7,303, excludes 0.
• DON'T resample the two models' listings separately. The same 3,788 gain with each model on its own resamples spread from -4,566 to 12,484.
• DON'T pick the model to keep from these intervals. Comparing fixed models on the test rows is a measurement, and choosing by it is choosing by the test score.

Hypothesis Tests and What a p-Value Says

What. A hypothesis test starts from a null hypothesis (the claim a test assumes in order to measure surprise, here that the two models are equally good) and asks how often, if it held, a gap at least as large as the observed one would turn up. That share is the p-value (the probability, if the null hypothesis held, of a result at least as far from it as the one observed). Why. An interval lists the gains a sample is consistent with. A test asks about one gain, 0: whether luck alone, with two equally good models, makes a gap this size on 600 listings rare. When. When a report must say whether a gap measured on held-out rows could be the luck of the draw, with the threshold for calling a gap rare set before the test runs.

The sign-flip test (a test that redraws the sign of every listing's difference at random) takes the null hypothesis in its strict form: the two models are interchangeable on every listing, so each difference is as likely to point one way as the other, and a gain with its signs redrawn is one the null hypothesis could have produced. rng.choice, which Lesson 2.2.1's Assignment used to draw labels, draws 10,000 rows of 600 signs, from a fresh generator so that these draws do not depend on the blocks above:

flips = np.random.default_rng(42)
signs = flips.choice([-1, 1], size=(10_000, 600))
# 10,000 gains the null hypothesis could produce
null_gains = (signs * gain).mean(axis=1)
null_size = np.abs(null_gains)
print(round(np.percentile(null_size, 95)))  # 3484
p_flip = (null_size >= abs(gain.mean())).mean()
print(round(p_flip, 3))                     # 0.033

signs * gain is Lesson 2.2.4's broadcasting: the 600 differences multiply into every row of signs, so each row is the gain list with its signs redrawn, and its mean is one gain the null hypothesis could produce. One redrawn gain in 20 reaches 3,484 or more in size, and the observed 3,788 is beyond that: 3.3 percent of the 10,000 reached it, so p = 0.033. The test counts sizes on both sides, -3,788 as well as 3,788, because before the test either model could have come out ahead.

The paired t-test (a test of whether paired differences could have a mean of 0, which measures their mean in standard errors) is ttest_rel in scipy.stats, handed the two error arrays in listing order:

from scipy import stats

gain_se = gain.std(ddof=1) / np.sqrt(len(gain))
t_by_hand = gain.mean() / gain_se
print(round(gain_se), round(t_by_hand, 3))  # 1771 2.139
result = stats.ttest_rel(linear, boosting)
print(round(result.statistic, 3))           # 2.139
print(round(result.pvalue, 3))              # 0.033

The gain's own standard error is 1,771, smaller than the linear MAE's 3,030, because the pairing takes out the difficulty the two models share, and the gain is 2.139 of them, the test's t statistic. SciPy turns 2.139 into a p-value with the t distribution (the spread a mean, measured in its own standard errors, would show across samples if the null hypothesis held), 0.033, the sign-flip test's figure to three decimals. Its null hypothesis is that the mean difference is 0, the sign-flip test's in a weaker form. The forest's gain, through both tests:

forest_null = (signs * forest_gain).mean(axis=1)
forest_p = (np.abs(forest_null) >= abs(forest_gain.mean())).mean()
forest_t = stats.ttest_rel(linear, forest)
# 0.095 0.092
print(round(forest_p, 3), round(forest_t.pvalue, 3))

The forest's gain gets p = 0.095 by sign flips and 0.092 by ttest_rel. A p below 0.05 is called statistically significant, and 0.05 is a widely shared convention, set before the test runs and never a property of the models: a line drawn after the results sits wherever they fell. Boosting's 0.033 falls below it and the forest's 0.092 does not, and that is all the line says.

What each p says is narrow. p = 0.033 says that if the two models were equally good, a gain at least 3,788 in size, either way, would turn up in about 3.3 percent of samples of 600 listings like these. It is computed by assuming the null hypothesis, so it cannot be the probability that the null hypothesis is true. The forest's 0.092 does not show the forest is no better: its paired interval holds gains up to 6,089 as well as 0, and a large p says only that these listings cannot tell. Neither p says how large a gain is, which is the interval's job.

DON'T:

• DON'T read p = 0.033 as a 3.3 percent chance that the two models are equally good. It is the share of redrawn gains at least 3,788 in size, computed by assuming they are.
• DON'T read the forest's p of 0.092 as proof that the forest is no better. Its paired interval runs from -393 to 6,089.
• DON'T read a small p as a large gain. Boosting's 0.033 sits beside a paired interval that starts at 365.

Many Comparisons

What. A test of a difference that is not there comes out below 0.05 about one time in 20. Run several such tests and the chance that at least one falls below 0.05 grows with their number, to 1 - 0.95 ** m for m independent tests. Why. A report that runs many comparisons and quotes the ones below 0.05 quotes luck along with any real gain, the way the search winner of Lesson 3.6.3 carried luck in its score. When. Whenever more than one test stands behind a result, across models, metrics, neighborhoods or settings, this lesson's two tests included.

ttest_1samp (the one-sample t-test, which asks whether one set of values could have a mean of 0), handed the paired differences, is ttest_rel by another route. Subtracting the gain's mean from every listing's difference keeps the listings' spread and removes the gain, so a test on a resample of the result tests a difference that is not there:

one_sample = stats.ttest_1samp(gain, 0)
print(round(one_sample.pvalue, 3))     # 0.033
centered = gain - gain.mean()          # the gain removed
print(round(abs(centered.mean()), 6))  # 0.0

Each batch below is 20 tests on 20 resamples of the centered differences, 1,000 batches in all, with axis=1 running one test per row:

tests = np.random.default_rng(42)
below = any_below = 0
for batch in range(1000):
    # 20 resamples of the centered differences, one test each
    rows = tests.integers(0, 600, size=(20, 600))
    p = stats.ttest_1samp(centered[rows], 0, axis=1).pvalue
    if batch == 0:
        first = np.sort(p)
    below += (p < 0.05).sum()
    any_below += (p < 0.05).any()
print(first[:3].round(3))
print(below / 20_000, any_below / 1000)
print(round(1 - 0.95 ** 20, 3))

[0.028 0.08  0.172]
0.0506 0.639
0.642

5.06 percent of the 20,000 tests fell below 0.05, the one in 20 the line allows, and 639 of the 1,000 batches held at least one, close to 1 - 0.95 ** 20, 0.642, the chance for 20 independent tests. The first batch's smallest p, 0.028, is below the boosting gain's 0.033, on a gain that is not there.

A report that runs many comparisons and quotes the one below 0.05 repeats the spent winner of Lesson 3.6.3: among many tests, luck alone hands out small p-values. This lesson ran two tests on the same 600 listings, boosting's and the forest's, and each carried the same one in 20 chance of a p below 0.05 with nothing to find.

DON'T:

• DON'T report the one comparison below 0.05 out of many. With no effect at all, 639 of 1,000 batches of 20 tests held at least one.
• DON'T read the smallest p of a batch as a finding. The first batch of 20 tests of a gain that is not there gave 0.028, below boosting's 0.033.

Lab Exercises

Lab Exercise 1: The Forest's Interval, and Where the Interval Falls Short

In forest_interval.py in ~/dev/bootcamp, importing errors from the Theory's held_out_errors.py, print the forest's test MAE, the MAEs of its six slices of 100 test listings, its standard error, the normal approximation's 95 percent interval, and the percentile interval and the standard deviation of 2,000 bootstrap MAEs drawn with np.random.default_rng(42).

Then add lines to the end of the Theory's coverage.py that print coverage of the linear errors at samples of 100 listings, the six largest linear errors, the size and the MAE of a market that is the linear errors without those six, and that market's coverage at samples of 600 and of 100. End with a comment saying where the interval's shortfall comes from.

Lab Exercise 2: Forest Against Boosting, and Batches of Every Size

In forest_or_boosting.py, importing errors, compare the forest with boosting on the same 600 test listings, each difference taken as the forest's error minus boosting's. Print the mean difference and how many listings boosting priced closer, then, from 2,000 resamples drawn with np.random.default_rng(42), the paired 95 percent interval and each model's own interval, then the sign-flip p from 10,000 rows of signs drawn with a fresh np.random.default_rng(42), and ttest_rel's p.

Then center the differences and, with one more fresh np.random.default_rng(42), run 1,000 batches each of 1, 5, 20 and 50 tests, printing for each size the share of batches holding at least one p below 0.05 beside 1 - 0.95 ** m. End with a comment on what a report may say about the two models.

import numpy as np
from held_out_errors import errors

forest = errors["forest"]
print(round(forest.mean()))   # 53802
print([round(mae) for mae in forest.reshape(6, 100).mean(axis=1)])
standard_error = forest.std(ddof=1) / np.sqrt(len(forest))
margin = 1.96 * standard_error
print(round(standard_error))  # 3024
print(round(forest.mean() - margin), round(forest.mean() + margin))
rng = np.random.default_rng(42)
picks = rng.integers(0, 600, size=(2000, 600))
boot_mae = forest[picks].mean(axis=1)
low, high = np.percentile(boot_mae, [2.5, 97.5])
print(round(low), round(high), round(boot_mae.std()))

[58188, 63270, 48811, 51307, 50326, 50908]
47875 59728
48402 60067 2992

market = errors["linear"]
print(coverage(market, 100))
print(np.sort(market)[-6:].round())
without_six = np.sort(market)[:-6]
print(len(without_six), round(without_six.mean()))
for n in [600, 100]:
    print(n, coverage(without_six, n))
# Short of 950 with every error, 927 at 600 listings and 882 at 100,
# and near it, 944 and 945, once the six largest are gone: the
# shortfall comes with the six listings missed by 274,431 to 986,058,
# and it is largest on the smaller samples.

(927, 810)
(882, 764)
[274431. 279776. 358208. 729700. 796591. 986058.]
594 51367
600 (944, 825)
100 (945, 846)

import numpy as np
from scipy import stats
from held_out_errors import errors

forest, boosting = errors["forest"], errors["boosting"]
gain = forest - boosting                     # above 0: boosting closer
print(round(gain.mean()), (gain > 0).sum())  # 1028 313
rng = np.random.default_rng(42)
picks = rng.integers(0, 600, size=(2000, 600))
for values in [gain, forest, boosting]:
    resampled = values[picks].mean(axis=1)
    print(np.percentile(resampled, [2.5, 97.5]).round())
flips = np.random.default_rng(42)
signs = flips.choice([-1, 1], size=(10_000, 600))
null_gains = (signs * gain).mean(axis=1)
p_flip = (np.abs(null_gains) >= abs(gain.mean())).mean()
print(round(p_flip, 3))                      # 0.386
# 0.382
print(round(stats.ttest_rel(forest, boosting).pvalue, 3))

[-1150.  3307.]
[48402. 60067.]
[47340. 59193.]

centered = gain - gain.mean()
tests = np.random.default_rng(42)
for m in [1, 5, 20, 50]:
    any_below = 0
    for _ in range(1000):
        rows = tests.integers(0, 600, size=(m, 600))
        p = stats.ttest_1samp(centered[rows], 0, axis=1).pvalue
        any_below += (p < 0.05).any()
    print(m, any_below / 1000, round(1 - 0.95 ** m, 3))
# Boosting priced 313 of the 600 listings closer, and its gain of
# 1,028 has a paired interval from -1,150 to 3,307 and a p of 0.382:
# these listings cannot tell the forest from boosting, and the report
# says that, never that the two are equally good.

1 0.049 0.05
5 0.196 0.226
20 0.636 0.642
50 0.919 0.923

Theory Assessment

Q1: Six Slices and One Standard Error

The linear pipeline's six slices of 100 test listings gave MAEs from 48,252 to 65,649, while its MAE over all 600 is 56,562 with a standard error of 3,030. Explain what the standard error estimates, why the slices spread far wider than 3,030 would suggest, and how many test listings with errors like these would give a standard error half as large. Then explain why the bootstrap's 3,008 and the formula's 3,030 answer one question.

Hint: Ask what each slice's MAE is a mean of, and where the count sits in the formula. For the last part, ask what each of the two numbers is the standard deviation of.

Q2: What the 95 Belongs To

A teammate writes in a report: "The linear pipeline's true MAE lies between 50,970 and 62,951 with probability 0.95, so the next test set's MAE will land in that range 95 times in 100." Using coverage.py's 927 and 810, explain what is wrong with each half of the sentence and what the 95 belongs to. The teammate adds that boosting's interval, 47,340 to 59,193, overlaps the linear pipeline's, so boosting's 3,788 lower MAE could be luck. Explain why the overlap cannot settle that, and which interval does.

Hint: Once an interval is computed, ask what is still random about whether it holds the true MAE. For the overlap, ask what each separate interval carries that a listing's own difference cancels.

Q3: A p of 0.092, and One Test in Twenty

Explain what boosting's p of 0.033 says, using the sign flips that produced it, and why it is not the probability that the two models are equally good. Then explain what the forest's p of 0.092 does and does not show about the forest. Last, a teammate tests the forest against the linear pipeline separately in each of 20 groups of listings, finds one group at p = 0.04, and reports it. Explain, from the Theory's batches, why that one p is weak evidence.

Hint: Ask what the flips assumed before they counted anything. For the 20 groups, ask how often a batch of 20 tests of a gain that is not there held one below 0.05.

Q4: Three Tests of One Gain

A teammate tests boosting's gain again. stats.ttest_ind, in scipy.stats, is the unpaired t-test: it treats its two arrays as two separate samples and compares their means against the spread within each array, so which value sits at which position plays no part. The teammate writes three_tests.py:

from scipy import stats
from held_out_errors import errors

linear, boosting = errors["linear"], errors["boosting"]
print(round((linear - boosting).mean()))
print(round(stats.ttest_rel(linear, boosting).pvalue, 3))
print(round(stats.ttest_ind(linear, boosting).pvalue, 3))
print(round(stats.ttest_rel(linear, boosting[::-1]).pvalue, 3))
print(round((linear - boosting[::-1]).mean()))

Read the script without running it. Say what each print shows, with this lesson's numbers where it printed them and a size where it did not, and explain why one gain gets three p-values, which of them answers whether boosting's gap over the linear pipeline on these 600 listings could be the luck of the draw, and what boosting[::-1] did to the pairs.

Hint: [::-1] reverses an array, so its position 0 holds the last test listing's error. Ask which call keeps each listing's two errors together, and recall the Theory's interval with the pairing taken away.

3788
0.033
0.377
0.374
3788

Assignment

Assignment: Three Models, Three Comparisons, One Sentence for the Report

The team's report ranks the three models by test MAE and says boosting beats the linear pipeline, p = 0.033. The reviewer points out that the report compared all three pairs on the same 600 listings, so three tests stand behind that sentence, and asks for all three with the count of tests taken into account.

Holm's method adjusts a set of p-values for the number of tests behind them. Sort the m p-values from smallest to largest, multiply the smallest by m, the next by m - 1 and so on down to the largest times 1, then walk up the sorted list keeping each adjusted value at least as large as the one before it, and cap every value at 1. A comparison is claimed at the 0.05 convention when its adjusted p is below 0.05, and the method is built to keep the chance that luck alone pushes any adjusted p below 0.05 at 0.05 or less, however many tests stand behind the result.

Write three_pairs.py. For the pairs linear and forest, linear and boosting, and forest and boosting, each difference taken as the first model's error minus the second's, build one DataFrame with a row per pair: the mean paired difference, how many listings the second model priced closer, the paired 95 percent interval from 2,000 resamples drawn with np.random.default_rng(42), the same resamples for every pair, ttest_rel's p, and the Holm-adjusted p. Then run the Theory's 1,000 batches of 20 tests of a gain that is not there again and print the share of batches holding a p below 0.05 before and after Holm's adjustment. End with comments saying which orderings of the three models the report may claim at the 0.05 convention, and the sentence the reviewer should write in place of "boosting beats the linear pipeline, p = 0.033".

Hints:

• report["p"].sort_values().index lists the pairs from the smallest p to the largest, and enumerate counts them from 0, so the pair at count k is multiplied by m - k. A dictionary of adjusted values, turned into a Series, lines up with the table by its labels.
• The smallest adjusted p is m times the smallest p, capped at 1, and no other adjusted p is smaller, so a batch holds an adjusted p below 0.05 exactly when its smallest p times 20 is below 0.05.

import numpy as np
import pandas as pd
from scipy import stats
from held_out_errors import errors

pairs = [("linear", "forest"), ("linear", "boosting"), ("forest", "boosting")]
rng = np.random.default_rng(42)
picks = rng.integers(0, 600, size=(2000, 600))
names, records = [], []
for first, second in pairs:
    gain = errors[first] - errors[second]  # above 0: second missed less
    low, high = np.percentile(gain[picks].mean(axis=1), [2.5, 97.5])
    names.append(f"{first}-{second}")
    records.append({
        "gain": round(gain.mean()),
        "closer": (gain > 0).sum(),
        "low": round(low),
        "high": round(high),
        "p": stats.ttest_rel(errors[first], errors[second]).pvalue,
    })
report = pd.DataFrame(records, index=names)

adjusted = {}
running = 0.0
for k, name in enumerate(report["p"].sort_values().index):
    # the k-th smallest p times the tests left, never below the last
    running = max(running, (len(report) - k) * report.loc[name, "p"])
    adjusted[name] = min(running, 1.0)
report["holm"] = pd.Series(adjusted)
print(report.round(3))

                 gain  closer   low  high      p   holm
linear-forest    2760     305  -393  6089  0.092  0.185
linear-boosting  3788     329   365  7303  0.033  0.099
forest-boosting  1028     313 -1150  3307  0.382  0.382

gain = errors["linear"] - errors["boosting"]
centered = gain - gain.mean()
tests = np.random.default_rng(42)
any_raw = any_holm = 0
for _ in range(1000):
    rows = tests.integers(0, 600, size=(20, 600))
    p = stats.ttest_1samp(centered[rows], 0, axis=1).pvalue
    any_raw += (p < 0.05).any()
    any_holm += p.min() * 20 < 0.05
print(any_raw / 1000, any_holm / 1000)  # 0.639 0.056
# Alone, boosting's gain over the linear pipeline has a paired
# interval that excludes 0 and a p below 0.05. Counted with the
# report's three comparisons its p is 0.099, and no ordering of the
# three models clears the 0.05 convention. On the Theory's batches
# of 20 tests of nothing, Holm cut the share of batches holding a p
# below 0.05 from 0.639 to 0.056, within the noise of 1,000 batches
# of the 0.05 it is built not to exceed.
# The reviewer writes: on these 600 test listings boosting missed by
# 3,788 less than the linear pipeline on average, 95 percent paired
# interval 365 to 7,303, p = 0.033 for that pair alone and 0.099
# with the report's three comparisons counted, so the test rows do
# not establish which of the three models prices best.

Additional Resources

Optional, and nothing later depends on either. Worth opening if you want another way to adjust many p-values at once, or the same comparison run on cross-validation folds.

• YouTube - "False Discovery Rates, FDR, clearly explained" by StatQuest with Josh Starmer. The Benjamini-Hochberg method, which adjusts the p-values of many tests at once to hold down the false discovery rate, the share of the results called significant that are false positives.
• Official docs - scikit-learn 1.9 example, "Statistical comparison of models using grid search". Models compared on the same cross-validation folds, why their fold scores are correlated, a paired t-test whose variance is corrected because those scores are not independent, a Bonferroni correction for comparing every pair, and a Bayesian alternative.

Summary

A test score turned out to belong as much to the 600 listings that happened to be held out as to the model, and the standard error and the bootstrap were two ways of asking those same listings how much of the number was theirs. Neither needed a new split or a refit. The errors held their own spread, and resampling them was enough to watch a score move by thousands with nothing about the model changed.

The interval's 95 belonged to the procedure from the start, which is why the coverage run counted a thousand intervals instead of reading one. A single range holds the true MAE or misses it, and the honest sentence about it names the procedure that built it. The same care decided the comparison of two models: a listing hard for one model is mostly hard for the other, and only a comparison that keeps each listing's two errors together sees past what they share.

A test answers a narrow question, whether luck alone, under an assumption stated in advance, makes a gap this size rare, and almost every misreading of it is a sentence written around the number. A p below 0.05 is not a probability about the models, a p above it is not a finding of equality, and among many tests luck alone hands out small p-values, which is why a report counts its tests before it quotes one.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

from scipy import stats

err = errors["linear"]                    # 600 absolute errors
se = err.std(ddof=1) / np.sqrt(len(err))  # standard error
rng = np.random.default_rng(42)
picks = rng.integers(0, 600, size=(2000, 600))
boot = err[picks].mean(axis=1)            # 2,000 resampled MAEs
np.percentile(boot, [2.5, 97.5])          # 95 percent interval
gain = err - errors["boosting"]           # paired differences
np.percentile(gain[picks].mean(axis=1), [2.5, 97.5])
stats.ttest_rel(err, errors["boosting"]).pvalue

signs = rng.choice([-1, 1], size=(10_000, 600))
null = np.abs((signs * gain).mean(axis=1))
p_flip = (null >= abs(gain.mean())).mean()  # sign-flip p