Hands-on lab — Day 118: Hypothesis Tests and Confidence Intervals

"""Exercise 1 -- two-sample z-test from scratch, checked against a hand
computation done with the standard library only.

Two fixed, non-random samples so the "hand" side of the comparison is
exact arithmetic anyone can redo with a calculator, not a random draw.
"""
import math
import statistics
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

from inference import phi, two_sample_z_test  # noqa: E402


def hand_two_sample_z(a: list[float], b: list[float]) -> tuple[float, float]:
    mean_a, var_a = statistics.mean(a), statistics.variance(a)
    mean_b, var_b = statistics.mean(b), statistics.variance(b)
    se = math.sqrt(var_a / len(a) + var_b / len(b))
    z = (mean_a - mean_b) / se
    p = 2.0 * (1.0 - phi(abs(z)))
    return z, p


def main() -> None:
    # Case 1: a clear, large separation.
    a = [50, 52, 49, 51, 53, 48, 50, 52, 51, 49]
    b = [54, 55, 53, 56, 54, 52, 55, 53, 54, 56]
    z_lib, p_lib = two_sample_z_test(a, b)
    z_hand, p_hand = hand_two_sample_z(a, b)
    print(f"Case 1 (clear separation): z_library={z_lib:.6f}  p_library={p_lib:.10e}")
    print(f"Case 1 (clear separation): z_hand={z_hand:.6f}  p_hand={p_hand:.10e}")
    assert abs(z_lib - z_hand) < 1e-9, "z from the library must match the hand computation"
    assert abs(p_lib - p_hand) < 1e-9, "p from the library must match the hand computation"
    assert p_lib < 0.001, "these two samples are obviously different -- p should be tiny"

    # Case 2: nearly identical samples -- p should be large (not significant).
    c = [50, 52, 49, 51, 53, 48, 50, 52, 51, 49]
    d = [51, 51, 50, 52, 52, 49, 51, 51, 52, 48]
    z_lib2, p_lib2 = two_sample_z_test(c, d)
    z_hand2, p_hand2 = hand_two_sample_z(c, d)
    print(f"Case 2 (near-identical): z_library={z_lib2:.6f}  p_library={p_lib2:.6f}")
    print(f"Case 2 (near-identical): z_hand={z_hand2:.6f}  p_hand={p_hand2:.6f}")
    assert abs(z_lib2 - z_hand2) < 1e-9
    assert abs(p_lib2 - p_hand2) < 1e-9
    assert p_lib2 > 0.30, "these two samples overlap heavily -- p should not be small"

    print("OK: two_sample_z_test agrees with an independent hand computation on both cases.")


if __name__ == "__main__":
    main()

"""Exercise 2 -- the centrepiece. What does "95% confidence" actually mean?

Build 10,000 nominal-95% confidence intervals from independent samples of a
population with a KNOWN true mean, and count how many of them actually
contain it. If the textbook claim is right, about 95% will -- not because
any one interval has a 95% chance of containing the fixed true mean (a
fixed number either is or is not in a fixed interval), but because the
PROCEDURE that builds the interval catches the true value 95% of the time
across repetition.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import confidence_interval_mean  # noqa: E402


def measure_coverage(rng: np.random.Generator, trials: int, n: int) -> float:
    hits = 0
    for _ in range(trials):
        sample = ds.normal_population(rng, n, ds.POP_MEAN, ds.POP_STD)
        lo, hi = confidence_interval_mean(sample, alpha=0.05)
        if lo <= ds.POP_MEAN <= hi:
            hits += 1
    return hits / trials


def main() -> None:
    rng = np.random.default_rng(42)
    coverage = measure_coverage(rng, ds.COVERAGE_TRIALS, ds.COVERAGE_SAMPLE_N)

    print(f"True population mean: {ds.POP_MEAN}")
    print(f"Nominal confidence level: {ds.COVERAGE_TARGET:.0%}")
    print(f"Intervals built: {ds.COVERAGE_TRIALS}")
    print(f"Measured coverage: {coverage:.4f}")
    print(f"Standard error of the measured coverage: {ds.COVERAGE_SE:.5f}")
    print(f"3-SE tolerance band: [{ds.COVERAGE_TARGET - ds.COVERAGE_TOLERANCE:.4f}, "
          f"{ds.COVERAGE_TARGET + ds.COVERAGE_TOLERANCE:.4f}]")

    deviation = abs(coverage - ds.COVERAGE_TARGET)
    assert deviation <= ds.COVERAGE_TOLERANCE, (
        f"measured coverage {coverage:.4f} is more than 3 SE "
        f"({ds.COVERAGE_TOLERANCE:.4f}) away from the nominal {ds.COVERAGE_TARGET:.2f}"
    )

    print(
        "OK: measured coverage is within three standard errors of the nominal "
        "95%. This IS the definition of a confidence interval -- not a claim "
        "about any one interval, a measured property of the procedure."
    )


if __name__ == "__main__":
    main()

"""Exercise 3 -- the test/interval duality.

A two-sided hypothesis test at level alpha rejects the null value exactly
when the (1 - alpha) confidence interval excludes it. Same underlying
machinery (the same z, the same standard error), presented two ways. This
is not a coincidence to be approximated -- it should hold EXACTLY, dataset
by dataset, because both are built from the identical z statistic.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import ci_excludes, confidence_interval_mean, one_sample_z_test_against_value  # noqa: E402


def main() -> None:
    rng = np.random.default_rng(7)
    n_datasets = 2000
    alpha = 0.05
    null_value = ds.POP_MEAN
    mismatches = 0
    rejections = 0

    for _ in range(n_datasets):
        n = int(rng.integers(15, 80))
        # Half the datasets are centred at the null value, half are shifted,
        # so both a "reject" and a "fail to reject" outcome actually occur.
        shift = 0.0 if rng.random() < 0.5 else rng.uniform(-6.0, 6.0)
        sample = ds.normal_population(rng, n, ds.POP_MEAN + shift, ds.POP_STD)

        _, p = one_sample_z_test_against_value(sample, null_value)
        rejects = p < alpha
        if rejects:
            rejections += 1

        interval = confidence_interval_mean(sample, alpha)
        excludes = ci_excludes(interval, null_value)

        if rejects != excludes:
            mismatches += 1

    print(f"Datasets checked: {n_datasets}")
    print(f"Test rejected the null in {rejections} of them")
    print(f"Test/interval disagreements: {mismatches}")

    assert mismatches == 0, (
        f"the test and the interval disagreed on {mismatches} of {n_datasets} "
        "datasets -- they should never disagree, since both come from the same z"
    )
    assert 0 < rejections < n_datasets, "the mix of shifted and unshifted datasets should produce both outcomes"

    print(
        "OK: across every dataset, rejecting at alpha and the interval "
        "excluding the null value agreed exactly -- zero mismatches."
    )


if __name__ == "__main__":
    main()

"""Exercise 4 -- the permutation test, from scratch.

No distributional assumption: shuffle the group labels, recompute the
statistic, repeat thousands of times, and read the p-value off how often a
shuffle produced something at least as extreme as what was actually
observed. Building this makes the meaning of "p-value" concrete in a way
a formula does not: it is a literal count of how surprising the real
arrangement was among all the ways the labels could have fallen.

Two cases:
  1. Two roughly normal populations at a moderate n, where the z-test's
     normal approximation is solid -- the two methods should land close.
  2. Two heavily right-skewed populations at a small n, where the normal
     approximation is shakier -- the two methods diverge more, and the
     permutation test needs no assumption about the population's shape to
     stay valid.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import permutation_test_diff_means, two_sample_z_test  # noqa: E402


def main() -> None:
    # Two independent generators, one per case, so neither case's random
    # draws depend on how many numbers the other case happened to consume.
    rng_case1 = np.random.default_rng(1)
    rng_case2 = np.random.default_rng(1)

    # Case 1: moderate n, normal populations -- the approximation should hold.
    a = ds.normal_population(rng_case1, 60, ds.POP_MEAN, ds.POP_STD)
    b = ds.normal_population(rng_case1, 60, ds.POP_B_MEAN, ds.POP_B_STD)
    z, p_z = two_sample_z_test(a, b)
    observed, p_perm = permutation_test_diff_means(a, b, 5000, rng_case1)
    print("Case 1 -- normal populations, n=60 each:")
    print(f"  observed difference in means: {observed:.4f}")
    print(f"  z-test p-value:          {p_z:.4f}")
    print(f"  permutation-test p-value: {p_perm:.4f}")
    diff_case1 = abs(p_z - p_perm)
    print(f"  |difference|: {diff_case1:.4f}")
    assert diff_case1 < 0.03, "at n=60 with normal populations, the two methods should agree closely"

    # Case 2: small n, right-skewed populations -- the normal approximation
    # is on shakier ground; the permutation test does not need it at all.
    c = ds.skewed_population(rng_case2, 8)
    d = ds.skewed_population(rng_case2, 8)
    z2, p_z2 = two_sample_z_test(c, d)
    observed2, p_perm2 = permutation_test_diff_means(c, d, 5000, rng_case2)
    print("\nCase 2 -- right-skewed populations, n=8 each:")
    print(f"  observed difference in means: {observed2:.4f}")
    print(f"  z-test p-value:          {p_z2:.4f}")
    print(f"  permutation-test p-value: {p_perm2:.4f}")
    diff_case2 = abs(p_z2 - p_perm2)
    print(f"  |difference|: {diff_case2:.4f}")
    # Both p-values must be sane probabilities, and the point of this case
    # is that they are allowed to diverge more than in Case 1 -- the
    # permutation test does not owe the normal approximation any agreement.
    assert 0.0 <= p_z2 <= 1.0 and 0.0 <= p_perm2 <= 1.0
    assert diff_case2 > diff_case1, (
        "the small-n, skewed case should show the two methods diverging "
        "more than the moderate-n, normal case did"
    )

    print(
        "\nOK: the permutation test agrees closely with the z-test where "
        "the normal approximation is solid, and diverges more where it is "
        "not -- while remaining a valid probability in both cases."
    )


if __name__ == "__main__":
    main()

"""Exercise 5 -- multiple comparisons and the Bonferroni correction.

A team that checks twenty independent metrics at alpha = 0.05 and calls
anything under 0.05 "significant" is not running one 5%-error test. It is
running a FAMILY of twenty tests, and the chance that at least one comes
back "significant" under pure noise is 1 - 0.95**20 -- an exact number,
derived below, confirmed by simulation, and then pulled back down with a
Bonferroni correction.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import bonferroni_alpha, z_critical_two_sided  # noqa: E402


def main() -> None:
    # --- The exact analytic value ---
    exact = 1 - (1 - ds.FWER_ALPHA) ** ds.FWER_TRIALS
    print(f"P(at least one false positive among {ds.FWER_TRIALS} independent "
          f"alpha={ds.FWER_ALPHA} tests) = 1 - (1-{ds.FWER_ALPHA})^{ds.FWER_TRIALS} "
          f"= {exact:.4f}")
    assert abs(exact - 0.6415) < 0.0001, f"expected 0.6415, got {exact:.4f}"

    # --- Confirm by simulation: many "families" of 20 independent tests
    # under a TRUE null (standard normal z-statistics), count how often at
    # least one exceeds the uncorrected critical value. ---
    rng = np.random.default_rng(42)
    z_crit = z_critical_two_sided(ds.FWER_ALPHA)
    at_least_one = 0
    for _ in range(ds.FWER_FAMILIES):
        zs = rng.standard_normal(ds.FWER_TRIALS)
        if np.any(np.abs(zs) > z_crit):
            at_least_one += 1
    simulated_fwer = at_least_one / ds.FWER_FAMILIES
    print(f"Simulated over {ds.FWER_FAMILIES} families: {simulated_fwer:.4f}")
    dev = abs(simulated_fwer - exact)
    assert dev <= ds.FWER_SIM_TOLERANCE, (
        f"simulated FWER {simulated_fwer:.4f} is too far from the exact "
        f"value {exact:.4f} (deviation {dev:.4f} > tolerance {ds.FWER_SIM_TOLERANCE})"
    )

    # --- Bonferroni: use alpha/m per test instead of alpha, and the
    # family-wise rate falls back down near the nominal alpha. ---
    corrected_alpha = bonferroni_alpha(ds.FWER_ALPHA, ds.FWER_TRIALS)
    print(f"Bonferroni-corrected per-test alpha: {ds.FWER_ALPHA}/{ds.FWER_TRIALS} = {corrected_alpha:.4f}")
    z_crit_corrected = z_critical_two_sided(corrected_alpha)
    at_least_one_corrected = 0
    for _ in range(ds.FWER_FAMILIES):
        zs = rng.standard_normal(ds.FWER_TRIALS)
        if np.any(np.abs(zs) > z_crit_corrected):
            at_least_one_corrected += 1
    simulated_bonferroni = at_least_one_corrected / ds.FWER_FAMILIES
    print(f"Simulated family-wise rate WITH Bonferroni: {simulated_bonferroni:.4f}")
    print(f"Analytic Bonferroni family-wise rate: {ds.BONFERRONI_EXPECTED:.4f}")
    dev_bonf = abs(simulated_bonferroni - ds.BONFERRONI_EXPECTED)
    assert dev_bonf <= ds.BONFERRONI_TOLERANCE, (
        f"Bonferroni-corrected simulated rate {simulated_bonferroni:.4f} is too far "
        f"from the analytic {ds.BONFERRONI_EXPECTED:.4f}"
    )
    assert abs(ds.BONFERRONI_EXPECTED - 0.0488) < 0.001, f"expected ~0.0488, got {ds.BONFERRONI_EXPECTED:.4f}"

    print(
        "\nOK: twenty uncorrected alpha=0.05 tests give a 64% chance of at "
        "least one false positive, confirmed by simulation. Bonferroni "
        "(alpha/m per test) pulls the family-wise rate back to about 4.9%. "
        "Most p-hacking is exactly this mechanism run silently, not fraud."
    )


if __name__ == "__main__":
    main()

"""Exercise 6 -- statistical power.

Power is P(reject H0 | H0 is false by a specific amount). It depends on
three things together: the true effect size, the sample size, and alpha.
"Found nothing" without a power figure is not evidence of absence -- it
might just mean the test never had a real chance to detect the effect
that was actually there.

This exercise computes power two ways: from the closed-form formula, and
by simulating the test itself thousands of times at a true effect size and
counting how often it actually rejects. The two must agree.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import power_two_sample_z, two_sample_z_test  # noqa: E402


def simulate_power(rng: np.random.Generator, effect: float, sigma: float, n: int, trials: int) -> float:
    rejections = 0
    for _ in range(trials):
        a = ds.normal_population(rng, n, ds.POP_MEAN, sigma)
        b = ds.normal_population(rng, n, ds.POP_MEAN + effect, sigma)
        _, p = two_sample_z_test(a, b)
        if p < 0.05:
            rejections += 1
    return rejections / trials


def main() -> None:
    sigma = ds.POP_STD

    # Power rises with n, at a fixed effect size.
    print(f"Power vs n, at a fixed effect of {ds.POWER_CHECK_EFFECT} (sigma={sigma}):")
    prev_power = -1.0
    for n in (10, 20, 50, 100, 200, 400):
        power = power_two_sample_z(ds.POWER_CHECK_EFFECT, sigma, n)
        print(f"  n={n:>4}: power={power:.4f}")
        assert power > prev_power, f"power should rise monotonically with n, failed at n={n}"
        prev_power = power

    # Power rises with effect size, at a fixed n.
    print(f"\nPower vs effect size, at a fixed n=100 (sigma={sigma}):")
    prev_power = -1.0
    for effect in (0.5, 1.0, 2.0, 4.0, 8.0):
        power = power_two_sample_z(effect, sigma, 100)
        print(f"  effect={effect:>4}: power={power:.4f}")
        assert power > prev_power, f"power should rise monotonically with effect size, failed at effect={effect}"
        prev_power = power

    # Check the formula against a direct simulation at one configuration.
    theoretical = power_two_sample_z(ds.POWER_CHECK_EFFECT, sigma, ds.POWER_CHECK_N)
    rng = np.random.default_rng(42)
    simulated = simulate_power(rng, ds.POWER_CHECK_EFFECT, sigma, ds.POWER_CHECK_N, ds.POWER_CHECK_TRIALS)
    print(
        f"\nConfiguration: effect={ds.POWER_CHECK_EFFECT}, n={ds.POWER_CHECK_N} per group, sigma={sigma}, alpha=0.05"
    )
    print(f"  theoretical power: {theoretical:.4f}")
    print(f"  simulated power ({ds.POWER_CHECK_TRIALS} trials): {simulated:.4f}")
    dev = abs(theoretical - simulated)
    assert dev <= ds.POWER_CHECK_TOLERANCE, (
        f"theoretical power {theoretical:.4f} and simulated power {simulated:.4f} "
        f"disagree by {dev:.4f}, more than the {ds.POWER_CHECK_TOLERANCE} tolerance"
    )

    print(
        "\nOK: power rises monotonically with both n and effect size, and "
        "the closed-form formula agrees with a direct simulation of the "
        "test itself."
    )


if __name__ == "__main__":
    main()

"""Exercise 7 -- effect size versus significance.

The same relative difference -- a fixed 0.5% shift in the population mean
-- is tested at a small sample size and at an enormous one. At small n it
is not statistically significant; at huge n it is, reliably. The
underlying effect size never changed. Only the amount of data collected
around it did. "Statistically significant" is a statement about whether
noise can be ruled out, not about whether an effect is large enough to
matter.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import two_sample_z_test  # noqa: E402


def main() -> None:
    mean, sigma = ds.POP_MEAN, ds.POP_STD
    effect = mean * ds.EFFECT_VS_N_RELATIVE_DIFF
    cohens_d = effect / sigma

    print(f"Population mean: {mean}, population std: {sigma}")
    print(f"Fixed relative difference: {ds.EFFECT_VS_N_RELATIVE_DIFF:.1%}")
    print(f"Absolute effect size (population): {effect:.4f}")
    print(f"Standardized effect size (Cohen's d = effect/sigma): {cohens_d:.5f}")

    rng_small = np.random.default_rng(42)
    a_small = ds.normal_population(rng_small, ds.EFFECT_VS_N_SMALL_N, mean, sigma)
    b_small = ds.normal_population(rng_small, ds.EFFECT_VS_N_SMALL_N, mean + effect, sigma)
    _, p_small = two_sample_z_test(a_small, b_small)

    rng_large = np.random.default_rng(42)
    a_large = ds.normal_population(rng_large, ds.EFFECT_VS_N_LARGE_N, mean, sigma)
    b_large = ds.normal_population(rng_large, ds.EFFECT_VS_N_LARGE_N, mean + effect, sigma)
    _, p_large = two_sample_z_test(a_large, b_large)

    print(f"\nn={ds.EFFECT_VS_N_SMALL_N} per group: p = {p_small:.4f}")
    print(f"n={ds.EFFECT_VS_N_LARGE_N} per group: p = {p_large:.6f}")

    assert p_small > 0.05, f"expected the small-n case to be NOT significant, got p={p_small:.4f}"
    assert p_large < 0.05, f"expected the huge-n case to be significant, got p={p_large:.6f}"

    # The effect size itself -- the population parameter the two samples
    # were built around -- is a single fixed number, untouched by n.
    cohens_d_small = effect / sigma
    cohens_d_large = effect / sigma
    assert cohens_d_small == cohens_d_large == cohens_d, "the effect size is a population parameter, independent of n"

    print(
        "\nOK: the exact same 0.5% relative difference is not significant "
        f"at n={ds.EFFECT_VS_N_SMALL_N} and is significant at n={ds.EFFECT_VS_N_LARGE_N}. "
        "The effect size (Cohen's d) is identical in both cases -- only the "
        "power to detect it changed. A huge n can make a trivial difference "
        "'significant'; significance alone never tells you whether the "
        "difference is big enough to matter."
    )


if __name__ == "__main__":
    main()

"""Exercise 8 -- peeking, the most common real sin.

A test's alpha (say 5%) is a promise about ONE look at ONE fixed sample
size. Checking the p-value after every batch of new data and stopping the
instant it dips below 0.05 is a different procedure entirely -- and its
real false-positive rate is far higher than 5%, even though every single
p-value computed along the way was calculated correctly.

This simulates many independent "experiments" under a population where
the null hypothesis is TRUE (mean exactly 0), each one collecting data in
batches of 10 and stopping at the first p < 0.05, and measures how often
that early stop happens purely from noise.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import one_sample_z_test_against_value  # noqa: E402


def run_experiment(rng: np.random.Generator) -> bool:
    """Returns True if peeking produced a false positive on this experiment."""
    data: list[float] = []
    for _ in range(ds.PEEK_MAX_BATCHES):
        data.extend(rng.normal(0.0, 1.0, ds.PEEK_BATCH_SIZE).tolist())
        arr = np.array(data)
        _, p = one_sample_z_test_against_value(arr, 0.0)
        if p < ds.PEEK_ALPHA:
            return True
    return False


def fixed_n_false_positive_rate(rng: np.random.Generator, n: int, trials: int) -> float:
    """The honest comparison: test ONCE at a fixed final sample size."""
    false_positives = 0
    for _ in range(trials):
        arr = rng.normal(0.0, 1.0, n)
        _, p = one_sample_z_test_against_value(arr, 0.0)
        if p < ds.PEEK_ALPHA:
            false_positives += 1
    return false_positives / trials


def main() -> None:
    rng = np.random.default_rng(42)
    peeked_false_positives = sum(run_experiment(rng) for _ in range(ds.PEEK_EXPERIMENTS))
    peeked_rate = peeked_false_positives / ds.PEEK_EXPERIMENTS

    final_n = ds.PEEK_BATCH_SIZE * ds.PEEK_MAX_BATCHES
    fixed_rate = fixed_n_false_positive_rate(rng, final_n, ds.PEEK_EXPERIMENTS)

    print(f"True null hypothesis (population mean is exactly 0), alpha={ds.PEEK_ALPHA}")
    print(f"Experiments simulated: {ds.PEEK_EXPERIMENTS}")
    print(f"Looks per experiment: {ds.PEEK_MAX_BATCHES} (every {ds.PEEK_BATCH_SIZE} observations, "
          f"up to n={final_n})")
    print(f"\nFalse-positive rate WITH peeking (stop at first p<0.05): {peeked_rate:.4f}")
    print(f"False-positive rate testing ONCE at n={final_n} (honest, no peeking): {fixed_rate:.4f}")
    print(f"Inflation factor: {peeked_rate / ds.PEEK_ALPHA:.2f}x nominal alpha")

    assert peeked_rate >= ds.PEEK_ALPHA * ds.PEEK_MIN_INFLATION_FACTOR, (
        f"peeking false-positive rate {peeked_rate:.4f} should be at least "
        f"{ds.PEEK_MIN_INFLATION_FACTOR}x alpha ({ds.PEEK_ALPHA * ds.PEEK_MIN_INFLATION_FACTOR:.4f})"
    )
    assert abs(fixed_rate - ds.PEEK_ALPHA) < 0.02, (
        f"the honest fixed-n rate {fixed_rate:.4f} should sit close to alpha={ds.PEEK_ALPHA}"
    )
    assert peeked_rate > fixed_rate, "peeking must produce a higher false-positive rate than the honest fixed-n test"

    print(
        "\nOK: testing repeatedly and stopping at the first p<0.05 inflates "
        "the true false-positive rate several times past the nominal alpha, "
        "even though every individual p-value was computed correctly and the "
        "null hypothesis was true the entire time. The fix is deciding the "
        "sample size in advance (or using a sequential-testing method built "
        "for repeated looks), not testing whenever the data happens to look "
        "promising."
    )


if __name__ == "__main__":
    main()

"""Exercise 9 -- bootstrap interval versus the normal-approximation interval.

Day 117 built the bootstrap from scratch: resample with replacement,
recompute the statistic, read the spread. Today's normal-approximation
interval (mean +/- z * standard_error) is a closed-form shortcut that is
valid under the same conditions the CLT needs. Where both are valid --
here, the mean of a reasonably-sized sample from a well-behaved
population -- they should agree closely, without needing to be identical.
"""
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).parent))

import numpy as np  # noqa: E402

import dataset as ds  # noqa: E402
from inference import bootstrap_ci, confidence_interval_mean  # noqa: E402


def main() -> None:
    rng = np.random.default_rng(42)
    n = 200
    sample = ds.normal_population(rng, n, ds.POP_MEAN, ds.POP_STD)

    normal_lo, normal_hi = confidence_interval_mean(sample, alpha=0.05)
    boot_lo, boot_hi = bootstrap_ci(sample, np.mean, ds.BOOTSTRAP_N_BOOT, 0.05, rng)

    normal_center = (normal_lo + normal_hi) / 2
    boot_center = (boot_lo + boot_hi) / 2
    normal_width = normal_hi - normal_lo
    boot_width = boot_hi - boot_lo
    se = sample.std(ddof=1) / np.sqrt(n)

    print(f"Sample: n={n}, mean={sample.mean():.4f}, std={sample.std(ddof=1):.4f}, SE={se:.4f}")
    print(f"Normal-approximation 95% CI:  [{normal_lo:.4f}, {normal_hi:.4f}]  width={normal_width:.4f}")
    print(f"Bootstrap ({ds.BOOTSTRAP_N_BOOT} resamples) 95% CI: [{boot_lo:.4f}, {boot_hi:.4f}]  "
          f"width={boot_width:.4f}")

    center_diff_in_se = abs(normal_center - boot_center) / se
    width_ratio = boot_width / normal_width
    print(f"Center difference: {center_diff_in_se:.4f} standard errors")
    print(f"Width ratio (bootstrap/normal): {width_ratio:.4f}")

    assert center_diff_in_se <= ds.BOOTSTRAP_CENTER_TOLERANCE_IN_SE, (
        f"centers differ by {center_diff_in_se:.4f} SE, more than the "
        f"{ds.BOOTSTRAP_CENTER_TOLERANCE_IN_SE} SE tolerance"
    )
    assert abs(width_ratio - 1.0) <= ds.BOOTSTRAP_WIDTH_RATIO_TOLERANCE, (
        f"width ratio {width_ratio:.4f} is more than "
        f"{ds.BOOTSTRAP_WIDTH_RATIO_TOLERANCE} away from 1.0"
    )

    print(
        "\nOK: the bootstrap interval and the normal-approximation interval "
        "agree closely where both are valid -- the bootstrap needed no "
        "formula for the standard error of the mean, and would work exactly "
        "the same way for a statistic that has no such formula."
    )


if __name__ == "__main__":
    main()

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

Both `examples/` and `starter/` contain modules called `inference` 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 ("inference", "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]

"""Population generators and tolerance constants shared by every exercise.

Every population here has a KNOWN true parameter, which is what makes the
coverage and duality exercises checkable at all: you cannot verify that a
confidence interval has 95% coverage without knowing the true value the
interval is supposed to be catching.

Tolerances below were sanity-checked by rerunning each exercise's core
simulation across seeds 1, 2, 3, 42 and 118 during development; the
comment beside each constant records the range actually observed.
"""
from __future__ import annotations

import numpy as np

# A population with a known mean and standard deviation, used for the
# z-test and coverage exercises. Deliberately not integers so that no
# accidental symmetry hides a bug.
POP_MEAN = 50.3
POP_STD = 12.7

# A second, shifted population -- same shape, different mean -- for the
# two-sample tests.
POP_B_MEAN = 53.1
POP_B_STD = 12.7

# A right-skewed population (exponential) for the permutation test's
# "where the normal approximation does not hold well at small n" arm.
SKEWED_SCALE = 8.0  # exponential scale (mean == scale)


def normal_population(rng: np.random.Generator, n: int, mean: float, std: float) -> np.ndarray:
    return rng.normal(loc=mean, scale=std, size=n)


def skewed_population(rng: np.random.Generator, n: int, scale: float = SKEWED_SCALE) -> np.ndarray:
    return rng.exponential(scale=scale, size=n)


# --- Tolerances, each derived from a standard error or a generous sanity
# band, not picked to make a single run pass. Every tolerance below was
# checked by rerunning the relevant simulation across seeds 1, 2, 3, 42
# and 118 during development; the observed range for that exact seed set
# is recorded next to each one. ---

# Exercise 2 (coverage): with n=300 per sample (large enough that the
# z-critical value is a good stand-in for the t-critical value -- at
# n=40 the true coverage undershoots 95% by about a point because of
# exactly that gap) and 10,000 trials of a Bernoulli(0.95) "did this
# interval cover?" indicator, the standard error of the measured coverage
# is sqrt(0.95 * 0.05 / 10000) = 0.00218, so three standard errors is
# 0.00654. Observed across seeds 1/2/3/42/118: 0.9468-0.9526, i.e. every
# seed landed within 0.0032 of 0.95, comfortably inside the 3-SE band.
COVERAGE_TARGET = 0.95
COVERAGE_TRIALS = 10_000
COVERAGE_SAMPLE_N = 300
COVERAGE_SE = (COVERAGE_TARGET * (1 - COVERAGE_TARGET) / COVERAGE_TRIALS) ** 0.5
COVERAGE_TOLERANCE = 3 * COVERAGE_SE

# Exercise 5 (multiple comparisons): the analytic family-wise error rate
# for 20 independent alpha=0.05 tests is exact: 1 - 0.95**20 = 0.641514...
FWER_TRIALS = 20
FWER_ALPHA = 0.05
FWER_EXACT = 1 - (1 - FWER_ALPHA) ** FWER_TRIALS  # 0.6415140775914581
# Simulated with FWER_FAMILIES families of 20 independent standard-normal
# z-statistics; SE of the simulated FWER at 20,000 families is
# sqrt(0.6415*0.3585/20000) = 0.0034, so 3 SE = 0.0102. Observed across
# five seeds: simulated FWER 0.6388-0.6436, always within 0.0027 of the
# exact value; the Bonferroni-corrected simulation (exact target 0.0488)
# came in 0.0479-0.0528, within 0.004.
FWER_FAMILIES = 20_000
FWER_SIM_TOLERANCE = 0.015
BONFERRONI_EXPECTED = 1 - (1 - FWER_ALPHA / FWER_TRIALS) ** FWER_TRIALS  # 0.048830...
BONFERRONI_TOLERANCE = 0.015

# Exercise 8 (peeking): under a true null, checking after every 10
# observations for up to PEEK_MAX_BATCHES looks and stopping at the first
# p < 0.05 is known to inflate the false-positive rate well past alpha --
# the more looks allowed, the worse it gets. This lab treats "far above
# 0.05" as at least PEEK_MIN_INFLATION_FACTOR times alpha. Observed across
# five seeds with 4,000 simulated experiments and 5 looks: the false-
# positive rate ranged 0.1668-0.1888, roughly 3.3x-3.8x alpha -- well
# clear of the 2x floor used as the assertion.
PEEK_ALPHA = 0.05
PEEK_EXPERIMENTS = 4_000
PEEK_BATCH_SIZE = 10
PEEK_MAX_BATCHES = 5  # up to 50 observations, 5 looks
PEEK_MIN_INFLATION_FACTOR = 2.0  # false-positive rate must be >= 2x alpha

# Exercise 9 (bootstrap vs normal CI): both are estimating the same
# interval for the mean of a normal population with n=200; they need not
# match to the decimal, but their centers should agree closely (in units
# of the normal interval's own standard error) and their widths should be
# within about 20% of each other. Observed across five seeds: center
# difference 0.003-0.053 standard errors (comfortably inside a 0.6-SE
# band), width ratio 0.972-1.023 (comfortably inside a 20% band).
BOOTSTRAP_N_BOOT = 5_000
BOOTSTRAP_WIDTH_RATIO_TOLERANCE = 0.20
BOOTSTRAP_CENTER_TOLERANCE_IN_SE = 0.6  # in units of the normal CI's own SE

# Exercise 6 (power): the closed-form power formula is checked against a
# simulated rejection rate at n=100, effect=2.8, sigma=POP_STD, 3,000
# simulated trials. SE of a simulated rate near 0.34 at 3,000 trials is
# sqrt(0.34*0.66/3000) = 0.0087, so 3 SE = 0.026. Observed across five
# seeds: simulated power 0.336-0.355 against a theoretical 0.3444, always
# within 0.011.
POWER_CHECK_EFFECT = 2.8
POWER_CHECK_N = 100
POWER_CHECK_TRIALS = 3_000
POWER_CHECK_TOLERANCE = 0.03

# Exercise 7 (effect size vs n): a fixed 0.5% relative difference in the
# population mean is tested at a small n and an enormous n. Observed
# across five seeds at n=30: p ranged 0.20-0.86 (never significant); at
# n=100000: p ranged 0.0 to 0.0008 (always significant at alpha=0.05).
EFFECT_VS_N_RELATIVE_DIFF = 0.005
EFFECT_VS_N_SMALL_N = 30
EFFECT_VS_N_LARGE_N = 100_000

"""Hand-rolled hypothesis-testing and confidence-interval machinery.

Every function here is built from `math.erf` and NumPy array arithmetic
only -- no `scipy.stats`, no `statsmodels`. That is the point of Day 118:
before trusting a library's `ttest_ind`, build the thing it computes.

All functions are pure (no hidden global state) and take an explicit
`numpy.random.Generator` wherever randomness is needed, per the
project's convention of never using an internally reseeded generator.
"""
from __future__ import annotations

import math

import numpy as np


def phi(z: float) -> float:
    """Standard normal CDF, computed from math.erf (no scipy needed)."""
    return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))


def p_from_z_two_sided(z: float) -> float:
    """Two-sided p-value for a standard normal test statistic z."""
    return 2.0 * (1.0 - phi(abs(z)))


def z_critical_two_sided(alpha: float) -> float:
    """The z such that P(|Z| > z) = alpha, found by bisection on phi.

    There is no closed form for the normal quantile function in terms of
    erf, so this inverts phi() numerically. 200 bisection steps on a
    [0, 10] bracket resolves z to well under 1e-9, far tighter than any
    tolerance this lab uses.
    """
    target = 1.0 - alpha / 2.0
    lo, hi = 0.0, 10.0
    for _ in range(200):
        mid = (lo + hi) / 2.0
        if phi(mid) < target:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2.0


def two_sample_z_test(a: np.ndarray, b: np.ndarray) -> tuple[float, float]:
    """Welch-style two-sample z-test (unpooled variance, large-n normal
    approximation to the sampling distribution of the difference in means).

    Returns (z, p_two_sided). This is the large-sample cousin of Welch's
    t-test (see the lesson's Tools section): it uses each sample's own
    variance rather than assuming the two populations share one variance,
    and it treats the standard error as known rather than estimated,
    which is accurate once each sample has on the order of 30+
    observations.
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    n_a, n_b = a.size, b.size
    mean_a, mean_b = a.mean(), b.mean()
    var_a, var_b = a.var(ddof=1), b.var(ddof=1)
    se = math.sqrt(var_a / n_a + var_b / n_b)
    z = (mean_a - mean_b) / se
    p = p_from_z_two_sided(z)
    return z, p


def confidence_interval_mean(sample: np.ndarray, alpha: float = 0.05) -> tuple[float, float]:
    """Normal-approximation (1 - alpha) confidence interval for a mean.

    mean +/- z_(alpha/2) * standard_error, where standard_error is the
    sample standard deviation divided by sqrt(n) -- Day 117's standard
    error, reused rather than re-derived.
    """
    sample = np.asarray(sample, dtype=float)
    n = sample.size
    mean = sample.mean()
    se = sample.std(ddof=1) / math.sqrt(n)
    z = z_critical_two_sided(alpha)
    return mean - z * se, mean + z * se


def ci_excludes(interval: tuple[float, float], value: float) -> bool:
    """True if `value` lies strictly outside the closed interval."""
    lo, hi = interval
    return value < lo or value > hi


def one_sample_z_test_against_value(sample: np.ndarray, null_value: float) -> tuple[float, float]:
    """One-sample z-test of H0: population mean == null_value."""
    sample = np.asarray(sample, dtype=float)
    n = sample.size
    se = sample.std(ddof=1) / math.sqrt(n)
    z = (sample.mean() - null_value) / se
    return z, p_from_z_two_sided(z)


def permutation_test_diff_means(
    a: np.ndarray, b: np.ndarray, n_perm: int, rng: np.random.Generator
) -> tuple[float, float]:
    """Two-sided permutation test for a difference in means.

    No distributional assumption: the null hypothesis is that the group
    label carries no information, so the labels are shuffled `n_perm`
    times, the difference in means is recomputed each time, and the
    p-value is the fraction of shuffles at least as extreme as what was
    actually observed (plus the observed arrangement itself, so the
    p-value can never read as exactly zero).

    Returns (observed_diff, p_two_sided).
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    n_a = a.size
    pooled = np.concatenate([a, b])
    observed = a.mean() - b.mean()
    count_as_extreme = 0
    for _ in range(n_perm):
        shuffled = rng.permutation(pooled)
        perm_diff = shuffled[:n_a].mean() - shuffled[n_a:].mean()
        if abs(perm_diff) >= abs(observed):
            count_as_extreme += 1
    p = (count_as_extreme + 1) / (n_perm + 1)
    return observed, p


def power_two_sample_z(effect: float, sigma: float, n_per_group: int, alpha: float = 0.05) -> float:
    """Power of the two-sample z-test above a true mean difference `effect`,
    with a common known standard deviation `sigma`, `n_per_group` in each
    arm, assuming both samples are the same size.

    Derivation: under H1 the test statistic Z = (Xbar_a - Xbar_b)/SE is
    Normal(effect/SE, 1), where SE = sigma * sqrt(2/n_per_group). The test
    rejects when |Z| > z_crit under H0, so power is the probability that a
    Normal(effect/SE, 1) variable falls outside [-z_crit, z_crit].
    """
    se = sigma * math.sqrt(2.0 / n_per_group)
    z_crit = z_critical_two_sided(alpha)
    shift = effect / se
    # P(Z > z_crit) + P(Z < -z_crit) for Z ~ Normal(shift, 1)
    return (1.0 - phi(z_crit - shift)) + phi(-z_crit - shift)


def bonferroni_alpha(alpha: float, m: int) -> float:
    """The per-test alpha that keeps the family-wise error rate at `alpha`
    across `m` independent tests, by the (conservative) Bonferroni bound.
    """
    return alpha / m


def bootstrap_ci(
    sample: np.ndarray,
    statistic,
    n_boot: int,
    alpha: float,
    rng: np.random.Generator,
) -> tuple[float, float]:
    """Percentile bootstrap confidence interval for an arbitrary statistic.

    Resample `sample` with replacement `n_boot` times, recompute
    `statistic` on each resample, and take the [alpha/2, 1 - alpha/2]
    percentiles of the resulting distribution. Built the same way as Day
    117's bootstrap, reused here for a statistic (a mean) that also has a
    closed-form normal-approximation interval, so the two can be checked
    against each other.
    """
    sample = np.asarray(sample, dtype=float)
    n = sample.size
    replicates = np.empty(n_boot)
    for i in range(n_boot):
        resample = rng.choice(sample, size=n, replace=True)
        replicates[i] = statistic(resample)
    lo = np.percentile(replicates, 100 * alpha / 2)
    hi = np.percentile(replicates, 100 * (1 - alpha / 2))
    return float(lo), float(hi)

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

Run from the lab directory:

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

import math

import numpy as np
import pytest

import dataset as D
import inference as I


# --------------------------------------------------------------------------
# phi / p_from_z_two_sided / z_critical_two_sided
# --------------------------------------------------------------------------


def test_phi_of_zero_is_one_half():
    assert I.phi(0.0) == pytest.approx(0.5)


def test_phi_matches_known_normal_table_values():
    assert I.phi(1.96) == pytest.approx(0.9750, abs=0.0001)
    assert I.phi(-1.96) == pytest.approx(0.0250, abs=0.0001)


def test_z_critical_two_sided_matches_known_constants():
    assert I.z_critical_two_sided(0.05) == pytest.approx(1.959964, abs=1e-5)
    assert I.z_critical_two_sided(0.01) == pytest.approx(2.575829, abs=1e-5)


def test_p_from_z_two_sided_round_trips_with_z_critical():
    z = I.z_critical_two_sided(0.05)
    assert I.p_from_z_two_sided(z) == pytest.approx(0.05, abs=1e-6)


# --------------------------------------------------------------------------
# two_sample_z_test
# --------------------------------------------------------------------------


def test_two_sample_z_test_matches_hand_computation():
    a = [50, 52, 49, 51, 53, 48, 50, 52, 51, 49]
    b = [54, 55, 53, 56, 54, 52, 55, 53, 54, 56]
    z, p = I.two_sample_z_test(a, b)
    import statistics

    mean_a, var_a = statistics.mean(a), statistics.variance(a)
    mean_b, var_b = statistics.mean(b), statistics.variance(b)
    se = math.sqrt(var_a / len(a) + var_b / len(b))
    z_hand = (mean_a - mean_b) / se
    assert z == pytest.approx(z_hand, abs=1e-9)
    assert p < 0.001


def test_two_sample_z_test_identical_samples_gives_p_near_one():
    a = [10.0, 20.0, 30.0, 40.0, 50.0]
    z, p = I.two_sample_z_test(a, a)
    assert z == pytest.approx(0.0, abs=1e-9)
    assert p == pytest.approx(1.0, abs=1e-9)


# --------------------------------------------------------------------------
# confidence_interval_mean / ci_excludes
# --------------------------------------------------------------------------


def test_confidence_interval_mean_is_centered_on_the_sample_mean():
    rng = np.random.default_rng(1)
    sample = D.normal_population(rng, 100, D.POP_MEAN, D.POP_STD)
    lo, hi = I.confidence_interval_mean(sample, alpha=0.05)
    assert (lo + hi) / 2 == pytest.approx(sample.mean(), abs=1e-9)
    assert lo < sample.mean() < hi


def test_wider_interval_for_smaller_alpha():
    rng = np.random.default_rng(2)
    sample = D.normal_population(rng, 100, D.POP_MEAN, D.POP_STD)
    lo95, hi95 = I.confidence_interval_mean(sample, alpha=0.05)
    lo99, hi99 = I.confidence_interval_mean(sample, alpha=0.01)
    assert (hi99 - lo99) > (hi95 - lo95)


def test_ci_excludes():
    assert I.ci_excludes((1.0, 2.0), 3.0) is True
    assert I.ci_excludes((1.0, 2.0), 1.5) is False
    assert I.ci_excludes((1.0, 2.0), 1.0) is False  # boundary counts as inside


def test_coverage_is_close_to_nominal():
    rng = np.random.default_rng(42)
    hits = 0
    trials = 2000
    for _ in range(trials):
        sample = D.normal_population(rng, D.COVERAGE_SAMPLE_N, D.POP_MEAN, D.POP_STD)
        lo, hi = I.confidence_interval_mean(sample, alpha=0.05)
        if lo <= D.POP_MEAN <= hi:
            hits += 1
    coverage = hits / trials
    se = math.sqrt(0.95 * 0.05 / trials)
    assert abs(coverage - 0.95) <= 3 * se


# --------------------------------------------------------------------------
# duality: test rejects <=> interval excludes the null value
# --------------------------------------------------------------------------


def test_duality_holds_exactly_across_many_datasets():
    rng = np.random.default_rng(9)
    for _ in range(200):
        n = int(rng.integers(15, 60))
        shift = 0.0 if rng.random() < 0.5 else rng.uniform(-6.0, 6.0)
        sample = D.normal_population(rng, n, D.POP_MEAN + shift, D.POP_STD)
        _, p = I.one_sample_z_test_against_value(sample, D.POP_MEAN)
        interval = I.confidence_interval_mean(sample, 0.05)
        assert (p < 0.05) == I.ci_excludes(interval, D.POP_MEAN)


# --------------------------------------------------------------------------
# permutation_test_diff_means
# --------------------------------------------------------------------------


def test_permutation_test_returns_valid_probability():
    rng = np.random.default_rng(3)
    a = D.normal_population(rng, 30, D.POP_MEAN, D.POP_STD)
    b = D.normal_population(rng, 30, D.POP_B_MEAN, D.POP_B_STD)
    _, p = I.permutation_test_diff_means(a, b, 500, rng)
    assert 0.0 <= p <= 1.0


def test_permutation_test_p_is_small_for_a_large_true_difference():
    rng = np.random.default_rng(4)
    a = D.normal_population(rng, 50, 0.0, 1.0)
    b = D.normal_population(rng, 50, 20.0, 1.0)
    _, p = I.permutation_test_diff_means(a, b, 500, rng)
    assert p < 0.01


# --------------------------------------------------------------------------
# power_two_sample_z
# --------------------------------------------------------------------------


def test_power_is_between_alpha_and_one():
    power = I.power_two_sample_z(effect=5.0, sigma=12.7, n_per_group=50)
    assert 0.05 < power < 1.0


def test_power_increases_with_n():
    p1 = I.power_two_sample_z(effect=3.0, sigma=12.7, n_per_group=20)
    p2 = I.power_two_sample_z(effect=3.0, sigma=12.7, n_per_group=200)
    assert p2 > p1


def test_power_increases_with_effect_size():
    p1 = I.power_two_sample_z(effect=1.0, sigma=12.7, n_per_group=100)
    p2 = I.power_two_sample_z(effect=10.0, sigma=12.7, n_per_group=100)
    assert p2 > p1


def test_power_at_zero_effect_equals_alpha():
    # With no true effect, "power" is just the false-positive rate: alpha.
    power = I.power_two_sample_z(effect=0.0, sigma=12.7, n_per_group=100, alpha=0.05)
    assert power == pytest.approx(0.05, abs=1e-6)


# --------------------------------------------------------------------------
# bonferroni_alpha
# --------------------------------------------------------------------------


def test_bonferroni_alpha_divides_by_m():
    assert I.bonferroni_alpha(0.05, 20) == pytest.approx(0.0025)


def test_exact_family_wise_error_rate_for_twenty_tests():
    exact = 1 - (1 - 0.05) ** 20
    assert exact == pytest.approx(0.6415, abs=0.0001)


def test_exact_bonferroni_corrected_rate_for_twenty_tests():
    corrected = 1 - (1 - 0.05 / 20) ** 20
    assert corrected == pytest.approx(0.0488, abs=0.0001)


# --------------------------------------------------------------------------
# bootstrap_ci
# --------------------------------------------------------------------------


def test_bootstrap_ci_contains_the_sample_mean_region():
    rng = np.random.default_rng(5)
    sample = D.normal_population(rng, 200, D.POP_MEAN, D.POP_STD)
    lo, hi = I.bootstrap_ci(sample, np.mean, 1000, 0.05, rng)
    assert lo < sample.mean() < hi


def test_bootstrap_ci_agrees_with_normal_ci_for_the_mean():
    rng = np.random.default_rng(6)
    sample = D.normal_population(rng, 200, D.POP_MEAN, D.POP_STD)
    normal_lo, normal_hi = I.confidence_interval_mean(sample, 0.05)
    boot_lo, boot_hi = I.bootstrap_ci(sample, np.mean, 3000, 0.05, rng)
    normal_center = (normal_lo + normal_hi) / 2
    boot_center = (boot_lo + boot_hi) / 2
    se = sample.std(ddof=1) / math.sqrt(200)
    assert abs(normal_center - boot_center) / se < 1.0
    width_ratio = (boot_hi - boot_lo) / (normal_hi - normal_lo)
    assert 0.7 < width_ratio < 1.3

lesson_id: D118
day: 118
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-118-hypothesis-tests-and-confidence-intervals
  - 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_two_sample_z_test.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_coverage.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_duality.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_permutation_test.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_multiple_comparisons.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_power.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_effect_size_vs_n.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 08_peeking.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 09_bootstrap_vs_normal_ci.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-19'
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 (captured directly, not through a pipeline). pytest examples -> 22 passed; pytest starter -> 1 passed, 15 skipped on an untouched checkout, and 16 passed against a fully solved copy of starter/ (verified by temporarily copying the reference inference.py into starter/, confirming all 16 tests passed, then restoring the blank skeleton -- the skip count after restoring was confirmed back to 15, and auto-discovering both suites with a bare `pytest -q` from the lab directory also reports 15 skipped, proving the two conftest.py import guards work for that invocation). All nine reference scripts exit 0 with every internal assertion holding, each ending with a line starting "OK:". 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 05 (multiple comparisons) with its expected exact family-wise error rate temporarily replaced with a deliberately wrong value (99.0 instead of 0.6415), confirms the run exits non-zero with the named AssertionError showing the real value, and does not touch the file on disk -- so the suite is demonstrated capable of failing rather than merely claimed to be. Separately during development, a real bug was introduced directly into inference.py''s two_sample_z_test (adding a constant offset to z) and the harness correctly reported 10 of 32 checks failing before the fix was reverted and a clean 32/32 run was reconfirmed. Four honesty notes from this run. FIRST: scipy and statsmodels are not installed in this environment. scipy.stats.ttest_ind and scipy.stats.norm.interval, and statsmodels.stats.multitest, 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 (measured coverage, permutation p-values, simulated family-wise error rates, simulated power, the peeking false-positive rate, and the bootstrap interval bounds) is a freshly measured number rather than a fixed literal, checked against a tolerance derived from a standard error or a generous, explicitly-derived sanity band rather than a value chosen to make the test pass. THIRD: every tolerance in dataset.py was checked by rerunning the exercise logic across five different seeds (1, 2, 3, 42, 118) 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 five runs, not as a formal guarantee no seed can ever violate it; the coverage exercise specifically uses n=300 per sample rather than a smaller n because at n=40 the true coverage of a normal-approximation interval undershoots 95% by roughly a point (a real, measured effect of using a normal rather than a t critical value at small n), which was observed directly during development and is recorded in troubleshooting.md rather than papered over. FOURTH: an explicit `pytest examples starter` invocation (both directories passed as two separate command-line arguments) was found during development to collide -- the starter suite picks up the reference inference.py and reports all tests passing instead of skipping -- while the two invocations this lab actually documents and tests (`pytest starter` alone, and bare `pytest` with no path argument, auto-discovering both directories from the lab root) both isolate correctly; this lab documents and tests only the latter two, and troubleshooting.md and run_tests.sh section 4 say so explicitly rather than silently avoiding the untested invocation. The exact family-wise error rate for twenty tests (0.6415) and its Bonferroni-corrected counterpart (0.0488) are closed-form arithmetic and 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 thousands of confidence intervals and permutation shuffles at once. |
| `pytest` | 9.1.1 | MIT | The reference suite (22 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 standard library's `math` module
(specifically `math.erf`) supplies every normal-distribution calculation --
no statistical package is needed for that part at all.

## 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`** would replace most of `examples/inference.py` with a
handful of function calls: `scipy.stats.ttest_ind(a, b, equal_var=False)`
runs Welch's t-test in one line, and `scipy.stats.norm.interval` builds a
confidence interval directly. 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. `inference.py`'s
`two_sample_z_test` and `confidence_interval_mean` are the exact ideas
`ttest_ind` and `norm.interval` implement, minus the engineering: a t
rather than a normal reference distribution for small samples, several
named alternatives, vectorised batch operation, and a stable public API.
Having written the from-scratch version, `scipy.stats.ttest_ind`'s
documentation reads as an implementation detail rather than a black box.

**`statsmodels`** is also not installed. Its `statsmodels.stats.multitest`
module implements Bonferroni, Holm and false-discovery-rate corrections in
one call each; this lab's Bonferroni correction (`bonferroni_alpha`) is
three lines, described further, with `statsmodels`, in the lesson's Tools
section.

## If you cannot install anything at all

You still need NumPy for this lab: every exercise draws random samples or
shuffles at a scale (thousands of permutations, thousands of simulated
confidence intervals) that the standard library's `random` module was not
designed to batch efficiently. If NumPy genuinely cannot be installed, the
*ideas* -- a p-value is P(data this extreme | null true), a confidence
interval's 95% is a property of the procedure not of any one interval, a
permutation test needs no distributional assumption -- 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 `inference.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 two-sample z-test, from scratch

`phi(z)`, `p_from_z_two_sided(z)`, `two_sample_z_test(a, b)`. Build the
standard normal CDF from `math.erf`, turn a z-statistic into a two-sided
p-value, and combine them into a two-sample test that uses each sample's
own variance. Check the result against a hand computation done with the
standard library's `statistics` module -- they must match to the ninth
decimal place, because both are exact arithmetic on the same numbers.

## 2. The critical value, and the centrepiece: coverage

`z_critical_two_sided(alpha)`, `confidence_interval_mean(sample, alpha)`.
There is no closed form for the inverse of `phi`, so find the z whose
two-sided tail probability is `alpha` by bisecting `phi` itself. Then build
a confidence interval as `mean +/- z * standard_error`. Build 10,000 such
intervals from a population with a KNOWN true mean and assert the fraction
that actually contain it lands within three standard errors of 0.95 --
this is what "95% confidence" measures.

## 3. `ci_excludes` and the duality

`ci_excludes(interval, value)`, `one_sample_z_test_against_value(sample,
null_value)`. A test at level alpha should reject the null value exactly
when the `(1 - alpha)` interval excludes it. Assert this holds with ZERO
mismatches across many datasets -- not approximately, exactly, because both
come from the same z.

## 4. The permutation test, from scratch

`permutation_test_diff_means(a, b, n_perm, rng)`. Shuffle the pooled group
labels, recompute the difference in means, repeat `n_perm` times, and count
how many shuffles were at least as extreme as what was actually observed.
No distributional assumption anywhere. Assert it agrees closely with the
z-test on a moderate-n normal case, and diverges more (while remaining a
valid probability) on a small-n, skewed case.

## 5. Multiple comparisons and Bonferroni

`bonferroni_alpha(alpha, m)`. Assert the exact family-wise false-positive
rate for 20 independent alpha=0.05 tests is `1 - 0.95**20 = 0.6415`,
confirm it by simulation, then assert the Bonferroni-corrected rate lands
near `0.0488`.

## 6. Power

`power_two_sample_z(effect, sigma, n_per_group, alpha)`. Derive the power
of the test above from the shifted normal distribution of the test
statistic under a true effect. Assert it rises monotonically with both `n`
and `effect`, and check the closed form against a direct simulation.

## 7. Effect size versus n

Nothing new to write here -- exercise 7 in
`examples/07_effect_size_vs_n.py` runs against whatever `two_sample_z_test`
you wrote for exercise 1. A tiny, fixed relative difference should not be
significant at a small n and should be significant at an enormous one, with
the effect size itself unchanged.

## 8. Peeking

Nothing new to write here either -- exercise 8 in `examples/08_peeking.py`
runs against `one_sample_z_test_against_value` from exercise 3. Checking
after every 10 observations and stopping at the first p < 0.05 should push
the true false-positive rate well above the nominal alpha, under a null
hypothesis that is true the entire time.

## 9. The bootstrap interval, versus the normal-approximation interval

`bootstrap_ci(sample, statistic, n_boot, alpha, rng)`. Resample with
replacement, recompute the statistic, take percentiles of the result.
Assert it agrees closely with `confidence_interval_mean` for the sample
mean -- both center and width -- where both are valid.

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

Both `examples/` and `starter/` contain modules called `inference` 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 ("inference", "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]

"""Population generators and tolerance constants shared by every exercise.

This file is complete and does not need editing -- your work for each
exercise lives in `inference.py`. Every population here has a KNOWN true
parameter, which is what makes the coverage and duality exercises checkable
at all.
"""
from __future__ import annotations

import numpy as np

POP_MEAN = 50.3
POP_STD = 12.7

POP_B_MEAN = 53.1
POP_B_STD = 12.7

SKEWED_SCALE = 8.0


def normal_population(rng: np.random.Generator, n: int, mean: float, std: float) -> np.ndarray:
    return rng.normal(loc=mean, scale=std, size=n)


def skewed_population(rng: np.random.Generator, n: int, scale: float = SKEWED_SCALE) -> np.ndarray:
    return rng.exponential(scale=scale, size=n)


COVERAGE_TARGET = 0.95
COVERAGE_TRIALS = 10_000
COVERAGE_SAMPLE_N = 300
COVERAGE_SE = (COVERAGE_TARGET * (1 - COVERAGE_TARGET) / COVERAGE_TRIALS) ** 0.5
COVERAGE_TOLERANCE = 3 * COVERAGE_SE

FWER_TRIALS = 20
FWER_ALPHA = 0.05
FWER_EXACT = 1 - (1 - FWER_ALPHA) ** FWER_TRIALS
FWER_FAMILIES = 20_000
FWER_SIM_TOLERANCE = 0.015
BONFERRONI_EXPECTED = 1 - (1 - FWER_ALPHA / FWER_TRIALS) ** FWER_TRIALS
BONFERRONI_TOLERANCE = 0.015

PEEK_ALPHA = 0.05
PEEK_EXPERIMENTS = 4_000
PEEK_BATCH_SIZE = 10
PEEK_MAX_BATCHES = 5
PEEK_MIN_INFLATION_FACTOR = 2.0

BOOTSTRAP_N_BOOT = 5_000
BOOTSTRAP_WIDTH_RATIO_TOLERANCE = 0.20
BOOTSTRAP_CENTER_TOLERANCE_IN_SE = 0.6

POWER_CHECK_EFFECT = 2.8
POWER_CHECK_N = 100
POWER_CHECK_TRIALS = 3_000
POWER_CHECK_TOLERANCE = 0.03

EFFECT_VS_N_RELATIVE_DIFF = 0.005
EFFECT_VS_N_SMALL_N = 30
EFFECT_VS_N_LARGE_N = 100_000

"""Hypothesis-testing and confidence-interval machinery -- 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

import numpy as np


def phi(z: float) -> float:
    """Exercise 1. Standard normal CDF, computed from math.erf. No scipy.

    Hint: math.erf(z / math.sqrt(2)) is the piece you need; phi(z) is
    0.5 * (1 + that).
    """
    return None


def p_from_z_two_sided(z: float) -> float:
    """Exercise 1. Two-sided p-value for a standard normal statistic z."""
    return None


def z_critical_two_sided(alpha: float) -> float:
    """Exercise 2. The z such that P(|Z| > z) = alpha, found by bisecting
    phi() on the interval [0, 10] until phi(mid) is within reach of
    1 - alpha/2. 200 iterations is far more than enough precision.
    """
    return None


def two_sample_z_test(a: np.ndarray, b: np.ndarray) -> tuple[float, float]:
    """Exercise 1. Two-sample z-test using each sample's own variance
    (unpooled): z = (mean_a - mean_b) / sqrt(var_a/n_a + var_b/n_b).
    Return (z, p_two_sided).
    """
    return None


def confidence_interval_mean(sample: np.ndarray, alpha: float = 0.05) -> tuple[float, float]:
    """Exercise 2 / 3. Normal-approximation (1 - alpha) confidence interval
    for a mean: mean +/- z_(alpha/2) * (sample_std / sqrt(n)).
    """
    return None


def ci_excludes(interval: tuple[float, float], value: float) -> bool:
    """Exercise 3. True if value lies strictly outside the closed interval."""
    return None


def one_sample_z_test_against_value(sample: np.ndarray, null_value: float) -> tuple[float, float]:
    """Exercise 3 / 8. One-sample z-test of H0: population mean == null_value."""
    return None


def permutation_test_diff_means(
    a: np.ndarray, b: np.ndarray, n_perm: int, rng: np.random.Generator
) -> tuple[float, float]:
    """Exercise 4. Shuffle the pooled labels n_perm times, recompute the
    difference in means each time, and count how many shuffles are at
    least as extreme (in absolute value) as the real observed difference.
    p = (count_at_least_as_extreme + 1) / (n_perm + 1) -- the "+1"s avoid a
    p-value of exactly zero. Return (observed_diff, p_two_sided).
    """
    return None


def power_two_sample_z(effect: float, sigma: float, n_per_group: int, alpha: float = 0.05) -> float:
    """Exercise 6. Power of the two-sample z-test given a true mean
    difference `effect`, common known std `sigma`, and `n_per_group` per
    arm. SE = sigma * sqrt(2/n_per_group); under H1 the test statistic is
    Normal(effect/SE, 1); power is the probability that variable falls
    outside [-z_crit, z_crit].
    """
    return None


def bonferroni_alpha(alpha: float, m: int) -> float:
    """Exercise 5. The per-test alpha that keeps the family-wise error rate
    at `alpha` across `m` independent tests."""
    return None


def bootstrap_ci(
    sample: np.ndarray,
    statistic,
    n_boot: int,
    alpha: float,
    rng: np.random.Generator,
) -> tuple[float, float]:
    """Exercise 9. Percentile bootstrap interval: resample `sample` with
    replacement `n_boot` times, recompute `statistic` on each resample, and
    take the [alpha/2, 1 - alpha/2] percentiles of the results.
    """
    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 inference as I


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_MEAN > 0


# --------------------------------------------------------------------------
# Exercise 1 -- phi, p_from_z_two_sided, two_sample_z_test
# --------------------------------------------------------------------------


def test_1_phi_of_zero_is_one_half():
    result = attempt(lambda: I.phi(0.0), "phi")
    assert result == pytest.approx(0.5), f"phi(0.0) should be 0.5, got {result}"


def test_1_p_from_z_two_sided_at_1_96_is_near_0_05():
    result = attempt(lambda: I.p_from_z_two_sided(1.959964), "p_from_z_two_sided")
    assert result == pytest.approx(0.05, abs=1e-4), f"expected ~0.05, got {result}"


def test_1_two_sample_z_test_matches_hand_computation():
    a = [50, 52, 49, 51, 53, 48, 50, 52, 51, 49]
    b = [54, 55, 53, 56, 54, 52, 55, 53, 54, 56]
    z, p = attempt(lambda: I.two_sample_z_test(a, b), "two_sample_z_test")
    import statistics

    mean_a, var_a = statistics.mean(a), statistics.variance(a)
    mean_b, var_b = statistics.mean(b), statistics.variance(b)
    se = math.sqrt(var_a / len(a) + var_b / len(b))
    z_hand = (mean_a - mean_b) / se
    assert z == pytest.approx(z_hand, abs=1e-6), f"expected z={z_hand}, got {z}"
    assert p < 0.001, f"these samples are clearly different -- expected p < 0.001, got {p}"


# --------------------------------------------------------------------------
# Exercise 2 -- z_critical_two_sided, confidence_interval_mean, coverage
# --------------------------------------------------------------------------


def test_2_z_critical_two_sided_matches_known_constant():
    result = attempt(lambda: I.z_critical_two_sided(0.05), "z_critical_two_sided")
    assert result == pytest.approx(1.959964, abs=1e-4), f"expected ~1.959964, got {result}"


def test_2_confidence_interval_mean_is_centered_on_the_sample_mean():
    rng = np.random.default_rng(1)
    sample = D.normal_population(rng, 100, D.POP_MEAN, D.POP_STD)
    lo, hi = attempt(lambda: I.confidence_interval_mean(sample, 0.05), "confidence_interval_mean")
    center = (lo + hi) / 2
    assert center == pytest.approx(sample.mean(), abs=1e-6), f"expected center {sample.mean()}, got {center}"


def test_2_coverage_is_close_to_nominal():
    rng = np.random.default_rng(42)

    def run():
        hits = 0
        trials = 2000
        for _ in range(trials):
            sample = D.normal_population(rng, D.COVERAGE_SAMPLE_N, D.POP_MEAN, D.POP_STD)
            lo, hi = I.confidence_interval_mean(sample, alpha=0.05)
            if lo <= D.POP_MEAN <= hi:
                hits += 1
        return hits / trials

    coverage = attempt(run, "confidence_interval_mean (coverage)")
    se = math.sqrt(0.95 * 0.05 / 2000)
    assert abs(coverage - 0.95) <= 3 * se, f"measured coverage {coverage} too far from 0.95"


# --------------------------------------------------------------------------
# Exercise 3 -- ci_excludes, one_sample_z_test_against_value, duality
# --------------------------------------------------------------------------


def test_3_ci_excludes():
    result = attempt(lambda: I.ci_excludes((1.0, 2.0), 3.0), "ci_excludes")
    assert result is True, f"3.0 is outside (1.0, 2.0) -- expected True, got {result}"
    result2 = attempt(lambda: I.ci_excludes((1.0, 2.0), 1.5), "ci_excludes")
    assert result2 is False, f"1.5 is inside (1.0, 2.0) -- expected False, got {result2}"


def test_3_duality_holds_exactly():
    rng = np.random.default_rng(9)

    def run():
        mismatches = 0
        for _ in range(200):
            n = int(rng.integers(15, 60))
            shift = 0.0 if rng.random() < 0.5 else rng.uniform(-6.0, 6.0)
            sample = D.normal_population(rng, n, D.POP_MEAN + shift, D.POP_STD)
            _, p = I.one_sample_z_test_against_value(sample, D.POP_MEAN)
            interval = I.confidence_interval_mean(sample, 0.05)
            if (p < 0.05) != I.ci_excludes(interval, D.POP_MEAN):
                mismatches += 1
        return mismatches

    mismatches = attempt(run, "one_sample_z_test_against_value / duality")
    assert mismatches == 0, f"expected zero test/interval disagreements, got {mismatches}"


# --------------------------------------------------------------------------
# Exercise 4 -- permutation_test_diff_means
# --------------------------------------------------------------------------


def test_4_permutation_test_returns_valid_probability():
    rng = np.random.default_rng(3)
    a = D.normal_population(rng, 30, D.POP_MEAN, D.POP_STD)
    b = D.normal_population(rng, 30, D.POP_B_MEAN, D.POP_B_STD)
    _, p = attempt(lambda: I.permutation_test_diff_means(a, b, 500, rng), "permutation_test_diff_means")
    assert 0.0 <= p <= 1.0, f"a p-value must be in [0, 1], got {p}"


def test_4_permutation_test_p_is_small_for_a_large_true_difference():
    rng = np.random.default_rng(4)
    a = D.normal_population(rng, 50, 0.0, 1.0)
    b = D.normal_population(rng, 50, 20.0, 1.0)
    _, p = attempt(lambda: I.permutation_test_diff_means(a, b, 500, rng), "permutation_test_diff_means")
    assert p < 0.01, f"a 20-sigma separation should be obviously significant, got p={p}"


# --------------------------------------------------------------------------
# Exercise 5 -- bonferroni_alpha
# --------------------------------------------------------------------------


def test_5_bonferroni_alpha_divides_by_m():
    result = attempt(lambda: I.bonferroni_alpha(0.05, 20), "bonferroni_alpha")
    assert result == pytest.approx(0.0025), f"expected 0.05/20=0.0025, got {result}"


# --------------------------------------------------------------------------
# Exercise 6 -- power_two_sample_z
# --------------------------------------------------------------------------


def test_6_power_increases_with_n():
    p1 = attempt(lambda: I.power_two_sample_z(3.0, 12.7, 20), "power_two_sample_z")
    p2 = attempt(lambda: I.power_two_sample_z(3.0, 12.7, 200), "power_two_sample_z")
    assert p2 > p1, f"power should rise with n, got power(20)={p1}, power(200)={p2}"


def test_6_power_at_zero_effect_equals_alpha():
    result = attempt(lambda: I.power_two_sample_z(0.0, 12.7, 100, alpha=0.05), "power_two_sample_z")
    assert result == pytest.approx(0.05, abs=1e-4), f"with no true effect, power should equal alpha, got {result}"


# --------------------------------------------------------------------------
# Exercise 9 -- bootstrap_ci
# --------------------------------------------------------------------------


def test_9_bootstrap_ci_contains_the_sample_mean():
    rng = np.random.default_rng(5)
    sample = D.normal_population(rng, 200, D.POP_MEAN, D.POP_STD)
    lo, hi = attempt(lambda: I.bootstrap_ci(sample, np.mean, 1000, 0.05, rng), "bootstrap_ci")
    assert lo < sample.mean() < hi, f"the sample mean {sample.mean()} should fall inside [{lo}, {hi}]"


def test_9_bootstrap_ci_agrees_with_normal_ci():
    rng = np.random.default_rng(6)
    sample = D.normal_population(rng, 200, D.POP_MEAN, D.POP_STD)
    normal_lo, normal_hi = attempt(
        lambda: I.confidence_interval_mean(sample, 0.05), "confidence_interval_mean"
    )
    boot_lo, boot_hi = attempt(lambda: I.bootstrap_ci(sample, np.mean, 3000, 0.05, rng), "bootstrap_ci")
    width_ratio = (boot_hi - boot_lo) / (normal_hi - normal_lo)
    assert 0.7 < width_ratio < 1.3, f"bootstrap and normal interval widths should roughly agree, ratio={width_ratio}"

#!/usr/bin/env bash
# Tests for the Day 118 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 two-sample z-test built from math.erf matches an independent hand
#     computation to nine decimal places;
#   * the centrepiece -- 10,000 nominal-95% confidence intervals actually
#     cover the known true mean about 95% of the time;
#   * the test/interval duality holds EXACTLY, with zero mismatches, across
#     hundreds of datasets;
#   * a from-scratch permutation test agrees closely with the z-test where
#     the normal approximation holds, and diverges more where it does not;
#   * twenty independent alpha=0.05 tests give a 64.15% chance of at least
#     one false positive -- exact arithmetic, confirmed by simulation -- and
#     Bonferroni pulls that back to about 4.9%;
#   * power rises with both n and effect size, and the closed-form formula
#     agrees with a direct simulation of the test itself;
#   * the same tiny relative difference is not significant at a small n and
#     is significant at an enormous one, with the effect size unchanged;
#   * checking after every 10 observations and stopping at the first
#     p < 0.05 inflates the true false-positive rate several times past
#     alpha, under a null hypothesis that never changed;
#   * a bootstrap interval agrees closely with the normal-approximation
#     interval where both are valid;
#   * 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 118 — Hypothesis Tests and Confidence Intervals"
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_two_sample_z_test 02_coverage 03_duality 04_permutation_test \
              05_multiple_comparisons 06_power 07_effect_size_vs_n 08_peeking \
              09_bootstrap_vs_normal_ci; 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
    *"OK:"*)
      check "${script}.py reports OK" "yes" ;;
    *) check "${script}.py reports OK" "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 18 ]; then
  check "the reference suite ran at least 18 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 18 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 `inference` and
# `dataset`, and pytest imports test files by putting their directory on
# sys.path -- so collecting both suites at once could 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: running pytest with NO explicit path (auto-discovering both
# directories from the lab root) must report the same skip count as running
# `pytest starter` alone.
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 "auto-discovering both suites 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 05 (multiple comparisons) with its expected exact
# family-wise error rate deliberately swapped for a wrong one, and asserts
# that the re-run reports the failure and exits non-zero.
if [ -z "${D118_SELF_TEST:-}" ]; then
  self_out="$(cd "${lab_dir}/examples" && D118_SELF_TEST=1 "${python_bin}" -c "
src = open('05_multiple_comparisons.py').read()
src = src.replace('abs(exact - 0.6415) < 0.0001', 'abs(exact - 99.0) < 0.0001')
g = {'__name__': '__main__', '__file__': '05_multiple_comparisons.py'}
exec(compile(src, '05_multiple_comparisons.py', 'exec'), g)
" 2>&1)"
  self_status=$?
  if [ "${self_status}" -ne 0 ]; then
    check "a deliberately wrong expectation makes script 05 exit non-zero (${self_status})" "yes"
  else
    check "a deliberately wrong expectation makes script 05 exit non-zero" "no"
  fi
  case "${self_out}" in
    *"AssertionError"*"0.6415"*)
      check "the failing assertion is named in the output with the real value" "yes" ;;
    *) check "the failing assertion is named in the output with the real value" "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 ]