Hands-on lab — Day 117: Sampling and the Central Limit Theorem

"""Exercise 1 -- the sampling distribution of the mean.

Draw 20,000 independent samples of size 40 from a skewed population, compute
each sample's mean, and look at the DISTRIBUTION of those 20,000 means -- not
any single one of them. Its own mean should sit close to the population mean,
and its own spread should sit close to sigma / sqrt(n).
"""

import numpy as np

import dataset as D
from sampling import population_mean_std, sampling_distribution, theoretical_standard_error

rng = np.random.default_rng(1)

pop_mean, pop_sigma = population_mean_std(D.SKEWED_POP)
means = sampling_distribution(D.SKEWED_POP, D.EX1_N, D.EX1_TRIALS, rng)

measured_mean = means.mean()
measured_se = means.std(ddof=1)
theoretical_se = theoretical_standard_error(pop_sigma, D.EX1_N)

# The standard error OF the measured mean itself, so "close" below is a
# statement with a real yardstick rather than a guessed number: with 20,000
# trials, the mean of the sampling distribution should sit within a few of
# these of the true population mean almost always.
se_of_measured_mean = measured_se / np.sqrt(D.EX1_TRIALS)
# The sample standard deviation of 20,000 draws has its own standard error,
# approximately sigma / sqrt(2 * trials) for a roughly-normal statistic.
se_of_measured_se = measured_se / np.sqrt(2 * D.EX1_TRIALS)

print(f"population mean = {pop_mean:.4f}, population sigma = {pop_sigma:.4f}")
print(f"n = {D.EX1_N}, trials = {D.EX1_TRIALS}")
print(f"measured mean of the sampling distribution = {measured_mean:.4f}")
print(f"measured standard error (std of the {D.EX1_TRIALS} sample means) = {measured_se:.4f}")
print(f"theoretical standard error (sigma / sqrt(n)) = {theoretical_se:.4f}")

mean_gap_in_ses = abs(measured_mean - pop_mean) / se_of_measured_mean
se_gap_in_ses = abs(measured_se - theoretical_se) / se_of_measured_se
print(f"gap between measured and population mean = {mean_gap_in_ses:.2f} standard errors")
print(f"gap between measured and theoretical SE = {se_gap_in_ses:.2f} standard errors")

assert mean_gap_in_ses < 3.0, "the sampling distribution's own mean drifted too far from the population mean"
assert se_gap_in_ses < 3.0, "the sampling distribution's own spread drifted too far from sigma / sqrt(n)"

print("01_sampling_distribution.py: every assertion held.")

"""Exercise 2 -- the sqrt(n) law.

Four sample sizes, each exactly 4x the one before it: 10, 40, 160, 640. Every
time n quadruples, the standard error of the mean should roughly HALVE, not
quarter -- that is the entire economics of measurement in one number. This
asserts the RATIOS between successive standard errors, not four hard-coded
values, so the check is honest about what simulation noise can and cannot
pin down.
"""

import math

import numpy as np

import dataset as D
from sampling import sampling_distribution

rng = np.random.default_rng(2)

standard_errors = {}
for n in D.SQRT_N_LAW_NS:
    means = sampling_distribution(D.SKEWED_POP, n, D.SQRT_N_LAW_TRIALS, rng)
    standard_errors[n] = means.std(ddof=1)
    print(f"n = {n:>4}  measured SE = {standard_errors[n]:.4f}")

print()
ratios = []
for smaller, larger in zip(D.SQRT_N_LAW_NS, D.SQRT_N_LAW_NS[1:]):
    ratio = standard_errors[smaller] / standard_errors[larger]
    ratios.append(ratio)
    print(f"SE(n={smaller}) / SE(n={larger}) = {ratio:.3f}  (expected: sqrt(4) = {math.sqrt(4):.3f})")

for (smaller, larger), ratio in zip(zip(D.SQRT_N_LAW_NS, D.SQRT_N_LAW_NS[1:]), ratios):
    gap = abs(ratio - 2.0)
    assert gap < D.SQRT_N_LAW_RATIO_TOLERANCE, (
        f"SE(n={smaller})/SE(n={larger}) = {ratio:.3f} strayed more than "
        f"{D.SQRT_N_LAW_RATIO_TOLERANCE} from the predicted 2.0"
    )

# The compounded law across all three quadruplings: from n=10 to n=640 is a
# 64x growth in sample size, so the standard error should have fallen by
# close to sqrt(64) = 8x, not 64x -- the mistake a 1/n law would predict.
overall_ratio = standard_errors[D.SQRT_N_LAW_NS[0]] / standard_errors[D.SQRT_N_LAW_NS[-1]]
print(f"\noverall SE(n=10) / SE(n=640) = {overall_ratio:.2f}  (expected: sqrt(64) = {math.sqrt(64):.2f})")
assert abs(overall_ratio - 8.0) < 1.5, "the compounded sqrt(n) law did not hold across the full range"
assert overall_ratio < 20.0, "a ratio anywhere near 64 would mean the error fell like 1/n, not 1/sqrt(n)"

print("02_the_sqrt_n_law.py: every assertion held.")

"""Exercise 3 -- the central limit theorem, measured rather than assumed.

The population itself is heavily right-skewed (an Exponential shape). As n
grows, the SKEWNESS of the sampling distribution of the mean should fall
toward zero, monotonically -- the population's lopsidedness washing out of
the statistic even though it never leaves the population itself.

Two more populations that look nothing like a bell -- a biased coin and a
lumpy two-spike distribution -- are reported alongside, to show the same
flattening is not a fact about this one population's shape.
"""

import numpy as np

import dataset as D
from sampling import population_mean_std, sampling_distribution, skewness

rng = np.random.default_rng(3)

pop_mean, pop_sigma = population_mean_std(D.SKEWED_POP)
pop_skew = skewness(D.SKEWED_POP)
print(f"population (Exponential-shaped): mean = {pop_mean:.3f}, skewness = {pop_skew:.3f}")
print()

skews = []
for n in D.SKEW_DEMO_NS:
    means = sampling_distribution(D.SKEWED_POP, n, D.SKEW_DEMO_TRIALS, rng)
    s = skewness(means)
    skews.append(s)
    print(f"n = {n:>4}  skewness of the sampling distribution = {s:.4f}")

print()
for smaller_n, larger_n, s_small, s_large in zip(
    D.SKEW_DEMO_NS, D.SKEW_DEMO_NS[1:], skews, skews[1:]
):
    assert s_large < s_small, (
        f"skewness did not decrease going from n={smaller_n} ({s_small:.4f}) "
        f"to n={larger_n} ({s_large:.4f})"
    )
print(f"skewness fell monotonically across n = {D.SKEW_DEMO_NS}: "
      + " > ".join(f"{s:.3f}" for s in skews))

# Two more non-bell populations, reported for the lesson's narrative. Both
# show the same qualitative flattening; neither is asserted to the same
# precision as the primary skewed population above, since a Bernoulli
# population's sampling-distribution skewness is much noisier at small n.
for label, population in (("biased coin", D.COIN_POP), ("two-spike", D.TWO_SPIKE_POP)):
    p_skew = skewness(population)
    small_n_skew = skewness(sampling_distribution(population, 5, D.SKEW_DEMO_TRIALS, rng))
    large_n_skew = skewness(sampling_distribution(population, 320, D.SKEW_DEMO_TRIALS, rng))
    print(
        f"{label}: population skewness = {p_skew:.3f}, "
        f"sampling-distribution skewness at n=5 -> {small_n_skew:.3f}, "
        f"at n=320 -> {large_n_skew:.3f}"
    )
    assert abs(large_n_skew) < abs(small_n_skew), (
        f"{label}: skewness at n=320 was not smaller in magnitude than at n=5"
    )

print("03_clt_from_a_skewed_population.py: every assertion held.")

"""Exercise 4 -- where the central limit theorem actually fails.

The Exponential distribution has a finite variance, so its sample mean's
spread should shrink by close to sqrt(100) = 10x when n grows from 10 to
1,000. The Cauchy distribution has NO defined mean or variance -- its tails
are too heavy -- and the mean of n Cauchy draws is itself standard-Cauchy
distributed, for every n. Averaging a thousand of them should be no better
than averaging ten.

The interquartile range (IQR) is used as the spread measure here rather than
the standard deviation, and that choice is load-bearing: a Cauchy sample's
standard deviation is not an estimate of anything, because the population
quantity it would estimate does not exist. The IQR depends only on the order
of the data, so it stays meaningful even when the mean and variance do not.
"""

import dataset as D
from sampling import cauchy_mean_iqr, exponential_mean_iqr
import numpy as np

rng = np.random.default_rng(4)

exp_iqr_small = exponential_mean_iqr(D.CAUCHY_DEMO_N_SMALL, D.CAUCHY_DEMO_TRIALS, rng, D.EXPONENTIAL_SCALE)
exp_iqr_large = exponential_mean_iqr(D.CAUCHY_DEMO_N_LARGE, D.CAUCHY_DEMO_TRIALS, rng, D.EXPONENTIAL_SCALE)
exp_ratio = exp_iqr_small / exp_iqr_large

cauchy_iqr_small = cauchy_mean_iqr(D.CAUCHY_DEMO_N_SMALL, D.CAUCHY_DEMO_TRIALS, rng)
cauchy_iqr_large = cauchy_mean_iqr(D.CAUCHY_DEMO_N_LARGE, D.CAUCHY_DEMO_TRIALS, rng)
cauchy_ratio = cauchy_iqr_small / cauchy_iqr_large

print(f"Exponential(scale={D.EXPONENTIAL_SCALE}): IQR of the mean at n={D.CAUCHY_DEMO_N_SMALL} = "
      f"{exp_iqr_small:.4f}, at n={D.CAUCHY_DEMO_N_LARGE} = {exp_iqr_large:.4f}")
print(f"  ratio = {exp_ratio:.2f}  (100x more data, expected shrink ~ sqrt(100) = 10x)")
print()
print(f"standard Cauchy: IQR of the mean at n={D.CAUCHY_DEMO_N_SMALL} = "
      f"{cauchy_iqr_small:.4f}, at n={D.CAUCHY_DEMO_N_LARGE} = {cauchy_iqr_large:.4f}")
print(f"  ratio = {cauchy_ratio:.2f}  (100x more data, expected shrink: NONE)")

assert exp_ratio > D.EXPONENTIAL_SHRINK_FLOOR, (
    f"the Exponential mean's spread shrank by only {exp_ratio:.2f}x, "
    f"expected close to 10x"
)
assert D.CAUCHY_NO_SHRINK_LOW < cauchy_ratio < D.CAUCHY_NO_SHRINK_HIGH, (
    f"the Cauchy mean's spread shrank by {cauchy_ratio:.2f}x -- it should have stayed "
    f"roughly flat, not shrunk toward the Exponential's ~10x"
)
assert cauchy_ratio < exp_ratio / 3.0, (
    "the Cauchy ratio was not clearly smaller than the Exponential ratio -- "
    "the whole point of this exercise is the contrast between the two"
)

print("04_the_cauchy_counterexample.py: every assertion held.")

"""Exercise 5 -- sampling bias is not sampling error, and n does not fix it.

An UNBIASED sampler draws from the whole population; a BIASED sampler draws
only from the half of the population strictly above the population's own
median, no matter how many draws it takes. Both errors -- mean absolute
distance from the true population mean -- are tracked as the sample size
grows 100x, from 30 to 3,000.

The unbiased sampler's error should shrink by close to 10x, exactly the
sqrt(n) law from exercise 2. The biased sampler's error should stay roughly
FLAT: more data buys it a more precise estimate of the wrong number, and the
mathematics does not know the difference between "confident" and "correct".
"""

import dataset as D
from sampling import biased_pool, mean_absolute_error, sampling_distribution
import numpy as np

rng = np.random.default_rng(5)

true_mean = float(D.SKEWED_POP.mean())
pool = biased_pool(D.SKEWED_POP)
print(f"true population mean = {true_mean:.4f}")
print(f"biased pool (values above the population median) mean = {pool.mean():.4f}  "
      f"-- this is what the biased sampler converges to, not the true mean")
print()

unbiased_small = mean_absolute_error(
    sampling_distribution(D.SKEWED_POP, D.BIAS_DEMO_N_SMALL, D.BIAS_DEMO_TRIALS, rng), true_mean
)
unbiased_large = mean_absolute_error(
    sampling_distribution(D.SKEWED_POP, D.BIAS_DEMO_N_LARGE, D.BIAS_DEMO_TRIALS, rng), true_mean
)
biased_small = mean_absolute_error(
    sampling_distribution(pool, D.BIAS_DEMO_N_SMALL, D.BIAS_DEMO_TRIALS, rng), true_mean
)
biased_large = mean_absolute_error(
    sampling_distribution(pool, D.BIAS_DEMO_N_LARGE, D.BIAS_DEMO_TRIALS, rng), true_mean
)

unbiased_ratio = unbiased_small / unbiased_large
biased_ratio = biased_small / biased_large

print(f"UNBIASED sampler: mean abs error at n={D.BIAS_DEMO_N_SMALL} = {unbiased_small:.4f}, "
      f"at n={D.BIAS_DEMO_N_LARGE} = {unbiased_large:.4f}  (ratio = {unbiased_ratio:.2f})")
print(f"BIASED sampler:   mean abs error at n={D.BIAS_DEMO_N_SMALL} = {biased_small:.4f}, "
      f"at n={D.BIAS_DEMO_N_LARGE} = {biased_large:.4f}  (ratio = {biased_ratio:.2f})")

assert unbiased_ratio > D.UNBIASED_SHRINK_FLOOR, (
    f"the unbiased sampler's error only shrank by {unbiased_ratio:.2f}x, expected close to 10x"
)
assert D.BIASED_FLAT_LOW < biased_ratio < D.BIASED_FLAT_HIGH, (
    f"the biased sampler's error changed by {biased_ratio:.2f}x -- it should have stayed roughly flat"
)
assert biased_large > 3.0 * unbiased_large, (
    "at the LARGE sample size the biased sampler's error should still dwarf the unbiased "
    "sampler's -- more data did not rescue it"
)

print("05_bias_does_not_shrink.py: every assertion held.")

"""Exercise 6 -- the bootstrap, from scratch.

Resample a single dataset with replacement, recompute a statistic on every
resample, and read the statistic's standard error straight off the spread of
the results -- no formula for the statistic's own sampling distribution
required.

First check it against the MEAN, where a formula (sigma_hat / sqrt(n)) does
exist, so the bootstrap can be judged against a known answer. Then apply the
exact same code to the MEDIAN, where no simple closed form exists, and check
its answer for sanity against the spread of medians computed from genuinely
fresh, independent samples of the same population -- the bootstrap's whole
reason for existing.
"""

import dataset as D
from sampling import bootstrap_standard_error
import numpy as np

rng = np.random.default_rng(6)

sample = rng.normal(loc=D.BOOTSTRAP_SAMPLE_MEAN, scale=D.BOOTSTRAP_SAMPLE_STD, size=D.BOOTSTRAP_SAMPLE_SIZE)
sigma_hat = sample.std(ddof=1)
theoretical_se = sigma_hat / np.sqrt(D.BOOTSTRAP_SAMPLE_SIZE)

boot_se_mean = bootstrap_standard_error(
    sample, lambda a: a.mean(axis=1), D.BOOTSTRAP_N_BOOT, rng
)
boot_se_median = bootstrap_standard_error(
    sample, lambda a: np.median(a, axis=1), D.BOOTSTRAP_N_BOOT, rng
)

# "Genuinely fresh samples", drawn straight from the same Normal population
# rather than resampled from the one dataset above -- this is the sanity
# check the median bootstrap is measured against, since no formula exists.
fresh = rng.normal(
    loc=D.BOOTSTRAP_SAMPLE_MEAN,
    scale=D.BOOTSTRAP_SAMPLE_STD,
    size=(D.FRESH_MEDIAN_REPLICATIONS, D.BOOTSTRAP_SAMPLE_SIZE),
)
fresh_median_se = np.median(fresh, axis=1).std(ddof=1)

relative_error_mean = abs(boot_se_mean - theoretical_se) / theoretical_se
median_ratio = boot_se_median / fresh_median_se

print(f"sample size = {D.BOOTSTRAP_SAMPLE_SIZE}, sigma_hat = {sigma_hat:.4f}")
print(f"theoretical SE of the mean (sigma_hat / sqrt(n)) = {theoretical_se:.4f}")
print(f"bootstrap SE of the mean ({D.BOOTSTRAP_N_BOOT} resamples) = {boot_se_mean:.4f}")
print(f"  relative error = {relative_error_mean:.3%}")
print()
print(f"bootstrap SE of the MEDIAN ({D.BOOTSTRAP_N_BOOT} resamples) = {boot_se_median:.4f}")
print(f"SE of the median from {D.FRESH_MEDIAN_REPLICATIONS} genuinely fresh samples = {fresh_median_se:.4f}")
print(f"  ratio (bootstrap / fresh) = {median_ratio:.2f}")

assert relative_error_mean < D.BOOTSTRAP_MEAN_RELATIVE_TOLERANCE, (
    f"the bootstrap SE of the mean disagreed with sigma_hat / sqrt(n) by "
    f"{relative_error_mean:.1%}, expected under {D.BOOTSTRAP_MEAN_RELATIVE_TOLERANCE:.0%}"
)
assert D.BOOTSTRAP_MEDIAN_SANITY_LOW < median_ratio < D.BOOTSTRAP_MEDIAN_SANITY_HIGH, (
    f"the bootstrap SE of the median ({boot_se_median:.4f}) is not within a sane "
    f"range of the fresh-sample benchmark ({fresh_median_se:.4f})"
)

print("06_bootstrap_from_scratch.py: every assertion held.")

"""Exercise 7 -- dependence inflates the true standard error, quietly.

Generate an autocorrelated AR(1) series -- each observation is 0.7 times the
previous one plus fresh noise, so consecutive observations are far from
independent. Measure the TRUE standard error of its sample mean the only
honest way: generate many independent replications of the whole series and
look at the spread of their means. Compare that against the NAIVE standard
error, sample_std / sqrt(n), which assumes independence the data does not
have.

The naive formula should understate the true standard error -- meaningfully,
not by a rounding error -- which is worse than an honest formula being
merely imprecise: it makes the analyst confident in exact proportion to how
wrong they are.
"""

import dataset as D
from sampling import ar1_series, naive_standard_error, true_standard_error_by_replication
import numpy as np

rng = np.random.default_rng(7)

true_se = true_standard_error_by_replication(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, D.AR1_REPLICATIONS, rng)

# The naive SE is reported as an average over many single series, since any
# one series' sample standard deviation is itself noisy.
naive_ses = [
    naive_standard_error(ar1_series(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, rng))
    for _ in range(500)
]
naive_se_avg = float(np.mean(naive_ses))

ratio = true_se / naive_se_avg

print(f"AR(1) series: n = {D.AR1_N}, phi = {D.AR1_PHI}, sigma = {D.AR1_SIGMA}")
print(f"TRUE standard error (from {D.AR1_REPLICATIONS} independent replications) = {true_se:.4f}")
print(f"NAIVE standard error (sample_std / sqrt(n), averaged over 500 series) = {naive_se_avg:.4f}")
print(f"ratio (true / naive) = {ratio:.2f}")
print()
print("An analyst who trusted the naive formula here would report a confidence "
      "interval that is too narrow -- not slightly, but by a factor that a "
      "single glance at the ratio above makes obvious.")

assert ratio > D.AR1_INFLATION_FLOOR, (
    f"the naive SE was not meaningfully smaller than the true SE (ratio = {ratio:.2f}, "
    f"expected above {D.AR1_INFLATION_FLOOR})"
)

print("07_dependence_inflates_se.py: every assertion held.")

"""Exercise 8 -- the evaluation-margin calculation.

Every accuracy number reported for a model on a held-out set is itself a
sample statistic -- a proportion, computed from a fixed number of examples --
and it carries a standard error exactly like any other. This exercise makes
the AI thread's claim a test: an accuracy of 91.4% measured on 500 examples
has a standard error of about 1.25 percentage points, so a 0.3-point
difference between two models on that same set is comfortably inside one
standard error of noise, not a demonstrated improvement.
"""

import dataset as D
from sampling import binomial_standard_error

se = binomial_standard_error(D.EVAL_ACCURACY, D.EVAL_N)
se_pct = se * 100.0

margin_in_se_units = D.EVAL_MARGIN_PCT / se_pct

print(f"accuracy = {D.EVAL_ACCURACY:.1%} on {D.EVAL_N} held-out examples")
print(f"standard error = sqrt(p_hat * (1 - p_hat) / n) = {se:.5f} = {se_pct:.3f} percentage points")
print()
print(f"a {D.EVAL_MARGIN_PCT} percentage-point difference between two models on this set "
      f"is {margin_in_se_units:.2f} standard errors -- ")
print("well inside one standard error, i.e. indistinguishable from noise on a set this size.")

assert abs(se_pct - D.EVAL_EXPECTED_SE_PCT) < D.EVAL_SE_TOLERANCE_PCT, (
    f"the binomial standard error came out to {se_pct:.3f} percentage points, "
    f"expected close to {D.EVAL_EXPECTED_SE_PCT}"
)
assert margin_in_se_units < 1.0, (
    f"a {D.EVAL_MARGIN_PCT}-point margin should sit well inside one standard error "
    f"({margin_in_se_units:.2f} SE), not outside it"
)

print("08_the_evaluation_margin.py: every assertion held.")

"""Exercise 9 -- reproducibility.

Every simulation in this lab depends on a seeded `numpy.random.Generator`.
The same seed must give bit-identical results -- otherwise none of the
"measured X" claims elsewhere in this lab mean anything, because a rerun
could silently produce a different number. A different seed should give a
DIFFERENT sampling distribution that still agrees with the first one within
a tolerance derived from the standard error, since both are estimating the
same underlying quantity.
"""

import numpy as np

import dataset as D
from sampling import sampling_distribution, theoretical_standard_error, population_mean_std

rng_a1 = np.random.default_rng(D.REPRO_SEED_A)
rng_a2 = np.random.default_rng(D.REPRO_SEED_A)
rng_b = np.random.default_rng(D.REPRO_SEED_B)

means_a1 = sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, rng_a1)
means_a2 = sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, rng_a2)
means_b = sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, rng_b)

identical = np.array_equal(means_a1, means_a2)
different = not np.array_equal(means_a1, means_b)

pop_mean, pop_sigma = population_mean_std(D.SKEWED_POP)
theoretical_se = theoretical_standard_error(pop_sigma, D.REPRO_N)
# The standard error of an estimate built from REPRO_TRIALS sample means.
se_of_mean_estimate = theoretical_se / np.sqrt(D.REPRO_TRIALS)

gap_a_b = abs(means_a1.mean() - means_b.mean())

print(f"seed {D.REPRO_SEED_A}, run 1: first three means = {means_a1[:3]}")
print(f"seed {D.REPRO_SEED_A}, run 2: first three means = {means_a2[:3]}")
print(f"seed {D.REPRO_SEED_B}, run 1: first three means = {means_b[:3]}")
print()
print(f"same seed produces bit-identical arrays: {identical}")
print(f"different seed produces a different array: {different}")
print(f"gap between seed {D.REPRO_SEED_A}'s and seed {D.REPRO_SEED_B}'s estimate of the mean "
      f"= {gap_a_b:.5f} ({gap_a_b / se_of_mean_estimate:.2f} standard errors)")

assert identical, "the same seed did not reproduce identical results"
assert different, "two different seeds produced identical results, which should not happen"
assert gap_a_b < 5.0 * se_of_mean_estimate, (
    "two different seeds' estimates of the population mean disagreed by more than "
    "5 standard errors -- they should agree within simulation noise"
)

print("09_reproducibility.py: every assertion held.")

"""Make this directory's own modules the ones its tests import.

Both `examples/` and `starter/` contain modules called `sampling` and
`dataset`, and pytest imports test files by putting their directory on
`sys.path`. Without this file, running `pytest` across both directories at
once would import whichever copy was seen first and reuse it for the other --
so the starter tests would silently pass against the reference solution
instead of skipping. That is a wrong answer with a green tick on it, which is
the worst kind.

So: put this directory first on the import path, and drop any already-imported
module of those names that came from somewhere else.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("sampling", "dataset", "answers"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""Every population, parameter and trial count used anywhere in this lab.

All data here is invented and deterministic: every population is generated
once, at import time, from `numpy.random.default_rng(DATASET_SEED)`. Nothing
is read from a file and nothing touches the network. If you want a different
population, change a parameter here and every script and test downstream
picks it up automatically.
"""

import numpy as np

DATASET_SEED = 117

_rng = np.random.default_rng(DATASET_SEED)

# ---------------------------------------------------------------------------
# Populations. POP_SIZE stands in for "the population" -- large enough that
# sampling from it with replacement behaves like sampling from the true
# infinite population to every decimal place this lab checks.
# ---------------------------------------------------------------------------

POP_SIZE = 200_000

# A skewed population: waiting-time shaped, heavily right-tailed. This is the
# population exercises 1, 2, 3 and 5 build the sampling distribution from.
SKEWED_SCALE = 3.0
SKEWED_POP = _rng.exponential(scale=SKEWED_SCALE, size=POP_SIZE)

# A biased coin: values are 0 or 1, P(1) = 0.2, so the population itself is
# skewed in the opposite direction from SKEWED_POP -- used in the lesson to
# show the CLT working on a population that is not just skewed but discrete
# with only two possible values.
COIN_P = 0.2
COIN_POP = (_rng.random(POP_SIZE) < COIN_P).astype(float)

# A lumpy, two-spike population: two clusters far apart with a little noise
# around each, and nothing in between -- about as far from bell-shaped as a
# population can look while still having finite variance. The two spikes
# carry UNEQUAL weight (80/20) so the population itself is genuinely skewed,
# not merely bimodal -- an equal-weight version is symmetric and has zero
# skewness already, which would make it useless for exercise 3's monotone-
# decrease check.
TWO_SPIKE_CENTERS = (-5.0, 5.0)
TWO_SPIKE_WEIGHTS = (0.8, 0.2)
TWO_SPIKE_NOISE = 0.4
_spike_choice = _rng.choice(TWO_SPIKE_CENTERS, size=POP_SIZE, p=TWO_SPIKE_WEIGHTS)
TWO_SPIKE_POP = _spike_choice + _rng.normal(0.0, TWO_SPIKE_NOISE, size=POP_SIZE)

# ---------------------------------------------------------------------------
# Exercise 1 -- the sampling distribution itself
# ---------------------------------------------------------------------------

EX1_N = 40
EX1_TRIALS = 20_000

# ---------------------------------------------------------------------------
# Exercise 2 -- the sqrt(n) law
# ---------------------------------------------------------------------------

SQRT_N_LAW_NS = (10, 40, 160, 640)
SQRT_N_LAW_TRIALS = 20_000
# Each successive n is exactly 4x the one before it, so the standard error
# should shrink by a factor of exactly sqrt(4) = 2 each step. Simulation noise
# keeps the measured ratio from landing on 2.000 exactly; this tolerance was
# checked across seeds 1, 2, 3, 42, 117 and 999 before being fixed here, and
# the worst observed ratio across those runs was 1.98.
SQRT_N_LAW_RATIO_TOLERANCE = 0.25

# ---------------------------------------------------------------------------
# Exercise 3 -- the CLT's skewness signature
# ---------------------------------------------------------------------------

SKEW_DEMO_NS = (2, 5, 20, 80, 320)
SKEW_DEMO_TRIALS = 50_000

# ---------------------------------------------------------------------------
# Exercise 4 -- the Cauchy counterexample
# ---------------------------------------------------------------------------

CAUCHY_DEMO_TRIALS = 20_000
CAUCHY_DEMO_N_SMALL = 10
CAUCHY_DEMO_N_LARGE = 1_000
EXPONENTIAL_SCALE = 1.0
# From n=10 to n=1000 the sample size grows 100x, so a finite-variance mean's
# spread should shrink by close to sqrt(100) = 10x. Checked across six seeds,
# the exponential ratio never fell below 9.7; 8.0 leaves comfortable margin
# while still ruling out "did not shrink".
EXPONENTIAL_SHRINK_FLOOR = 8.0
# The Cauchy mean's spread should NOT shrink at all -- checked across the same
# six seeds, the ratio stayed within 0.98-1.03. A factor-of-3 band in either
# direction is generous and still miles from the exponential's ~10x.
CAUCHY_NO_SHRINK_LOW = 1.0 / 3.0
CAUCHY_NO_SHRINK_HIGH = 3.0

# ---------------------------------------------------------------------------
# Exercise 5 -- bias does not shrink
# ---------------------------------------------------------------------------

BIAS_DEMO_TRIALS = 5_000
BIAS_DEMO_N_SMALL = 30
BIAS_DEMO_N_LARGE = 3_000
# Same 100x growth in n as the Cauchy exercise. The unbiased sampler's mean
# absolute error should shrink by close to 10x; checked across six seeds the
# ratio stayed in 9.8-10.3, so a floor of 7 is comfortable.
UNBIASED_SHRINK_FLOOR = 7.0
# The biased sampler's error should stay roughly flat. Checked across six
# seeds the ratio (error at n=30 / error at n=3000) stayed within 0.99-1.00;
# a band of 0.4-2.5 is generous and still clearly distinguishes "flat" from
# the unbiased sampler's ~10x drop.
BIASED_FLAT_LOW = 0.4
BIASED_FLAT_HIGH = 2.5

# ---------------------------------------------------------------------------
# Exercise 6 -- the bootstrap, from scratch
# ---------------------------------------------------------------------------

BOOTSTRAP_SAMPLE_SIZE = 200
BOOTSTRAP_SAMPLE_MEAN = 50.0
BOOTSTRAP_SAMPLE_STD = 10.0
BOOTSTRAP_N_BOOT = 5_000
# The bootstrap standard error of the MEAN has a closed form to check against
# (sigma_hat / sqrt(n)), so this tolerance can be tight: checked across six
# seeds, the relative error never exceeded 0.7%.
BOOTSTRAP_MEAN_RELATIVE_TOLERANCE = 0.15
# The MEDIAN has no closed form. "Sane" here means: within a generous factor
# of the spread of medians computed from genuinely fresh independent samples
# of the population. Checked across six seeds, the ratio of the two spreads
# stayed within 0.57-1.11; a factor of 3 in either direction leaves room for
# the small-sample noise inherent in estimating a spread of a spread.
BOOTSTRAP_MEDIAN_SANITY_LOW = 1.0 / 3.0
BOOTSTRAP_MEDIAN_SANITY_HIGH = 3.0
FRESH_MEDIAN_REPLICATIONS = 2_000

# ---------------------------------------------------------------------------
# Exercise 7 -- dependence inflates the true standard error
# ---------------------------------------------------------------------------

AR1_N = 200
AR1_PHI = 0.7
AR1_SIGMA = 1.0
AR1_REPLICATIONS = 3_000
# Checked across six seeds, true_se / naive_se stayed within 2.34-2.46 for
# this (phi, n) combination, so a floor of 1.5 leaves a wide, safe margin
# while still requiring a real, meaningful understatement.
AR1_INFLATION_FLOOR = 1.5

# ---------------------------------------------------------------------------
# Exercise 8 -- the evaluation-margin calculation
# ---------------------------------------------------------------------------

EVAL_ACCURACY = 0.914
EVAL_N = 500
EVAL_EXPECTED_SE_PCT = 1.25  # percentage points, to 2 decimal places
EVAL_SE_TOLERANCE_PCT = 0.05
EVAL_MARGIN_PCT = 0.3

# ---------------------------------------------------------------------------
# Exercise 9 -- reproducibility
# ---------------------------------------------------------------------------

REPRO_N = 40
REPRO_TRIALS = 2_000
REPRO_SEED_A = 7
REPRO_SEED_B = 8

"""Sampling distributions, the standard error, and the bootstrap -- from
scratch.

Every function here works from first principles: `numpy.random.Generator`
supplies uniform and named draws, and everything else -- the sampling
distribution itself, its standard error, its skewness, the bootstrap, the
naive-versus-true standard error under dependence -- is built on top of that,
not imported from a library that already does it.
"""

from __future__ import annotations

import math
from collections.abc import Callable

import numpy as np


def sampling_distribution(
    population: np.ndarray, n: int, trials: int, rng: np.random.Generator
) -> np.ndarray:
    """Draw `trials` independent samples of size `n`, with replacement, from
    `population`, and return the array of `trials` sample means.

    This is the central object of the whole lesson: not a single statistic,
    but the distribution *of* a statistic, built by literally repeating the
    experiment.
    """
    idx = rng.integers(0, population.shape[0], size=(trials, n))
    return population[idx].mean(axis=1)


def population_mean_std(population: np.ndarray) -> tuple[float, float]:
    """The population mean and the population standard deviation (divided by
    n, not n-1 -- there is no estimation happening here, the population is
    fully known)."""
    return float(population.mean()), float(population.std(ddof=0))


def theoretical_standard_error(sigma: float, n: int) -> float:
    """The standard error of the sample mean: sigma / sqrt(n)."""
    return sigma / math.sqrt(n)


def skewness(values: np.ndarray) -> float:
    """The sample skewness g1 = m3 / m2^1.5, where mk is the k-th central
    moment. Zero for a symmetric distribution (including the Normal),
    positive for a right-tailed one."""
    values = np.asarray(values, dtype=float)
    centered = values - values.mean()
    m2 = np.mean(centered**2)
    m3 = np.mean(centered**3)
    return float(m3 / m2**1.5)


def iqr(values: np.ndarray) -> float:
    """The interquartile range: a spread measure that stays meaningful even
    when the standard deviation does not, because it depends only on the
    order of the data, not on its moments. A Cauchy-distributed sample has an
    undefined population variance, so its sample standard deviation is not
    an estimate of anything -- the IQR is not affected by that at all."""
    values = np.asarray(values, dtype=float)
    q75, q25 = np.percentile(values, [75, 25])
    return float(q75 - q25)


def exponential_mean_iqr(
    n: int, trials: int, rng: np.random.Generator, scale: float = 1.0
) -> float:
    """The IQR of the sampling distribution of the mean, for `trials`
    samples of size `n` drawn directly from Exponential(scale)."""
    draws = rng.exponential(scale=scale, size=(trials, n))
    return iqr(draws.mean(axis=1))


def cauchy_mean_iqr(n: int, trials: int, rng: np.random.Generator) -> float:
    """The IQR of the sampling distribution of the mean, for `trials`
    samples of size `n` drawn from the standard Cauchy distribution.

    The Cauchy distribution has no defined mean or variance -- its tails are
    too heavy for either integral to converge. Averaging n Cauchy draws
    produces a value that is ITSELF standard-Cauchy distributed, for every n,
    which this function's caller checks by comparing this IQR at two very
    different values of n and finding them the same.
    """
    draws = rng.standard_cauchy(size=(trials, n))
    return iqr(draws.mean(axis=1))


def mean_absolute_error(estimates: np.ndarray, truth: float) -> float:
    """The average absolute distance between a set of estimates and the
    truth they are estimating -- a single number combining bias and
    variance, which is exactly why it is the right thing to track across
    exercise 5."""
    return float(np.mean(np.abs(np.asarray(estimates, dtype=float) - truth)))


def biased_pool(population: np.ndarray) -> np.ndarray:
    """The subset of the population strictly above its own median -- a
    sampling frame that can never produce a draw from the lower half, no
    matter how many draws it makes."""
    threshold = np.median(population)
    return population[population > threshold]


def bootstrap_replicates(
    data: np.ndarray,
    statistic: Callable[[np.ndarray], np.ndarray],
    n_boot: int,
    rng: np.random.Generator,
) -> np.ndarray:
    """Resample `data` with replacement `n_boot` times (each resample the
    same size as `data`), apply `statistic` to every resample, and return the
    `n_boot` results.

    `statistic` must accept a 2-D array of shape (n_boot, len(data)) and
    return a 1-D array of shape (n_boot,) -- `lambda a: a.mean(axis=1)` or
    `lambda a: numpy.median(a, axis=1)` are the two used in this lab.
    """
    data = np.asarray(data, dtype=float)
    n = data.shape[0]
    idx = rng.integers(0, n, size=(n_boot, n))
    resamples = data[idx]
    return np.asarray(statistic(resamples), dtype=float)


def bootstrap_standard_error(
    data: np.ndarray,
    statistic: Callable[[np.ndarray], np.ndarray],
    n_boot: int,
    rng: np.random.Generator,
) -> float:
    """The bootstrap estimate of a statistic's standard error: the standard
    deviation of the statistic computed across the bootstrap replicates.
    No formula for the statistic's sampling distribution is used or needed."""
    replicates = bootstrap_replicates(data, statistic, n_boot, rng)
    return float(replicates.std(ddof=1))


def ar1_series(n: int, phi: float, sigma: float, rng: np.random.Generator) -> np.ndarray:
    """A length-n AR(1) series: x[0] ~ Normal(0, sigma), and each later value
    is `phi * previous + innovation`, with the innovation's variance chosen
    so the series' own marginal variance stays sigma^2 throughout (the
    stationary AR(1) variance is innovation_variance / (1 - phi^2))."""
    x = np.empty(n, dtype=float)
    x[0] = rng.normal(0.0, sigma)
    innovation_sigma = sigma * math.sqrt(1.0 - phi**2)
    for t in range(1, n):
        x[t] = phi * x[t - 1] + rng.normal(0.0, innovation_sigma)
    return x


def naive_standard_error(series: np.ndarray) -> float:
    """The textbook formula, sample_std / sqrt(n), applied as if the
    observations were independent -- which is exactly the assumption an
    autocorrelated series violates."""
    series = np.asarray(series, dtype=float)
    n = series.shape[0]
    return float(series.std(ddof=1) / math.sqrt(n))


def true_standard_error_by_replication(
    n: int, phi: float, sigma: float, replications: int, rng: np.random.Generator
) -> float:
    """The actual standard error of the sample mean of an AR(1) series,
    measured the only way that does not assume independence: generate many
    independent length-n series, compute each one's sample mean, and take
    the standard deviation of THOSE means."""
    means = np.array(
        [ar1_series(n, phi, sigma, rng).mean() for _ in range(replications)]
    )
    return float(means.std(ddof=1))


def binomial_standard_error(phat: float, n: int) -> float:
    """The standard error of a sample proportion (equivalently, an accuracy
    measured on n examples): sqrt(phat * (1 - phat) / n)."""
    return math.sqrt(phat * (1.0 - phat) / n)

"""The reference pytest suite: every function in `sampling.py`, checked
against real values from real seeded runs.

Run from the lab directory:

    .venv/bin/pytest examples -q
"""

import math

import numpy as np
import pytest

import dataset as D
import sampling as S


# --------------------------------------------------------------------------
# sampling_distribution / theoretical_standard_error
# --------------------------------------------------------------------------


def test_sampling_distribution_returns_one_mean_per_trial():
    rng = np.random.default_rng(101)
    means = S.sampling_distribution(D.SKEWED_POP, n=20, trials=500, rng=rng)
    assert means.shape == (500,)


def test_sampling_distribution_mean_is_close_to_population_mean():
    rng = np.random.default_rng(102)
    pop_mean, pop_sigma = S.population_mean_std(D.SKEWED_POP)
    n, trials = 50, 20_000
    means = S.sampling_distribution(D.SKEWED_POP, n, trials, rng)
    se_of_mean = S.theoretical_standard_error(pop_sigma, n) / math.sqrt(trials)
    assert abs(means.mean() - pop_mean) < 4.0 * se_of_mean


def test_theoretical_standard_error_scales_as_inverse_sqrt_n():
    assert S.theoretical_standard_error(10.0, 100) == pytest.approx(1.0)
    assert S.theoretical_standard_error(10.0, 400) == pytest.approx(0.5)


# --------------------------------------------------------------------------
# skewness / iqr
# --------------------------------------------------------------------------


def test_skewness_of_a_symmetric_sample_is_near_zero():
    rng = np.random.default_rng(103)
    symmetric = rng.normal(0.0, 1.0, size=200_000)
    assert abs(S.skewness(symmetric)) < 0.02


def test_skewness_of_an_exponential_population_is_positive_and_near_two():
    # The Exponential distribution's population skewness is exactly 2,
    # regardless of scale. This checks the from-scratch estimator against
    # that known constant on a large sample.
    assert S.skewness(D.SKEWED_POP) == pytest.approx(2.0, abs=0.05)


def test_iqr_of_a_known_uniform_sample():
    values = np.arange(0, 101)  # 0..100, IQR should be exactly 50
    assert S.iqr(values) == pytest.approx(50.0)


def test_iqr_ignores_a_single_extreme_outlier():
    values = np.concatenate([np.arange(0, 99), [1_000_000.0]])
    assert S.iqr(values) < 100.0  # unaffected by the one wild value


# --------------------------------------------------------------------------
# exponential_mean_iqr / cauchy_mean_iqr
# --------------------------------------------------------------------------


def test_exponential_mean_iqr_shrinks_with_more_data():
    rng = np.random.default_rng(104)
    small = S.exponential_mean_iqr(10, 5_000, rng, scale=1.0)
    large = S.exponential_mean_iqr(1_000, 5_000, rng, scale=1.0)
    assert small / large > 5.0


def test_cauchy_mean_iqr_does_not_shrink_with_more_data():
    rng = np.random.default_rng(105)
    small = S.cauchy_mean_iqr(10, 5_000, rng)
    large = S.cauchy_mean_iqr(1_000, 5_000, rng)
    assert 0.3 < small / large < 3.0


# --------------------------------------------------------------------------
# mean_absolute_error / biased_pool
# --------------------------------------------------------------------------


def test_mean_absolute_error_is_zero_for_exact_estimates():
    assert S.mean_absolute_error(np.array([5.0, 5.0, 5.0]), 5.0) == 0.0


def test_mean_absolute_error_is_positive_for_biased_estimates():
    assert S.mean_absolute_error(np.array([6.0, 7.0, 8.0]), 5.0) == pytest.approx(2.0)


def test_biased_pool_contains_only_values_above_the_median():
    values = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
    pool = S.biased_pool(values)
    assert (pool > np.median(values)).all()
    assert len(pool) < len(values)


# --------------------------------------------------------------------------
# bootstrap_replicates / bootstrap_standard_error
# --------------------------------------------------------------------------


def test_bootstrap_replicates_shape():
    rng = np.random.default_rng(106)
    data = np.arange(50, dtype=float)
    reps = S.bootstrap_replicates(data, lambda a: a.mean(axis=1), 300, rng)
    assert reps.shape == (300,)


def test_bootstrap_se_of_the_mean_agrees_with_the_formula():
    rng = np.random.default_rng(107)
    data = rng.normal(0.0, 5.0, size=300)
    sigma_hat = data.std(ddof=1)
    theoretical = sigma_hat / math.sqrt(len(data))
    boot_se = S.bootstrap_standard_error(data, lambda a: a.mean(axis=1), 4_000, rng)
    assert abs(boot_se - theoretical) / theoretical < 0.15


# --------------------------------------------------------------------------
# ar1_series / naive_standard_error / true_standard_error_by_replication
# --------------------------------------------------------------------------


def test_ar1_series_has_the_right_length():
    rng = np.random.default_rng(108)
    series = S.ar1_series(n=50, phi=0.5, sigma=1.0, rng=rng)
    assert series.shape == (50,)


def test_ar1_series_is_independent_when_phi_is_zero():
    rng = np.random.default_rng(109)
    naive = np.mean(
        [S.naive_standard_error(S.ar1_series(300, 0.0, 1.0, rng)) for _ in range(200)]
    )
    true_se = S.true_standard_error_by_replication(300, 0.0, 1.0, 600, rng)
    # With phi = 0, there is no dependence, so naive and true should agree
    # closely -- this is the control case for exercise 7's main result.
    assert abs(naive - true_se) / true_se < 0.25


def test_dependence_makes_naive_se_understate_the_true_se():
    rng = np.random.default_rng(110)
    naive = np.mean(
        [S.naive_standard_error(S.ar1_series(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, rng)) for _ in range(200)]
    )
    true_se = S.true_standard_error_by_replication(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, 1_500, rng)
    assert true_se > naive


# --------------------------------------------------------------------------
# binomial_standard_error
# --------------------------------------------------------------------------


def test_binomial_standard_error_matches_the_brief():
    se_pct = S.binomial_standard_error(0.914, 500) * 100.0
    assert se_pct == pytest.approx(1.25, abs=0.05)


def test_binomial_standard_error_is_largest_at_p_one_half():
    se_half = S.binomial_standard_error(0.5, 100)
    se_extreme = S.binomial_standard_error(0.99, 100)
    assert se_half > se_extreme

lesson_id: D117
day: 117
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-117-sampling-and-the-central-limit-theorem
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy; print(numpy.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_sampling_distribution.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_the_sqrt_n_law.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_clt_from_a_skewed_population.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_the_cauchy_counterexample.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_bias_does_not_shrink.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_bootstrap_from_scratch.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_dependence_inflates_se.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 08_the_evaluation_margin.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 09_reproducibility.py && cd ..'
  - .venv/bin/pytest examples -q -p no:cacheprovider
  - .venv/bin/pytest starter -q -p no:cacheprovider
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - "find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {} +"
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 35
last_executed: '2026-08-17'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 -- bash tests/run_tests.sh -> 32 checks, 0 failure(s), exit 0; pytest examples -> 19 passed; pytest starter -> 1 passed, 12 skipped on an untouched checkout, and 13 passed against a fully solved copy of starter/ (verified by temporarily copying the reference sampling.py into starter/, confirming all 13 tests passed, then restoring the blank skeleton -- the skip count after restoring was confirmed back to 12, and collecting both suites together also reports 12 skipped, proving the two conftest.py import guards work). All nine reference scripts exit 0 with every internal assertion holding. Everything was run through a real lab-local .venv created by the documented setup commands, not through an authoring environment; pip install used the network exactly once, as documented. Section 5 of the harness re-runs script 08 (the evaluation-margin calculation) with its expected standard error temporarily replaced with a deliberately wrong value (99.0 instead of 1.25 percentage points), confirms the run exits non-zero with the named AssertionError showing both the wrong expectation and the correctly-computed value, and does not touch the real dataset.py file on disk -- so the suite is demonstrated to be capable of failing rather than merely claimed to be. Three honesty notes from this run. FIRST: scipy, pandas and matplotlib are not installed in this environment. scipy.stats (bootstrap, sem) and pandas'' DataFrame.sample are described from their public documentation in the lesson''s Tools section and explicitly marked as not run here; no output attributed to either anywhere in this lab or its lesson was actually produced by them. SECOND: every sampled figure in this lab (every standard error, skewness value, IQR, bootstrap replicate spread and AR(1) replication result) is a freshly measured number rather than a fixed literal, checked against a tolerance derived from a standard error, an inverse-square-root ratio, or a generous sanity band rather than a value chosen to make the test pass; the specific figures will differ slightly on another machine or NumPy version, as documented in expected-output/FIELDS.md. THIRD: every tolerance in dataset.py (the sqrt(n)-law ratio band, the Cauchy/Exponential shrink floors, the bias-flat band, the bootstrap-median sanity band, and the AR(1) inflation floor) was checked by rerunning the exercise logic across six different seeds (1, 2, 3, 42, 117, 999) during development before being fixed in that file, and the specific ranges observed across those seeds are recorded in the comments beside each tolerance -- this is reported as observed behaviour across six runs, not as a formal guarantee no seed can ever violate it. The binomial standard error (1.254 percentage points for 91.4% accuracy on 500 examples) is exact arithmetic from a closed-form formula and is identical on any correct implementation, anywhere.'

# What is installed, why, and what it costs

Two packages, both free and open source, both installed into a lab-local
virtual environment that `rm -rf .venv` completely undoes.

| Package | Version pinned | Licence | What this lab uses it for |
| --- | --- | --- | --- |
| `numpy` | 2.5.2 | BSD 3-Clause | `numpy.random.default_rng` for every seeded draw, and vectorised array operations for building sampling distributions of tens of thousands of trials at once. |
| `pytest` | 9.1.1 | MIT | The reference suite (19 tests) and your running score in `starter/`. |

There is no paid tier of anything in this lab, no account, no key and no
signup, personally or commercially.

## The one time the network is needed

```bash
.venv/bin/pip install -r requirements/requirements.txt
```

That is the only command in the lab that opens a connection. Section 6 of
`tests/run_tests.sh` greps every source file in `examples/` and `starter/`
to prove that nothing else does.

## What is deliberately *not* installed

**`scipy.stats`** does the two heaviest pieces of this lab's work for you:
`scipy.stats.bootstrap` runs exercise 6 in a few lines, with several
confidence-interval methods to choose from, and `scipy.stats.sem` computes
exercise 1's standard error of the mean directly from a sample. Neither is
installed here, and **no output from scipy is reproduced anywhere** in
this lab or its lesson. The lesson's Tools section describes both from
their public documentation.

That is not a limitation to apologise for. The bootstrap you write in
`starter/sampling.py` -- resample with replacement, recompute the
statistic, read the standard error off the spread -- is the exact idea
`scipy.stats.bootstrap` implements. The difference between the two is
engineering: bias-corrected-and-accelerated confidence intervals, batched
vectorisation, and a stable public API. Having written the fifteen-line
version, you will read `scipy.stats.bootstrap`'s documentation differently.

**`pandas`** is also not installed. `DataFrame.sample` does a version of
this lab's population sampling on tabular data, and is described from its
documentation in the lesson's Tools section, not run here.

**`matplotlib`** is not installed either. Every figure this lesson
describes -- the sampling distribution's shape, the Cauchy-versus-
Exponential contrast -- is described in words and numbers rather than
plotted; the two hand-authored SVG diagrams in the lesson carry the visual
argument instead.

## If you cannot install anything at all

You still need NumPy for this lab, unlike some earlier days in this
section: every exercise draws random samples, and the standard library's
`random` module does not offer the same vectorised batch sampling this lab
relies on to draw tens of thousands of trials in a single call without
writing a Python-level loop for each one. If NumPy genuinely cannot be
installed, the *ideas* -- a statistic has a distribution, the standard
error shrinks as `1/sqrt(n)`, bias does not shrink, the bootstrap resamples
with replacement -- can still be worked through by hand on a small dataset
with `random.Random`, but this lab's exercises and tests are not written
against that path.

numpy==2.5.2
pytest==9.1.1

# The nine exercises

Work through these in order, in `sampling.py`. Check yourself as you go:

```bash
.venv/bin/pytest starter -q
```

Unattempted work reports as **skipped**, never failed. Wrong work **fails**
with your answer printed beside the correct one.

## 1. The sampling distribution itself

`sampling_distribution(population, n, trials, rng)` and
`population_mean_std(population)`. Draw `trials` independent samples of size
`n` from `population`, with replacement, and return the array of sample
means. Assert its own mean is close to the population mean, and its own
standard deviation close to `sigma / sqrt(n)`, both within three standard
errors — a statistic has a distribution too, and this is what it looks like.

## 2. The standard error and the sqrt(n) law

`theoretical_standard_error(sigma, n)`. For `n` in `{10, 40, 160, 640}` --
each exactly 4x the one before -- assert the measured standard error roughly
**halves** each time, not quarters. Assert the *ratios* between successive
standard errors, not four hard-coded values.

## 3. The CLT from a skewed population

`skewness(values)`. Measure the skewness of the sampling distribution of the
mean for an increasingly large `n`, drawn from a heavily right-skewed
population. Assert it falls monotonically toward zero as `n` grows -- the
population's lopsidedness washes out of the statistic.

## 4. The Cauchy counterexample

`iqr(values)`, `exponential_mean_iqr(n, trials, rng, scale)`,
`cauchy_mean_iqr(n, trials, rng)`. For an Exponential population, assert the
spread of the sample mean shrinks by roughly the expected factor from n=10
to n=1000. For a **Cauchy** population, assert it does **not** shrink. Use
the IQR rather than the standard deviation -- a Cauchy sample's standard
deviation is not an estimate of anything, because the population variance it
would estimate does not exist.

## 5. Bias does not shrink

`mean_absolute_error(estimates, truth)`, `biased_pool(population)`. Build a
sampler that can only draw from the half of the population above its own
median. Assert its error stays roughly constant as `n` grows by a factor of
100, while an unbiased sampler's error shrinks by roughly ten -- more data
buys a biased sampler a more *precise* wrong answer, never a correct one.

## 6. The bootstrap, from scratch

`bootstrap_replicates(data, statistic, n_boot, rng)`,
`bootstrap_standard_error(data, statistic, n_boot, rng)`. Resample a dataset
with replacement, recompute a statistic on every resample, and read its
standard error off the spread of the results. Assert it agrees with
`sigma_hat / sqrt(n)` for the mean, then apply the same code to the
**median**, where no simple formula exists, and check the result for
sanity against the spread of medians from genuinely fresh samples.

## 7. Dependence inflates the true standard error

`ar1_series(n, phi, sigma, rng)`, `naive_standard_error(series)`,
`true_standard_error_by_replication(n, phi, sigma, replications, rng)`.
Generate an autocorrelated series, measure its TRUE standard error by
replication, and assert the naive `sample_std / sqrt(n)` formula understates
it meaningfully.

## 8. The evaluation-margin calculation

`binomial_standard_error(phat, n)`. For an accuracy of 91.4% on 500
examples, assert the standard error is about 1.25 percentage points, and
that a 0.3-point difference between two models is well inside one standard
error -- noise, not a demonstrated improvement.

## 9. Reproducibility

Nothing new to write here -- exercise 9 in `examples/09_reproducibility.py`
runs against whatever `sampling_distribution` you wrote for exercise 1. The
same seed must give identical results; a different seed must give different
results that still agree within tolerance. If exercise 1 is correct, this
one follows for free -- and if it does not hold, that is a sign exercise 1's
random-index construction is not using `rng` correctly.

"""Make this directory's own modules the ones its tests import.

Both `examples/` and `starter/` contain modules called `sampling` and
`dataset`, and pytest imports test files by putting their directory on
`sys.path`. Without this file, running `pytest` across both directories at
once would import whichever copy was seen first and reuse it for the other --
so these starter tests would silently pass against the reference solution
instead of skipping. That is a wrong answer with a green tick on it, which is
the worst kind.

So: put this directory first on the import path, and drop any already-imported
module of those names that came from somewhere else.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("sampling", "dataset", "answers"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""Every population, parameter and trial count used anywhere in this lab.

This file is complete and does not need editing -- your work for each
exercise lives in `sampling.py`. All data here is invented and deterministic:
every population is generated once, at import time, from
`numpy.random.default_rng(DATASET_SEED)`. Nothing is read from a file and
nothing touches the network.
"""

import numpy as np

DATASET_SEED = 117

_rng = np.random.default_rng(DATASET_SEED)

# ---------------------------------------------------------------------------
# Populations. POP_SIZE stands in for "the population" -- large enough that
# sampling from it with replacement behaves like sampling from the true
# infinite population to every decimal place this lab checks.
# ---------------------------------------------------------------------------

POP_SIZE = 200_000

# A skewed population: waiting-time shaped, heavily right-tailed. This is the
# population exercises 1, 2, 3 and 5 build the sampling distribution from.
SKEWED_SCALE = 3.0
SKEWED_POP = _rng.exponential(scale=SKEWED_SCALE, size=POP_SIZE)

# A biased coin: values are 0 or 1, P(1) = 0.2, so the population itself is
# skewed in the opposite direction from SKEWED_POP -- used in the lesson to
# show the CLT working on a population that is not just skewed but discrete
# with only two possible values.
COIN_P = 0.2
COIN_POP = (_rng.random(POP_SIZE) < COIN_P).astype(float)

# A lumpy, two-spike population: two clusters far apart with a little noise
# around each, and nothing in between -- about as far from bell-shaped as a
# population can look while still having finite variance. The two spikes
# carry UNEQUAL weight (80/20) so the population itself is genuinely skewed,
# not merely bimodal.
TWO_SPIKE_CENTERS = (-5.0, 5.0)
TWO_SPIKE_WEIGHTS = (0.8, 0.2)
TWO_SPIKE_NOISE = 0.4
_spike_choice = _rng.choice(TWO_SPIKE_CENTERS, size=POP_SIZE, p=TWO_SPIKE_WEIGHTS)
TWO_SPIKE_POP = _spike_choice + _rng.normal(0.0, TWO_SPIKE_NOISE, size=POP_SIZE)

# ---------------------------------------------------------------------------
# Exercise 1 -- the sampling distribution itself
# ---------------------------------------------------------------------------

EX1_N = 40
EX1_TRIALS = 20_000

# ---------------------------------------------------------------------------
# Exercise 2 -- the sqrt(n) law
# ---------------------------------------------------------------------------

SQRT_N_LAW_NS = (10, 40, 160, 640)
SQRT_N_LAW_TRIALS = 20_000
SQRT_N_LAW_RATIO_TOLERANCE = 0.25

# ---------------------------------------------------------------------------
# Exercise 3 -- the CLT's skewness signature
# ---------------------------------------------------------------------------

SKEW_DEMO_NS = (2, 5, 20, 80, 320)
SKEW_DEMO_TRIALS = 50_000

# ---------------------------------------------------------------------------
# Exercise 4 -- the Cauchy counterexample
# ---------------------------------------------------------------------------

CAUCHY_DEMO_TRIALS = 20_000
CAUCHY_DEMO_N_SMALL = 10
CAUCHY_DEMO_N_LARGE = 1_000
EXPONENTIAL_SCALE = 1.0
EXPONENTIAL_SHRINK_FLOOR = 8.0
CAUCHY_NO_SHRINK_LOW = 1.0 / 3.0
CAUCHY_NO_SHRINK_HIGH = 3.0

# ---------------------------------------------------------------------------
# Exercise 5 -- bias does not shrink
# ---------------------------------------------------------------------------

BIAS_DEMO_TRIALS = 5_000
BIAS_DEMO_N_SMALL = 30
BIAS_DEMO_N_LARGE = 3_000
UNBIASED_SHRINK_FLOOR = 7.0
BIASED_FLAT_LOW = 0.4
BIASED_FLAT_HIGH = 2.5

# ---------------------------------------------------------------------------
# Exercise 6 -- the bootstrap, from scratch
# ---------------------------------------------------------------------------

BOOTSTRAP_SAMPLE_SIZE = 200
BOOTSTRAP_SAMPLE_MEAN = 50.0
BOOTSTRAP_SAMPLE_STD = 10.0
BOOTSTRAP_N_BOOT = 5_000
BOOTSTRAP_MEAN_RELATIVE_TOLERANCE = 0.15
BOOTSTRAP_MEDIAN_SANITY_LOW = 1.0 / 3.0
BOOTSTRAP_MEDIAN_SANITY_HIGH = 3.0
FRESH_MEDIAN_REPLICATIONS = 2_000

# ---------------------------------------------------------------------------
# Exercise 7 -- dependence inflates the true standard error
# ---------------------------------------------------------------------------

AR1_N = 200
AR1_PHI = 0.7
AR1_SIGMA = 1.0
AR1_REPLICATIONS = 3_000
AR1_INFLATION_FLOOR = 1.5

# ---------------------------------------------------------------------------
# Exercise 8 -- the evaluation-margin calculation
# ---------------------------------------------------------------------------

EVAL_ACCURACY = 0.914
EVAL_N = 500
EVAL_EXPECTED_SE_PCT = 1.25  # percentage points, to 2 decimal places
EVAL_SE_TOLERANCE_PCT = 0.05
EVAL_MARGIN_PCT = 0.3

# ---------------------------------------------------------------------------
# Exercise 9 -- reproducibility
# ---------------------------------------------------------------------------

REPRO_N = 40
REPRO_TRIALS = 2_000
REPRO_SEED_A = 7
REPRO_SEED_B = 8

"""Sampling distributions, the standard error, and the bootstrap -- from
scratch.

Work through `starter/00_brief.md` in order, filling in the functions below.
Check your progress with:

    .venv/bin/pytest starter -q

Unattempted work reports as skipped, never failed. Every function currently
returns None -- replace the body, not just the return statement.
"""

from __future__ import annotations

import math
from collections.abc import Callable

import numpy as np


def sampling_distribution(
    population: np.ndarray, n: int, trials: int, rng: np.random.Generator
) -> np.ndarray:
    """Exercise 1. Draw `trials` independent samples of size `n`, with
    replacement, from `population`, and return the array of `trials` sample
    means.

    Hint: `rng.integers(0, population.shape[0], size=(trials, n))` gives you
    every sample's indices in one call; index `population` with that array
    and take `.mean(axis=1)`.
    """
    return None


def population_mean_std(population: np.ndarray) -> tuple[float, float]:
    """The population mean and the population standard deviation (ddof=0,
    since the whole population is known -- there is no estimation here)."""
    return None


def theoretical_standard_error(sigma: float, n: int) -> float:
    """Exercise 1 / 2. The standard error of the sample mean: sigma / sqrt(n)."""
    return None


def skewness(values: np.ndarray) -> float:
    """Exercise 3. The sample skewness g1 = m3 / m2^1.5, where mk is the k-th
    central moment: mk = mean((x - mean(x))**k)."""
    return None


def iqr(values: np.ndarray) -> float:
    """Exercise 4. The interquartile range: the 75th percentile minus the
    25th. Use `numpy.percentile`."""
    return None


def exponential_mean_iqr(
    n: int, trials: int, rng: np.random.Generator, scale: float = 1.0
) -> float:
    """Exercise 4. The IQR of the sampling distribution of the mean, for
    `trials` samples of size `n` drawn directly from Exponential(scale) via
    `rng.exponential`."""
    return None


def cauchy_mean_iqr(n: int, trials: int, rng: np.random.Generator) -> float:
    """Exercise 4. The same idea, drawing from `rng.standard_cauchy` instead."""
    return None


def mean_absolute_error(estimates: np.ndarray, truth: float) -> float:
    """Exercise 5. The average absolute distance between a set of estimates
    and the truth they are estimating."""
    return None


def biased_pool(population: np.ndarray) -> np.ndarray:
    """Exercise 5. The subset of `population` strictly above its own median.
    Use `numpy.median` and boolean indexing."""
    return None


def bootstrap_replicates(
    data: np.ndarray,
    statistic: Callable[[np.ndarray], np.ndarray],
    n_boot: int,
    rng: np.random.Generator,
) -> np.ndarray:
    """Exercise 6. Resample `data` with replacement `n_boot` times (each
    resample the same size as `data`), apply `statistic` to every resample,
    and return the `n_boot` results.

    Hint: build an index array of shape (n_boot, len(data)) with
    `rng.integers`, index `data` with it to get all the resamples at once,
    and pass that 2-D array straight to `statistic`.
    """
    return None


def bootstrap_standard_error(
    data: np.ndarray,
    statistic: Callable[[np.ndarray], np.ndarray],
    n_boot: int,
    rng: np.random.Generator,
) -> float:
    """Exercise 6. The standard deviation (ddof=1) of `bootstrap_replicates`."""
    return None


def ar1_series(n: int, phi: float, sigma: float, rng: np.random.Generator) -> np.ndarray:
    """Exercise 7. A length-n AR(1) series: x[0] ~ Normal(0, sigma), and each
    later value is `phi * previous + innovation`, where the innovation's own
    standard deviation is `sigma * sqrt(1 - phi**2)` -- chosen so the
    series' marginal variance stays sigma**2 throughout."""
    return None


def naive_standard_error(series: np.ndarray) -> float:
    """Exercise 7. sample_std(ddof=1) / sqrt(n) -- the textbook formula,
    applied as if the observations were independent."""
    return None


def true_standard_error_by_replication(
    n: int, phi: float, sigma: float, replications: int, rng: np.random.Generator
) -> float:
    """Exercise 7. Generate `replications` independent length-n AR(1) series,
    compute each one's sample mean, and return the standard deviation
    (ddof=1) of those means."""
    return None


def binomial_standard_error(phat: float, n: int) -> float:
    """Exercise 8. sqrt(phat * (1 - phat) / n)."""
    return None

"""Your running score. Unattempted work SKIPS; wrong work FAILS with both
values.

Run from the lab directory:

    .venv/bin/pytest starter -q

On an untouched checkout this reports one pass and everything else skipped.
A skip means "not attempted". A failure means "attempted and wrong", and the
message shows your answer next to the real one so you can see the gap rather
than guess at it.
"""

import math

import numpy as np
import pytest

import dataset as D
import sampling as S


def attempt(fn, what):
    """Call something that may not be written yet, and skip if it is not."""
    try:
        result = fn()
    except (TypeError, AttributeError, NotImplementedError):
        pytest.skip(f"not attempted yet: {what}")
    if result is None:
        pytest.skip(f"not attempted yet: {what}")
    return result


def test_the_suite_itself_runs():
    """One test that always passes, so a green run is distinguishable from a
    collection error that quietly ran nothing at all."""
    assert D.POP_SIZE > 0


# --------------------------------------------------------------------------
# Exercise 1 -- the sampling distribution itself
# --------------------------------------------------------------------------


def test_1_sampling_distribution_has_one_mean_per_trial():
    rng = np.random.default_rng(201)
    means = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, n=20, trials=300, rng=rng),
        "sampling_distribution",
    )
    assert means.shape == (300,)


def test_1_sampling_distribution_mean_is_close_to_population_mean():
    rng = np.random.default_rng(202)
    pop_mean, pop_sigma = attempt(
        lambda: S.population_mean_std(D.SKEWED_POP), "population_mean_std"
    )
    n, trials = 50, 20_000
    means = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, n, trials, rng), "sampling_distribution"
    )
    se_theory = attempt(
        lambda: S.theoretical_standard_error(pop_sigma, n), "theoretical_standard_error"
    )
    se_of_mean = se_theory / math.sqrt(trials)
    assert abs(means.mean() - pop_mean) < 4.0 * se_of_mean


# --------------------------------------------------------------------------
# Exercise 2 -- the sqrt(n) law
# --------------------------------------------------------------------------


def test_2_standard_error_halves_when_n_quadruples():
    rng = np.random.default_rng(203)
    ses = {}
    for n in D.SQRT_N_LAW_NS:
        means = attempt(
            lambda n=n: S.sampling_distribution(D.SKEWED_POP, n, D.SQRT_N_LAW_TRIALS, rng),
            "sampling_distribution",
        )
        ses[n] = means.std(ddof=1)
    for smaller, larger in zip(D.SQRT_N_LAW_NS, D.SQRT_N_LAW_NS[1:]):
        ratio = ses[smaller] / ses[larger]
        assert abs(ratio - 2.0) < D.SQRT_N_LAW_RATIO_TOLERANCE, (
            f"SE(n={smaller})/SE(n={larger}) = {ratio:.3f}, expected near 2.0"
        )


# --------------------------------------------------------------------------
# Exercise 3 -- CLT from a skewed population
# --------------------------------------------------------------------------


def test_3_skewness_of_the_sampling_distribution_decreases_monotonically():
    rng = np.random.default_rng(204)
    skews = []
    for n in D.SKEW_DEMO_NS:
        means = attempt(
            lambda n=n: S.sampling_distribution(D.SKEWED_POP, n, D.SKEW_DEMO_TRIALS, rng),
            "sampling_distribution",
        )
        skews.append(attempt(lambda means=means: S.skewness(means), "skewness"))
    for a, b in zip(skews, skews[1:]):
        assert b < a, f"skewness did not decrease: {skews}"


# --------------------------------------------------------------------------
# Exercise 4 -- the Cauchy counterexample
# --------------------------------------------------------------------------


def test_4_exponential_mean_shrinks_but_cauchy_mean_does_not():
    rng = np.random.default_rng(205)
    exp_small = attempt(
        lambda: S.exponential_mean_iqr(D.CAUCHY_DEMO_N_SMALL, D.CAUCHY_DEMO_TRIALS, rng, D.EXPONENTIAL_SCALE),
        "exponential_mean_iqr",
    )
    exp_large = attempt(
        lambda: S.exponential_mean_iqr(D.CAUCHY_DEMO_N_LARGE, D.CAUCHY_DEMO_TRIALS, rng, D.EXPONENTIAL_SCALE),
        "exponential_mean_iqr",
    )
    cauchy_small = attempt(
        lambda: S.cauchy_mean_iqr(D.CAUCHY_DEMO_N_SMALL, D.CAUCHY_DEMO_TRIALS, rng), "cauchy_mean_iqr"
    )
    cauchy_large = attempt(
        lambda: S.cauchy_mean_iqr(D.CAUCHY_DEMO_N_LARGE, D.CAUCHY_DEMO_TRIALS, rng), "cauchy_mean_iqr"
    )
    assert exp_small / exp_large > D.EXPONENTIAL_SHRINK_FLOOR
    assert D.CAUCHY_NO_SHRINK_LOW < cauchy_small / cauchy_large < D.CAUCHY_NO_SHRINK_HIGH


# --------------------------------------------------------------------------
# Exercise 5 -- bias does not shrink
# --------------------------------------------------------------------------


def test_5_biased_sampler_error_stays_flat_while_unbiased_shrinks():
    rng = np.random.default_rng(206)
    true_mean = float(D.SKEWED_POP.mean())
    pool = attempt(lambda: S.biased_pool(D.SKEWED_POP), "biased_pool")

    def err(population, n):
        means = S.sampling_distribution(population, n, D.BIAS_DEMO_TRIALS, rng)
        return attempt(lambda: S.mean_absolute_error(means, true_mean), "mean_absolute_error")

    unbiased_small = err(D.SKEWED_POP, D.BIAS_DEMO_N_SMALL)
    unbiased_large = err(D.SKEWED_POP, D.BIAS_DEMO_N_LARGE)
    biased_small = err(pool, D.BIAS_DEMO_N_SMALL)
    biased_large = err(pool, D.BIAS_DEMO_N_LARGE)

    assert unbiased_small / unbiased_large > D.UNBIASED_SHRINK_FLOOR
    assert D.BIASED_FLAT_LOW < biased_small / biased_large < D.BIASED_FLAT_HIGH


# --------------------------------------------------------------------------
# Exercise 6 -- the bootstrap, from scratch
# --------------------------------------------------------------------------


def test_6_bootstrap_se_of_mean_agrees_with_formula():
    rng = np.random.default_rng(207)
    data = rng.normal(D.BOOTSTRAP_SAMPLE_MEAN, D.BOOTSTRAP_SAMPLE_STD, D.BOOTSTRAP_SAMPLE_SIZE)
    theoretical = data.std(ddof=1) / math.sqrt(len(data))
    boot_se = attempt(
        lambda: S.bootstrap_standard_error(data, lambda a: a.mean(axis=1), D.BOOTSTRAP_N_BOOT, rng),
        "bootstrap_standard_error",
    )
    assert abs(boot_se - theoretical) / theoretical < D.BOOTSTRAP_MEAN_RELATIVE_TOLERANCE


def test_6_bootstrap_se_of_median_is_sane():
    rng = np.random.default_rng(208)
    data = rng.normal(D.BOOTSTRAP_SAMPLE_MEAN, D.BOOTSTRAP_SAMPLE_STD, D.BOOTSTRAP_SAMPLE_SIZE)
    boot_se_median = attempt(
        lambda: S.bootstrap_standard_error(data, lambda a: np.median(a, axis=1), D.BOOTSTRAP_N_BOOT, rng),
        "bootstrap_standard_error",
    )
    fresh = rng.normal(
        D.BOOTSTRAP_SAMPLE_MEAN, D.BOOTSTRAP_SAMPLE_STD, size=(D.FRESH_MEDIAN_REPLICATIONS, D.BOOTSTRAP_SAMPLE_SIZE)
    )
    fresh_se = np.median(fresh, axis=1).std(ddof=1)
    ratio = boot_se_median / fresh_se
    assert D.BOOTSTRAP_MEDIAN_SANITY_LOW < ratio < D.BOOTSTRAP_MEDIAN_SANITY_HIGH


# --------------------------------------------------------------------------
# Exercise 7 -- dependence inflates the true standard error
# --------------------------------------------------------------------------


def test_7_naive_se_understates_the_true_se_under_dependence():
    rng = np.random.default_rng(209)
    naive_ses = [
        attempt(
            lambda: S.naive_standard_error(S.ar1_series(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, rng)),
            "naive_standard_error / ar1_series",
        )
        for _ in range(300)
    ]
    naive_avg = float(np.mean(naive_ses))
    true_se = attempt(
        lambda: S.true_standard_error_by_replication(D.AR1_N, D.AR1_PHI, D.AR1_SIGMA, D.AR1_REPLICATIONS, rng),
        "true_standard_error_by_replication",
    )
    assert true_se / naive_avg > D.AR1_INFLATION_FLOOR


# --------------------------------------------------------------------------
# Exercise 8 -- the evaluation-margin calculation
# --------------------------------------------------------------------------


def test_8_binomial_se_matches_the_brief_and_the_margin_is_noise():
    se = attempt(lambda: S.binomial_standard_error(D.EVAL_ACCURACY, D.EVAL_N), "binomial_standard_error")
    se_pct = se * 100.0
    assert abs(se_pct - D.EVAL_EXPECTED_SE_PCT) < D.EVAL_SE_TOLERANCE_PCT
    assert D.EVAL_MARGIN_PCT / se_pct < 1.0


# --------------------------------------------------------------------------
# Exercise 9 -- reproducibility (depends only on exercise 1)
# --------------------------------------------------------------------------


def test_9_same_seed_reproduces_identical_results():
    a1 = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, np.random.default_rng(D.REPRO_SEED_A)),
        "sampling_distribution",
    )
    a2 = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, np.random.default_rng(D.REPRO_SEED_A)),
        "sampling_distribution",
    )
    assert np.array_equal(a1, a2)


def test_9_different_seed_gives_a_different_but_compatible_result():
    a = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, np.random.default_rng(D.REPRO_SEED_A)),
        "sampling_distribution",
    )
    b = attempt(
        lambda: S.sampling_distribution(D.SKEWED_POP, D.REPRO_N, D.REPRO_TRIALS, np.random.default_rng(D.REPRO_SEED_B)),
        "sampling_distribution",
    )
    assert not np.array_equal(a, b)
    pop_mean, pop_sigma = attempt(lambda: S.population_mean_std(D.SKEWED_POP), "population_mean_std")
    se_theory = attempt(lambda: S.theoretical_standard_error(pop_sigma, D.REPRO_N), "theoretical_standard_error")
    se_of_mean = se_theory / math.sqrt(D.REPRO_TRIALS)
    assert abs(a.mean() - b.mean()) < 5.0 * se_of_mean

#!/usr/bin/env bash
# Tests for the Day 117 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running code and reading real
# values, never by reading source:
#
#   * a statistic has a distribution -- the sampling distribution of the mean
#     has its own mean close to the population mean and its own spread close
#     to sigma / sqrt(n);
#   * the sqrt(n) law -- quadrupling n roughly halves the standard error,
#     asserted as ratios across four sample sizes, not four fixed values;
#   * the central limit theorem -- the skewness of the sampling distribution
#     of the mean, built from a heavily skewed population, falls
#     monotonically toward zero as n grows;
#   * the Cauchy counterexample -- an Exponential population's sample mean
#     tightens by roughly 10x from n=10 to n=1000, while a Cauchy
#     population's sample mean does not tighten at all;
#   * sampling bias is not sampling error -- a biased sampling frame's error
#     stays flat as n grows 100x, while an unbiased sampler's error falls by
#     roughly ten;
#   * the bootstrap, from scratch, agrees with sigma_hat / sqrt(n) for the
#     mean and gives a sane answer for the median, where no formula exists;
#   * dependence quietly inflates the true standard error beyond what the
#     naive sigma_hat / sqrt(n) formula reports;
#   * the evaluation-margin calculation -- a 0.3-point accuracy difference on
#     500 examples sits well inside one binomial standard error;
#   * the same seed reproduces identical results, and a different seed gives
#     a different but statistically compatible one;
#   * nothing is left behind on disk.
#
# Everything after the one-time install runs offline. Nothing binds a port,
# nothing writes outside the lab, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

lab_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"

# Bytecode left by an EARLIER command is not this run's litter. The README
# documents `pytest starter -q`, and running it writes .pyc files that would
# then fail the cleanliness check at the end of this script -- failing the
# reader for following the instructions. Clearing them here makes that final
# check measure what it claims to: what THIS run left behind. `.venv` is
# untouched, because the packages' own bytecode is theirs, not ours.
find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -exec rm -rf {} + 2>/dev/null || true
find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -exec rm -rf {} + 2>/dev/null || true

failures=0
checks=0

check() {
  local label="$1" ok="$2"
  checks=$((checks + 1))
  if [ "${ok}" = "yes" ]; then
    echo "  ok: ${label}"
  else
    echo "  FAIL: ${label}"
    failures=$((failures + 1))
  fi
}

check_eq() {
  # check_eq <label> <expected> <actual>
  if [ "$2" = "$3" ]; then
    check "$1" "yes"
  else
    check "$1 (expected [$2], got [$3])" "no"
  fi
}

# Resolve pytest: an explicit override, then this lab's .venv, then PATH.
# Fails loudly with instructions rather than silently skipping checks.
resolve_tool() {
  local tool="$1" override="$2"
  if [ -n "${override}" ] && [ -x "${override}" ]; then echo "${override}"; return 0; fi
  if [ -x "${lab_dir}/.venv/bin/${tool}" ]; then echo "${lab_dir}/.venv/bin/${tool}"; return 0; fi
  if command -v "${tool}" >/dev/null 2>&1; then command -v "${tool}"; return 0; fi
  return 1
}

pytest_bin="$(resolve_tool pytest "${PYTEST:-}")" || {
  echo "FAIL: pytest not found." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  echo "  Or point this suite at an existing pytest:" >&2
  echo "    PYTEST=/path/to/pytest bash tests/run_tests.sh" >&2
  exit 1
}

python_bin="$(dirname "${pytest_bin}")/python3"
if [ ! -x "${python_bin}" ]; then
  python_bin="$(command -v python3 || true)"
fi
if [ -z "${python_bin}" ]; then
  echo "FAIL: python3 not found on PATH." >&2
  exit 1
fi

if ! "${python_bin}" -c "import numpy" >/dev/null 2>&1; then
  echo "FAIL: numpy is not importable from ${python_bin}." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  exit 1
fi

echo "Day 117 — Sampling and the Central Limit Theorem"
echo

# --------------------------------------------------------------------------
echo "1. The tools and the versions this lab was written against"
# --------------------------------------------------------------------------

versions="$("${python_bin}" - <<'PY'
import platform
import sys
from importlib.metadata import version

print(f"python   {platform.python_version()}")
for name in ("numpy", "pytest"):
    print(f"{name:<8} {version(name)}")
print(f"platform {platform.platform()}")
print(f"exe      {sys.executable.rsplit('/', 3)[-1]}")
PY
)"
echo "${versions}" | sed 's/^/  /'

pinned_numpy="$(grep -E '^numpy==' "${lab_dir}/requirements/requirements.txt" | cut -d= -f3)"
installed_numpy="$("${python_bin}" -c "from importlib.metadata import version; print(version('numpy'))")"
check_eq "installed numpy matches requirements.txt" "${pinned_numpy}" "${installed_numpy}"

major="$("${python_bin}" -c "import numpy; print(numpy.__version__.split('.')[0])")"
check_eq "numpy is version 2 or later" "2" "${major}"

# --------------------------------------------------------------------------
echo
echo "2. Every reference script runs and every assertion inside it holds"
# --------------------------------------------------------------------------

for script in 01_sampling_distribution 02_the_sqrt_n_law 03_clt_from_a_skewed_population \
              04_the_cauchy_counterexample 05_bias_does_not_shrink 06_bootstrap_from_scratch \
              07_dependence_inflates_se 08_the_evaluation_margin 09_reproducibility; do
  out="$(cd "${lab_dir}/examples" && "${python_bin}" "${script}.py" 2>&1)"
  status=$?
  if [ "${status}" -ne 0 ]; then
    check "${script}.py exits 0" "no"
    echo "${out}" | tail -5 | sed 's/^/      /'
  else
    check "${script}.py exits 0" "yes"
  fi
  case "${out}" in
    *"${script}.py: every assertion held."*)
      check "${script}.py reports every assertion held" "yes" ;;
    *) check "${script}.py reports every assertion held" "no" ;;
  esac
done

# --------------------------------------------------------------------------
echo
echo "3. The reference pytest suite: real values, real exceptions"
# --------------------------------------------------------------------------

ref_out="$(cd "${lab_dir}" && "${pytest_bin}" examples -q -p no:cacheprovider 2>&1)"
ref_status=$?
echo "${ref_out}" | tail -3 | sed 's/^/  /'
if [ "${ref_status}" -eq 0 ]; then
  check "pytest examples exits 0" "yes"
else
  check "pytest examples exits 0" "no"
fi
case "${ref_out}" in
  *" failed"*) check "no test in the reference suite failed" "no" ;;
  *)           check "no test in the reference suite failed" "yes" ;;
esac
ref_passed="$(printf '%s\n' "${ref_out}" | grep -o '[0-9][0-9]* passed' | head -1 | cut -d' ' -f1)"
if [ "${ref_passed:-0}" -ge 15 ]; then
  check "the reference suite ran at least 15 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 15 tests (ran ${ref_passed:-0})" "no"
fi

# --------------------------------------------------------------------------
echo
echo "4. The starter suite skips unattempted work instead of failing it"
# --------------------------------------------------------------------------

start_out="$(cd "${lab_dir}" && "${pytest_bin}" starter -q -p no:cacheprovider 2>&1)"
start_status=$?
echo "${start_out}" | tail -3 | sed 's/^/  /'
if [ "${start_status}" -eq 0 ]; then
  check "pytest starter exits 0 on an untouched checkout" "yes"
else
  check "pytest starter exits 0 on an untouched checkout" "no"
fi
case "${start_out}" in
  *" failed"*) check "the starter suite reports no failures" "no" ;;
  *)           check "the starter suite reports no failures" "yes" ;;
esac
case "${start_out}" in
  *skipped*) check "unwritten exercises are reported as skipped, not passed" "yes" ;;
  *) check "unwritten exercises are reported as skipped, not passed" "no" ;;
esac

# The import guard. Both directories contain modules called `sampling` and
# `dataset`, and pytest imports test files by putting their directory on
# sys.path -- so collecting both suites at once would otherwise let the
# starter tests import the REFERENCE solution and report unwritten exercises
# as passing. Each directory's conftest.py prevents that. This check proves
# it still does: across both suites, the skip count must be unchanged.
both_out="$(cd "${lab_dir}" && "${pytest_bin}" -q -p no:cacheprovider 2>&1)"
start_skipped="$(printf '%s\n' "${start_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
both_skipped="$(printf '%s\n' "${both_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
check_eq "collecting both suites at once does not turn skips into passes" \
  "${start_skipped:-none}" "${both_skipped:-none}"

# --------------------------------------------------------------------------
echo
echo "5. The harness can actually fail"
# --------------------------------------------------------------------------

# A green test suite proves nothing until you have watched it go red. This
# section re-runs script 08 (the evaluation margin) with its expected
# standard error deliberately swapped for a wrong one, and asserts that the
# re-run reports the failure and exits non-zero. If this section passes,
# section 2 is not decorative.
if [ -z "${D117_SELF_TEST:-}" ]; then
  self_out="$(cd "${lab_dir}/examples" && D117_SELF_TEST=1 "${python_bin}" -c "
import dataset as D
D.EVAL_EXPECTED_SE_PCT = 99.0  # a deliberately wrong expectation
exec(open('08_the_evaluation_margin.py').read())
" 2>&1)"
  self_status=$?
  if [ "${self_status}" -ne 0 ]; then
    check "a deliberately wrong expectation makes script 08 exit non-zero (${self_status})" "yes"
  else
    check "a deliberately wrong expectation makes script 08 exit non-zero" "no"
  fi
  case "${self_out}" in
    *"AssertionError"*"the binomial standard error came out to"*)
      check "the failing assertion is named in the output with both values" "yes" ;;
    *) check "the failing assertion is named in the output with both values" "no" ;;
  esac
else
  echo "  (self-test run: section 5 does not recurse)"
fi

# --------------------------------------------------------------------------
echo
echo "6. Nothing was left behind"
# --------------------------------------------------------------------------

# `.venv` is pruned from both searches below. The virtual environment ships
# NumPy's and pytest's own precompiled bytecode -- hundreds of __pycache__
# directories that came with the packages and have nothing to do with whether
# THIS lab tidied up after itself. Searching them would report a failure the
# reader cannot fix and did not cause. Everything the lab itself writes lives
# outside `.venv`, which is exactly what these two checks look at.

if find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -print -quit 2>/dev/null | grep -q .; then
  check "no __pycache__ directory left by the lab's own code" "no"
else
  check "no __pycache__ directory left by the lab's own code" "yes"
fi

if find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -print -quit 2>/dev/null | grep -q .; then
  check "no .pytest_cache directory left under the lab" "no"
else
  check "no .pytest_cache directory left under the lab" "yes"
fi

if grep -rqE 'urlopen|requests\.|socket\.|http://|https://' \
     "${lab_dir}/examples" "${lab_dir}/starter" 2>/dev/null; then
  check "no lab source opens a network connection" "no"
else
  check "no lab source opens a network connection" "yes"
fi

echo
echo "${checks} checks, ${failures} failure(s)."
[ "${failures}" -eq 0 ]