Hands-on lab — Day 116: Descriptive Statistics That Don’t Lie

"""Exercise 1: mean, median and mode, from scratch, checked against the
`statistics` module on the same inputs."""

import statistics as st

import dataset as D
import descriptive as F


def main() -> None:
    m = F.mean(D.ODD_LIST)
    ref_mean = st.fmean(D.ODD_LIST)
    print(f"odd list         : {D.ODD_LIST}")
    print(f"  mean (ours)     = {m}")
    print(f"  mean (stdlib)   = {ref_mean}")
    assert m == ref_mean

    med = F.median(D.ODD_LIST)
    ref_med = st.median(D.ODD_LIST)
    print(f"  median (ours)   = {med}")
    print(f"  median (stdlib) = {ref_med}")
    assert med == ref_med

    mo = F.modes(D.ODD_LIST)
    print(f"  modes           = {mo}")
    assert mo == [7]

    print(f"even list         : {D.EVEN_LIST}")
    med_even = F.median(D.EVEN_LIST)
    ref_med_even = st.median(D.EVEN_LIST)
    print(f"  median (ours)   = {med_even}")
    print(f"  median (stdlib) = {ref_med_even}")
    assert med_even == ref_med_even == 6.0  # average of 4.0 and 8.0

    print(f"multimodal list   : {D.MULTIMODAL_LIST}")
    mo_multi = F.modes(D.MULTIMODAL_LIST)
    ref_multi = sorted(st.multimode(D.MULTIMODAL_LIST))
    print(f"  modes (ours)    = {mo_multi}")
    print(f"  modes (stdlib)  = {ref_multi}")
    assert mo_multi == ref_multi == [3, 8]
    # statistics.mode() (singular) does NOT raise on ties -- it silently
    # returns whichever tied value it saw first, which is exactly the
    # "average income" trap in miniature: a single-number summary hiding
    # that there were two equally valid answers.
    single = st.mode(D.MULTIMODAL_LIST)
    print(f"  statistics.mode() picks just one: {single} (silently drops {mo_multi[1] if single == mo_multi[0] else mo_multi[0]})")
    assert single in mo_multi

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


if __name__ == "__main__":
    main()

"""Exercise 2: the breakdown point. Corrupt exactly one salary out of nine
and watch what happens to the mean versus the median."""

import dataset as D
import descriptive as F


def main() -> None:
    print(f"salaries           : {D.SALARY_LIST}")
    print(f"corrupted value    : {D.CORRUPTED_SALARY:,}")

    mean_before, mean_after = F.breakdown_point_mean(D.SALARY_LIST, D.CORRUPTED_SALARY)
    print(f"mean before        = {mean_before:,.2f}")
    print(f"mean after         = {mean_after:,.2f}")
    mean_shift = mean_after - mean_before
    print(f"mean moved by      = {mean_shift:,.2f}")
    assert mean_shift > D.BREAKDOWN_MEAN_SHIFT_FLOOR

    median_before, median_after = F.breakdown_point_median(D.SALARY_LIST, D.CORRUPTED_SALARY)
    print(f"median before      = {median_before:,.2f}")
    print(f"median after       = {median_after:,.2f}")
    median_shift = median_after - median_before
    print(f"median moved by    = {median_shift:,.2f}")
    # Exact equality is the right assertion: one corrupted value out of nine
    # cannot move the median AT ALL, because the median only cares about the
    # RANK of the middle value, and the corrupted value (however extreme)
    # is still the single largest value, occupying the same rank position
    # (9th of 9) whether it is $60,000 or $10,000,000.
    assert median_shift == 0.0

    print(
        "One value out of nine dragged the mean by "
        f"${mean_shift:,.0f} and the median by exactly $0."
    )
    print("02_breakdown_point.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 3: Bessel's correction, measured rather than asserted.

Draw many samples of size 5 from a population of known variance. The
divide-by-n estimator should be biased LOW by the factor (n-1)/n on
average; the divide-by-(n-1) estimator should be unbiased, within a few
standard errors of its own sampling mean.
"""

import numpy as np

import dataset as D
import simulate as S


def main() -> None:
    rng = np.random.default_rng(D.BESSEL_SEED)
    biased, unbiased = S.bessel_trial_variances(
        rng,
        D.BESSEL_POPULATION_MEAN,
        D.BESSEL_POPULATION_SIGMA,
        D.BESSEL_SAMPLE_SIZE,
        D.BESSEL_TRIALS,
    )

    print(f"true population variance        = {D.BESSEL_TRUE_VARIANCE}")
    print(f"sample size (n)                 = {D.BESSEL_SAMPLE_SIZE}")
    print(f"trials                          = {D.BESSEL_TRIALS:,}")

    mean_biased = float(biased.mean())
    ratio_biased = mean_biased / D.BESSEL_TRUE_VARIANCE
    print(f"mean of divide-by-n estimator   = {mean_biased:.4f}")
    print(f"  ratio to true variance        = {ratio_biased:.4f}")
    print(f"  predicted ratio (n-1)/n       = {D.BESSEL_EXPECTED_BIAS_FACTOR:.4f}")
    assert abs(ratio_biased - D.BESSEL_EXPECTED_BIAS_FACTOR) < D.BESSEL_BIAS_FACTOR_TOLERANCE

    mean_unbiased = float(unbiased.mean())
    se_unbiased = float(unbiased.std(ddof=1)) / (D.BESSEL_TRIALS**0.5)
    deviations_in_se = abs(mean_unbiased - D.BESSEL_TRUE_VARIANCE) / se_unbiased
    print(f"mean of divide-by-(n-1) estimator = {mean_unbiased:.4f}")
    print(f"  standard error of that mean     = {se_unbiased:.4f}")
    print(f"  distance from truth, in SEs     = {deviations_in_se:.2f}")
    assert deviations_in_se < D.BESSEL_UNBIASED_SE_TOLERANCE

    print(
        f"Dividing by n underestimates the true variance by a factor of "
        f"~{1 - ratio_biased:.1%}. Dividing by n-1 lands within "
        f"{deviations_in_se:.2f} standard errors of the truth."
    )
    print("03_bessel_correction.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 4: the 75th percentile of one small array, under several of
NumPy's `method=` conventions. The assertion is that "the percentile" is
not a well-defined number: at least two conventions disagree."""

import dataset as D
import descriptive as F


def main() -> None:
    print(f"array   : {D.PERCENTILE_ARRAY}")
    print(f"target  : {D.PERCENTILE_TARGET}th percentile")

    results: dict[str, float] = {}
    for method in D.PERCENTILE_METHODS:
        value = F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, method)
        results[method] = value
        print(f"  method={method:<16} -> {value}")

    distinct = sorted(set(results.values()))
    print(f"distinct values across {len(D.PERCENTILE_METHODS)} conventions: {distinct}")
    assert len(distinct) >= 2, "expected at least two conventions to disagree"

    # The default ('linear') is the convention pandas' DataFrame.describe()
    # also uses -- worth naming explicitly, since it is the one most people
    # get without ever choosing it.
    print(f"default ('linear') result: {results['linear']}")
    assert results["linear"] == 8.25

    # 'lower' and 'higher' are not just close -- they land on two different
    # ACTUAL DATA POINTS, one full step apart.
    assert results["lower"] != results["higher"]
    print(f"'lower' picks an actual data point: {results['lower']}")
    print(f"'higher' picks a DIFFERENT actual data point: {results['higher']}")

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


if __name__ == "__main__":
    main()

"""Exercise 5: Pearson measures LINEAR association only; Spearman measures
MONOTONE association. A perfect parabola fools Pearson; a perfect (but
non-linear) monotone curve does not fool Spearman."""

import dataset as D
import descriptive as F


def main() -> None:
    print(f"parabola x = {D.PARABOLA_X}")
    print(f"parabola y = x^2 = {D.PARABOLA_Y}")
    pear_parabola = F.pearson(D.PARABOLA_X, D.PARABOLA_Y)
    print(f"  Pearson correlation  = {pear_parabola}")
    assert abs(pear_parabola) < D.PARABOLA_PEARSON_TOLERANCE
    print("  -> essentially zero, despite y being EXACTLY determined by x")

    print()
    print(f"monotone x = {D.MONOTONE_X}")
    print(f"monotone y = x^3 = {D.MONOTONE_Y}")
    pear_monotone = F.pearson(D.MONOTONE_X, D.MONOTONE_Y)
    spear_monotone = F.spearman(D.MONOTONE_X, D.MONOTONE_Y)
    print(f"  Pearson correlation  = {pear_monotone:.6f}  (strong, but not perfect)")
    print(f"  Spearman correlation = {spear_monotone:.6f}  (exactly 1.0)")
    assert spear_monotone == 1.0
    assert pear_monotone < 1.0  # the cubic is not a straight line

    print()
    print(
        "Same shape of evidence, two different questions: Pearson asks "
        "'how well does a straight line fit', Spearman asks 'does y always "
        "increase when x does'. The parabola answers the first question "
        "with 'not at all' and would answer the second with 'no' too, "
        "since it isn't monotone. The cubic answers the first with "
        "'pretty well' and the second with 'perfectly'."
    )
    print("05_pearson_vs_spearman.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 6: Anscombe's quartet. Four datasets that agree on every
familiar summary statistic to the documented precision -- and three shape
diagnostics that separate them completely, each for a different structural
reason.

Anscombe, F. J. (1973). "Graphs in Statistical Analysis."
The American Statistician, 27(1), 17-21.
"""

import dataset as D
import descriptive as F


def main() -> None:
    summaries = {name: F.anscombe_summary(x, y) for name, (x, y) in D.ANSCOMBE_SETS.items()}

    print(f"{'set':<5}{'mean x':>9}{'mean y':>9}{'var x':>9}{'var y':>9}{'r':>9}{'slope':>9}")
    for name, s in summaries.items():
        print(
            f"{name:<5}{s['mean_x']:>9.2f}{s['mean_y']:>9.2f}{s['var_x']:>9.2f}"
            f"{s['var_y']:>9.2f}{s['correlation']:>9.4f}{s['slope']:>9.4f}"
        )

    # All four sets agree to the documented precision on every classic
    # summary statistic.
    dec = D.ANSCOMBE_AGREEMENT_DECIMALS
    reference = summaries["I"]
    for name, s in summaries.items():
        assert round(s["mean_x"], dec) == round(reference["mean_x"], dec), name
        assert round(s["mean_y"], dec) == round(reference["mean_y"], dec), name
        assert round(s["var_x"], dec) == round(reference["var_x"], dec), name
        assert round(s["var_y"], dec) == round(reference["var_y"], dec), name
        assert round(s["correlation"], 1) == round(reference["correlation"], 1), name
        assert round(s["slope"], 1) == round(reference["slope"], 1), name

    print()
    print("Every classic summary statistic agrees. Three diagnostics those")
    print("summaries cannot see now tell the four sets apart:")
    print()

    shapes = {name: F.shape_statistics(x, y) for name, (x, y) in D.ANSCOMBE_SETS.items()}
    print(f"{'set':<5}{'max leverage':>14}{'outlier ratio':>15}{'sign changes':>14}")
    for name, s in shapes.items():
        print(
            f"{name:<5}{s['max_leverage']:>14.3f}{s['outlier_ratio']:>15.3f}"
            f"{s['residual_sign_changes']:>14.0f}"
        )

    # Set IV: the x-values are identical except for one, so ONE point
    # controls the entire slope -- its leverage swamps everyone else's,
    # while sets I, II and III (sharing the same x column) have identical,
    # unremarkable leverage.
    assert shapes["IV"]["max_leverage"] > 3.0 * shapes["I"]["max_leverage"]
    assert shapes["I"]["max_leverage"] == shapes["II"]["max_leverage"] == shapes["III"]["max_leverage"]

    # Set III: a perfect line plus one outlier -- that outlier's residual
    # dwarfs the combined residuals of every other point.
    assert shapes["III"]["outlier_ratio"] > 2.0 * shapes["I"]["outlier_ratio"]

    # Set II: a perfect parabola -- a straight line fit to it produces a
    # smooth, systematic curve of residuals that changes sign rarely,
    # unlike set I's honestly scattered noise, which flips sign often.
    assert shapes["II"]["residual_sign_changes"] < shapes["I"]["residual_sign_changes"]

    print()
    print(
        "Set IV's one non-repeated x-value carries most of the leverage. "
        "Set III's one outlier carries most of the residual. Set II's "
        "residuals trace a smooth curve instead of scattering. Set I is "
        "unremarkable on all three -- which is exactly what a genuinely "
        "linear relationship with honest noise should look like."
    )
    print("06_anscombes_quartet.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 7: Simpson's paradox. Treatment A beats treatment B in EVERY
subgroup, and treatment B beats treatment A overall -- from the same
table, both directions verified by direct arithmetic."""

import dataset as D
import descriptive as F


def main() -> None:
    a_easy_rate = F.success_rate(*D.TREATMENT_A_EASY)
    a_hard_rate = F.success_rate(*D.TREATMENT_A_HARD)
    b_easy_rate = F.success_rate(*D.TREATMENT_B_EASY)
    b_hard_rate = F.success_rate(*D.TREATMENT_B_HARD)

    print("                 easy subgroup      hard subgroup")
    print(
        f"treatment A   {D.TREATMENT_A_EASY[0]}/{D.TREATMENT_A_EASY[1]} = {a_easy_rate:.1%}"
        f"         {D.TREATMENT_A_HARD[0]}/{D.TREATMENT_A_HARD[1]} = {a_hard_rate:.1%}"
    )
    print(
        f"treatment B   {D.TREATMENT_B_EASY[0]}/{D.TREATMENT_B_EASY[1]} = {b_easy_rate:.1%}"
        f"        {D.TREATMENT_B_HARD[0]}/{D.TREATMENT_B_HARD[1]} = {b_hard_rate:.1%}"
    )

    # Direction one: A wins BOTH subgroups.
    assert a_easy_rate > b_easy_rate
    assert a_hard_rate > b_hard_rate
    print()
    print(f"A beats B in the easy subgroup:  {a_easy_rate:.1%} > {b_easy_rate:.1%}")
    print(f"A beats B in the hard subgroup:  {a_hard_rate:.1%} > {b_hard_rate:.1%}")

    a_total = F.combined_rate(D.TREATMENT_A_EASY, D.TREATMENT_A_HARD)
    b_total = F.combined_rate(D.TREATMENT_B_EASY, D.TREATMENT_B_HARD)
    a_n = D.TREATMENT_A_EASY[1] + D.TREATMENT_A_HARD[1]
    b_n = D.TREATMENT_B_EASY[1] + D.TREATMENT_B_HARD[1]
    print()
    print(f"treatment A overall: {a_total:.1%}  (out of {a_n} trials)")
    print(f"treatment B overall: {b_total:.1%}  (out of {b_n} trials)")

    # Direction two: B wins OVERALL. Both directions, same table.
    assert b_total > a_total

    print()
    print(
        "A won every subgroup and lost overall. The mechanism: A took "
        f"{D.TREATMENT_A_HARD[1]} of its {a_n} trials in the hard subgroup "
        f"(where success is rare for everyone), while B took only "
        f"{D.TREATMENT_B_HARD[1]} of its {b_n} trials there. The overall "
        "rate is a WEIGHTED average, and the weights -- not the treatment "
        "-- are what flipped the ranking."
    )
    print("07_simpsons_paradox.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 8: robust spread under contamination. A few percent of extreme
values inflate the standard deviation dramatically and barely move the
median absolute deviation."""

import numpy as np

import dataset as D
import descriptive as F
import simulate as S


def main() -> None:
    rng = np.random.default_rng(D.CONTAMINATION_SEED)
    clean, contaminated = S.contaminated_sample(
        rng,
        D.CONTAMINATION_BASE_MEAN,
        D.CONTAMINATION_BASE_SIGMA,
        D.CONTAMINATION_BASE_N,
        D.CONTAMINATION_OUTLIERS,
    )
    contamination_fraction = len(D.CONTAMINATION_OUTLIERS) / len(contaminated)
    print(f"clean sample size        = {len(clean)}")
    print(f"outliers added           = {list(D.CONTAMINATION_OUTLIERS)}")
    print(f"contamination fraction   = {contamination_fraction:.1%}")

    std_clean = float(np.std(clean, ddof=1))
    std_contam = float(np.std(contaminated, ddof=1))
    std_multiplier = std_contam / std_clean
    print(f"standard deviation, clean        = {std_clean:.3f}")
    print(f"standard deviation, contaminated = {std_contam:.3f}")
    print(f"  multiplier                     = {std_multiplier:.2f}x")
    assert std_multiplier > D.CONTAMINATION_STD_MULTIPLIER_FLOOR

    mad_clean = F.median_absolute_deviation(clean)
    mad_contam = F.median_absolute_deviation(contaminated)
    mad_multiplier = mad_contam / mad_clean
    print(f"median absolute deviation, clean        = {mad_clean:.3f}")
    print(f"median absolute deviation, contaminated = {mad_contam:.3f}")
    print(f"  multiplier                            = {mad_multiplier:.2f}x")
    assert mad_multiplier < D.CONTAMINATION_MAD_MULTIPLIER_CEILING

    print()
    print(
        f"{contamination_fraction:.0%} contamination inflated the standard "
        f"deviation by {std_multiplier:.1f}x and the MAD by only "
        f"{mad_multiplier:.2f}x."
    )
    print("08_robust_spread_under_contamination.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 9: standardisation and z-scores. Standardising gives mean 0 and
standard deviation 1 by construction, and it does NOT change the Pearson
correlation between two variables."""

import numpy as np

import dataset as D
import descriptive as F


def main() -> None:
    rng = np.random.default_rng(D.STANDARDIZATION_SEED)
    x = rng.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
    noise = rng.normal(0.0, D.STANDARDIZATION_Y_NOISE_SIGMA, D.STANDARDIZATION_N)
    y = D.STANDARDIZATION_Y_SLOPE * x + noise

    zx = F.zscores(x)
    zy = F.zscores(y)

    mean_zx = sum(zx) / len(zx)
    std_zx = (sum((v - mean_zx) ** 2 for v in zx) / len(zx)) ** 0.5
    print(f"standardised x: mean = {mean_zx:.2e}, std = {std_zx:.10f}")
    assert abs(mean_zx) < D.STANDARDIZATION_MEAN_TOLERANCE
    assert abs(std_zx - 1.0) < D.STANDARDIZATION_STD_TOLERANCE

    r_original = F.pearson(x, y)
    r_standardised = F.pearson(zx, zy)
    diff = abs(r_original - r_standardised)
    print(f"Pearson correlation, original x,y      = {r_original:.10f}")
    print(f"Pearson correlation, standardised x,y  = {r_standardised:.10f}")
    print(f"  difference                           = {diff:.2e}")
    assert diff < D.STANDARDIZATION_CORRELATION_TOLERANCE

    print()
    print(
        "Standardising rescaled every value's units, but it moved nothing "
        "relative to anything else -- correlation, which is exactly a "
        "statement about relative structure, could not have changed."
    )
    print("09_standardization.py: every assertion held.")


if __name__ == "__main__":
    main()

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

Both `examples/` and `starter/` contain modules called `dataset`,
`descriptive`, `simulate` and `answers`, 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 `descriptive` was
seen first and then reuse it for the other suite -- 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 ("dataset", "descriptive", "simulate", "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 dataset and tolerance this lab compares against.

Read this file. Nothing here is tuned to make a test pass: the salary list
is small and hand-checkable, Anscombe's quartet is the published 1973
dataset (Anscombe, "Graphs in Statistical Analysis", The American
Statistician, 1973), and the Bessel-correction tolerance is derived from the
standard error of the simulated mean, with the arithmetic written out beside
it, never chosen by running a test and loosening the number until it
passed.
"""

import math
from fractions import Fraction

# --------------------------------------------------------------------------
# Exercise 1: mean, median, mode
# --------------------------------------------------------------------------

#: An odd-length list with a single clear mode (7 appears three times).
ODD_LIST: tuple[int, ...] = (2, 4, 4, 7, 7, 7, 9, 12, 15)

#: An even-length list, so the median must average the two middle values.
EVEN_LIST: tuple[float, ...] = (1.0, 3.0, 4.0, 8.0, 10.0, 12.0)

#: A multimodal list: 3 and 8 are each the most frequent value.
MULTIMODAL_LIST: tuple[int, ...] = (3, 3, 5, 6, 8, 8, 9)

# --------------------------------------------------------------------------
# Exercise 2: the breakdown point
# --------------------------------------------------------------------------

#: Nine ordinary salaries, in dollars. Small enough to check the median of
#: by eye: sorted, the middle (5th of 9) value is 50000.
SALARY_LIST: tuple[int, ...] = (
    42000,
    45000,
    47000,
    48000,
    50000,
    52000,
    55000,
    58000,
    60000,
)

#: Replace the single largest salary with one wildly corrupted value. One
#: value out of nine changing is "one value out of many" -- nowhere near the
#: 50% the median can absorb before it moves at all.
CORRUPTED_SALARY: int = 10_000_000

#: The mean must move by at least this many dollars for the demonstration to
#: count as "dragged anywhere at all" -- chosen as a small fraction of the
#: actual shift (over $1,000,000), so the assertion is not a coin flip.
BREAKDOWN_MEAN_SHIFT_FLOOR: float = 500_000.0

# --------------------------------------------------------------------------
# Exercise 3: Bessel's correction, measured
# --------------------------------------------------------------------------

#: The population this exercise draws from: mean 50, standard deviation 10,
#: so the true variance is exactly 100. Both are known exactly because the
#: population is synthetic, which is what makes the bias measurable at all
#: -- there is a ground truth to compare the estimators against.
BESSEL_POPULATION_MEAN: float = 50.0
BESSEL_POPULATION_SIGMA: float = 10.0
BESSEL_TRUE_VARIANCE: float = BESSEL_POPULATION_SIGMA**2

#: Small samples, where the fitted-mean bias is largest relative to n.
BESSEL_SAMPLE_SIZE: int = 5

#: Many repeated samples, so the estimators' average is itself a precise
#: number rather than one noisy draw.
BESSEL_TRIALS: int = 20_000

BESSEL_SEED: int = 116

#: The textbook claim this exercise measures: dividing by n underestimates
#: the true variance, on average, by exactly the factor (n-1)/n.
BESSEL_EXPECTED_BIAS_FACTOR: float = (BESSEL_SAMPLE_SIZE - 1) / BESSEL_SAMPLE_SIZE

#: How close the measured biased-estimator ratio must land to (n-1)/n. This
#: is a Monte Carlo measurement, not exact arithmetic, so the tolerance is
#: generous but not vacuous -- wide enough to pass on any honest run of
#: BESSEL_TRIALS trials, narrow enough that a genuinely broken implementation
#: (dividing by the wrong thing, or not at all) fails it.
BESSEL_BIAS_FACTOR_TOLERANCE: float = 0.02

#: How many standard errors of the *unbiased* estimator's own sampling mean
#: it is allowed to sit from the true variance. The unbiased estimator is
#: correct only on average, not on every run, so this must be generous
#: enough to pass reliably -- 4 standard errors covers essentially all honest
#: runs (about 1 in 16,000 fails purely by chance under a normal
#: approximation).
BESSEL_UNBIASED_SE_TOLERANCE: float = 4.0

# --------------------------------------------------------------------------
# Exercise 4: percentile ambiguity
# --------------------------------------------------------------------------

#: A small, deliberately awkward array: 8 values, so the 75th percentile
#: falls between two of them under every interpolation convention, and the
#: conventions genuinely disagree about where.
PERCENTILE_ARRAY: tuple[int, ...] = (1, 2, 3, 4, 6, 8, 9, 15)

PERCENTILE_TARGET: float = 75.0

#: A representative slice of NumPy's nine documented `method=` conventions
#: for `numpy.percentile` (NumPy >= 1.22). `linear` is the default and the
#: one pandas' `DataFrame.describe()` also uses.
PERCENTILE_METHODS: tuple[str, ...] = (
    "linear",
    "lower",
    "higher",
    "nearest",
    "midpoint",
    "weibull",
    "median_unbiased",
    "normal_unbiased",
    "hazen",
)

# --------------------------------------------------------------------------
# Exercise 5: Pearson versus Spearman
# --------------------------------------------------------------------------

#: A perfect, symmetric parabola: y = x^2 over a symmetric range of x. Every
#: unit increase in x on the left half is mirrored by an equal decrease in y
#: on the right half's counterpart, so the *linear* trend that Pearson
#: measures cancels out exactly, even though y is perfectly determined by x.
PARABOLA_X: tuple[int, ...] = tuple(range(-5, 6))
PARABOLA_Y: tuple[int, ...] = tuple(v**2 for v in PARABOLA_X)

#: The tolerance "essentially zero" is checked against -- not equality to 0,
#: because that would be asserting a specific floating-point outcome rather
#: than the mathematical fact.
PARABOLA_PEARSON_TOLERANCE: float = 1e-9

#: A monotone but non-linear relationship: y = x^3. Every increase in x
#: produces an increase in y, so the *rank order* is perfectly preserved,
#: even though the relationship is not a straight line.
MONOTONE_X: tuple[int, ...] = tuple(range(-5, 6))
MONOTONE_Y: tuple[int, ...] = tuple(v**3 for v in MONOTONE_X)

# --------------------------------------------------------------------------
# Exercise 6: Anscombe's quartet
# --------------------------------------------------------------------------
# The published 1973 dataset, reproduced exactly. Anscombe, F. J. (1973).
# "Graphs in Statistical Analysis." The American Statistician, 27(1), 17-21.

ANSCOMBE_X_I: tuple[float, ...] = (10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0)
ANSCOMBE_Y_I: tuple[float, ...] = (
    8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68,
)

ANSCOMBE_X_II: tuple[float, ...] = ANSCOMBE_X_I
ANSCOMBE_Y_II: tuple[float, ...] = (
    9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74,
)

ANSCOMBE_X_III: tuple[float, ...] = ANSCOMBE_X_I
ANSCOMBE_Y_III: tuple[float, ...] = (
    7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73,
)

ANSCOMBE_X_IV: tuple[float, ...] = (8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 19.0, 8.0, 8.0, 8.0)
ANSCOMBE_Y_IV: tuple[float, ...] = (
    6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89,
)

ANSCOMBE_SETS: dict[str, tuple[tuple[float, ...], tuple[float, ...]]] = {
    "I": (ANSCOMBE_X_I, ANSCOMBE_Y_I),
    "II": (ANSCOMBE_X_II, ANSCOMBE_Y_II),
    "III": (ANSCOMBE_X_III, ANSCOMBE_Y_III),
    "IV": (ANSCOMBE_X_IV, ANSCOMBE_Y_IV),
}

#: How closely the four sets must agree on mean x, mean y, variance x,
#: variance y, correlation and regression slope -- rounded to these many
#: decimal places, as Anscombe's original table reports them (2 decimal
#: places for the summary statistics themselves).
ANSCOMBE_AGREEMENT_DECIMALS: int = 1

# --------------------------------------------------------------------------
# Exercise 7: Simpson's paradox
# --------------------------------------------------------------------------
# The smallest integer table that shows the paradox: treatment A wins both
# subgroups, treatment B wins overall, purely because A's trials are
# concentrated in the harder subgroup and B's in the easier one.

#: (successes, trials) for treatment A, by subgroup.
TREATMENT_A_EASY: tuple[int, int] = (1, 1)     # 100%
TREATMENT_A_HARD: tuple[int, int] = (9, 90)    # 10%

#: (successes, trials) for treatment B, by subgroup.
TREATMENT_B_EASY: tuple[int, int] = (9, 10)    # 90%
TREATMENT_B_HARD: tuple[int, int] = (0, 1)     # 0%

# --------------------------------------------------------------------------
# Exercise 8: robust spread under contamination
# --------------------------------------------------------------------------

CONTAMINATION_SEED: int = 116
CONTAMINATION_BASE_MEAN: float = 100.0
CONTAMINATION_BASE_SIGMA: float = 5.0
CONTAMINATION_BASE_N: int = 97

#: Three extreme values added to 97 clean ones -- 3 out of 100, i.e. 3%
#: contamination.
CONTAMINATION_OUTLIERS: tuple[float, ...] = (500.0, 520.0, 480.0)

#: The standard deviation must inflate by at least this multiplier for the
#: demonstration to count as "much more" than the MAD's shift.
CONTAMINATION_STD_MULTIPLIER_FLOOR: float = 5.0

#: The MAD must stay below this multiplier of its clean value -- "barely
#: moves".
CONTAMINATION_MAD_MULTIPLIER_CEILING: float = 1.5

# --------------------------------------------------------------------------
# Exercise 9: standardisation and z-scores
# --------------------------------------------------------------------------

STANDARDIZATION_SEED: int = 116
STANDARDIZATION_N: int = 30
STANDARDIZATION_X_MEAN: float = 50.0
STANDARDIZATION_X_SIGMA: float = 10.0
STANDARDIZATION_Y_SLOPE: float = 2.0
STANDARDIZATION_Y_NOISE_SIGMA: float = 5.0

STANDARDIZATION_MEAN_TOLERANCE: float = 1e-9
STANDARDIZATION_STD_TOLERANCE: float = 1e-9
STANDARDIZATION_CORRELATION_TOLERANCE: float = 1e-9


def standard_error_of_mean(sample_variance: float, n: int) -> float:
    """The standard error of a sample mean estimated from n draws."""
    return math.sqrt(sample_variance / n)

"""Exercises 1, 2, 4, 5, 6, 7 and 9: statistics computed from scratch, and
checked against a second, independent way of getting the same number.

Every function that has an exact answer avoids floating-point surprises
where it can (the mean and mode use exact arithmetic on the inputs given
here); every function that depends on an interpolation convention (the
percentile) makes that convention an explicit argument rather than hiding a
default.
"""

import math
from collections import Counter
from typing import Sequence


# ---------------------------------------------------------------------------
# Exercise 1: mean, median, mode, from scratch
# ---------------------------------------------------------------------------


def mean(values: Sequence[float]) -> float:
    """The arithmetic mean: the sum divided by the count."""
    values = list(values)
    return sum(values) / len(values)


def median(values: Sequence[float]) -> float:
    """The middle value once sorted; the average of the two middle values
    when the count is even."""
    ordered = sorted(values)
    n = len(ordered)
    mid = n // 2
    if n % 2 == 1:
        return float(ordered[mid])
    return (ordered[mid - 1] + ordered[mid]) / 2.0


def modes(values: Sequence[float]) -> list[float]:
    """Every value that occurs the maximum number of times. A list, not a
    single value, because a distribution can be multimodal."""
    counts = Counter(values)
    top = max(counts.values())
    return sorted(v for v, c in counts.items() if c == top)


# ---------------------------------------------------------------------------
# Exercise 2: the breakdown point
# ---------------------------------------------------------------------------


def breakdown_point_mean(values: Sequence[float], corrupted_value: float) -> tuple[float, float]:
    """Replace the largest value with `corrupted_value`; return
    (mean before, mean after)."""
    ordered = sorted(values)
    before = mean(ordered)
    corrupted = ordered[:-1] + [corrupted_value]
    after = mean(corrupted)
    return before, after


def breakdown_point_median(values: Sequence[float], corrupted_value: float) -> tuple[float, float]:
    """The same replacement, but tracking the median instead of the mean."""
    ordered = sorted(values)
    before = median(ordered)
    corrupted = ordered[:-1] + [corrupted_value]
    after = median(corrupted)
    return before, after


# ---------------------------------------------------------------------------
# Exercise 4: percentile ambiguity -- deliberately NOT resolved to one
# function. The point of this exercise is that "the" 75th percentile does
# not exist; this helper just calls NumPy with an explicit method, so every
# call site is forced to say which convention it means.
# ---------------------------------------------------------------------------


def percentile_under(values: Sequence[float], target: float, method: str) -> float:
    """`numpy.percentile` under one named interpolation convention."""
    import numpy as np

    return float(np.percentile(np.asarray(values, dtype=float), target, method=method))


# ---------------------------------------------------------------------------
# Exercise 5: Pearson versus Spearman
# ---------------------------------------------------------------------------


def pearson(x: Sequence[float], y: Sequence[float]) -> float:
    """Pearson's r: the standardized covariance, measuring LINEAR
    association only."""
    import statistics as st

    return st.correlation(list(x), list(y))


def _rank(values: Sequence[float]) -> list[float]:
    """Fractional (average) ranks, so tied values share the mean of the
    rank positions they occupy -- the standard convention Spearman's
    correlation uses."""
    order = sorted(range(len(values)), key=lambda i: values[i])
    ranks = [0.0] * len(values)
    i = 0
    while i < len(order):
        j = i
        while j + 1 < len(order) and values[order[j + 1]] == values[order[i]]:
            j += 1
        average_rank = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[order[k]] = average_rank
        i = j + 1
    return ranks


def spearman(x: Sequence[float], y: Sequence[float]) -> float:
    """Spearman's rank correlation: Pearson's r computed on the RANKS of x
    and y, measuring monotone association regardless of shape."""
    return pearson(_rank(list(x)), _rank(list(y)))


# ---------------------------------------------------------------------------
# Exercise 6: Anscombe's quartet
# ---------------------------------------------------------------------------


def anscombe_summary(x: Sequence[float], y: Sequence[float]) -> dict[str, float]:
    """The five classic summary statistics: mean x, mean y, variance x,
    variance y, Pearson correlation, and the fitted regression slope."""
    import statistics as st

    x, y = list(x), list(y)
    slope, intercept = st.linear_regression(x, y)
    return {
        "mean_x": st.fmean(x),
        "mean_y": st.fmean(y),
        "var_x": st.variance(x),
        "var_y": st.variance(y),
        "correlation": st.correlation(x, y),
        "slope": slope,
        "intercept": intercept,
    }


def shape_statistics(x: Sequence[float], y: Sequence[float]) -> dict[str, float]:
    """Three diagnostics the five classic summary numbers do NOT capture --
    each one, on its own, is unremarkable for an ordinary linear
    relationship, and each one is exactly what makes one Anscombe set
    different from the other three, for three different structural
    reasons.

    - ``max_leverage``: how much a single x-value alone (before y is even
      considered) could determine the fitted line. Depends only on the x
      values, so it is identical for any two datasets that share an x
      column -- and dramatically different for a dataset whose x values
      are not spread out the same way.
    - ``outlier_ratio``: the largest single residual, divided by the sum
      of every OTHER residual's magnitude. Large when one point's
      deviation from the fitted line dwarfs everyone else's combined.
    - ``residual_sign_changes``: walking the residuals in x-order, how
      many times the sign flips. Scattered, honestly linear noise flips
      sign often; one smooth systematic curve (a fitted line failing to
      follow a genuine bend) flips rarely.
    """
    import statistics as st

    x, y = list(x), list(y)
    n = len(x)
    slope, intercept = st.linear_regression(x, y)
    residuals = [yi - (slope * xi + intercept) for xi, yi in zip(x, y)]

    mean_x = st.fmean(x)
    ss_x = sum((xi - mean_x) ** 2 for xi in x)
    leverages = [1.0 / n + (xi - mean_x) ** 2 / ss_x for xi in x]

    abs_residuals = sorted((abs(r) for r in residuals), reverse=True)
    largest, rest = abs_residuals[0], sum(abs_residuals[1:])
    outlier_ratio = largest / rest if rest > 0 else math.inf

    order = sorted(range(n), key=lambda i: x[i])
    ordered_residuals = [residuals[i] for i in order]
    sign_changes = sum(
        1
        for i in range(n - 1)
        if ordered_residuals[i] * ordered_residuals[i + 1] < 0
    )

    return {
        "max_leverage": max(leverages),
        "outlier_ratio": outlier_ratio,
        "residual_sign_changes": float(sign_changes),
    }


# ---------------------------------------------------------------------------
# Exercise 7: Simpson's paradox
# ---------------------------------------------------------------------------


def success_rate(successes: int, trials: int) -> float:
    return successes / trials


def combined_rate(*subgroups: tuple[int, int]) -> float:
    """The overall success rate across several (successes, trials)
    subgroups -- NOT the average of the subgroup rates, but the total
    successes over the total trials, which is what "overall rate" means and
    exactly where Simpson's paradox hides."""
    total_successes = sum(s for s, _ in subgroups)
    total_trials = sum(t for _, t in subgroups)
    return total_successes / total_trials


# ---------------------------------------------------------------------------
# Exercise 8 helper: median absolute deviation (the simulation itself lives
# in simulate.py, since it draws random contamination)
# ---------------------------------------------------------------------------


def median_absolute_deviation(values: Sequence[float]) -> float:
    """MAD: the median of the absolute deviations from the median. A
    robust measure of spread -- corrupting a small fraction of the data
    moves it far less than it moves the standard deviation."""
    values = list(values)
    m = median(values)
    return median([abs(v - m) for v in values])


def population_std(values: Sequence[float]) -> float:
    """The (ddof=0) standard deviation, used for the clean/contaminated
    comparison in exercise 8 -- either ddof gives the same qualitative
    story, this lab's tests use the sample (ddof=1) version via NumPy for
    consistency with `numpy.std(ddof=1)`."""
    import statistics as st

    return st.pstdev(values)


# ---------------------------------------------------------------------------
# Exercise 9: standardisation and z-scores
# ---------------------------------------------------------------------------


def zscores(values: Sequence[float]) -> list[float]:
    """(x - mean) / population standard deviation, for every value."""
    import statistics as st

    values = list(values)
    m = st.fmean(values)
    s = st.pstdev(values)
    return [(v - m) / s for v in values]

"""Exercises 3 and 8: the two places this lab draws random numbers.

Both use `numpy.random.default_rng(seed)`, an independent generator object
rather than the legacy global-state `numpy.random.seed()` -- so the same
seed gives byte-identical results regardless of what else has run.
"""

import numpy as np


# ---------------------------------------------------------------------------
# Exercise 3: Bessel's correction, measured by simulation
# ---------------------------------------------------------------------------


def bessel_trial_variances(
    rng: np.random.Generator,
    population_mean: float,
    population_sigma: float,
    sample_size: int,
    trials: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Draw `trials` independent samples of size `sample_size` from a
    normal population, and compute both the biased (divide-by-n) and
    unbiased (divide-by-(n-1)) sample variance for every one of them.

    Returns (biased_variances, unbiased_variances), one value per trial.
    """
    samples = rng.normal(
        loc=population_mean, scale=population_sigma, size=(trials, sample_size)
    )
    sample_means = samples.mean(axis=1, keepdims=True)
    squared_deviations = (samples - sample_means) ** 2
    sum_sq = squared_deviations.sum(axis=1)
    biased = sum_sq / sample_size
    unbiased = sum_sq / (sample_size - 1)
    return biased, unbiased


# ---------------------------------------------------------------------------
# Exercise 8: robust spread under contamination
# ---------------------------------------------------------------------------


def contaminated_sample(
    rng: np.random.Generator,
    clean_mean: float,
    clean_sigma: float,
    clean_n: int,
    outliers: tuple[float, ...],
) -> tuple[np.ndarray, np.ndarray]:
    """A clean sample from a normal distribution, and the same sample with
    a handful of extreme values appended.

    Returns (clean, contaminated).
    """
    clean = rng.normal(loc=clean_mean, scale=clean_sigma, size=clean_n)
    contaminated = np.concatenate([clean, np.asarray(outliers, dtype=float)])
    return clean, contaminated

"""The reference pytest suite: real values, real exceptions, run from
inside `examples/` (see `conftest.py`)."""

import statistics as st

import numpy as np
import pytest

import dataset as D
import descriptive as F
import simulate as S


# ---------------------------------------------------------------------------
# Exercise 1: mean, median, mode
# ---------------------------------------------------------------------------


def test_mean_matches_statistics_module():
    assert F.mean(D.ODD_LIST) == st.fmean(D.ODD_LIST)


def test_median_odd_length_matches_statistics_module():
    assert F.median(D.ODD_LIST) == st.median(D.ODD_LIST)


def test_median_even_length_averages_the_two_middle_values():
    assert F.median(D.EVEN_LIST) == st.median(D.EVEN_LIST) == 6.0


def test_mode_single_peak():
    assert F.modes(D.ODD_LIST) == [7]


def test_mode_multimodal_returns_every_tied_value():
    assert F.modes(D.MULTIMODAL_LIST) == sorted(st.multimode(D.MULTIMODAL_LIST)) == [3, 8]


def test_statistics_mode_singular_silently_drops_a_tied_value():
    single = st.mode(D.MULTIMODAL_LIST)
    assert single in F.modes(D.MULTIMODAL_LIST)
    assert len(F.modes(D.MULTIMODAL_LIST)) > 1


# ---------------------------------------------------------------------------
# Exercise 2: the breakdown point
# ---------------------------------------------------------------------------


def test_mean_breakdown_point_is_dragged_far():
    before, after = F.breakdown_point_mean(D.SALARY_LIST, D.CORRUPTED_SALARY)
    assert after - before > D.BREAKDOWN_MEAN_SHIFT_FLOOR


def test_median_breakdown_point_does_not_move_at_all():
    before, after = F.breakdown_point_median(D.SALARY_LIST, D.CORRUPTED_SALARY)
    assert after == before  # exact equality: the median's rank did not change


def test_mean_and_median_agree_on_the_uncorrupted_data():
    ordered = sorted(D.SALARY_LIST)
    assert F.mean(ordered) == pytest.approx(50777.78, abs=0.01)
    assert F.median(ordered) == 50000


# ---------------------------------------------------------------------------
# Exercise 3: Bessel's correction, measured
# ---------------------------------------------------------------------------


def test_bessel_divide_by_n_estimator_is_biased_low_by_n_minus_1_over_n():
    rng = np.random.default_rng(D.BESSEL_SEED)
    biased, _ = S.bessel_trial_variances(
        rng,
        D.BESSEL_POPULATION_MEAN,
        D.BESSEL_POPULATION_SIGMA,
        D.BESSEL_SAMPLE_SIZE,
        D.BESSEL_TRIALS,
    )
    ratio = float(biased.mean()) / D.BESSEL_TRUE_VARIANCE
    assert ratio == pytest.approx(D.BESSEL_EXPECTED_BIAS_FACTOR, abs=D.BESSEL_BIAS_FACTOR_TOLERANCE)


def test_bessel_divide_by_n_minus_1_estimator_is_unbiased_within_tolerance():
    rng = np.random.default_rng(D.BESSEL_SEED)
    _, unbiased = S.bessel_trial_variances(
        rng,
        D.BESSEL_POPULATION_MEAN,
        D.BESSEL_POPULATION_SIGMA,
        D.BESSEL_SAMPLE_SIZE,
        D.BESSEL_TRIALS,
    )
    mean_unbiased = float(unbiased.mean())
    se = float(unbiased.std(ddof=1)) / (D.BESSEL_TRIALS**0.5)
    assert abs(mean_unbiased - D.BESSEL_TRUE_VARIANCE) < D.BESSEL_UNBIASED_SE_TOLERANCE * se


def test_bessel_unbiased_estimator_is_closer_to_truth_than_biased_on_average():
    rng = np.random.default_rng(D.BESSEL_SEED)
    biased, unbiased = S.bessel_trial_variances(
        rng,
        D.BESSEL_POPULATION_MEAN,
        D.BESSEL_POPULATION_SIGMA,
        D.BESSEL_SAMPLE_SIZE,
        D.BESSEL_TRIALS,
    )
    assert abs(float(unbiased.mean()) - D.BESSEL_TRUE_VARIANCE) < abs(
        float(biased.mean()) - D.BESSEL_TRUE_VARIANCE
    )


# ---------------------------------------------------------------------------
# Exercise 4: percentile ambiguity
# ---------------------------------------------------------------------------


def test_percentile_conventions_disagree():
    values = {
        method: F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, method)
        for method in D.PERCENTILE_METHODS
    }
    assert len(set(values.values())) >= 2


def test_percentile_default_linear_method_matches_the_documented_value():
    assert F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, "linear") == 8.25


def test_percentile_lower_and_higher_land_on_different_real_data_points():
    lower = F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, "lower")
    higher = F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, "higher")
    assert lower != higher
    assert lower in D.PERCENTILE_ARRAY
    assert higher in D.PERCENTILE_ARRAY


# ---------------------------------------------------------------------------
# Exercise 5: Pearson versus Spearman
# ---------------------------------------------------------------------------


def test_pearson_on_a_symmetric_parabola_is_essentially_zero():
    r = F.pearson(D.PARABOLA_X, D.PARABOLA_Y)
    assert abs(r) < D.PARABOLA_PEARSON_TOLERANCE


def test_spearman_on_a_monotone_cubic_is_exactly_one():
    assert F.spearman(D.MONOTONE_X, D.MONOTONE_Y) == 1.0


def test_pearson_on_the_same_cubic_is_strong_but_not_perfect():
    r = F.pearson(D.MONOTONE_X, D.MONOTONE_Y)
    assert 0.0 < r < 1.0


# ---------------------------------------------------------------------------
# Exercise 6: Anscombe's quartet
# ---------------------------------------------------------------------------


@pytest.mark.parametrize("name", ["I", "II", "III", "IV"])
def test_anscombe_summaries_agree_to_documented_precision(name):
    reference = F.anscombe_summary(*D.ANSCOMBE_SETS["I"])
    summary = F.anscombe_summary(*D.ANSCOMBE_SETS[name])
    dec = D.ANSCOMBE_AGREEMENT_DECIMALS
    assert round(summary["mean_x"], dec) == round(reference["mean_x"], dec)
    assert round(summary["mean_y"], dec) == round(reference["mean_y"], dec)
    assert round(summary["var_x"], dec) == round(reference["var_x"], dec)
    assert round(summary["var_y"], dec) == round(reference["var_y"], dec)
    assert round(summary["correlation"], 1) == round(reference["correlation"], 1)
    assert round(summary["slope"], 1) == round(reference["slope"], 1)


def test_anscombe_set_iv_has_dramatically_higher_leverage():
    shapes = {name: F.shape_statistics(x, y) for name, (x, y) in D.ANSCOMBE_SETS.items()}
    assert shapes["IV"]["max_leverage"] > 3.0 * shapes["I"]["max_leverage"]
    assert shapes["I"]["max_leverage"] == shapes["II"]["max_leverage"] == shapes["III"]["max_leverage"]


def test_anscombe_set_iii_has_a_dominant_outlier_residual():
    shapes = {name: F.shape_statistics(x, y) for name, (x, y) in D.ANSCOMBE_SETS.items()}
    assert shapes["III"]["outlier_ratio"] > 2.0 * shapes["I"]["outlier_ratio"]


def test_anscombe_set_ii_residuals_change_sign_less_often_than_set_i():
    shapes = {name: F.shape_statistics(x, y) for name, (x, y) in D.ANSCOMBE_SETS.items()}
    assert shapes["II"]["residual_sign_changes"] < shapes["I"]["residual_sign_changes"]


# ---------------------------------------------------------------------------
# Exercise 7: Simpson's paradox
# ---------------------------------------------------------------------------


def test_treatment_a_wins_every_subgroup():
    assert F.success_rate(*D.TREATMENT_A_EASY) > F.success_rate(*D.TREATMENT_B_EASY)
    assert F.success_rate(*D.TREATMENT_A_HARD) > F.success_rate(*D.TREATMENT_B_HARD)


def test_treatment_b_wins_overall():
    a_total = F.combined_rate(D.TREATMENT_A_EASY, D.TREATMENT_A_HARD)
    b_total = F.combined_rate(D.TREATMENT_B_EASY, D.TREATMENT_B_HARD)
    assert b_total > a_total


def test_overall_rate_is_the_pooled_total_not_the_average_of_subgroup_rates():
    a_total = F.combined_rate(D.TREATMENT_A_EASY, D.TREATMENT_A_HARD)
    naive_average = (
        F.success_rate(*D.TREATMENT_A_EASY) + F.success_rate(*D.TREATMENT_A_HARD)
    ) / 2
    # The pooled rate is dominated by the much larger hard subgroup, so it
    # sits far closer to the hard-subgroup rate than a naive 50/50 average
    # of the two subgroup rates would.
    assert abs(a_total - F.success_rate(*D.TREATMENT_A_HARD)) < abs(a_total - naive_average)


# ---------------------------------------------------------------------------
# Exercise 8: robust spread under contamination
# ---------------------------------------------------------------------------


def test_contamination_inflates_standard_deviation_a_lot():
    rng = np.random.default_rng(D.CONTAMINATION_SEED)
    clean, contaminated = S.contaminated_sample(
        rng,
        D.CONTAMINATION_BASE_MEAN,
        D.CONTAMINATION_BASE_SIGMA,
        D.CONTAMINATION_BASE_N,
        D.CONTAMINATION_OUTLIERS,
    )
    std_clean = float(np.std(clean, ddof=1))
    std_contam = float(np.std(contaminated, ddof=1))
    assert std_contam / std_clean > D.CONTAMINATION_STD_MULTIPLIER_FLOOR


def test_contamination_barely_moves_the_median_absolute_deviation():
    rng = np.random.default_rng(D.CONTAMINATION_SEED)
    clean, contaminated = S.contaminated_sample(
        rng,
        D.CONTAMINATION_BASE_MEAN,
        D.CONTAMINATION_BASE_SIGMA,
        D.CONTAMINATION_BASE_N,
        D.CONTAMINATION_OUTLIERS,
    )
    mad_clean = F.median_absolute_deviation(clean)
    mad_contam = F.median_absolute_deviation(contaminated)
    assert mad_contam / mad_clean < D.CONTAMINATION_MAD_MULTIPLIER_CEILING


# ---------------------------------------------------------------------------
# Exercise 9: standardisation
# ---------------------------------------------------------------------------


def test_standardized_sample_has_mean_zero_and_std_one():
    rng = np.random.default_rng(D.STANDARDIZATION_SEED)
    x = rng.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
    zx = F.zscores(x)
    mean_zx = sum(zx) / len(zx)
    std_zx = (sum((v - mean_zx) ** 2 for v in zx) / len(zx)) ** 0.5
    assert abs(mean_zx) < D.STANDARDIZATION_MEAN_TOLERANCE
    assert abs(std_zx - 1.0) < D.STANDARDIZATION_STD_TOLERANCE


def test_standardizing_does_not_change_pearson_correlation():
    rng = np.random.default_rng(D.STANDARDIZATION_SEED)
    x = rng.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
    noise = rng.normal(0.0, D.STANDARDIZATION_Y_NOISE_SIGMA, D.STANDARDIZATION_N)
    y = D.STANDARDIZATION_Y_SLOPE * x + noise
    zx, zy = F.zscores(x), F.zscores(y)
    r_original = F.pearson(x, y)
    r_standardized = F.pearson(zx, zy)
    assert abs(r_original - r_standardized) < D.STANDARDIZATION_CORRELATION_TOLERANCE

lesson_id: D116
day: 116
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-116-descriptive-statistics-that-dont-lie
  - 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_mean_median_mode.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_breakdown_point.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_bessel_correction.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_percentile_ambiguity.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_pearson_vs_spearman.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_anscombes_quartet.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_simpsons_paradox.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 08_robust_spread_under_contamination.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 09_standardization.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 -> 55 checks, 0 failure(s), exit 0; pytest examples -> 32 passed; pytest starter -> 1 passed, 38 skipped on an untouched checkout, and 39 passed against a fully solved copy of starter/ (verified by temporarily copying the reference descriptive.py and simulate.py into starter/ and filling in every answers.py prediction, then restoring the blank skeletons -- the skip counts before and after that restore were confirmed identical: 1 passed, 38 skipped both times). 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 6 of the harness re-runs the reference pytest suite with one assertion (the breakdown-point median-does-not-move check) temporarily flipped to the wrong direction, confirms the run exits non-zero with exactly one failure named in the output, and restores the original file -- so the suite is demonstrated to be capable of failing rather than merely claimed to be. Three honesty notes from this run. FIRST: scipy and pandas are not installed in this environment. Both 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: exercises 3, 8 and 9 report real, freshly measured numbers from numpy.random.default_rng(116) rather than fixed literals -- on this run, dividing by n underestimated the true variance (100.0) by a measured ratio of 0.8007 against a predicted (n-1)/n of 0.8000, and the divide-by-(n-1) estimator landed 0.17 standard errors from the truth; 3% contamination inflated the standard deviation by a measured 15.10x and the median absolute deviation by a measured 1.00x; standardising left the Pearson correlation unchanged to a measured difference of 3.33e-16. These exact figures will differ slightly on another machine or NumPy version; the tests assert against fixed floors, ceilings and tolerances rather than these specific digits, as documented in expected-output/FIELDS.md. THIRD: Anscombe''s quartet (exercise 6) is the published 1973 dataset (Anscombe, "Graphs in Statistical Analysis," The American Statistician, 27(1), 17-21) reproduced exactly, not generated for this lab; the Simpson''s-paradox table (exercise 7) is an invented illustrative example, not real clinical or experimental data, and is described as such in security.md and the lesson.'

# 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.percentile` under nine interpolation conventions (exercise 4), and the `numpy.random.Generator` built by `default_rng(seed)` used for the Bessel-correction simulation (exercise 3) and the contamination sample (exercise 8). |
| `pytest` | 9.1.1 | MIT | The reference suite 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 5 of
`tests/run_tests.sh` greps every source file in `examples/` and `starter/`
to prove that nothing else does.

## What most of this lab does not need NumPy for at all

Exercises 1, 2, 5, 6, 7 and 9 use only the standard library —
`statistics`, `collections.Counter`, and plain arithmetic. Only exercise 4
(the percentile conventions) genuinely needs `numpy.percentile`'s
`method=` argument, which the standard library has no equivalent for.
Exercises 3 and 8 need `numpy.random.Generator` for fast, reproducible
simulation, though Python's own `random.Random(seed)` could stand in at
the cost of a slower Python-level loop instead of a vectorised call —
`troubleshooting.md` shows the substitution.

## What is deliberately *not* installed

`pandas` and `scipy.stats` are **not installed in this environment, and no
output from either is reproduced anywhere** in this lab or its lesson.
`pandas.DataFrame.describe()` and `scipy.stats` are both described from
their public documentation in the lesson's Tools section, and both are
explicitly marked as not run here.

That is not a limitation to apologise for. Every statistic this lab
computes is either exact arithmetic from the standard library, an explicit
call to a named NumPy convention, or a simulation you run yourself with
`numpy.random.Generator` — nothing here depends on `pandas` or `scipy` to
be correct.

numpy==2.5.2
pytest==9.1.1

# The nine exercises

Work through these in order. Predict the answer to each `answers.py`
question *before* running anything — the breakdown-point question and the
percentile-agreement question only catch you if you commit to a guess
first.

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. Mean, median, mode (`descriptive.py`)

Write `mean()`, `median()` (handling both odd and even lengths) and
`modes()` (returning every tied value, not just one) from scratch. Assert
against the `statistics` module on the same inputs.

## 2. The breakdown point (`descriptive.py`)

`breakdown_point_mean()` and `breakdown_point_median()`. Replace the salary
list's largest value with an absurd one and watch the mean move by over a
million dollars while the median does not move at all — exact equality is
the right assertion for the median.

## 3. Bessel's correction, measured (`simulate.py`)

`bessel_trial_variances()` draws many small samples from a population of
known variance and computes both the divide-by-n and divide-by-(n-1)
estimators for each one. The divide-by-n estimator should average out
biased low by the factor `(n-1)/n`; the divide-by-(n-1) estimator should
land within a few standard errors of the truth.

## 4. Percentile ambiguity (`descriptive.py`)

`percentile_under()` wraps `numpy.percentile` with an explicit `method=`.
Call it under several conventions on the same small array and confirm that
at least two of them genuinely disagree about the 75th percentile.

## 5. Pearson versus Spearman (`descriptive.py`)

`pearson()` and `spearman()`. A perfect symmetric parabola should give a
Pearson correlation essentially zero; a perfect monotone cubic should give
a Spearman correlation of exactly 1.0.

## 6. Anscombe's quartet (`descriptive.py`)

`anscombe_summary()` computes the five classic statistics (mean x, mean y,
variance x, variance y, correlation, slope) for each of the four published
sets — they should all agree. `shape_statistics()` computes three
diagnostics the classic five cannot see, and those should tell the four
sets apart.

## 7. Simpson's paradox (`descriptive.py`)

`success_rate()` and `combined_rate()`. Confirm treatment A beats treatment
B in *both* subgroups of the smallest table that shows the paradox, and
that treatment B still wins *overall* — both directions, from the same
four numbers.

## 8. Robust spread under contamination (`descriptive.py`, `simulate.py`)

`contaminated_sample()` adds a handful of extreme values to a clean sample.
`median_absolute_deviation()` measures spread the way the standard
deviation does, but robustly. Confirm the standard deviation inflates by a
large multiplier while the MAD barely moves.

## 9. Standardisation (`descriptive.py`)

`zscores()`. Confirm the standardised sample has mean 0 and standard
deviation 1, and that standardising does not change the Pearson correlation
between two variables.

"""Exercises 1 through 9 -- seventeen predictions.

Replace each `None` with the value you think is correct. A `None` is a skip,
not a failure: `pytest starter -q` counts only what you have attempted. When
you are wrong it prints both your answer and the real one, so a wrong guess
is worth more than a blank.

Predict BEFORE you run anything. The two that catch almost everyone are the
breakdown-point question (exercise 2) and the percentile-agreement question
(exercise 4) -- and they only catch you if you commit to a guess first.

Every answer is a number or a Python bool.
"""

ANSWERS: dict[str, object] = {
    # ----------------------------------------------------------------------
    # Exercise 1 -- mean, median, mode
    # ----------------------------------------------------------------------
    # 1.1 The mean of (2, 4, 4, 7, 7, 7, 9, 12, 15), as a decimal.
    "odd_list_mean": None,
    # 1.2 The median of (1.0, 3.0, 4.0, 8.0, 10.0, 12.0) -- an even-length
    #     list, so this is an average of two values.
    "even_list_median": None,
    # ----------------------------------------------------------------------
    # Exercise 2 -- the breakdown point
    # ----------------------------------------------------------------------
    # 2.1 Replace the salary list's largest value ($60,000) with
    #     $10,000,000. Does the MEAN move by more than $500,000?
    "mean_breaks_down": None,
    # 2.2 Under the same corruption, does the MEDIAN move at all
    #     (True/False)?
    "median_breaks_down": None,
    # ----------------------------------------------------------------------
    # Exercise 3 -- Bessel's correction
    # ----------------------------------------------------------------------
    # 3.1 Divide by n instead of n-1: is the resulting variance estimate
    #     biased HIGH, LOW, or unbiased? Answer "high", "low", or
    #     "unbiased".
    "divide_by_n_bias_direction": None,
    # ----------------------------------------------------------------------
    # Exercise 4 -- percentile ambiguity
    # ----------------------------------------------------------------------
    # 4.1 Across NumPy's nine `method=` conventions, do at least two of them
    #     disagree on the 75th percentile of (1, 2, 3, 4, 6, 8, 9, 15)?
    "percentile_methods_disagree": None,
    # 4.2 The default ('linear') method's answer, as a decimal.
    "percentile_linear_value": None,
    # ----------------------------------------------------------------------
    # Exercise 5 -- Pearson versus Spearman
    # ----------------------------------------------------------------------
    # 5.1 Pearson correlation of a symmetric parabola (y = x^2 over a
    #     symmetric range of x): close to 1, close to 0, or close to -1?
    #     Answer "close_to_zero", "close_to_one", or "close_to_negative_one".
    "parabola_pearson_magnitude": None,
    # 5.2 Spearman correlation of a monotone cubic (y = x^3): exactly what
    #     number?
    "monotone_spearman_value": None,
    # ----------------------------------------------------------------------
    # Exercise 6 -- Anscombe's quartet
    # ----------------------------------------------------------------------
    # 6.1 Do all four Anscombe sets share the same mean of x (to one
    #     decimal place)?
    "anscombe_means_agree": None,
    # 6.2 Does set IV have dramatically higher "leverage" on its one
    #     non-repeated x-value than set I?
    "anscombe_set_iv_leverage_dominant": None,
    # ----------------------------------------------------------------------
    # Exercise 7 -- Simpson's paradox
    # ----------------------------------------------------------------------
    # 7.1 Treatment A beats treatment B in BOTH subgroups. Does treatment B
    #     still win OVERALL?
    "simpson_b_wins_overall": None,
    # ----------------------------------------------------------------------
    # Exercise 8 -- robust spread under contamination
    # ----------------------------------------------------------------------
    # 8.1 3% contamination: does the standard deviation inflate by more
    #     than 5x?
    "contamination_inflates_std": None,
    # 8.2 Under the same contamination, does the median absolute deviation
    #     stay under 1.5x its clean value?
    "contamination_mad_stable": None,
    # ----------------------------------------------------------------------
    # Exercise 9 -- standardisation
    # ----------------------------------------------------------------------
    # 9.1 After standardising a sample, its mean is (approximately) what
    #     number?
    "standardized_mean": None,
    # 9.2 After standardising a sample, its standard deviation is
    #     (approximately) what number?
    "standardized_std": None,
    # 9.3 Does standardising CHANGE the Pearson correlation between two
    #     variables?
    "standardizing_changes_correlation": None,
}

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

Both `examples/` and `starter/` contain modules called `dataset`,
`descriptive`, `simulate` and `answers`, 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 `descriptive` was
seen first and then reuse it for the other suite -- 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 ("dataset", "descriptive", "simulate", "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 dataset and tolerance this lab compares against.

Read this file. Nothing here is tuned to make a test pass: the salary list
is small and hand-checkable, Anscombe's quartet is the published 1973
dataset (Anscombe, "Graphs in Statistical Analysis", The American
Statistician, 1973), and the Bessel-correction tolerance is derived from the
standard error of the simulated mean, with the arithmetic written out beside
it, never chosen by running a test and loosening the number until it
passed.
"""

import math
from fractions import Fraction

# --------------------------------------------------------------------------
# Exercise 1: mean, median, mode
# --------------------------------------------------------------------------

#: An odd-length list with a single clear mode (7 appears three times).
ODD_LIST: tuple[int, ...] = (2, 4, 4, 7, 7, 7, 9, 12, 15)

#: An even-length list, so the median must average the two middle values.
EVEN_LIST: tuple[float, ...] = (1.0, 3.0, 4.0, 8.0, 10.0, 12.0)

#: A multimodal list: 3 and 8 are each the most frequent value.
MULTIMODAL_LIST: tuple[int, ...] = (3, 3, 5, 6, 8, 8, 9)

# --------------------------------------------------------------------------
# Exercise 2: the breakdown point
# --------------------------------------------------------------------------

#: Nine ordinary salaries, in dollars. Small enough to check the median of
#: by eye: sorted, the middle (5th of 9) value is 50000.
SALARY_LIST: tuple[int, ...] = (
    42000,
    45000,
    47000,
    48000,
    50000,
    52000,
    55000,
    58000,
    60000,
)

#: Replace the single largest salary with one wildly corrupted value. One
#: value out of nine changing is "one value out of many" -- nowhere near the
#: 50% the median can absorb before it moves at all.
CORRUPTED_SALARY: int = 10_000_000

#: The mean must move by at least this many dollars for the demonstration to
#: count as "dragged anywhere at all" -- chosen as a small fraction of the
#: actual shift (over $1,000,000), so the assertion is not a coin flip.
BREAKDOWN_MEAN_SHIFT_FLOOR: float = 500_000.0

# --------------------------------------------------------------------------
# Exercise 3: Bessel's correction, measured
# --------------------------------------------------------------------------

#: The population this exercise draws from: mean 50, standard deviation 10,
#: so the true variance is exactly 100. Both are known exactly because the
#: population is synthetic, which is what makes the bias measurable at all
#: -- there is a ground truth to compare the estimators against.
BESSEL_POPULATION_MEAN: float = 50.0
BESSEL_POPULATION_SIGMA: float = 10.0
BESSEL_TRUE_VARIANCE: float = BESSEL_POPULATION_SIGMA**2

#: Small samples, where the fitted-mean bias is largest relative to n.
BESSEL_SAMPLE_SIZE: int = 5

#: Many repeated samples, so the estimators' average is itself a precise
#: number rather than one noisy draw.
BESSEL_TRIALS: int = 20_000

BESSEL_SEED: int = 116

#: The textbook claim this exercise measures: dividing by n underestimates
#: the true variance, on average, by exactly the factor (n-1)/n.
BESSEL_EXPECTED_BIAS_FACTOR: float = (BESSEL_SAMPLE_SIZE - 1) / BESSEL_SAMPLE_SIZE

#: How close the measured biased-estimator ratio must land to (n-1)/n. This
#: is a Monte Carlo measurement, not exact arithmetic, so the tolerance is
#: generous but not vacuous -- wide enough to pass on any honest run of
#: BESSEL_TRIALS trials, narrow enough that a genuinely broken implementation
#: (dividing by the wrong thing, or not at all) fails it.
BESSEL_BIAS_FACTOR_TOLERANCE: float = 0.02

#: How many standard errors of the *unbiased* estimator's own sampling mean
#: it is allowed to sit from the true variance. The unbiased estimator is
#: correct only on average, not on every run, so this must be generous
#: enough to pass reliably -- 4 standard errors covers essentially all honest
#: runs (about 1 in 16,000 fails purely by chance under a normal
#: approximation).
BESSEL_UNBIASED_SE_TOLERANCE: float = 4.0

# --------------------------------------------------------------------------
# Exercise 4: percentile ambiguity
# --------------------------------------------------------------------------

#: A small, deliberately awkward array: 8 values, so the 75th percentile
#: falls between two of them under every interpolation convention, and the
#: conventions genuinely disagree about where.
PERCENTILE_ARRAY: tuple[int, ...] = (1, 2, 3, 4, 6, 8, 9, 15)

PERCENTILE_TARGET: float = 75.0

#: A representative slice of NumPy's nine documented `method=` conventions
#: for `numpy.percentile` (NumPy >= 1.22). `linear` is the default and the
#: one pandas' `DataFrame.describe()` also uses.
PERCENTILE_METHODS: tuple[str, ...] = (
    "linear",
    "lower",
    "higher",
    "nearest",
    "midpoint",
    "weibull",
    "median_unbiased",
    "normal_unbiased",
    "hazen",
)

# --------------------------------------------------------------------------
# Exercise 5: Pearson versus Spearman
# --------------------------------------------------------------------------

#: A perfect, symmetric parabola: y = x^2 over a symmetric range of x. Every
#: unit increase in x on the left half is mirrored by an equal decrease in y
#: on the right half's counterpart, so the *linear* trend that Pearson
#: measures cancels out exactly, even though y is perfectly determined by x.
PARABOLA_X: tuple[int, ...] = tuple(range(-5, 6))
PARABOLA_Y: tuple[int, ...] = tuple(v**2 for v in PARABOLA_X)

#: The tolerance "essentially zero" is checked against -- not equality to 0,
#: because that would be asserting a specific floating-point outcome rather
#: than the mathematical fact.
PARABOLA_PEARSON_TOLERANCE: float = 1e-9

#: A monotone but non-linear relationship: y = x^3. Every increase in x
#: produces an increase in y, so the *rank order* is perfectly preserved,
#: even though the relationship is not a straight line.
MONOTONE_X: tuple[int, ...] = tuple(range(-5, 6))
MONOTONE_Y: tuple[int, ...] = tuple(v**3 for v in MONOTONE_X)

# --------------------------------------------------------------------------
# Exercise 6: Anscombe's quartet
# --------------------------------------------------------------------------
# The published 1973 dataset, reproduced exactly. Anscombe, F. J. (1973).
# "Graphs in Statistical Analysis." The American Statistician, 27(1), 17-21.

ANSCOMBE_X_I: tuple[float, ...] = (10.0, 8.0, 13.0, 9.0, 11.0, 14.0, 6.0, 4.0, 12.0, 7.0, 5.0)
ANSCOMBE_Y_I: tuple[float, ...] = (
    8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68,
)

ANSCOMBE_X_II: tuple[float, ...] = ANSCOMBE_X_I
ANSCOMBE_Y_II: tuple[float, ...] = (
    9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74,
)

ANSCOMBE_X_III: tuple[float, ...] = ANSCOMBE_X_I
ANSCOMBE_Y_III: tuple[float, ...] = (
    7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73,
)

ANSCOMBE_X_IV: tuple[float, ...] = (8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 19.0, 8.0, 8.0, 8.0)
ANSCOMBE_Y_IV: tuple[float, ...] = (
    6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89,
)

ANSCOMBE_SETS: dict[str, tuple[tuple[float, ...], tuple[float, ...]]] = {
    "I": (ANSCOMBE_X_I, ANSCOMBE_Y_I),
    "II": (ANSCOMBE_X_II, ANSCOMBE_Y_II),
    "III": (ANSCOMBE_X_III, ANSCOMBE_Y_III),
    "IV": (ANSCOMBE_X_IV, ANSCOMBE_Y_IV),
}

#: How closely the four sets must agree on mean x, mean y, variance x,
#: variance y, correlation and regression slope -- rounded to these many
#: decimal places, as Anscombe's original table reports them (2 decimal
#: places for the summary statistics themselves).
ANSCOMBE_AGREEMENT_DECIMALS: int = 1

# --------------------------------------------------------------------------
# Exercise 7: Simpson's paradox
# --------------------------------------------------------------------------
# The smallest integer table that shows the paradox: treatment A wins both
# subgroups, treatment B wins overall, purely because A's trials are
# concentrated in the harder subgroup and B's in the easier one.

#: (successes, trials) for treatment A, by subgroup.
TREATMENT_A_EASY: tuple[int, int] = (1, 1)     # 100%
TREATMENT_A_HARD: tuple[int, int] = (9, 90)    # 10%

#: (successes, trials) for treatment B, by subgroup.
TREATMENT_B_EASY: tuple[int, int] = (9, 10)    # 90%
TREATMENT_B_HARD: tuple[int, int] = (0, 1)     # 0%

# --------------------------------------------------------------------------
# Exercise 8: robust spread under contamination
# --------------------------------------------------------------------------

CONTAMINATION_SEED: int = 116
CONTAMINATION_BASE_MEAN: float = 100.0
CONTAMINATION_BASE_SIGMA: float = 5.0
CONTAMINATION_BASE_N: int = 97

#: Three extreme values added to 97 clean ones -- 3 out of 100, i.e. 3%
#: contamination.
CONTAMINATION_OUTLIERS: tuple[float, ...] = (500.0, 520.0, 480.0)

#: The standard deviation must inflate by at least this multiplier for the
#: demonstration to count as "much more" than the MAD's shift.
CONTAMINATION_STD_MULTIPLIER_FLOOR: float = 5.0

#: The MAD must stay below this multiplier of its clean value -- "barely
#: moves".
CONTAMINATION_MAD_MULTIPLIER_CEILING: float = 1.5

# --------------------------------------------------------------------------
# Exercise 9: standardisation and z-scores
# --------------------------------------------------------------------------

STANDARDIZATION_SEED: int = 116
STANDARDIZATION_N: int = 30
STANDARDIZATION_X_MEAN: float = 50.0
STANDARDIZATION_X_SIGMA: float = 10.0
STANDARDIZATION_Y_SLOPE: float = 2.0
STANDARDIZATION_Y_NOISE_SIGMA: float = 5.0

STANDARDIZATION_MEAN_TOLERANCE: float = 1e-9
STANDARDIZATION_STD_TOLERANCE: float = 1e-9
STANDARDIZATION_CORRELATION_TOLERANCE: float = 1e-9


def standard_error_of_mean(sample_variance: float, n: int) -> float:
    """The standard error of a sample mean estimated from n draws."""
    return math.sqrt(sample_variance / n)

"""Exercises 1, 2, 4, 5, 6, 7 and 9. Write the functions below.

Every function that has an exact answer should avoid floating-point
surprises where possible; every function that depends on an interpolation
convention (the percentile) should make that convention an explicit
argument rather than hiding a default.

Read `00_brief.md` for the exercise-by-exercise instructions, and
`dataset.py` for the exact inputs each exercise is checked against.
"""

from collections import Counter
from typing import Sequence


# ---------------------------------------------------------------------------
# Exercise 1: mean, median, mode, from scratch
# ---------------------------------------------------------------------------


def mean(values: Sequence[float]) -> float:
    """The arithmetic mean: the sum divided by the count."""
    raise NotImplementedError


def median(values: Sequence[float]) -> float:
    """The middle value once sorted; the average of the two middle values
    when the count is even."""
    raise NotImplementedError


def modes(values: Sequence[float]) -> list[float]:
    """Every value that occurs the maximum number of times, sorted."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 2: the breakdown point
# ---------------------------------------------------------------------------


def breakdown_point_mean(values: Sequence[float], corrupted_value: float) -> tuple[float, float]:
    """Replace the largest value with `corrupted_value`; return
    (mean before, mean after)."""
    raise NotImplementedError


def breakdown_point_median(values: Sequence[float], corrupted_value: float) -> tuple[float, float]:
    """The same replacement, but tracking the median instead of the mean."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 4: percentile ambiguity
# ---------------------------------------------------------------------------


def percentile_under(values: Sequence[float], target: float, method: str) -> float:
    """`numpy.percentile` under one named interpolation convention."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 5: Pearson versus Spearman
# ---------------------------------------------------------------------------


def pearson(x: Sequence[float], y: Sequence[float]) -> float:
    """Pearson's r: the standardized covariance, measuring LINEAR
    association only. `statistics.correlation` does this correctly."""
    raise NotImplementedError


def spearman(x: Sequence[float], y: Sequence[float]) -> float:
    """Spearman's rank correlation: Pearson's r computed on the RANKS of x
    and y. Rank ties by their AVERAGE rank position."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 6: Anscombe's quartet
# ---------------------------------------------------------------------------


def anscombe_summary(x: Sequence[float], y: Sequence[float]) -> dict[str, float]:
    """mean_x, mean_y, var_x, var_y, correlation, slope, intercept."""
    raise NotImplementedError


def shape_statistics(x: Sequence[float], y: Sequence[float]) -> dict[str, float]:
    """max_leverage, outlier_ratio, residual_sign_changes -- see
    `examples/descriptive.py` for the full explanation of each, once you
    have tried this yourself."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 7: Simpson's paradox
# ---------------------------------------------------------------------------


def success_rate(successes: int, trials: int) -> float:
    raise NotImplementedError


def combined_rate(*subgroups: tuple[int, int]) -> float:
    """The overall success rate across several (successes, trials)
    subgroups -- total successes over total trials, NOT the average of the
    subgroup rates."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 8 helper: median absolute deviation
# ---------------------------------------------------------------------------


def median_absolute_deviation(values: Sequence[float]) -> float:
    """The median of the absolute deviations from the median."""
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 9: standardisation and z-scores
# ---------------------------------------------------------------------------


def zscores(values: Sequence[float]) -> list[float]:
    """(x - mean) / population standard deviation, for every value."""
    raise NotImplementedError

"""Exercises 3 and 8: the two places this lab draws random numbers.

Use `numpy.random.default_rng(seed)` -- an independent generator object,
never the legacy global-state `numpy.random.seed()`.
"""

import numpy as np


def bessel_trial_variances(
    rng: np.random.Generator,
    population_mean: float,
    population_sigma: float,
    sample_size: int,
    trials: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Draw `trials` independent samples of size `sample_size` from a
    normal population, and compute both the biased (divide-by-n) and
    unbiased (divide-by-(n-1)) sample variance for every one of them.

    Returns (biased_variances, unbiased_variances), one value per trial.
    """
    raise NotImplementedError


def contaminated_sample(
    rng: np.random.Generator,
    clean_mean: float,
    clean_sigma: float,
    clean_n: int,
    outliers: tuple[float, ...],
) -> tuple[np.ndarray, np.ndarray]:
    """A clean sample from a normal distribution, and the same sample with
    a handful of extreme values appended.

    Returns (clean, contaminated).
    """
    raise NotImplementedError

"""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 numpy as np
import pytest

import answers
import dataset as D
import descriptive as F
import simulate as S


def need(value, what):
    if value is None:
        pytest.skip(f"not attempted yet: {what}")
    return value


def attempt(fn, what):
    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 close(got, want, tol, what):
    assert abs(float(got) - float(want)) < tol, (
        f"{what}: your answer {got!r}, expected {want!r} "
        f"(difference {abs(float(got) - float(want)):.3e}, tolerance {tol:g})"
    )


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 len(D.SALARY_LIST) == 9


# ---------------------------------------------------------------------------
# Exercise 1 -- mean, median, mode
# ---------------------------------------------------------------------------


def test_1_mean_matches_statistics_module():
    got = attempt(lambda: F.mean(D.ODD_LIST), "mean")
    import statistics as st

    assert got == st.fmean(D.ODD_LIST)


def test_1_median_odd_length():
    got = attempt(lambda: F.median(D.ODD_LIST), "median")
    assert got == 7.0


def test_1_median_even_length_averages_the_middle_two():
    got = attempt(lambda: F.median(D.EVEN_LIST), "median")
    assert got == 6.0


def test_1_modes_finds_the_single_peak():
    got = attempt(lambda: F.modes(D.ODD_LIST), "modes")
    assert got == [7]


def test_1_modes_finds_every_tied_value():
    got = attempt(lambda: F.modes(D.MULTIMODAL_LIST), "modes")
    assert got == [3, 8]


def test_1_prediction_odd_list_mean():
    predicted = need(answers.ANSWERS["odd_list_mean"], "odd_list_mean prediction")
    close(predicted, 7.444444444444445, 1e-6, "odd_list_mean")


def test_1_prediction_even_list_median():
    predicted = need(answers.ANSWERS["even_list_median"], "even_list_median prediction")
    close(predicted, 6.0, 1e-9, "even_list_median")


# ---------------------------------------------------------------------------
# Exercise 2 -- the breakdown point
# ---------------------------------------------------------------------------


def test_2_mean_moves_far():
    before, after = attempt(
        lambda: F.breakdown_point_mean(D.SALARY_LIST, D.CORRUPTED_SALARY),
        "breakdown_point_mean",
    )
    assert after - before > D.BREAKDOWN_MEAN_SHIFT_FLOOR


def test_2_median_does_not_move_at_all():
    before, after = attempt(
        lambda: F.breakdown_point_median(D.SALARY_LIST, D.CORRUPTED_SALARY),
        "breakdown_point_median",
    )
    assert after == before


def test_2_prediction_mean_breaks_down():
    predicted = need(answers.ANSWERS["mean_breaks_down"], "mean_breaks_down prediction")
    assert predicted is True


def test_2_prediction_median_breaks_down():
    predicted = need(answers.ANSWERS["median_breaks_down"], "median_breaks_down prediction")
    assert predicted is False


# ---------------------------------------------------------------------------
# Exercise 3 -- Bessel's correction
# ---------------------------------------------------------------------------


def test_3_divide_by_n_is_biased_low():
    rng = np.random.default_rng(D.BESSEL_SEED)
    biased, _ = attempt(
        lambda: S.bessel_trial_variances(
            rng,
            D.BESSEL_POPULATION_MEAN,
            D.BESSEL_POPULATION_SIGMA,
            D.BESSEL_SAMPLE_SIZE,
            D.BESSEL_TRIALS,
        ),
        "bessel_trial_variances",
    )
    ratio = float(biased.mean()) / D.BESSEL_TRUE_VARIANCE
    assert ratio < 1.0, "dividing by n should UNDERestimate the true variance"
    assert abs(ratio - D.BESSEL_EXPECTED_BIAS_FACTOR) < D.BESSEL_BIAS_FACTOR_TOLERANCE


def test_3_divide_by_n_minus_1_is_unbiased():
    rng = np.random.default_rng(D.BESSEL_SEED)
    _, unbiased = attempt(
        lambda: S.bessel_trial_variances(
            rng,
            D.BESSEL_POPULATION_MEAN,
            D.BESSEL_POPULATION_SIGMA,
            D.BESSEL_SAMPLE_SIZE,
            D.BESSEL_TRIALS,
        ),
        "bessel_trial_variances",
    )
    mean_unbiased = float(unbiased.mean())
    se = float(unbiased.std(ddof=1)) / (D.BESSEL_TRIALS**0.5)
    assert abs(mean_unbiased - D.BESSEL_TRUE_VARIANCE) < D.BESSEL_UNBIASED_SE_TOLERANCE * se


def test_3_prediction_bias_direction():
    predicted = need(
        answers.ANSWERS["divide_by_n_bias_direction"], "divide_by_n_bias_direction prediction"
    )
    assert predicted == "low"


# ---------------------------------------------------------------------------
# Exercise 4 -- percentile ambiguity
# ---------------------------------------------------------------------------


def test_4_conventions_disagree():
    values = {
        method: attempt(
            lambda method=method: F.percentile_under(
                D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, method
            ),
            "percentile_under",
        )
        for method in D.PERCENTILE_METHODS
    }
    assert len(set(values.values())) >= 2


def test_4_default_linear_matches_documented_value():
    got = attempt(
        lambda: F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, "linear"),
        "percentile_under",
    )
    assert got == 8.25


def test_4_prediction_disagree():
    predicted = need(
        answers.ANSWERS["percentile_methods_disagree"], "percentile_methods_disagree prediction"
    )
    assert predicted is True


def test_4_prediction_linear_value():
    predicted = need(
        answers.ANSWERS["percentile_linear_value"], "percentile_linear_value prediction"
    )
    close(predicted, 8.25, 1e-9, "percentile_linear_value")


# ---------------------------------------------------------------------------
# Exercise 5 -- Pearson versus Spearman
# ---------------------------------------------------------------------------


def test_5_pearson_on_parabola_is_essentially_zero():
    got = attempt(lambda: F.pearson(D.PARABOLA_X, D.PARABOLA_Y), "pearson")
    assert abs(got) < D.PARABOLA_PEARSON_TOLERANCE


def test_5_spearman_on_monotone_cubic_is_exactly_one():
    got = attempt(lambda: F.spearman(D.MONOTONE_X, D.MONOTONE_Y), "spearman")
    assert got == 1.0


def test_5_prediction_parabola_pearson():
    predicted = need(
        answers.ANSWERS["parabola_pearson_magnitude"], "parabola_pearson_magnitude prediction"
    )
    assert predicted == "close_to_zero"


def test_5_prediction_monotone_spearman():
    predicted = need(
        answers.ANSWERS["monotone_spearman_value"], "monotone_spearman_value prediction"
    )
    close(predicted, 1.0, 1e-9, "monotone_spearman_value")


# ---------------------------------------------------------------------------
# Exercise 6 -- Anscombe's quartet
# ---------------------------------------------------------------------------


def test_6_all_four_sets_agree_on_summaries():
    reference = attempt(
        lambda: F.anscombe_summary(*D.ANSCOMBE_SETS["I"]), "anscombe_summary"
    )
    for name in ("II", "III", "IV"):
        s = attempt(lambda name=name: F.anscombe_summary(*D.ANSCOMBE_SETS[name]), "anscombe_summary")
        assert round(s["mean_x"], 1) == round(reference["mean_x"], 1)
        assert round(s["correlation"], 1) == round(reference["correlation"], 1)


def test_6_shape_statistics_separate_set_iv():
    shapes = {
        name: attempt(lambda name=name: F.shape_statistics(*D.ANSCOMBE_SETS[name]), "shape_statistics")
        for name in D.ANSCOMBE_SETS
    }
    assert shapes["IV"]["max_leverage"] > 3.0 * shapes["I"]["max_leverage"]


def test_6_prediction_means_agree():
    predicted = need(answers.ANSWERS["anscombe_means_agree"], "anscombe_means_agree prediction")
    assert predicted is True


def test_6_prediction_set_iv_leverage():
    predicted = need(
        answers.ANSWERS["anscombe_set_iv_leverage_dominant"],
        "anscombe_set_iv_leverage_dominant prediction",
    )
    assert predicted is True


# ---------------------------------------------------------------------------
# Exercise 7 -- Simpson's paradox
# ---------------------------------------------------------------------------


def test_7_a_wins_every_subgroup():
    easy_a = attempt(lambda: F.success_rate(*D.TREATMENT_A_EASY), "success_rate")
    easy_b = attempt(lambda: F.success_rate(*D.TREATMENT_B_EASY), "success_rate")
    hard_a = attempt(lambda: F.success_rate(*D.TREATMENT_A_HARD), "success_rate")
    hard_b = attempt(lambda: F.success_rate(*D.TREATMENT_B_HARD), "success_rate")
    assert easy_a > easy_b
    assert hard_a > hard_b


def test_7_b_wins_overall():
    a_total = attempt(
        lambda: F.combined_rate(D.TREATMENT_A_EASY, D.TREATMENT_A_HARD), "combined_rate"
    )
    b_total = attempt(
        lambda: F.combined_rate(D.TREATMENT_B_EASY, D.TREATMENT_B_HARD), "combined_rate"
    )
    assert b_total > a_total


def test_7_prediction_b_wins_overall():
    predicted = need(answers.ANSWERS["simpson_b_wins_overall"], "simpson_b_wins_overall prediction")
    assert predicted is True


# ---------------------------------------------------------------------------
# Exercise 8 -- robust spread under contamination
# ---------------------------------------------------------------------------


def test_8_std_inflates_a_lot():
    rng = np.random.default_rng(D.CONTAMINATION_SEED)
    clean, contaminated = attempt(
        lambda: S.contaminated_sample(
            rng,
            D.CONTAMINATION_BASE_MEAN,
            D.CONTAMINATION_BASE_SIGMA,
            D.CONTAMINATION_BASE_N,
            D.CONTAMINATION_OUTLIERS,
        ),
        "contaminated_sample",
    )
    std_clean = float(np.std(clean, ddof=1))
    std_contam = float(np.std(contaminated, ddof=1))
    assert std_contam / std_clean > D.CONTAMINATION_STD_MULTIPLIER_FLOOR


def test_8_mad_barely_moves():
    rng = np.random.default_rng(D.CONTAMINATION_SEED)
    clean, contaminated = attempt(
        lambda: S.contaminated_sample(
            rng,
            D.CONTAMINATION_BASE_MEAN,
            D.CONTAMINATION_BASE_SIGMA,
            D.CONTAMINATION_BASE_N,
            D.CONTAMINATION_OUTLIERS,
        ),
        "contaminated_sample",
    )
    mad_clean = attempt(lambda: F.median_absolute_deviation(clean), "median_absolute_deviation")
    mad_contam = attempt(
        lambda: F.median_absolute_deviation(contaminated), "median_absolute_deviation"
    )
    assert mad_contam / mad_clean < D.CONTAMINATION_MAD_MULTIPLIER_CEILING


def test_8_prediction_std_inflates():
    predicted = need(
        answers.ANSWERS["contamination_inflates_std"], "contamination_inflates_std prediction"
    )
    assert predicted is True


def test_8_prediction_mad_stable():
    predicted = need(
        answers.ANSWERS["contamination_mad_stable"], "contamination_mad_stable prediction"
    )
    assert predicted is True


# ---------------------------------------------------------------------------
# Exercise 9 -- standardisation
# ---------------------------------------------------------------------------


def test_9_standardized_mean_and_std():
    rng = np.random.default_rng(D.STANDARDIZATION_SEED)
    x = rng.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
    zx = attempt(lambda: F.zscores(x), "zscores")
    mean_zx = sum(zx) / len(zx)
    std_zx = (sum((v - mean_zx) ** 2 for v in zx) / len(zx)) ** 0.5
    assert abs(mean_zx) < D.STANDARDIZATION_MEAN_TOLERANCE
    assert abs(std_zx - 1.0) < D.STANDARDIZATION_STD_TOLERANCE


def test_9_correlation_unchanged_by_standardizing():
    rng = np.random.default_rng(D.STANDARDIZATION_SEED)
    x = rng.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
    noise = rng.normal(0.0, D.STANDARDIZATION_Y_NOISE_SIGMA, D.STANDARDIZATION_N)
    y = D.STANDARDIZATION_Y_SLOPE * x + noise
    zx = attempt(lambda: F.zscores(x), "zscores")
    zy = attempt(lambda: F.zscores(y), "zscores")
    r_original = attempt(lambda: F.pearson(x, y), "pearson")
    r_standardized = attempt(lambda: F.pearson(zx, zy), "pearson")
    assert abs(r_original - r_standardized) < D.STANDARDIZATION_CORRELATION_TOLERANCE


def test_9_prediction_mean():
    predicted = need(answers.ANSWERS["standardized_mean"], "standardized_mean prediction")
    close(predicted, 0.0, 1e-6, "standardized_mean")


def test_9_prediction_std():
    predicted = need(answers.ANSWERS["standardized_std"], "standardized_std prediction")
    close(predicted, 1.0, 1e-6, "standardized_std")


def test_9_prediction_correlation_unchanged():
    predicted = need(
        answers.ANSWERS["standardizing_changes_correlation"],
        "standardizing_changes_correlation prediction",
    )
    assert predicted is False

#!/usr/bin/env bash
# Tests for the Day 116 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:
#
#   * mean, median and mode from scratch match the `statistics` module;
#   * one corrupted salary out of nine drags the mean by over a million
#     dollars and does not move the median at all -- exact equality;
#   * dividing by n underestimates the true variance on average by exactly
#     the factor (n-1)/n; dividing by n-1 lands within a few standard
#     errors of the truth, measured by simulation over 20,000 trials;
#   * NumPy's percentile conventions genuinely disagree about the 75th
#     percentile of the same eight numbers;
#   * Pearson on a perfect parabola is essentially zero; Spearman on a
#     monotone cubic is exactly 1.0;
#   * Anscombe's published 1973 quartet agrees on every classic summary
#     statistic to the documented precision, and three diagnostics those
#     summaries cannot see tell the four sets apart;
#   * treatment A beats treatment B in every subgroup of the smallest
#     table that shows Simpson's paradox, and treatment B wins overall --
#     both directions, from the same table;
#   * 3% contamination inflates the standard deviation by more than 5x and
#     moves the median absolute deviation by less than 1.5x;
#   * standardising gives mean 0 and standard deviation 1, and leaves the
#     Pearson correlation between two variables unchanged;
#   * 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 116 — Statistics That Don't Lie"
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_mean_median_mode 02_breakdown_point 03_bessel_correction \
              04_percentile_ambiguity 05_pearson_vs_spearman 06_anscombes_quartet \
              07_simpsons_paradox 08_robust_spread_under_contamination \
              09_standardization; 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 25 ]; then
  check "the reference suite ran at least 25 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 25 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 `dataset`,
# `descriptive`, `simulate` and `answers`, 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 lesson's claims, checked one value at a time"
# --------------------------------------------------------------------------

facts="$(cd "${lab_dir}/examples" && "${python_bin}" - <<'PY'
import statistics as st

import numpy as np

import dataset as D
import descriptive as F
import simulate as S

print("odd_list_mean", F.mean(D.ODD_LIST))
print("odd_list_median", F.median(D.ODD_LIST))
print("multimodal_modes", F.modes(D.MULTIMODAL_LIST))

mean_before, mean_after = F.breakdown_point_mean(D.SALARY_LIST, D.CORRUPTED_SALARY)
median_before, median_after = F.breakdown_point_median(D.SALARY_LIST, D.CORRUPTED_SALARY)
print("mean_shift_large", (mean_after - mean_before) > D.BREAKDOWN_MEAN_SHIFT_FLOOR)
print("median_shift_zero", median_after == median_before)

rng = np.random.default_rng(D.BESSEL_SEED)
biased, unbiased = S.bessel_trial_variances(
    rng, D.BESSEL_POPULATION_MEAN, D.BESSEL_POPULATION_SIGMA,
    D.BESSEL_SAMPLE_SIZE, D.BESSEL_TRIALS,
)
ratio_biased = float(biased.mean()) / D.BESSEL_TRUE_VARIANCE
print("bessel_biased_ratio", round(ratio_biased, 4))
print("bessel_biased_near_n_minus_1_over_n", abs(ratio_biased - D.BESSEL_EXPECTED_BIAS_FACTOR) < D.BESSEL_BIAS_FACTOR_TOLERANCE)
se_unbiased = float(unbiased.std(ddof=1)) / (D.BESSEL_TRIALS ** 0.5)
print("bessel_unbiased_within_tolerance", abs(float(unbiased.mean()) - D.BESSEL_TRUE_VARIANCE) < D.BESSEL_UNBIASED_SE_TOLERANCE * se_unbiased)

pvals = {m: F.percentile_under(D.PERCENTILE_ARRAY, D.PERCENTILE_TARGET, m) for m in D.PERCENTILE_METHODS}
print("percentile_distinct_count", len(set(pvals.values())))
print("percentile_linear", pvals["linear"])
print("percentile_lower", pvals["lower"])
print("percentile_higher", pvals["higher"])

pear_parabola = F.pearson(D.PARABOLA_X, D.PARABOLA_Y)
print("parabola_pearson_near_zero", abs(pear_parabola) < D.PARABOLA_PEARSON_TOLERANCE)
spear_monotone = F.spearman(D.MONOTONE_X, D.MONOTONE_Y)
print("monotone_spearman_is_one", spear_monotone == 1.0)

anscombe_ref = F.anscombe_summary(*D.ANSCOMBE_SETS["I"])
agree = all(
    round(F.anscombe_summary(*D.ANSCOMBE_SETS[name])["mean_x"], 1) == round(anscombe_ref["mean_x"], 1)
    and round(F.anscombe_summary(*D.ANSCOMBE_SETS[name])["correlation"], 1) == round(anscombe_ref["correlation"], 1)
    for name in D.ANSCOMBE_SETS
)
print("anscombe_summaries_agree", agree)
shapes = {name: F.shape_statistics(*D.ANSCOMBE_SETS[name]) for name in D.ANSCOMBE_SETS}
print("anscombe_set_iv_leverage_dominant", shapes["IV"]["max_leverage"] > 3.0 * shapes["I"]["max_leverage"])
print("anscombe_set_iii_outlier_dominant", shapes["III"]["outlier_ratio"] > 2.0 * shapes["I"]["outlier_ratio"])

a_easy = F.success_rate(*D.TREATMENT_A_EASY)
a_hard = F.success_rate(*D.TREATMENT_A_HARD)
b_easy = F.success_rate(*D.TREATMENT_B_EASY)
b_hard = F.success_rate(*D.TREATMENT_B_HARD)
a_total = F.combined_rate(D.TREATMENT_A_EASY, D.TREATMENT_A_HARD)
b_total = F.combined_rate(D.TREATMENT_B_EASY, D.TREATMENT_B_HARD)
print("simpson_a_wins_both_subgroups", a_easy > b_easy and a_hard > b_hard)
print("simpson_b_wins_overall", b_total > a_total)

rng2 = np.random.default_rng(D.CONTAMINATION_SEED)
clean, contaminated = S.contaminated_sample(
    rng2, D.CONTAMINATION_BASE_MEAN, D.CONTAMINATION_BASE_SIGMA,
    D.CONTAMINATION_BASE_N, D.CONTAMINATION_OUTLIERS,
)
std_clean = float(np.std(clean, ddof=1))
std_contam = float(np.std(contaminated, ddof=1))
mad_clean = F.median_absolute_deviation(clean)
mad_contam = F.median_absolute_deviation(contaminated)
print("contamination_std_multiplier", round(std_contam / std_clean, 2))
print("contamination_mad_multiplier", round(mad_contam / mad_clean, 2))
print("contamination_std_inflates", (std_contam / std_clean) > D.CONTAMINATION_STD_MULTIPLIER_FLOOR)
print("contamination_mad_stable", (mad_contam / mad_clean) < D.CONTAMINATION_MAD_MULTIPLIER_CEILING)

rng3 = np.random.default_rng(D.STANDARDIZATION_SEED)
x = rng3.normal(D.STANDARDIZATION_X_MEAN, D.STANDARDIZATION_X_SIGMA, D.STANDARDIZATION_N)
noise = rng3.normal(0.0, D.STANDARDIZATION_Y_NOISE_SIGMA, D.STANDARDIZATION_N)
y = D.STANDARDIZATION_Y_SLOPE * x + noise
zx, zy = F.zscores(x), F.zscores(y)
mean_zx = sum(zx) / len(zx)
std_zx = (sum((v - mean_zx) ** 2 for v in zx) / len(zx)) ** 0.5
print("standardized_mean_near_zero", abs(mean_zx) < D.STANDARDIZATION_MEAN_TOLERANCE)
print("standardized_std_near_one", abs(std_zx - 1.0) < D.STANDARDIZATION_STD_TOLERANCE)
r_orig = F.pearson(x, y)
r_std = F.pearson(zx, zy)
print("standardizing_preserves_correlation", abs(r_orig - r_std) < D.STANDARDIZATION_CORRELATION_TOLERANCE)
PY
)"

get() { printf '%s\n' "${facts}" | grep "^$1 " | cut -d' ' -f2-; }

check_eq "the odd list's mean matches the worked figure" "7.444444444444445" "$(get odd_list_mean)"
check_eq "the odd list's median is 7.0" "7.0" "$(get odd_list_median)"
check_eq "the multimodal list has modes [3, 8]" "[3, 8]" "$(get multimodal_modes)"
check_eq "one corrupted salary drags the mean by over the stated floor" "True" "$(get mean_shift_large)"
check_eq "the same corruption leaves the median exactly unchanged" "True" "$(get median_shift_zero)"
check_eq "dividing by n is biased low, near (n-1)/n" "True" "$(get bessel_biased_near_n_minus_1_over_n)"
check_eq "dividing by n-1 lands within tolerance of the true variance" "True" "$(get bessel_unbiased_within_tolerance)"
check_eq "at least 2 percentile conventions disagree" "True" "$([ "$(get percentile_distinct_count)" -ge 2 ] && echo True || echo False)"
check_eq "the default ('linear') 75th percentile is 8.25" "8.25" "$(get percentile_linear)"
check_eq "'lower' and 'higher' land on different real data points" "True" "$([ "$(get percentile_lower)" != "$(get percentile_higher)" ] && echo True || echo False)"
check_eq "Pearson on the symmetric parabola is essentially zero" "True" "$(get parabola_pearson_near_zero)"
check_eq "Spearman on the monotone cubic is exactly 1.0" "True" "$(get monotone_spearman_is_one)"
check_eq "all four Anscombe sets agree on the classic summaries" "True" "$(get anscombe_summaries_agree)"
check_eq "set IV's leverage dramatically dominates set I's" "True" "$(get anscombe_set_iv_leverage_dominant)"
check_eq "set III's outlier residual dramatically dominates set I's" "True" "$(get anscombe_set_iii_outlier_dominant)"
check_eq "treatment A wins both Simpson's-paradox subgroups" "True" "$(get simpson_a_wins_both_subgroups)"
check_eq "treatment B still wins overall" "True" "$(get simpson_b_wins_overall)"
check_eq "3% contamination inflates the standard deviation past the floor" "True" "$(get contamination_std_inflates)"
check_eq "the same contamination leaves the MAD under the ceiling" "True" "$(get contamination_mad_stable)"
check_eq "the standardized sample's mean is (numerically) zero" "True" "$(get standardized_mean_near_zero)"
check_eq "the standardized sample's standard deviation is (numerically) one" "True" "$(get standardized_std_near_one)"
check_eq "standardising leaves the Pearson correlation unchanged" "True" "$(get standardizing_preserves_correlation)"

echo "  contamination std multiplier measured at: $(get contamination_std_multiplier)x"
echo "  contamination MAD multiplier measured at: $(get contamination_mad_multiplier)x"

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

# A green test suite proves nothing until you have watched it go red. This
# section re-runs the reference pytest suite with one assertion deliberately
# broken -- the breakdown-point median check flipped from "equal" to
# "not equal" -- and asserts that the run reports the failure and exits
# non-zero.
sentinel_file="${lab_dir}/examples/test_reference.py"
if grep -q "def test_median_breakdown_point_does_not_move_at_all" "${sentinel_file}"; then
  backup="$(mktemp)"
  cp "${sentinel_file}" "${backup}"
  python3 - "${sentinel_file}" <<'PY'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert after == before  # exact equality: the median's rank did not change"
replacement = "assert after != before  # DELIBERATELY WRONG: self-test only"
assert needle in text, "sentinel assertion not found"
open(path, "w").write(text.replace(needle, replacement, 1))
PY
  self_out="$(cd "${lab_dir}" && "${pytest_bin}" examples -q -p no:cacheprovider 2>&1)"
  self_status=$?
  cp "${backup}" "${sentinel_file}"
  rm -f "${backup}"
  find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -exec rm -rf {} + 2>/dev/null || true

  if [ "${self_status}" -ne 0 ]; then
    check "a deliberately broken assertion makes the suite exit non-zero (${self_status})" "yes"
  else
    check "a deliberately broken assertion makes the suite exit non-zero" "no"
  fi
  case "${self_out}" in
    *"test_median_breakdown_point_does_not_move_at_all"*"failed"*|*"FAILED"*"test_median_breakdown_point_does_not_move_at_all"*)
      check "the failing test is named in the output" "yes" ;;
    *) check "the failing test is named in the output" "no" ;;
  esac
  case "${self_out}" in
    *"1 failed"*) check "exactly one test failed" "yes" ;;
    *) check "exactly one test failed" "no" ;;
  esac
else
  check "sentinel assertion found in examples/test_reference.py" "no"
fi

# --------------------------------------------------------------------------
echo
echo "7. 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.

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 ]