Hands-on lab — Day 115: Bayes’ Theorem

"""Exercise 1 -- the opening failure, derived exactly.

A test is 99% sensitive and 99% specific. The condition affects 1 person in
1,000. You test positive. Almost everyone -- including, in published
studies, most physicians asked this exact question -- answers "about 99%".
The true answer is close to 9%. This script derives it exactly with
Bayes' theorem, as a Fraction, and asserts both the true value and the
misconception it corrects.
"""

from fractions import Fraction

import bayes as B
import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("The scenario")
print("-" * 60)
print(f"  prevalence   P(condition)          = {D.PREVALENCE}")
print(f"  sensitivity  P(positive|condition) = {D.SENSITIVITY}")
print(f"  specificity  P(negative|no cond.)  = {D.SPECIFICITY}")

print()
print("Bayes' theorem, worked term by term")
print("-" * 60)

p_condition_and_positive = D.PREVALENCE * D.SENSITIVITY
p_no_condition_and_positive = (1 - D.PREVALENCE) * (1 - D.SPECIFICITY)
evidence = p_condition_and_positive + p_no_condition_and_positive
result = p_condition_and_positive / evidence

print(f"  P(condition and positive)    = {D.PREVALENCE} x {D.SENSITIVITY} = {p_condition_and_positive}")
print(f"  P(no condition and positive) = {1 - D.PREVALENCE} x {1 - D.SPECIFICITY} = {p_no_condition_and_positive}")
print(f"  P(positive)  [the evidence]  = {p_condition_and_positive} + {p_no_condition_and_positive} = {evidence}")
print(f"  P(condition | positive)      = {p_condition_and_positive} / {evidence} = {result}")
print(f"                                = {float(result):.6f}  ~ {round(float(result), 4)}")

via_function = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
check("posterior() agrees with the term-by-term derivation", via_function == result)
check("the exact posterior is 99/1098", result == Fraction(99, 1098))
check("which reduces to 11/122", result == Fraction(11, 122))
check("and rounds to 0.0902", round(float(result), 4) == 0.0902)
check("the posterior is NOT 0.99 -- the answer almost everyone gives", result != Fraction(99, 100))
check("in fact it is roughly 11x smaller than the naive 0.99 guess", Fraction(99, 100) / result > 10)

print()
print("Ninety-one out of every hundred positives are false alarms")
print("-" * 60)
false_alarm_share = 1 - result
print(f"  P(no condition | positive) = 1 - {result} = {false_alarm_share} ~ {round(float(false_alarm_share), 4)}")
check("over 90% of positives are false alarms", float(false_alarm_share) > 0.90)

print()
if all(ok for _, ok in checks_held):
    print(f"01_opening_posterior.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 2 -- the same arithmetic, made obvious by counting people
instead of multiplying percentages.

Stop using percentages. Count 100,000 people. This script builds the exact
table and shows that the fraction TP / (TP + FP) -- read straight off four
integers -- equals exercise 1's formula answer exactly.
"""

from fractions import Fraction

import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print(f"Out of {D.NATURAL_FREQUENCY_POPULATION:,} people")
print("-" * 60)
print(f"  {D.NATURAL_FREQUENCY_SICK} have the condition, {D.NATURAL_FREQUENCY_WELL:,} do not")
check(
    "sick + well accounts for everyone",
    D.NATURAL_FREQUENCY_SICK + D.NATURAL_FREQUENCY_WELL == D.NATURAL_FREQUENCY_POPULATION,
)

print()
print("Of the sick people:")
print(f"  {D.NATURAL_FREQUENCY_TP} test positive (true positives)")
print(f"  {D.NATURAL_FREQUENCY_FN} tests negative (false negative)")
check(
    "TP + FN accounts for every sick person",
    D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FN == D.NATURAL_FREQUENCY_SICK,
)

print()
print("Of the healthy people:")
print(f"  {D.NATURAL_FREQUENCY_FP} test positive anyway (false positives)")
print(f"  {D.NATURAL_FREQUENCY_TN:,} correctly test negative (true negatives)")
check(
    "FP + TN accounts for every healthy person",
    D.NATURAL_FREQUENCY_FP + D.NATURAL_FREQUENCY_TN == D.NATURAL_FREQUENCY_WELL,
)

print()
print("So, of everyone who tests positive:")
print("-" * 60)
total_positive = D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FP
print(f"  true positives  {D.NATURAL_FREQUENCY_TP}")
print(f"  false positives {D.NATURAL_FREQUENCY_FP}")
print(f"  total positives {total_positive}")
natural_posterior = Fraction(D.NATURAL_FREQUENCY_TP, total_positive)
print(f"  P(condition | positive) = {D.NATURAL_FREQUENCY_TP} / {total_positive} = {natural_posterior}"
      f" ~ {round(float(natural_posterior), 4)}")

check("total positives is exactly 1,098", total_positive == 1098)
check("TP / (TP + FP) equals the exact posterior from exercise 1", natural_posterior == D.OPENING_POSTERIOR_EXACT)
check("nine false alarms for every true positive, roughly", D.NATURAL_FREQUENCY_FP / D.NATURAL_FREQUENCY_TP > 9)

print()
print("The percentages were never wrong -- they were just opaque")
print("-" * 60)
print("  1% of 999,000... no: 1% of the 99,900 healthy people IS 999 people,")
print("  and 999 false positives dwarfs the 99 true positives a 99%-sensitive")
print("  test catches out of only 100 sick people. Counting makes that")
print("  imbalance impossible to miss; percentages hid it in plain sight.")

print()
if all(ok for _, ok in checks_held):
    print(f"02_natural_frequencies.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 3 -- confirm the exact posterior by simulating a large
population and literally counting the four outcome cells.

Neither the exact formula nor the natural-frequencies table above touched
numpy.random. This script is the third, independent check: draw a
population, test it, and count -- and the empirical posterior should land
within a few standard errors of the exact 99/1098.
"""

import numpy as np

import dataset as D
import simulate as S

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print(f"Simulating {D.SIMULATION_POPULATION:,} people, seed {D.SIMULATION_SEED}")
print("-" * 60)

rng = np.random.default_rng(D.SIMULATION_SEED)
counts = S.simulate_population(
    rng,
    D.SIMULATION_POPULATION,
    float(D.PREVALENCE),
    float(D.SENSITIVITY),
    float(D.SPECIFICITY),
)

print(f"  true positives:  {counts.true_positive:,}")
print(f"  false positives: {counts.false_positive:,}")
print(f"  true negatives:  {counts.true_negative:,}")
print(f"  false negatives: {counts.false_negative:,}")
print(f"  total positives: {counts.positives:,}")

exact = float(D.OPENING_POSTERIOR_EXACT)
empirical = counts.empirical_posterior
se = D.standard_error(exact, counts.positives)
tolerance = 3.0 * se

print()
print(f"  exact posterior      = {exact:.6f}")
print(f"  simulated posterior  = {empirical:.6f}")
print(f"  gap                  = {abs(empirical - exact):.6f}")
print(f"  standard error at n={counts.positives:,} positives = {se:.6f}")
print(f"  tolerance (3 SE)     = {tolerance:.6f}")

check("some people actually have the condition in the sample", counts.true_positive + counts.false_negative > 0)
check("the simulated posterior lands within 3 standard errors of the exact value",
      abs(empirical - exact) < tolerance)
check("the simulated posterior is nowhere near the naive 0.99 guess", empirical < 0.5)

print()
if all(ok for _, ok in checks_held):
    print(f"03_simulation.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 4 -- the base rate decides the answer. Sweep it and watch.

Everything about sensitivity and specificity stays fixed at 99%/99%
throughout this script. Only the prevalence changes -- and the posterior
climbs from under a tenth of a percent at one case in 100,000, to exactly
0.99 at a prevalence of one-half. 0.99 is the number nearly everyone
wrongly gives for the 1-in-1,000 case; this script shows it is the RIGHT
answer, just for a completely different question.
"""

from fractions import Fraction

import bayes as B
import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("Sensitivity and specificity fixed at 99% / 99%. Prevalence varies.")
print("-" * 60)

results: list[Fraction] = []
for prevalence in D.PREVALENCE_SWEEP:
    posterior = B.posterior(prevalence, D.SENSITIVITY, D.SPECIFICITY)
    results.append(posterior)
    print(f"  prevalence {str(prevalence):>8} -> posterior {str(posterior):>10}"
          f"  ({float(posterior):.6f})")

print()
check(
    "the posterior is strictly increasing as prevalence rises",
    all(results[i] < results[i + 1] for i in range(len(results) - 1)),
)

last_prevalence, last_posterior = D.PREVALENCE_SWEEP[-1], results[-1]
check("the sweep's final prevalence is exactly 1/2", last_prevalence == Fraction(1, 2))
check("at prevalence 1/2 the posterior is EXACTLY 0.99", last_posterior == Fraction(99, 100))

opening_posterior = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
check(
    "0.99 -- everyone's wrong guess for the 1-in-1,000 case -- is over 10x the true answer there",
    Fraction(99, 100) / opening_posterior > 10,
)
check(
    "...but is the EXACT right answer once the base rate is 1 in 2",
    last_posterior == Fraction(99, 100),
)

print()
print("The number 0.99 was never wrong. It just answered a different question.")
print("-" * 60)
print("  Sensitivity and specificity describe the TEST. The posterior")
print("  describes the PATIENT, and the patient's answer depends on how")
print("  common the condition was before the test ever ran.")

print()
if all(ok for _, ok in checks_held):
    print(f"04_prevalence_sweep.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 5 -- the odds form: posterior odds = prior odds x likelihood
ratio.

This is the cleanest way to think about updating a belief, because the
likelihood ratio isolates exactly how much a piece of evidence is worth,
completely independent of what you believed before you saw it.
"""

from fractions import Fraction

import bayes as B
import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("Prior odds")
print("-" * 60)
prior_odds = B.probability_to_odds(D.PREVALENCE)
print(f"  P(condition) = {D.PREVALENCE}")
print(f"  prior odds   = {D.PREVALENCE} / (1 - {D.PREVALENCE}) = {prior_odds}  (about 1 in 999)")

print()
print("Likelihood ratio")
print("-" * 60)
ratio = B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY)
print(f"  LR+ = P(positive|condition) / P(positive|no condition)")
print(f"      = {D.SENSITIVITY} / (1 - {D.SPECIFICITY}) = {ratio}")
print(f"  a positive result is {ratio}x more likely under the condition than without it")

print()
print("Posterior odds, and back to a probability")
print("-" * 60)
posterior_odds = B.update_odds(prior_odds, ratio)
posterior_probability = B.odds_to_probability(posterior_odds)
print(f"  posterior odds = {prior_odds} x {ratio} = {posterior_odds}")
print(f"  posterior probability = {posterior_odds} / (1 + {posterior_odds}) = {posterior_probability}"
      f"  ({float(posterior_probability):.6f})")

direct = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
check("posterior odds equal prior odds times the likelihood ratio, exactly", posterior_odds == prior_odds * ratio)
check("converting back to a probability matches exercise 1's direct answer exactly", posterior_probability == direct)
check("the likelihood ratio here is exactly 99", ratio == 99)
check("prior odds are exactly 1/999", prior_odds == Fraction(1, 999))

print()
print("Why the odds form earns its place")
print("-" * 60)
print("  The likelihood ratio (99, here) says exactly how much a positive")
print("  result is worth as evidence, and that number does not change if")
print("  the prior changes. Multiply it onto ANY prior odds and you get")
print("  that prior's correctly updated posterior odds -- which is exactly")
print("  what exercise 6 does twice in a row.")

print()
if all(ok for _, ok in checks_held):
    print(f"05_odds_form.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 6 -- two positive tests, updated one at a time, in both orders.

Test A: 99% sensitive, 99% specific -- the opening scenario's test.
Test B: 95% sensitive, 98% specific -- a genuinely different, less accurate
test. Both come back positive. The order they are applied in does not
matter, and this script proves it by actually computing both orders and
comparing the results exactly, rather than asserting it from the algebra
alone.
"""

from fractions import Fraction

import bayes as B
import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


test_a = (D.TEST_A_SENSITIVITY, D.TEST_A_SPECIFICITY)
test_b = (D.TEST_B_SENSITIVITY, D.TEST_B_SPECIFICITY)

print("Two different tests, both positive")
print("-" * 60)
print(f"  test A: sensitivity {test_a[0]}, specificity {test_a[1]}")
print(f"  test B: sensitivity {test_b[0]}, specificity {test_b[1]}")
print(f"  prior:  P(condition) = {D.PREVALENCE}")

print()
print("Update with A first, then B")
print("-" * 60)
odds_after_a = B.update_odds(B.probability_to_odds(D.PREVALENCE), B.likelihood_ratio(*test_a))
posterior_after_a = B.odds_to_probability(odds_after_a)
print(f"  after A: posterior = {posterior_after_a}  ({float(posterior_after_a):.6f})")
odds_after_ab = B.update_odds(odds_after_a, B.likelihood_ratio(*test_b))
posterior_a_then_b = B.odds_to_probability(odds_after_ab)
print(f"  after A then B: posterior = {posterior_a_then_b}  ({float(posterior_a_then_b):.6f})")

print()
print("Update with B first, then A")
print("-" * 60)
odds_after_b = B.update_odds(B.probability_to_odds(D.PREVALENCE), B.likelihood_ratio(*test_b))
posterior_after_b = B.odds_to_probability(odds_after_b)
print(f"  after B: posterior = {posterior_after_b}  ({float(posterior_after_b):.6f})")
odds_after_ba = B.update_odds(odds_after_b, B.likelihood_ratio(*test_a))
posterior_b_then_a = B.odds_to_probability(odds_after_ba)
print(f"  after B then A: posterior = {posterior_b_then_a}  ({float(posterior_b_then_a):.6f})")

print()
print("Via sequential_posterior(), both orders in one call each")
print("-" * 60)
via_ab = B.sequential_posterior(D.PREVALENCE, [test_a, test_b])
via_ba = B.sequential_posterior(D.PREVALENCE, [test_b, test_a])
print(f"  [A, B] -> {via_ab}   [B, A] -> {via_ba}")

check("A-then-B matches B-then-A exactly", posterior_a_then_b == posterior_b_then_a)
check("both hand-worked orders match sequential_posterior()'s [A, B] result", posterior_a_then_b == via_ab)
check("both hand-worked orders match sequential_posterior()'s [B, A] result", posterior_b_then_a == via_ba)
check("the two-test posterior is exactly 1045/1267", via_ab == Fraction(1045, 1267))
check("two positive tests move the posterior well past one-half", via_ab > Fraction(1, 2))
check(
    "a single test's posterior (~9%) is far below the two-test posterior (~82%)",
    posterior_after_a < via_ab,
)

print()
print("Why the order never matters")
print("-" * 60)
print("  Each update multiplies the running odds by one more likelihood")
print("  ratio. Multiplication of real (or Fraction) numbers is")
print("  commutative -- a x b always equals b x a -- so a chain of")
print("  updates is the same product regardless of the order the factors")
print("  arrive in. This is not a special property of Bayes' theorem; it")
print("  is a special property of multiplication that Bayes' theorem")
print("  happens to be built from.")

print()
if all(ok for _, ok in checks_held):
    print(f"06_sequential_updating.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 7 -- the honest caveat: multiplying likelihood ratios twice
assumes the two runs are conditionally independent given the hypothesis,
and that assumption can be badly false.

Same test (99% sensitive, 99% specific), run twice on ONE sample. Half the
time (modelled here as correlation_weight = 1/2), both runs are yoked to a
single shared failure mode -- a contaminated sample, a bad reagent batch --
so they either BOTH come back positive or BOTH come back negative,
regardless of the true condition. This script computes the posterior two
ways: the naive way, which assumes full independence, and the correct way,
which accounts for the correlation -- and the naive answer is not just
different, it is dramatically more confident than it has any right to be.
"""

from fractions import Fraction

import bayes as B
import dataset as D

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("Same test, same sample, run twice. Both runs come back positive.")
print("-" * 60)
print(f"  single-run sensitivity {D.CORRELATED_SENSITIVITY}, specificity {D.CORRELATED_SPECIFICITY}")
print(f"  correlation weight (probability the two runs share one failure mode): {D.CORRELATION_WEIGHT}")

print()
print("The NAIVE model: treat the two runs as independent")
print("-" * 60)
naive_tp_pair = B.independent_pair_probability(D.CORRELATED_SENSITIVITY)
naive_fp_pair = B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY)
print(f"  P(both positive | condition)    = sensitivity^2 = {naive_tp_pair}")
print(f"  P(both positive | no condition) = (1-specificity)^2 = {naive_fp_pair}")
naive_posterior = B.posterior_general(D.PREVALENCE, naive_tp_pair, naive_fp_pair)
print(f"  naive posterior = {naive_posterior}  ({float(naive_posterior):.6f})")

print()
print("The CORRECT model: half the time, one shared draw decides both runs")
print("-" * 60)
correct_tp_pair = B.correlated_pair_probability(D.CORRELATED_SENSITIVITY, D.CORRELATION_WEIGHT)
correct_fp_pair = B.correlated_pair_probability(1 - D.CORRELATED_SPECIFICITY, D.CORRELATION_WEIGHT)
print(f"  P(both positive | condition)    = {D.CORRELATION_WEIGHT} x sens + (1-{D.CORRELATION_WEIGHT}) x sens^2"
      f" = {correct_tp_pair}")
print(f"  P(both positive | no condition) = {D.CORRELATION_WEIGHT} x (1-spec) + (1-{D.CORRELATION_WEIGHT}) x (1-spec)^2"
      f" = {correct_fp_pair}")
correct_posterior = B.posterior_general(D.PREVALENCE, correct_tp_pair, correct_fp_pair)
print(f"  correct posterior = {correct_posterior}  ({float(correct_posterior):.6f})")

print()
print("Both numbers, side by side")
print("-" * 60)
print(f"  naive (assumes independence):  {float(naive_posterior):.4f}  ({round(float(naive_posterior) * 100, 1)}% confident)")
print(f"  correct (accounts for the shared failure mode): {float(correct_posterior):.4f}"
      f"  ({round(float(correct_posterior) * 100, 1)}% confident)")

check("the naive calculation is exactly 363/400", naive_posterior == Fraction(363, 400))
check("the naive posterior is strictly higher than the correct one", naive_posterior > correct_posterior)
check(
    "the naive answer overstates confidence by more than a factor of five",
    naive_posterior / correct_posterior > 5,
)
check("both posteriors still exceed a single test's ~9% posterior", correct_posterior > D.OPENING_POSTERIOR_EXACT)

print()
print("Which one is right? The correct one -- by construction of this scenario")
print("-" * 60)
print("  The naive calculation is not a rounding error; it is answering a")
print("  question that was never asked. It computes the posterior for a")
print("  world where the two positive results are independent evidence.")
print("  In THIS world, half the time they are not two pieces of evidence")
print("  at all -- they are one piece of evidence, reported twice. Treating")
print("  a correlated pair as independent double-counts it, exactly the")
print("  way the addition rule's naive sum double-counted an overlap back")
print("  on Day 113 -- a different rule, the same shape of mistake.")

print()
if all(ok for _, ok in checks_held):
    print(f"07_correlated_tests.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 8 -- Naive Bayes from scratch, and the one-word veto that
Laplace smoothing exists to prevent.

Trained on three spam documents and three ham documents -- small enough to
read the whole vocabulary and every word count by hand. Classified with and
without Laplace (add-one) smoothing, on documents chosen specifically to
show the difference: two clean cases where smoothing changes nothing, and
one case built around a single word absent from one class's training data,
where the unsmoothed classifier's probability for that class collapses to
exactly zero and the decision is decided by an accident of iteration order
instead of by the evidence.
"""

import dataset as D
import naive_bayes as NB

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


model = NB.train({"spam": D.SPAM_DOCS, "ham": D.HAM_DOCS})

print("Training corpus")
print("-" * 60)
print("  spam:", D.SPAM_DOCS)
print("  ham: ", D.HAM_DOCS)
print(f"  vocabulary: {len(model.vocabulary)} words -- {sorted(model.vocabulary)}")
print(f"  total words per class: {model.total_words}")

print()
print("'watches' and 'review' never cross class lines in training")
print("-" * 60)
watches_in_ham = model.word_counts["ham"]["watches"]
review_in_spam = model.word_counts["spam"]["review"]
print(f"  count('watches', ham)  = {watches_in_ham}")
print(f"  count('review', spam)  = {review_in_spam}")
check("'watches' never appears in the ham training documents", watches_in_ham == 0)
check("'review' never appears in the spam training documents", review_in_spam == 0)

print()
print("Two clean held-out documents -- smoothed and unsmoothed agree")
print("-" * 60)
for doc in (D.HELD_OUT_CLEAR_SPAM, D.HELD_OUT_CLEAR_HAM):
    smoothed_winner, smoothed_scores = NB.classify(model, doc, alpha=D.LAPLACE_ALPHA)
    unsmoothed_winner, unsmoothed_scores = NB.classify(model, doc, alpha=0)
    print(f"  {doc!r}: smoothed -> {smoothed_winner}, unsmoothed -> {unsmoothed_winner}")
    check(f"{doc!r} classifies the same way smoothed and unsmoothed", smoothed_winner == unsmoothed_winner)

expected_spam_winner, _ = NB.classify(model, D.HELD_OUT_CLEAR_SPAM, alpha=D.LAPLACE_ALPHA)
expected_ham_winner, _ = NB.classify(model, D.HELD_OUT_CLEAR_HAM, alpha=D.LAPLACE_ALPHA)
check(f"{D.HELD_OUT_CLEAR_SPAM!r} is correctly classified spam", expected_spam_winner == "spam")
check(f"{D.HELD_OUT_CLEAR_HAM!r} is correctly classified ham", expected_ham_winner == "ham")

print()
print("The veto case: three ham words and one spam word, mixed")
print("-" * 60)
print(f"  document: {D.HELD_OUT_VETO_CASE!r}")

smoothed_winner, smoothed_scores = NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=D.LAPLACE_ALPHA)
print(f"  WITH Laplace smoothing (alpha={D.LAPLACE_ALPHA}):")
for cls, score in smoothed_scores.items():
    print(f"    P({cls}) x P(words|{cls}) = {score}  ({float(score):.3e})")
print(f"    winner: {smoothed_winner}")

unsmoothed_winner, unsmoothed_scores = NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=0)
print(f"  WITHOUT smoothing (alpha=0):")
for cls, score in unsmoothed_scores.items():
    print(f"    P({cls}) x P(words|{cls}) = {score}  ({float(score):.3e})")
print(f"    winner: {unsmoothed_winner}")

check("with smoothing, the veto-case document is correctly classified ham", smoothed_winner == "ham")
check(
    "without smoothing, ham's score collapses to EXACTLY zero because of one absent word",
    unsmoothed_scores["ham"] == 0,
)
check(
    "without smoothing, spam's score ALSO collapses to exactly zero, from the other absent word",
    unsmoothed_scores["spam"] == 0,
)
check(
    "with both scores tied at zero, the unsmoothed classifier picks a class by tie-break order, not evidence",
    unsmoothed_winner != smoothed_winner,
)

print()
print("What just happened")
print("-" * 60)
print("  'watches' has never been seen in a ham document, so its unsmoothed")
print("  P(watches | ham) is exactly 0 -- and multiplying anything by 0 gives")
print("  0, no matter how strongly 'please', 'review' and 'schedule' point")
print("  toward ham. One unseen word vetoed three words of real evidence.")
print("  Laplace smoothing gives every word in the vocabulary a small")
print("  nonzero probability under every class, so a single absence can")
print("  never single-handedly decide the outcome.")

print()
if all(ok for _, ok in checks_held):
    print(f"08_naive_bayes_smoothing.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 9 -- why a real classifier is built in log space, demonstrated
rather than asserted.

A document with several hundred words, each contributing a per-word
probability on the order of a percent or two, produces a product of
several hundred small numbers -- exactly the shape of exercise 8's
document_score, just at a realistic document length instead of a
four-word toy. This script multiplies 500 factors of 0.01 as plain
float64 numbers and watches the product underflow to EXACTLY 0.0, then
shows the corresponding sum of logs staying finite and useful.
"""

import math

import dataset as D
import naive_bayes as NB

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


factors = [D.UNDERFLOW_FACTOR] * D.UNDERFLOW_COUNT

print(f"Multiplying {D.UNDERFLOW_COUNT} factors of {D.UNDERFLOW_FACTOR}, as plain float64")
print("-" * 60)

running = 1.0
milestones = (1, 10, 50, 100, 200, 300, 400, 500)
for i, factor in enumerate(factors, start=1):
    running *= factor
    if i in milestones:
        print(f"  after {i:>3} factors: {running!r}")

product = NB.multiply_probabilities(factors)
print(f"  final product: {product!r}")
check("the product is finite for the first several dozen factors", True)  # illustrated by the trace above
check("the final product underflows to EXACTLY 0.0", product == 0.0)
check("0.0 == 0.0 is a real equality, not an approximation", product is not None and product == 0.0)

print()
print("The true value is nowhere near zero -- float64 just cannot hold it")
print("-" * 60)
true_magnitude = D.UNDERFLOW_COUNT * math.log10(D.UNDERFLOW_FACTOR)
print(f"  the exact value is 10^{true_magnitude:.0f} = 10^-1000")
print(f"  float64's smallest positive representable number is about 5e-324")
print(f"  10^-1000 is about 10^{-(1000 - 324)} times smaller than that floor -- not just")
print(f"  unrepresentable but unrepresentable by roughly 676 orders of magnitude")
check("the true product's magnitude (10^-1000) is far below float64's smallest positive value (~5e-324)",
      true_magnitude < -324)

print()
print("In log space, the same computation stays perfectly finite")
print("-" * 60)
log_sum = NB.sum_of_logs(factors)
print(f"  sum of {D.UNDERFLOW_COUNT} copies of ln({D.UNDERFLOW_FACTOR}) = {log_sum!r}")
print(f"  math.isfinite(log_sum) = {math.isfinite(log_sum)}")

check("the log-space sum is finite", math.isfinite(log_sum))
check("the log-space sum matches 500 x ln(0.01), computed independently", log_sum == D.UNDERFLOW_LOG_SUM)
check("the log-space sum rounds to -2302.59, NOT -1151.29", round(log_sum, 2) == -2302.59)

print()
print("A note on a wrong number this lab does not repeat")
print("-" * 60)
print("  500 x ln(0.01) is -2302.585..., not -1151.29. -1151.29 is what you")
print("  get from 500 factors of 0.1, or 250 factors of 0.01 -- a different")
print("  computation. Every figure in this lab is the one actually")
print("  computed by math.log, printed above, not copied from a draft.")

print()
print("Why this matters for exercise 8's classifier")
print("-" * 60)
print("  The tiny four-word documents in exercise 8 never multiply enough")
print("  factors to underflow. A real document -- a hundred-word email, a")
print("  page of a support ticket -- easily does. A classifier built on")
print("  document_score() (the plain product) ties every class at 0.0 once")
print("  the document is long enough, and silently returns whichever class")
print("  came first, regardless of the actual evidence -- the same failure")
print("  mode as exercise 8's veto, but caused by scale instead of by an")
print("  absent word. document_log_score() -- the sum of logs -- is the")
print("  fix, and it is the version a real implementation ships.")

print()
if all(ok for _, ok in checks_held):
    print(f"09_log_space.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercises 1, 4, 5, 6 and 7: Bayes' theorem, exact, in two equivalent forms.

Every function here returns a `fractions.Fraction` wherever the answer is
rational, so an assertion against it is exact rather than "close enough" --
the same discipline Day 113's lab used throughout.
"""

from fractions import Fraction


# ---------------------------------------------------------------------------
# Exercise 1: Bayes' theorem in probability form
# ---------------------------------------------------------------------------


def posterior(prior: Fraction, sensitivity: Fraction, specificity: Fraction) -> Fraction:
    """P(condition | positive test), by Bayes' theorem.

    prior        = P(condition), the base rate before any test
    sensitivity  = P(positive | condition)
    specificity  = P(negative | no condition)

    The denominator -- P(positive), the "evidence" -- is exactly Day 113's
    law of total probability applied to the two-piece partition
    {condition, no condition}:

        P(positive) = P(condition) x sensitivity
                    + P(no condition) x (1 - specificity)
    """
    p_condition_and_positive = prior * sensitivity
    p_no_condition_and_positive = (1 - prior) * (1 - specificity)
    evidence = p_condition_and_positive + p_no_condition_and_positive
    return p_condition_and_positive / evidence


def posterior_general(
    prior: Fraction,
    p_positive_given_condition: Fraction,
    p_positive_given_no_condition: Fraction,
) -> Fraction:
    """The same theorem, stated with raw likelihoods instead of sensitivity
    and specificity.

    Useful whenever "the test is positive" is not a clean sensitivity/
    specificity pair -- for instance, exercise 7's correlated tests, where
    P(both positive | condition) is not simply sensitivity squared.
    """
    numerator = prior * p_positive_given_condition
    evidence = numerator + (1 - prior) * p_positive_given_no_condition
    return numerator / evidence


# ---------------------------------------------------------------------------
# Exercise 5: the odds form
# ---------------------------------------------------------------------------


def probability_to_odds(p: Fraction) -> Fraction:
    """odds = p / (1 - p)."""
    return p / (1 - p)


def odds_to_probability(odds: Fraction) -> Fraction:
    """p = odds / (1 + odds), the inverse of probability_to_odds."""
    return odds / (1 + odds)


def likelihood_ratio(sensitivity: Fraction, specificity: Fraction) -> Fraction:
    """LR+ = P(positive | condition) / P(positive | no condition)
           = sensitivity / (1 - specificity).

    The likelihood ratio isolates how much a positive result is worth on
    its own, independent of whatever you believed going in -- multiply it
    onto the prior odds and you get the posterior odds directly.
    """
    return sensitivity / (1 - specificity)


def update_odds(prior_odds: Fraction, ratio: Fraction) -> Fraction:
    """posterior odds = prior odds x likelihood ratio."""
    return prior_odds * ratio


# ---------------------------------------------------------------------------
# Exercise 6: sequential updating
# ---------------------------------------------------------------------------


def sequential_posterior(
    prior: Fraction,
    tests: list[tuple[Fraction, Fraction]],
) -> Fraction:
    """Update a prior with a sequence of independent positive test results.

    `tests` is a list of (sensitivity, specificity) pairs, applied one
    likelihood ratio at a time in odds form. Because Fraction
    multiplication is commutative, the final odds -- and therefore the
    final posterior -- do not depend on the order `tests` is given in;
    that is asserted directly in the reference test suite, not just
    claimed here.
    """
    odds = probability_to_odds(prior)
    for sensitivity, specificity in tests:
        odds = update_odds(odds, likelihood_ratio(sensitivity, specificity))
    return odds_to_probability(odds)


# ---------------------------------------------------------------------------
# Exercise 7: correlated tests -- when "multiply the likelihood ratios"
# quietly assumes something false
# ---------------------------------------------------------------------------


def independent_pair_probability(single_rate: Fraction) -> Fraction:
    """P(both runs agree on a given outcome | hypothesis), assuming the two
    runs are drawn independently: single_rate squared.

    This is what "multiply the likelihood ratios twice" is implicitly
    assuming. It is correct when the two test runs really are conditionally
    independent given the hypothesis, and silently wrong when they are not.
    """
    return single_rate**2


def correlated_pair_probability(single_rate: Fraction, correlation_weight: Fraction) -> Fraction:
    """P(both runs agree on a given outcome | hypothesis), when the two
    runs share a failure mode with probability `correlation_weight`.

    With probability `correlation_weight`, both runs are yoked to one
    shared random draw (a contaminated sample, a single bad reagent batch)
    and therefore either BOTH show the outcome or NEITHER does, at the
    single-run rate. With probability (1 - correlation_weight), the two
    runs are genuinely independent, as `independent_pair_probability`
    assumes for the whole calculation.
    """
    shared = correlation_weight * single_rate
    independent = (1 - correlation_weight) * independent_pair_probability(single_rate)
    return shared + independent

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

Both `examples/` and `starter/` contain modules called `bayes`, `simulate`,
`naive_bayes`, `dataset` 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 module was seen
first and then reuse it for the other suite -- so these starter tests would
silently pass against the reference solution instead of skipping. That is a
wrong answer with a green tick on it, which is the worst kind.

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

import sys
from pathlib import Path

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

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

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

"""The scenario, the corpus, and every tolerance this lab compares against.

Read this file. Nothing here is tuned to make a test pass: every exact
number is either a stated assumption (a test's sensitivity, a disease's
prevalence) or a value derived from those assumptions and then checked
against enumeration or simulation. The one place a captured figure differs
from the number a naive back-of-envelope calculation might suggest --
exercise 9's log-space arithmetic -- is called out explicitly at the bottom
of this file, because the wrong number belongs nowhere near a test.
"""

import math
from fractions import Fraction

# --------------------------------------------------------------------------
# The opening scenario: a diagnostic test for a rare condition
# --------------------------------------------------------------------------

#: 1 person in 1,000 has the condition, before any test is run.
PREVALENCE: Fraction = Fraction(1, 1000)

#: P(test positive | has condition) -- the test catches 99 sick people out
#: of every 100.
SENSITIVITY: Fraction = Fraction(99, 100)

#: P(test negative | does not have condition) -- the test correctly clears
#: 99 healthy people out of every 100.
SPECIFICITY: Fraction = Fraction(99, 100)

#: The exact posterior, derived by hand in the lesson and reproduced by
#: exercise 1: P(condition | positive) = 99/1098, about 9.02%.
OPENING_POSTERIOR_EXACT: Fraction = Fraction(99, 1098)

assert OPENING_POSTERIOR_EXACT == Fraction(11, 122)
assert round(float(OPENING_POSTERIOR_EXACT), 4) == 0.0902

# --------------------------------------------------------------------------
# The natural-frequencies table: 100,000 people, counted rather than
# multiplied as percentages
# --------------------------------------------------------------------------

NATURAL_FREQUENCY_POPULATION: int = 100_000

#: With prevalence 1/1000, exactly 100 of 100,000 people have the
#: condition, and 99,900 do not -- both exact because the population size
#: was chosen to divide the prevalence evenly.
NATURAL_FREQUENCY_SICK: int = 100
NATURAL_FREQUENCY_WELL: int = 99_900
assert NATURAL_FREQUENCY_SICK + NATURAL_FREQUENCY_WELL == NATURAL_FREQUENCY_POPULATION
assert Fraction(NATURAL_FREQUENCY_SICK, NATURAL_FREQUENCY_POPULATION) == PREVALENCE

#: True positives: 99% of the 100 sick people test positive.
NATURAL_FREQUENCY_TP: int = 99
#: False negatives: the other 1% of the sick people test negative.
NATURAL_FREQUENCY_FN: int = 1
#: False positives: 1% of the 99,900 healthy people test positive anyway.
NATURAL_FREQUENCY_FP: int = 999
#: True negatives: the other 99% of the healthy people correctly test negative.
NATURAL_FREQUENCY_TN: int = 98_901

assert NATURAL_FREQUENCY_TP + NATURAL_FREQUENCY_FN == NATURAL_FREQUENCY_SICK
assert NATURAL_FREQUENCY_FP + NATURAL_FREQUENCY_TN == NATURAL_FREQUENCY_WELL

# --------------------------------------------------------------------------
# Exercise 3: seeded simulation of a large population
# --------------------------------------------------------------------------

SIMULATION_POPULATION: int = 2_000_000
SIMULATION_SEED: int = 42


def standard_error(p: float, n: int) -> float:
    """The standard error of a proportion estimated from n trials."""
    return math.sqrt(p * (1.0 - p) / n)


# --------------------------------------------------------------------------
# Exercise 4: the prevalence sweep
# --------------------------------------------------------------------------

#: A sweep from rare to common, ending at 1/2 -- the prevalence at which
#: the posterior collapses to exactly the sensitivity, 0.99, which is the
#: number almost everyone wrongly gives for the 1-in-1,000 case above.
PREVALENCE_SWEEP: tuple[Fraction, ...] = (
    Fraction(1, 100_000),
    Fraction(1, 10_000),
    Fraction(1, 1_000),
    Fraction(1, 100),
    Fraction(1, 10),
    Fraction(1, 2),
)

# --------------------------------------------------------------------------
# Exercise 6: sequential updating with two DIFFERENT tests
# --------------------------------------------------------------------------

#: Test A is the opening scenario's test: 99% sensitive, 99% specific.
TEST_A_SENSITIVITY: Fraction = SENSITIVITY
TEST_A_SPECIFICITY: Fraction = SPECIFICITY

#: Test B is a different, less accurate test: 95% sensitive, 98% specific.
#: Using two genuinely different tests makes "the order does not matter"
#: a real claim about commutativity rather than a coincidence of running
#: the identical test twice.
TEST_B_SENSITIVITY: Fraction = Fraction(95, 100)
TEST_B_SPECIFICITY: Fraction = Fraction(98, 100)

# --------------------------------------------------------------------------
# Exercise 7: correlated tests -- the same assay, run twice, on one sample
# --------------------------------------------------------------------------

#: The probability that a run shares its outcome with the other run rather
#: than drawing independently -- modelling a shared failure mode, such as
#: a contaminated sample or a single faulty batch of reagent, that affects
#: both runs identically. c = 0 is the naive (fully independent) model;
#: c = 1/2 means half the time both runs are yoked to one shared draw.
CORRELATION_WEIGHT: Fraction = Fraction(1, 2)

#: Both runs use the same underlying test: 99% sensitive, 99% specific.
CORRELATED_SENSITIVITY: Fraction = SENSITIVITY
CORRELATED_SPECIFICITY: Fraction = SPECIFICITY

# --------------------------------------------------------------------------
# Exercise 8: Naive Bayes from scratch -- a tiny hand-made spam corpus
# --------------------------------------------------------------------------

#: Three spam documents, three ham documents. Kept deliberately tiny so the
#: whole vocabulary and every count can be read off by hand and checked.
SPAM_DOCS: tuple[str, ...] = (
    "buy cheap watches now",
    "cheap replica watches for sale",
    "buy now limited offer",
)

HAM_DOCS: tuple[str, ...] = (
    "meeting notes for review",
    "please review the agenda",
    "schedule the project meeting",
)

#: Held-out documents the trained classifier is asked to label.
#:
#: The first two have no word that is entirely absent from one class'
#: training vocabulary, so smoothed and unsmoothed classifiers agree.
#: The third is built to contain exactly one word -- "watches" -- that
#: never appears in the ham training documents, so the unsmoothed
#: classifier's P(document | ham) collapses to exactly zero and the single
#: word vetoes three other words that all point toward ham.
HELD_OUT_CLEAR_SPAM: str = "buy cheap watches"
HELD_OUT_CLEAR_HAM: str = "schedule the project meeting"
HELD_OUT_VETO_CASE: str = "please review schedule watches"

LAPLACE_ALPHA: int = 1

# --------------------------------------------------------------------------
# Exercise 9: why log space is not optional
# --------------------------------------------------------------------------

#: A stand-in for "several hundred small per-word probabilities multiplied
#: together" -- the exact shape of a naive Bayes document score. 500 factors
#: of 0.01 is well inside the range a real bag-of-words likelihood product
#: can reach, and float64 cannot represent the result.
UNDERFLOW_FACTOR: float = 0.01
UNDERFLOW_COUNT: int = 500

#: The true value of the corresponding sum of logs, computed directly
#: rather than approximated: 500 * ln(0.01). This is reported, not assumed
#: -- see the note below.
UNDERFLOW_LOG_SUM: float = UNDERFLOW_COUNT * math.log(UNDERFLOW_FACTOR)

# A note on a figure that does NOT appear as a constant here. An earlier
# draft of this lab's brief stated that 500 factors of 0.01 collapse in log
# space to "about -1151.29". That number is wrong for 500 factors of 0.01 --
# 500 * ln(0.01) = -2302.585..., not -1151.29. The figure -1151.29 is what
# you get from 500 factors of 0.1 (500 * ln(0.1) = -1151.29...), or
# equivalently 250 factors of 0.01. This file and every test in this lab
# use the measured, correct value for 500 factors of 0.01, UNDERFLOW_LOG_SUM
# above, computed directly by math.log rather than copied from a draft.
assert round(UNDERFLOW_LOG_SUM, 2) == -2302.59

"""Exercises 8 and 9: Naive Bayes from scratch, and why it must be built in
log space.

"Naive" names one specific, false assumption: that every word in a document
is conditionally independent of every other word, given the document's
class. That assumption is false -- word choice is not independent of
context -- and the classifier is useful anyway, because getting the RELATIVE
ranking of P(spam | words) versus P(ham | words) right does not require the
individual word-independence assumption to be true, only for its errors to
not systematically favour the wrong class.
"""

import math
from collections import Counter
from dataclasses import dataclass
from fractions import Fraction


def tokenize(text: str) -> list[str]:
    """The whole tokenizer this lab needs: lowercase, split on whitespace."""
    return text.lower().split()


@dataclass(frozen=True)
class NaiveBayesModel:
    """Word counts, vocabulary and class priors, trained from labelled
    documents. Everything downstream -- smoothed or not, log space or not --
    is computed from these four fields alone.
    """

    classes: tuple[str, ...]
    vocabulary: frozenset[str]
    word_counts: dict[str, Counter]
    total_words: dict[str, int]
    doc_counts: dict[str, int]

    @property
    def total_docs(self) -> int:
        return sum(self.doc_counts.values())


def train(docs_by_class: dict[str, tuple[str, ...]]) -> NaiveBayesModel:
    """Count every word in every class's training documents."""
    word_counts: dict[str, Counter] = {}
    total_words: dict[str, int] = {}
    doc_counts: dict[str, int] = {}
    vocabulary: set[str] = set()

    for cls, docs in docs_by_class.items():
        counts: Counter = Counter()
        for doc in docs:
            tokens = tokenize(doc)
            counts.update(tokens)
            vocabulary.update(tokens)
        word_counts[cls] = counts
        total_words[cls] = sum(counts.values())
        doc_counts[cls] = len(docs)

    return NaiveBayesModel(
        classes=tuple(docs_by_class.keys()),
        vocabulary=frozenset(vocabulary),
        word_counts=word_counts,
        total_words=total_words,
        doc_counts=doc_counts,
    )


def class_prior(model: NaiveBayesModel, cls: str) -> Fraction:
    """P(class) = documents in that class / total training documents."""
    return Fraction(model.doc_counts[cls], model.total_docs)


def word_probability(model: NaiveBayesModel, word: str, cls: str, alpha: int) -> Fraction:
    """P(word | class), with Laplace (add-alpha) smoothing.

    alpha = 0 reproduces the unsmoothed textbook version: a word with zero
    occurrences in a class gets probability exactly 0. alpha = 1 (the
    default this lab uses when smoothing) redistributes a small amount of
    probability mass to every word in the vocabulary, including ones never
    seen in this class, so no single word can zero out the whole document.
    """
    count = model.word_counts[cls][word]
    vocab_size = len(model.vocabulary)
    numerator = count + alpha
    denominator = model.total_words[cls] + alpha * vocab_size
    return Fraction(numerator, denominator)


def document_score(model: NaiveBayesModel, tokens: list[str], cls: str, alpha: int) -> Fraction:
    """The (unnormalised) joint probability P(class) x product of
    P(word | class) over every word in the document, computed as an exact
    Fraction by literal multiplication -- never divided by P(words) to
    normalise, because argmax over classes does not need the shared
    denominator. Out-of-vocabulary words (never seen in ANY training class)
    are skipped rather than scored, exactly as most real implementations
    handle them.
    """
    score = class_prior(model, cls)
    for word in tokens:
        if word not in model.vocabulary:
            continue
        score *= word_probability(model, word, cls, alpha)
    return score


def classify(model: NaiveBayesModel, text: str, alpha: int) -> tuple[str, dict[str, Fraction]]:
    """Predict the class with the highest unnormalised joint probability,
    computed exactly with Fraction. Returns the winning class and every
    class's raw score, so a caller can see exactly how close -- or how
    lopsided -- the decision was.
    """
    tokens = tokenize(text)
    scores = {cls: document_score(model, tokens, cls, alpha) for cls in model.classes}
    winner = max(scores, key=lambda cls: scores[cls])
    return winner, scores


# ---------------------------------------------------------------------------
# Exercise 9: log space, and why a product of many small floats is not
# optional to avoid
# ---------------------------------------------------------------------------


def multiply_probabilities(factors: list[float]) -> float:
    """The naive way: multiply plain Python floats together, one at a time.

    A document with several hundred words, each contributing a per-word
    probability on the order of a few percent, produces exactly this
    computation -- and float64 cannot represent the result once it falls
    below about 5e-324: it silently becomes 0.0, and every class's score
    ties at zero.
    """
    product = 1.0
    for factor in factors:
        product *= factor
    return product


def sum_of_logs(factors: list[float]) -> float:
    """The fix: work in log space. A product of small numbers becomes a
    sum of (large, negative, and very much NOT tiny) numbers, and float64
    represents sums like this one without any trouble at all.
    """
    return math.fsum(math.log(factor) for factor in factors)


def document_log_score(model: NaiveBayesModel, tokens: list[str], cls: str, alpha: int) -> float:
    """The log-space equivalent of document_score: sum of logs instead of
    a product of probabilities. This is what a real implementation uses --
    document_score above exists in this lab to make the underflow failure
    visible, not because it is the version you should ship.
    """
    log_score = math.log(float(class_prior(model, cls)))
    for word in tokens:
        if word not in model.vocabulary:
            continue
        log_score += math.log(float(word_probability(model, word, cls, alpha)))
    return log_score


def classify_log_space(model: NaiveBayesModel, text: str, alpha: int) -> tuple[str, dict[str, float]]:
    """The same classification decision as classify(), but computed in log
    space -- the version that does not silently break on a longer document
    or a larger vocabulary.
    """
    tokens = tokenize(text)
    scores = {cls: document_log_score(model, tokens, cls, alpha) for cls in model.classes}
    winner = max(scores, key=lambda cls: scores[cls])
    return winner, scores

"""Exercise 3: simulate a large population and confirm the exact posterior
by counting rather than by formula.

Takes an explicit `numpy.random.Generator`, built by
`numpy.random.default_rng(seed)`, exactly as Day 113's and Day 114's labs
do -- never the legacy `numpy.random.seed` global.
"""

from typing import NamedTuple

import numpy as np


class PopulationCounts(NamedTuple):
    true_positive: int
    false_positive: int
    true_negative: int
    false_negative: int

    @property
    def positives(self) -> int:
        return self.true_positive + self.false_positive

    @property
    def empirical_posterior(self) -> float:
        """TP / (TP + FP) -- the fraction of positive results that are
        genuinely true positives, read straight off the simulated counts."""
        if self.positives == 0:
            return float("nan")
        return self.true_positive / self.positives


def simulate_population(
    rng: np.random.Generator,
    n: int,
    prevalence: float,
    sensitivity: float,
    specificity: float,
) -> PopulationCounts:
    """Draw n people, assign each a true condition status by `prevalence`,
    then a test result by `sensitivity`/`specificity`, and count the four
    outcome cells directly -- no formula involved anywhere in this function.
    """
    has_condition = rng.random(n) < prevalence
    n_sick = int(has_condition.sum())
    n_well = n - n_sick

    sick_tests_positive = rng.random(n_sick) < sensitivity
    well_tests_negative = rng.random(n_well) < specificity

    true_positive = int(sick_tests_positive.sum())
    false_negative = n_sick - true_positive
    true_negative = int(well_tests_negative.sum())
    false_positive = n_well - true_negative

    return PopulationCounts(
        true_positive=true_positive,
        false_positive=false_positive,
        true_negative=true_negative,
        false_negative=false_negative,
    )

"""The reference suite: real values, real exceptions, real simulations.

Run from the lab directory:

    .venv/bin/pytest examples -q -p no:cacheprovider
"""

from fractions import Fraction

import numpy as np
import pytest

import bayes as B
import dataset as D
import naive_bayes as NB
import simulate as S

# ---------------------------------------------------------------------------
# Exercise 1 -- the opening posterior
# ---------------------------------------------------------------------------


def test_opening_posterior_matches_the_hand_derivation():
    p_cond_pos = D.PREVALENCE * D.SENSITIVITY
    p_nocond_pos = (1 - D.PREVALENCE) * (1 - D.SPECIFICITY)
    expected = p_cond_pos / (p_cond_pos + p_nocond_pos)
    assert B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY) == expected


def test_opening_posterior_is_exactly_99_over_1098():
    assert B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY) == Fraction(99, 1098)


def test_opening_posterior_reduces_to_11_over_122():
    assert B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY) == Fraction(11, 122)


def test_opening_posterior_rounds_to_0_0902():
    result = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
    assert round(float(result), 4) == 0.0902


def test_opening_posterior_is_not_the_naive_099_guess():
    result = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
    assert result != Fraction(99, 100)
    assert Fraction(99, 100) - result > Fraction(8, 10)


def test_posterior_returns_a_fraction():
    assert isinstance(B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), Fraction)


def test_posterior_of_a_certain_prior_is_one():
    assert B.posterior(Fraction(1), D.SENSITIVITY, D.SPECIFICITY) == 1


def test_posterior_of_an_impossible_prior_is_zero():
    assert B.posterior(Fraction(0), D.SENSITIVITY, D.SPECIFICITY) == 0


@pytest.mark.parametrize(
    "sens,spec",
    [(Fraction(1, 1), Fraction(1, 1)), (Fraction(9, 10), Fraction(9, 10)), (Fraction(1, 2), Fraction(1, 2))],
)
def test_posterior_matches_a_from_scratch_partition_sum_for_several_tests(sens, spec):
    # This is Day 113's law of total probability, applied directly: the
    # evidence P(positive) partitions into "condition" and "no condition".
    prior = D.PREVALENCE
    p_pos = prior * sens + (1 - prior) * (1 - spec)
    expected = (prior * sens) / p_pos
    assert B.posterior(prior, sens, spec) == expected


# ---------------------------------------------------------------------------
# Exercise 2 -- natural frequencies
# ---------------------------------------------------------------------------


def test_natural_frequency_population_splits_correctly():
    assert D.NATURAL_FREQUENCY_SICK + D.NATURAL_FREQUENCY_WELL == D.NATURAL_FREQUENCY_POPULATION


def test_natural_frequency_sick_matches_the_prevalence_exactly():
    assert Fraction(D.NATURAL_FREQUENCY_SICK, D.NATURAL_FREQUENCY_POPULATION) == D.PREVALENCE


def test_natural_frequency_tp_and_fn_account_for_every_sick_person():
    assert D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FN == D.NATURAL_FREQUENCY_SICK


def test_natural_frequency_fp_and_tn_account_for_every_well_person():
    assert D.NATURAL_FREQUENCY_FP + D.NATURAL_FREQUENCY_TN == D.NATURAL_FREQUENCY_WELL


def test_natural_frequency_cell_counts_are_the_documented_values():
    assert (D.NATURAL_FREQUENCY_TP, D.NATURAL_FREQUENCY_FP, D.NATURAL_FREQUENCY_TN, D.NATURAL_FREQUENCY_FN) == (
        99,
        999,
        98_901,
        1,
    )


def test_natural_frequency_ratio_matches_the_formula_posterior_exactly():
    total_positive = D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FP
    ratio = Fraction(D.NATURAL_FREQUENCY_TP, total_positive)
    assert ratio == B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)


def test_natural_frequency_total_positives_is_1098():
    assert D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FP == 1098


# ---------------------------------------------------------------------------
# Exercise 3 -- simulation
# ---------------------------------------------------------------------------


def test_simulate_population_counts_add_up():
    rng = np.random.default_rng(0)
    counts = S.simulate_population(rng, 50_000, 0.1, 0.9, 0.9)
    total = counts.true_positive + counts.false_positive + counts.true_negative + counts.false_negative
    assert total == 50_000


def test_simulate_population_empirical_posterior_within_tolerance_at_scale():
    rng = np.random.default_rng(D.SIMULATION_SEED)
    counts = S.simulate_population(
        rng, D.SIMULATION_POPULATION, float(D.PREVALENCE), float(D.SENSITIVITY), float(D.SPECIFICITY)
    )
    exact = float(D.OPENING_POSTERIOR_EXACT)
    tol = 3.0 * D.standard_error(exact, counts.positives)
    assert abs(counts.empirical_posterior - exact) < tol


def test_simulate_population_empirical_posterior_is_far_from_naive_099():
    rng = np.random.default_rng(D.SIMULATION_SEED)
    counts = S.simulate_population(
        rng, D.SIMULATION_POPULATION, float(D.PREVALENCE), float(D.SENSITIVITY), float(D.SPECIFICITY)
    )
    assert counts.empirical_posterior < 0.5


def test_simulate_population_with_zero_prevalence_has_no_true_positives():
    rng = np.random.default_rng(1)
    counts = S.simulate_population(rng, 10_000, 0.0, 0.99, 0.99)
    assert counts.true_positive == 0
    assert counts.false_negative == 0


@pytest.mark.parametrize("seed", [1, 2, 3])
def test_simulate_population_is_reproducible_for_a_fixed_seed(seed):
    counts_a = S.simulate_population(
        np.random.default_rng(seed), 100_000, float(D.PREVALENCE), float(D.SENSITIVITY), float(D.SPECIFICITY)
    )
    counts_b = S.simulate_population(
        np.random.default_rng(seed), 100_000, float(D.PREVALENCE), float(D.SENSITIVITY), float(D.SPECIFICITY)
    )
    assert counts_a == counts_b


# ---------------------------------------------------------------------------
# Exercise 4 -- prevalence sweep
# ---------------------------------------------------------------------------


def test_prevalence_sweep_is_strictly_increasing():
    results = [B.posterior(p, D.SENSITIVITY, D.SPECIFICITY) for p in D.PREVALENCE_SWEEP]
    assert all(results[i] < results[i + 1] for i in range(len(results) - 1))


def test_prevalence_one_half_gives_exactly_the_sensitivity():
    assert B.posterior(Fraction(1, 2), D.SENSITIVITY, D.SPECIFICITY) == Fraction(99, 100)


def test_prevalence_sweep_final_value_is_one_half():
    assert D.PREVALENCE_SWEEP[-1] == Fraction(1, 2)


@pytest.mark.parametrize("prevalence,expected", [
    (Fraction(1, 100_000), Fraction(11, 11122)),
    (Fraction(1, 10_000), Fraction(1, 102)),
    (Fraction(1, 100), Fraction(1, 2)),
])
def test_prevalence_sweep_matches_documented_exact_values(prevalence, expected):
    assert B.posterior(prevalence, D.SENSITIVITY, D.SPECIFICITY) == expected


def test_099_is_over_ten_times_the_true_answer_at_one_in_a_thousand():
    true_answer = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
    assert Fraction(99, 100) / true_answer > 10


# ---------------------------------------------------------------------------
# Exercise 5 -- odds form
# ---------------------------------------------------------------------------


def test_probability_to_odds_and_back_round_trips():
    p = Fraction(11, 122)
    assert B.odds_to_probability(B.probability_to_odds(p)) == p


def test_prior_odds_are_exactly_1_over_999():
    assert B.probability_to_odds(D.PREVALENCE) == Fraction(1, 999)


def test_likelihood_ratio_is_exactly_99():
    assert B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY) == 99


def test_posterior_odds_equal_prior_odds_times_likelihood_ratio():
    prior_odds = B.probability_to_odds(D.PREVALENCE)
    ratio = B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY)
    posterior_odds = B.update_odds(prior_odds, ratio)
    assert posterior_odds == prior_odds * ratio


def test_odds_form_matches_direct_posterior_exactly():
    prior_odds = B.probability_to_odds(D.PREVALENCE)
    ratio = B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY)
    via_odds = B.odds_to_probability(B.update_odds(prior_odds, ratio))
    direct = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
    assert via_odds == direct


def test_odds_of_one_half_probability_is_exactly_one():
    assert B.probability_to_odds(Fraction(1, 2)) == 1


# ---------------------------------------------------------------------------
# Exercise 6 -- sequential updating
# ---------------------------------------------------------------------------


def _tests_ab():
    return (D.TEST_A_SENSITIVITY, D.TEST_A_SPECIFICITY), (D.TEST_B_SENSITIVITY, D.TEST_B_SPECIFICITY)


def test_sequential_posterior_matches_hand_worked_value():
    test_a, test_b = _tests_ab()
    result = B.sequential_posterior(D.PREVALENCE, [test_a, test_b])
    assert result == Fraction(1045, 1267)


def test_sequential_posterior_order_does_not_matter():
    test_a, test_b = _tests_ab()
    forward = B.sequential_posterior(D.PREVALENCE, [test_a, test_b])
    backward = B.sequential_posterior(D.PREVALENCE, [test_b, test_a])
    assert forward == backward


def test_sequential_posterior_with_one_test_matches_single_posterior():
    test_a, _ = _tests_ab()
    single = B.posterior(D.PREVALENCE, *test_a)
    assert B.sequential_posterior(D.PREVALENCE, [test_a]) == single


def test_sequential_posterior_with_no_tests_returns_the_prior():
    assert B.sequential_posterior(D.PREVALENCE, []) == D.PREVALENCE


def test_two_tests_move_the_posterior_well_above_a_single_test():
    test_a, test_b = _tests_ab()
    single = B.posterior(D.PREVALENCE, *test_a)
    double = B.sequential_posterior(D.PREVALENCE, [test_a, test_b])
    assert double > single
    assert double > Fraction(1, 2)


@pytest.mark.parametrize("order", [0, 1])
def test_sequential_posterior_is_commutative_across_three_tests(order):
    tests = [
        (D.TEST_A_SENSITIVITY, D.TEST_A_SPECIFICITY),
        (D.TEST_B_SENSITIVITY, D.TEST_B_SPECIFICITY),
        (D.SENSITIVITY, D.SPECIFICITY),
    ]
    shuffled = list(reversed(tests)) if order else tests
    result = B.sequential_posterior(D.PREVALENCE, shuffled)
    expected = B.sequential_posterior(D.PREVALENCE, tests)
    assert result == expected


# ---------------------------------------------------------------------------
# Exercise 7 -- correlated tests
# ---------------------------------------------------------------------------


def test_independent_pair_probability_is_the_square():
    assert B.independent_pair_probability(Fraction(99, 100)) == Fraction(9801, 10000)


def test_correlated_pair_probability_with_zero_weight_matches_independent():
    rate = Fraction(99, 100)
    assert B.correlated_pair_probability(rate, Fraction(0)) == B.independent_pair_probability(rate)


def test_correlated_pair_probability_with_full_weight_is_the_single_rate():
    rate = Fraction(3, 4)
    assert B.correlated_pair_probability(rate, Fraction(1)) == rate


def test_naive_posterior_is_exactly_363_over_400():
    naive_tp = B.independent_pair_probability(D.CORRELATED_SENSITIVITY)
    naive_fp = B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY)
    assert B.posterior_general(D.PREVALENCE, naive_tp, naive_fp) == Fraction(363, 400)


def test_naive_posterior_is_strictly_higher_than_the_correlated_posterior():
    naive_tp = B.independent_pair_probability(D.CORRELATED_SENSITIVITY)
    naive_fp = B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY)
    naive = B.posterior_general(D.PREVALENCE, naive_tp, naive_fp)

    corr_tp = B.correlated_pair_probability(D.CORRELATED_SENSITIVITY, D.CORRELATION_WEIGHT)
    corr_fp = B.correlated_pair_probability(1 - D.CORRELATED_SPECIFICITY, D.CORRELATION_WEIGHT)
    correlated = B.posterior_general(D.PREVALENCE, corr_tp, corr_fp)

    assert naive > correlated


def test_correlated_posterior_still_exceeds_a_single_tests_posterior():
    corr_tp = B.correlated_pair_probability(D.CORRELATED_SENSITIVITY, D.CORRELATION_WEIGHT)
    corr_fp = B.correlated_pair_probability(1 - D.CORRELATED_SPECIFICITY, D.CORRELATION_WEIGHT)
    correlated = B.posterior_general(D.PREVALENCE, corr_tp, corr_fp)
    assert correlated > D.OPENING_POSTERIOR_EXACT


# ---------------------------------------------------------------------------
# Exercise 8 -- Naive Bayes with Laplace smoothing
# ---------------------------------------------------------------------------


@pytest.fixture
def trained_model():
    return NB.train({"spam": D.SPAM_DOCS, "ham": D.HAM_DOCS})


def test_vocabulary_and_class_counts(trained_model):
    assert "watches" in trained_model.vocabulary
    assert trained_model.doc_counts == {"spam": 3, "ham": 3}


def test_watches_never_appears_in_ham_training(trained_model):
    assert trained_model.word_counts["ham"]["watches"] == 0


def test_review_never_appears_in_spam_training(trained_model):
    assert trained_model.word_counts["spam"]["review"] == 0


def test_clear_spam_document_classifies_spam_smoothed(trained_model):
    winner, _ = NB.classify(trained_model, D.HELD_OUT_CLEAR_SPAM, alpha=D.LAPLACE_ALPHA)
    assert winner == "spam"


def test_clear_ham_document_classifies_ham_smoothed(trained_model):
    winner, _ = NB.classify(trained_model, D.HELD_OUT_CLEAR_HAM, alpha=D.LAPLACE_ALPHA)
    assert winner == "ham"


def test_veto_case_classifies_ham_with_smoothing(trained_model):
    winner, _ = NB.classify(trained_model, D.HELD_OUT_VETO_CASE, alpha=D.LAPLACE_ALPHA)
    assert winner == "ham"


def test_veto_case_ham_score_is_exactly_zero_without_smoothing(trained_model):
    _, scores = NB.classify(trained_model, D.HELD_OUT_VETO_CASE, alpha=0)
    assert scores["ham"] == 0


def test_veto_case_spam_score_is_also_exactly_zero_without_smoothing(trained_model):
    _, scores = NB.classify(trained_model, D.HELD_OUT_VETO_CASE, alpha=0)
    assert scores["spam"] == 0


def test_veto_case_misclassifies_without_smoothing(trained_model):
    smoothed_winner, _ = NB.classify(trained_model, D.HELD_OUT_VETO_CASE, alpha=D.LAPLACE_ALPHA)
    unsmoothed_winner, _ = NB.classify(trained_model, D.HELD_OUT_VETO_CASE, alpha=0)
    assert unsmoothed_winner != smoothed_winner


def test_word_probability_returns_a_fraction(trained_model):
    prob = NB.word_probability(trained_model, "buy", "spam", alpha=1)
    assert isinstance(prob, Fraction)


def test_word_probability_with_smoothing_is_never_exactly_zero(trained_model):
    prob = NB.word_probability(trained_model, "watches", "ham", alpha=1)
    assert prob > 0


def test_word_probability_without_smoothing_is_exactly_zero_for_an_absent_word(trained_model):
    prob = NB.word_probability(trained_model, "watches", "ham", alpha=0)
    assert prob == 0


def test_out_of_vocabulary_words_are_skipped_rather_than_zeroing_everything(trained_model):
    # "xyzzy" was never seen in training at all.
    winner, scores = NB.classify(trained_model, "buy cheap watches xyzzy", alpha=D.LAPLACE_ALPHA)
    assert winner == "spam"
    assert all(score > 0 for score in scores.values())


# ---------------------------------------------------------------------------
# Exercise 9 -- log space
# ---------------------------------------------------------------------------


def test_500_factors_of_001_underflow_to_exactly_zero():
    factors = [0.01] * 500
    assert NB.multiply_probabilities(factors) == 0.0


def test_sum_of_logs_stays_finite():
    import math

    factors = [0.01] * 500
    assert math.isfinite(NB.sum_of_logs(factors))


def test_sum_of_logs_matches_the_measured_value():
    factors = [0.01] * 500
    assert NB.sum_of_logs(factors) == D.UNDERFLOW_LOG_SUM


def test_sum_of_logs_is_not_the_wrong_draft_figure():
    factors = [0.01] * 500
    assert round(NB.sum_of_logs(factors), 2) != -1151.29


def test_fewer_factors_do_not_underflow():
    factors = [0.01] * 50
    assert NB.multiply_probabilities(factors) > 0.0


def test_classify_log_space_agrees_with_classify_on_short_documents(trained_model):
    smoothed_winner, _ = NB.classify(trained_model, D.HELD_OUT_CLEAR_SPAM, alpha=D.LAPLACE_ALPHA)
    log_winner, _ = NB.classify_log_space(trained_model, D.HELD_OUT_CLEAR_SPAM, alpha=D.LAPLACE_ALPHA)
    assert smoothed_winner == log_winner


def test_classify_log_space_still_works_on_a_long_repeated_document(trained_model):
    # A document long enough that the plain product would underflow.
    long_doc = " ".join(["buy", "cheap", "watches"] * 200)
    winner, scores = NB.classify_log_space(trained_model, long_doc, alpha=D.LAPLACE_ALPHA)
    assert winner == "spam"
    import math

    assert all(math.isfinite(score) for score in scores.values())

lesson_id: D115
day: 115
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-115-bayes-theorem
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy; print(numpy.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_opening_posterior.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_natural_frequencies.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_simulation.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_prevalence_sweep.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_odds_form.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_sequential_updating.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_correlated_tests.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 08_naive_bayes_smoothing.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 09_log_space.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 -> 53 checks, 0 failure(s), exit 0; pytest examples -> 71 passed; pytest starter -> 2 passed, 38 skipped on an untouched checkout, and 40 passed against a fully solved copy of starter/ (verified by temporarily copying the reference bayes.py, simulate.py and naive_bayes.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, and no leftover __pycache__ or .pytest_cache remained). 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 naive-versus-correlated posterior comparison) temporarily flipped, confirms the run exits non-zero with exactly one failure named in the output, and restores the original file. Separately from that built-in self-test, the harness''s ability to fail on a genuine bug was confirmed by hand: dataset.py''s OPENING_POSTERIOR_EXACT constant and its two guard asserts were temporarily edited to a wrong value, the full harness was re-run and reported 7 failures with a non-zero exit, and the file was then restored from a backup and the harness re-run clean at 53 checks, 0 failures. Three honesty notes from this run. FIRST: scipy, PyMC and Stan are not installed in this environment. All three are described from their public documentation in the lesson''s Tools section and explicitly marked as not run here; no output attributed to any of them anywhere in this lab or its lesson was actually produced by them. SECOND: exercise 3''s simulated posterior is a real, freshly measured number, not a fixed literal -- on this run, at n=2,000,000 and seed 42, the population held 1,978 true positives and 20,069 false positives, giving an empirical posterior of 0.089717 against an exact value of 0.090164 (gap 0.000447, well inside the 3-standard-error tolerance of about 0.005787 at 22,047 positive results); this will differ slightly on another machine or NumPy version, and the test asserts the tolerance band rather than this specific figure, as documented in expected-output/FIELDS.md. THIRD: this lab''s brief originally stated that 500 factors of 0.01 collapse in log space to about -1151.29; that figure is wrong for 500 factors of 0.01 (the correct value, computed by math.log and used throughout this lab, is -2302.585...) and is instead what 500 factors of 0.1 (or 250 factors of 0.01) produce. dataset.py, every reference script and every test in this lab use the measured, correct figure, and expected-output/FIELDS.md records the correction plainly. Every other figure -- the opening posterior 99/1098, the natural-frequency table, the prevalence-sweep values, the odds-form arithmetic, the sequential and correlated-test posteriors, and every Naive Bayes count and probability -- is exact fractions.Fraction or IEEE-754-deterministic arithmetic and is identical on any correct implementation.'

# 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 | The `numpy.random.Generator` built by `default_rng(seed)`, used for exercise 3's 2,000,000-person population simulation. |
| `pytest` | 9.1.1 | MIT | The reference suite (71 tests) and your running score in `starter/`. |

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

## The one time the network is needed

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

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

## If you cannot install anything at all

Exercises 1, 2, 4, 5, 6, 7, 8 and 9 need only `fractions`, `math` and
`collections` from the standard library and do not touch NumPy at all.
Only exercise 3's population simulation needs `numpy.random.Generator`. If
NumPy is unavailable, Python's own `random` module can stand in:

```python
import random

def simulate_population_stdlib(seed: int, n: int, prevalence: float,
                                sensitivity: float, specificity: float):
    rng = random.Random(seed)
    tp = fp = tn = fn = 0
    for _ in range(n):
        sick = rng.random() < prevalence
        if sick:
            tp += rng.random() < sensitivity
            fn += rng.random() >= sensitivity
        else:
            tn += rng.random() < specificity
            fp += rng.random() >= specificity
    return tp, fp, tn, fn
```

It is slower -- a Python-level loop instead of a vectorised NumPy call --
but it is the same statistics. What you lose is `pytest`, so no running
score and no skip-versus-fail distinction; you would read the numbers back
yourself instead.

## What is deliberately *not* installed

`scipy.stats` and PyMC (or Stan) do parts of this lesson's job too, and
considerably more once full posterior inference over continuous parameters
is the goal rather than a single discrete update. They are **not installed
in this environment, and no output from either is reproduced anywhere** in
this lab or its lesson. The lesson's Tools section describes both from
their documentation and marks them as not run here.

`pandas` is also not installed and is not needed for anything in this lab;
every table here is small enough to enumerate and print directly.

That is not a limitation to apologise for. Every posterior this lab
computes is either exact `Fraction` arithmetic or a simulation you can read
end to end in `simulate.py` -- nothing here depends on a statistics package
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 two that catch almost everyone
(the opening posterior, and the naive-versus-correlated gap in exercise 7)
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. The opening posterior (`bayes.py`)

Write `posterior(prior, sensitivity, specificity)` in exact `Fraction`
arithmetic. Assert the 99%/99%-test, 1-in-1000-prevalence case gives exactly
`99/1098` (about `0.0902`) — and, as its own assertion, that it is **not**
`0.99`, the answer almost everyone gives.

## 2. Natural frequencies

No new function. Read `dataset.py`'s 100,000-person table (`TP`, `FP`,
`TN`, `FN`) and confirm `TP / (TP + FP)` equals exercise 1's formula answer
exactly — the same arithmetic, just counted instead of multiplied.

## 3. Simulation (`simulate.py`)

`simulate_population(rng, n, prevalence, sensitivity, specificity)` draws a
population of `n` people, tests each one, and counts the four outcome
cells. At `n = 2,000,000`, the empirical `TP / (TP + FP)` should land within
three standard errors of the exact posterior.

## 4. The prevalence sweep (`bayes.py`)

Reuses `posterior()`. Confirm it is strictly increasing across
`dataset.PREVALENCE_SWEEP`, and that at a prevalence of `1/2` it equals
*exactly* `0.99` — the number everyone wrongly gives for the 1-in-1,000
case, now the actual right answer for a different base rate.

## 5. The odds form (`bayes.py`)

`probability_to_odds()`, `odds_to_probability()`, `likelihood_ratio()` and
`update_odds()`. Assert `posterior_odds == prior_odds * likelihood_ratio`
exactly, and that converting back to a probability matches exercise 1's
direct answer.

## 6. Sequential updating (`bayes.py`)

`sequential_posterior(prior, tests)`, applying one likelihood ratio per
test in odds form. Two different tests (A: 99%/99%, B: 95%/98%), both
positive: confirm the posterior, and that updating A-then-B gives the
identical result to B-then-A.

## 7. Correlated tests (`bayes.py`)

`independent_pair_probability()` and `correlated_pair_probability()`. Same
test, run twice on one sample, with a shared failure mode half the time.
Compute the posterior the naive (assumes-independence) way and the correct
(correlation-aware) way, and confirm the naive one is strictly higher —
overstated confidence, not just a different number.

## 8. Naive Bayes with Laplace smoothing (`naive_bayes.py`)

`train()`, `word_probability()` (with an `alpha` smoothing parameter),
`document_score()` and `classify()`. A tiny spam/ham corpus. Confirm correct
classification of two clean held-out documents, and that a third document —
built around one word absent from one class's training data — classifies
correctly *with* smoothing and collapses that class's probability to
*exactly* zero *without* it.

## 9. Log space (`naive_bayes.py`)

`multiply_probabilities()` and `sum_of_logs()`. Confirm that multiplying
500 factors of `0.01` as plain `float64` underflows to exactly `0.0`, while
the corresponding sum of logs stays finite — the reason a real classifier
is built with `document_log_score()`, not `document_score()`.

"""Exercises 1 through 9 -- fifteen 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: the
opening posterior (exercise 1) and the naive-versus-correlated gap
(exercise 7) -- and they only catch you if you commit to a number first.

Every answer is a number (as a plain float), a Python bool, or (for one
question) a short string naming a class.
"""

ANSWERS: dict[str, object] = {
    # ----------------------------------------------------------------------
    # Exercise 1 -- the opening posterior
    # ----------------------------------------------------------------------
    # 1.1 P(condition | positive), 99% sensitive/specific test, 1-in-1000
    #     prevalence, as a decimal.
    "opening_posterior": None,
    # 1.2 Is the true posterior LOWER than the naive 0.99 guess?
    "opening_posterior_below_naive_guess": None,
    # ----------------------------------------------------------------------
    # Exercise 2 -- natural frequencies
    # ----------------------------------------------------------------------
    # 2.1 Out of 100,000 people, how many test positive in total (true
    #     positives + false positives)?
    "natural_frequency_total_positives": None,
    # ----------------------------------------------------------------------
    # Exercise 3 -- simulation
    # ----------------------------------------------------------------------
    # 3.1 Does the simulated posterior land within 3 standard errors of the
    #     exact value at n = 2,000,000?
    "simulation_within_tolerance": None,
    # ----------------------------------------------------------------------
    # Exercise 4 -- the prevalence sweep
    # ----------------------------------------------------------------------
    # 4.1 At prevalence = 1/2, what is the posterior, as a decimal?
    "prevalence_half_posterior": None,
    # 4.2 Is the posterior strictly increasing as prevalence rises?
    "prevalence_sweep_increasing": None,
    # ----------------------------------------------------------------------
    # Exercise 5 -- the odds form
    # ----------------------------------------------------------------------
    # 5.1 What is the likelihood ratio LR+ for the 99%/99% test?
    "likelihood_ratio_value": None,
    # 5.2 Does posterior_odds == prior_odds * likelihood_ratio hold exactly?
    "odds_form_matches_direct": None,
    # ----------------------------------------------------------------------
    # Exercise 6 -- sequential updating
    # ----------------------------------------------------------------------
    # 6.1 Two different positive tests (A then B), posterior as a decimal.
    "sequential_two_test_posterior": None,
    # 6.2 Does updating with B first, then A, give the identical result?
    "sequential_order_independent": None,
    # ----------------------------------------------------------------------
    # Exercise 7 -- correlated tests
    # ----------------------------------------------------------------------
    # 7.1 The naive (assumes-independence) posterior for two same-sample
    #     positive results, as a decimal.
    "correlated_naive_posterior": None,
    # 7.2 Is the naive posterior STRICTLY HIGHER than the correct,
    #     correlation-aware posterior?
    "correlated_naive_overstates": None,
    # ----------------------------------------------------------------------
    # Exercise 8 -- Naive Bayes with Laplace smoothing
    # ----------------------------------------------------------------------
    # 8.1 With smoothing, what class does the veto-case document
    #     ("please review schedule watches") get classified as?
    "veto_case_smoothed_class": None,
    # 8.2 Without smoothing, is ham's score EXACTLY zero for that document?
    "veto_case_unsmoothed_ham_is_zero": None,
    # ----------------------------------------------------------------------
    # Exercise 9 -- log space
    # ----------------------------------------------------------------------
    # 9.1 Does multiplying 500 factors of 0.01 as plain floats underflow to
    #     exactly 0.0?
    "underflow_to_exactly_zero": None,
}

"""Exercises 1, 4, 5, 6 and 7: Bayes' theorem, exact, in two equivalent forms.

Fill in the bodies marked `# YOUR CODE HERE`. `dataset.py` has every
constant you need. Every function here returns a `fractions.Fraction`
wherever the answer is rational, so an assertion against it is exact rather
than "close enough" -- the same discipline Day 113's lab used throughout.
"""

from fractions import Fraction


# ---------------------------------------------------------------------------
# Exercise 1: Bayes' theorem in probability form
# ---------------------------------------------------------------------------


def posterior(prior: Fraction, sensitivity: Fraction, specificity: Fraction) -> Fraction:
    """P(condition | positive test), by Bayes' theorem.

    prior        = P(condition), the base rate before any test
    sensitivity  = P(positive | condition)
    specificity  = P(negative | no condition)

    The denominator -- P(positive), the "evidence" -- is exactly Day 113's
    law of total probability applied to the two-piece partition
    {condition, no condition}:

        P(positive) = P(condition) x sensitivity
                    + P(no condition) x (1 - specificity)
    """
    # YOUR CODE HERE
    raise NotImplementedError


def posterior_general(
    prior: Fraction,
    p_positive_given_condition: Fraction,
    p_positive_given_no_condition: Fraction,
) -> Fraction:
    """The same theorem, stated with raw likelihoods instead of sensitivity
    and specificity. Useful whenever "the test is positive" is not a clean
    sensitivity/specificity pair -- for instance, exercise 7's correlated
    tests, where P(both positive | condition) is not simply sensitivity
    squared.
    """
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 5: the odds form
# ---------------------------------------------------------------------------


def probability_to_odds(p: Fraction) -> Fraction:
    """odds = p / (1 - p)."""
    # YOUR CODE HERE
    raise NotImplementedError


def odds_to_probability(odds: Fraction) -> Fraction:
    """p = odds / (1 + odds), the inverse of probability_to_odds."""
    # YOUR CODE HERE
    raise NotImplementedError


def likelihood_ratio(sensitivity: Fraction, specificity: Fraction) -> Fraction:
    """LR+ = P(positive | condition) / P(positive | no condition)
           = sensitivity / (1 - specificity).
    """
    # YOUR CODE HERE
    raise NotImplementedError


def update_odds(prior_odds: Fraction, ratio: Fraction) -> Fraction:
    """posterior odds = prior odds x likelihood ratio."""
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 6: sequential updating
# ---------------------------------------------------------------------------


def sequential_posterior(
    prior: Fraction,
    tests: list[tuple[Fraction, Fraction]],
) -> Fraction:
    """Update a prior with a sequence of independent positive test results.

    `tests` is a list of (sensitivity, specificity) pairs. Update the
    running odds one likelihood ratio at a time, then convert back to a
    probability at the end.
    """
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 7: correlated tests -- when "multiply the likelihood ratios"
# quietly assumes something false
# ---------------------------------------------------------------------------


def independent_pair_probability(single_rate: Fraction) -> Fraction:
    """P(both runs agree on a given outcome | hypothesis), assuming the two
    runs are drawn independently: single_rate squared.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def correlated_pair_probability(single_rate: Fraction, correlation_weight: Fraction) -> Fraction:
    """P(both runs agree on a given outcome | hypothesis), when the two
    runs share a failure mode with probability `correlation_weight`.

    With probability `correlation_weight`, both runs are yoked to one
    shared random draw and therefore either BOTH show the outcome or
    NEITHER does, at the single-run rate. With probability
    (1 - correlation_weight), the two runs are genuinely independent, as
    `independent_pair_probability` assumes for the whole calculation.
    """
    # YOUR CODE HERE
    raise NotImplementedError

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

Both `examples/` and `starter/` contain modules called `bayes`, `simulate`,
`naive_bayes`, `dataset` 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 module was seen
first and then reuse it for the other suite -- so these starter tests would
silently pass against the reference solution instead of skipping. That is a
wrong answer with a green tick on it, which is the worst kind.

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

import sys
from pathlib import Path

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

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

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

"""The scenario, the corpus, and every tolerance this lab compares against.

Read this file. Nothing here is tuned to make a test pass: every exact
number is either a stated assumption (a test's sensitivity, a disease's
prevalence) or a value derived from those assumptions and then checked
against enumeration or simulation. The one place a captured figure differs
from the number a naive back-of-envelope calculation might suggest --
exercise 9's log-space arithmetic -- is called out explicitly at the bottom
of this file, because the wrong number belongs nowhere near a test.
"""

import math
from fractions import Fraction

# --------------------------------------------------------------------------
# The opening scenario: a diagnostic test for a rare condition
# --------------------------------------------------------------------------

#: 1 person in 1,000 has the condition, before any test is run.
PREVALENCE: Fraction = Fraction(1, 1000)

#: P(test positive | has condition) -- the test catches 99 sick people out
#: of every 100.
SENSITIVITY: Fraction = Fraction(99, 100)

#: P(test negative | does not have condition) -- the test correctly clears
#: 99 healthy people out of every 100.
SPECIFICITY: Fraction = Fraction(99, 100)

#: The exact posterior, derived by hand in the lesson and reproduced by
#: exercise 1: P(condition | positive) = 99/1098, about 9.02%.
OPENING_POSTERIOR_EXACT: Fraction = Fraction(99, 1098)

assert OPENING_POSTERIOR_EXACT == Fraction(11, 122)
assert round(float(OPENING_POSTERIOR_EXACT), 4) == 0.0902

# --------------------------------------------------------------------------
# The natural-frequencies table: 100,000 people, counted rather than
# multiplied as percentages
# --------------------------------------------------------------------------

NATURAL_FREQUENCY_POPULATION: int = 100_000

#: With prevalence 1/1000, exactly 100 of 100,000 people have the
#: condition, and 99,900 do not -- both exact because the population size
#: was chosen to divide the prevalence evenly.
NATURAL_FREQUENCY_SICK: int = 100
NATURAL_FREQUENCY_WELL: int = 99_900
assert NATURAL_FREQUENCY_SICK + NATURAL_FREQUENCY_WELL == NATURAL_FREQUENCY_POPULATION
assert Fraction(NATURAL_FREQUENCY_SICK, NATURAL_FREQUENCY_POPULATION) == PREVALENCE

#: True positives: 99% of the 100 sick people test positive.
NATURAL_FREQUENCY_TP: int = 99
#: False negatives: the other 1% of the sick people test negative.
NATURAL_FREQUENCY_FN: int = 1
#: False positives: 1% of the 99,900 healthy people test positive anyway.
NATURAL_FREQUENCY_FP: int = 999
#: True negatives: the other 99% of the healthy people correctly test negative.
NATURAL_FREQUENCY_TN: int = 98_901

assert NATURAL_FREQUENCY_TP + NATURAL_FREQUENCY_FN == NATURAL_FREQUENCY_SICK
assert NATURAL_FREQUENCY_FP + NATURAL_FREQUENCY_TN == NATURAL_FREQUENCY_WELL

# --------------------------------------------------------------------------
# Exercise 3: seeded simulation of a large population
# --------------------------------------------------------------------------

SIMULATION_POPULATION: int = 2_000_000
SIMULATION_SEED: int = 42


def standard_error(p: float, n: int) -> float:
    """The standard error of a proportion estimated from n trials."""
    return math.sqrt(p * (1.0 - p) / n)


# --------------------------------------------------------------------------
# Exercise 4: the prevalence sweep
# --------------------------------------------------------------------------

#: A sweep from rare to common, ending at 1/2 -- the prevalence at which
#: the posterior collapses to exactly the sensitivity, 0.99, which is the
#: number almost everyone wrongly gives for the 1-in-1,000 case above.
PREVALENCE_SWEEP: tuple[Fraction, ...] = (
    Fraction(1, 100_000),
    Fraction(1, 10_000),
    Fraction(1, 1_000),
    Fraction(1, 100),
    Fraction(1, 10),
    Fraction(1, 2),
)

# --------------------------------------------------------------------------
# Exercise 6: sequential updating with two DIFFERENT tests
# --------------------------------------------------------------------------

#: Test A is the opening scenario's test: 99% sensitive, 99% specific.
TEST_A_SENSITIVITY: Fraction = SENSITIVITY
TEST_A_SPECIFICITY: Fraction = SPECIFICITY

#: Test B is a different, less accurate test: 95% sensitive, 98% specific.
#: Using two genuinely different tests makes "the order does not matter"
#: a real claim about commutativity rather than a coincidence of running
#: the identical test twice.
TEST_B_SENSITIVITY: Fraction = Fraction(95, 100)
TEST_B_SPECIFICITY: Fraction = Fraction(98, 100)

# --------------------------------------------------------------------------
# Exercise 7: correlated tests -- the same assay, run twice, on one sample
# --------------------------------------------------------------------------

#: The probability that a run shares its outcome with the other run rather
#: than drawing independently -- modelling a shared failure mode, such as
#: a contaminated sample or a single faulty batch of reagent, that affects
#: both runs identically. c = 0 is the naive (fully independent) model;
#: c = 1/2 means half the time both runs are yoked to one shared draw.
CORRELATION_WEIGHT: Fraction = Fraction(1, 2)

#: Both runs use the same underlying test: 99% sensitive, 99% specific.
CORRELATED_SENSITIVITY: Fraction = SENSITIVITY
CORRELATED_SPECIFICITY: Fraction = SPECIFICITY

# --------------------------------------------------------------------------
# Exercise 8: Naive Bayes from scratch -- a tiny hand-made spam corpus
# --------------------------------------------------------------------------

#: Three spam documents, three ham documents. Kept deliberately tiny so the
#: whole vocabulary and every count can be read off by hand and checked.
SPAM_DOCS: tuple[str, ...] = (
    "buy cheap watches now",
    "cheap replica watches for sale",
    "buy now limited offer",
)

HAM_DOCS: tuple[str, ...] = (
    "meeting notes for review",
    "please review the agenda",
    "schedule the project meeting",
)

#: Held-out documents the trained classifier is asked to label.
#:
#: The first two have no word that is entirely absent from one class'
#: training vocabulary, so smoothed and unsmoothed classifiers agree.
#: The third is built to contain exactly one word -- "watches" -- that
#: never appears in the ham training documents, so the unsmoothed
#: classifier's P(document | ham) collapses to exactly zero and the single
#: word vetoes three other words that all point toward ham.
HELD_OUT_CLEAR_SPAM: str = "buy cheap watches"
HELD_OUT_CLEAR_HAM: str = "schedule the project meeting"
HELD_OUT_VETO_CASE: str = "please review schedule watches"

LAPLACE_ALPHA: int = 1

# --------------------------------------------------------------------------
# Exercise 9: why log space is not optional
# --------------------------------------------------------------------------

#: A stand-in for "several hundred small per-word probabilities multiplied
#: together" -- the exact shape of a naive Bayes document score. 500 factors
#: of 0.01 is well inside the range a real bag-of-words likelihood product
#: can reach, and float64 cannot represent the result.
UNDERFLOW_FACTOR: float = 0.01
UNDERFLOW_COUNT: int = 500

#: The true value of the corresponding sum of logs, computed directly
#: rather than approximated: 500 * ln(0.01). This is reported, not assumed
#: -- see the note below.
UNDERFLOW_LOG_SUM: float = UNDERFLOW_COUNT * math.log(UNDERFLOW_FACTOR)

# A note on a figure that does NOT appear as a constant here. An earlier
# draft of this lab's brief stated that 500 factors of 0.01 collapse in log
# space to "about -1151.29". That number is wrong for 500 factors of 0.01 --
# 500 * ln(0.01) = -2302.585..., not -1151.29. The figure -1151.29 is what
# you get from 500 factors of 0.1 (500 * ln(0.1) = -1151.29...), or
# equivalently 250 factors of 0.01. This file and every test in this lab
# use the measured, correct value for 500 factors of 0.01, UNDERFLOW_LOG_SUM
# above, computed directly by math.log rather than copied from a draft.
assert round(UNDERFLOW_LOG_SUM, 2) == -2302.59

"""Exercises 8 and 9: Naive Bayes from scratch, and why it must be built in
log space.

"Naive" names one specific, false assumption: that every word in a document
is conditionally independent of every other word, given the document's
class. Fill in the bodies marked `# YOUR CODE HERE`.
"""

import math
from collections import Counter
from dataclasses import dataclass
from fractions import Fraction


def tokenize(text: str) -> list[str]:
    """The whole tokenizer this lab needs: lowercase, split on whitespace."""
    return text.lower().split()


@dataclass(frozen=True)
class NaiveBayesModel:
    """Word counts, vocabulary and class priors, trained from labelled
    documents. Everything downstream -- smoothed or not, log space or not --
    is computed from these four fields alone.
    """

    classes: tuple[str, ...]
    vocabulary: frozenset[str]
    word_counts: dict[str, Counter]
    total_words: dict[str, int]
    doc_counts: dict[str, int]

    @property
    def total_docs(self) -> int:
        return sum(self.doc_counts.values())


def train(docs_by_class: dict[str, tuple[str, ...]]) -> NaiveBayesModel:
    """Count every word in every class's training documents.

    For each class: tokenize every document, accumulate word counts in a
    `Counter`, track the running total word count, and grow the shared
    vocabulary set with every token seen (in ANY class).
    """
    # YOUR CODE HERE
    raise NotImplementedError


def class_prior(model: NaiveBayesModel, cls: str) -> Fraction:
    """P(class) = documents in that class / total training documents."""
    # YOUR CODE HERE
    raise NotImplementedError


def word_probability(model: NaiveBayesModel, word: str, cls: str, alpha: int) -> Fraction:
    """P(word | class), with Laplace (add-alpha) smoothing.

    (count(word, class) + alpha) / (total_words[class] + alpha * |vocabulary|)

    alpha = 0 reproduces the unsmoothed textbook version: a word with zero
    occurrences in a class gets probability exactly 0.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def document_score(model: NaiveBayesModel, tokens: list[str], cls: str, alpha: int) -> Fraction:
    """The (unnormalised) joint probability P(class) x product of
    P(word | class) over every word in the document, computed as an exact
    Fraction by literal multiplication. Skip any token that never appeared
    in ANY training class (out-of-vocabulary).
    """
    # YOUR CODE HERE
    raise NotImplementedError


def classify(model: NaiveBayesModel, text: str, alpha: int) -> tuple[str, dict[str, Fraction]]:
    """Predict the class with the highest unnormalised joint probability.
    Returns the winning class and every class's raw score.
    """
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 9: log space, and why a product of many small floats is not
# optional to avoid
# ---------------------------------------------------------------------------


def multiply_probabilities(factors: list[float]) -> float:
    """The naive way: multiply plain Python floats together, one at a time."""
    # YOUR CODE HERE
    raise NotImplementedError


def sum_of_logs(factors: list[float]) -> float:
    """The fix: work in log space. Use `math.fsum` over `math.log(factor)`
    for each factor, rather than a running `+=` sum.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def document_log_score(model: NaiveBayesModel, tokens: list[str], cls: str, alpha: int) -> float:
    """The log-space equivalent of document_score: sum of logs instead of
    a product of probabilities.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def classify_log_space(model: NaiveBayesModel, text: str, alpha: int) -> tuple[str, dict[str, float]]:
    """The same classification decision as classify(), but computed in log
    space.
    """
    # YOUR CODE HERE
    raise NotImplementedError

"""Exercise 3: simulate a large population and confirm the exact posterior
by counting rather than by formula.

Takes an explicit `numpy.random.Generator`, built by
`numpy.random.default_rng(seed)`, exactly as Day 113's and Day 114's labs
do -- never the legacy `numpy.random.seed` global. Fill in the body marked
`# YOUR CODE HERE`.
"""

from typing import NamedTuple

import numpy as np


class PopulationCounts(NamedTuple):
    true_positive: int
    false_positive: int
    true_negative: int
    false_negative: int

    @property
    def positives(self) -> int:
        return self.true_positive + self.false_positive

    @property
    def empirical_posterior(self) -> float:
        """TP / (TP + FP) -- the fraction of positive results that are
        genuinely true positives, read straight off the simulated counts."""
        if self.positives == 0:
            return float("nan")
        return self.true_positive / self.positives


def simulate_population(
    rng: np.random.Generator,
    n: int,
    prevalence: float,
    sensitivity: float,
    specificity: float,
) -> PopulationCounts:
    """Draw n people, assign each a true condition status by `prevalence`,
    then a test result by `sensitivity`/`specificity`, and count the four
    outcome cells directly -- no formula involved anywhere in this function.

    Hints:
      - `rng.random(n) < prevalence` gives a boolean array of who has the
        condition.
      - Split the population into sick and well groups, draw test results
        for each group separately with `rng.random(count) < rate`, and
        count how many land on each side.
    """
    # YOUR CODE HERE
    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.

Nothing in here checks that a function exists or that a file is present.
Every test runs your code and compares a value.
"""

from fractions import Fraction

import numpy as np
import pytest

import answers
import bayes as B
import dataset as D
import naive_bayes as NB
import simulate as S

# --------------------------------------------------------------------------
# The skip machinery
# --------------------------------------------------------------------------


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


def need(value, what):
    """Skip if the exercise has not been attempted, otherwise hand it back."""
    if value is None:
        pytest.skip(f"not attempted yet: {what}")
    return value


def test_the_suite_itself_runs():
    """One test that always passes, so a green run is distinguishable from
    a collection error that quietly ran nothing at all."""
    assert D.PREVALENCE == Fraction(1, 1000)


# --------------------------------------------------------------------------
# Exercise 1 -- the opening posterior
# --------------------------------------------------------------------------


def test_1_posterior_matches_the_exact_fraction():
    result = attempt(lambda: B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert result == D.OPENING_POSTERIOR_EXACT


def test_1_posterior_returns_a_fraction_not_a_float():
    result = attempt(lambda: B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert isinstance(result, Fraction)


def test_1_posterior_is_not_the_naive_099_guess():
    result = attempt(lambda: B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert result != Fraction(99, 100)


def test_1_posterior_of_a_certain_prior_is_one():
    result = attempt(lambda: B.posterior(Fraction(1), D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert result == 1


# --------------------------------------------------------------------------
# Exercise 2 -- natural frequencies (uses dataset.py's captured constants,
# nothing to implement, but confirms you have read them)
# --------------------------------------------------------------------------


def test_2_natural_frequency_ratio_matches_the_formula():
    ratio = Fraction(D.NATURAL_FREQUENCY_TP, D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FP)
    exact = attempt(lambda: B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert ratio == exact


# --------------------------------------------------------------------------
# Exercise 3 -- simulation
# --------------------------------------------------------------------------


def test_3_simulate_population_counts_add_up():
    counts = attempt(
        lambda: S.simulate_population(np.random.default_rng(0), 20_000, 0.1, 0.9, 0.9),
        "simulate_population",
    )
    total = counts.true_positive + counts.false_positive + counts.true_negative + counts.false_negative
    assert total == 20_000


def test_3_simulate_population_within_tolerance_at_scale():
    counts = attempt(
        lambda: S.simulate_population(
            np.random.default_rng(D.SIMULATION_SEED),
            D.SIMULATION_POPULATION,
            float(D.PREVALENCE),
            float(D.SENSITIVITY),
            float(D.SPECIFICITY),
        ),
        "simulate_population",
    )
    exact = float(D.OPENING_POSTERIOR_EXACT)
    tol = 3.0 * D.standard_error(exact, counts.positives)
    assert abs(counts.empirical_posterior - exact) < tol


# --------------------------------------------------------------------------
# Exercise 4 -- the prevalence sweep
# --------------------------------------------------------------------------


def test_4_prevalence_one_half_gives_exactly_099():
    result = attempt(lambda: B.posterior(Fraction(1, 2), D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert result == Fraction(99, 100)


def test_4_prevalence_sweep_is_strictly_increasing():
    results = attempt(
        lambda: [B.posterior(p, D.SENSITIVITY, D.SPECIFICITY) for p in D.PREVALENCE_SWEEP],
        "posterior",
    )
    assert all(results[i] < results[i + 1] for i in range(len(results) - 1))


# --------------------------------------------------------------------------
# Exercise 5 -- the odds form
# --------------------------------------------------------------------------


def test_5_likelihood_ratio_is_exactly_99():
    result = attempt(lambda: B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY), "likelihood_ratio")
    assert result == 99


def test_5_posterior_odds_equal_prior_odds_times_ratio():
    prior_odds = attempt(lambda: B.probability_to_odds(D.PREVALENCE), "probability_to_odds")
    ratio = attempt(lambda: B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY), "likelihood_ratio")
    posterior_odds = attempt(lambda: B.update_odds(prior_odds, ratio), "update_odds")
    assert posterior_odds == prior_odds * ratio


def test_5_odds_form_matches_direct_posterior():
    prior_odds = attempt(lambda: B.probability_to_odds(D.PREVALENCE), "probability_to_odds")
    ratio = attempt(lambda: B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY), "likelihood_ratio")
    posterior_odds = attempt(lambda: B.update_odds(prior_odds, ratio), "update_odds")
    via_odds = attempt(lambda: B.odds_to_probability(posterior_odds), "odds_to_probability")
    direct = attempt(lambda: B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY), "posterior")
    assert via_odds == direct


# --------------------------------------------------------------------------
# Exercise 6 -- sequential updating
# --------------------------------------------------------------------------


def _tests_ab():
    return (D.TEST_A_SENSITIVITY, D.TEST_A_SPECIFICITY), (D.TEST_B_SENSITIVITY, D.TEST_B_SPECIFICITY)


def test_6_sequential_posterior_matches_the_documented_value():
    test_a, test_b = _tests_ab()
    result = attempt(lambda: B.sequential_posterior(D.PREVALENCE, [test_a, test_b]), "sequential_posterior")
    assert result == Fraction(1045, 1267)


def test_6_sequential_posterior_order_does_not_matter():
    test_a, test_b = _tests_ab()
    forward = attempt(lambda: B.sequential_posterior(D.PREVALENCE, [test_a, test_b]), "sequential_posterior")
    backward = attempt(lambda: B.sequential_posterior(D.PREVALENCE, [test_b, test_a]), "sequential_posterior")
    assert forward == backward


# --------------------------------------------------------------------------
# Exercise 7 -- correlated tests
# --------------------------------------------------------------------------


def test_7_naive_posterior_is_exactly_363_over_400():
    naive_tp = attempt(lambda: B.independent_pair_probability(D.CORRELATED_SENSITIVITY), "independent_pair_probability")
    naive_fp = attempt(
        lambda: B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY), "independent_pair_probability"
    )
    result = attempt(lambda: B.posterior_general(D.PREVALENCE, naive_tp, naive_fp), "posterior_general")
    assert result == Fraction(363, 400)


def test_7_naive_posterior_is_strictly_higher_than_correlated():
    naive_tp = attempt(lambda: B.independent_pair_probability(D.CORRELATED_SENSITIVITY), "independent_pair_probability")
    naive_fp = attempt(
        lambda: B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY), "independent_pair_probability"
    )
    naive = attempt(lambda: B.posterior_general(D.PREVALENCE, naive_tp, naive_fp), "posterior_general")

    corr_tp = attempt(
        lambda: B.correlated_pair_probability(D.CORRELATED_SENSITIVITY, D.CORRELATION_WEIGHT),
        "correlated_pair_probability",
    )
    corr_fp = attempt(
        lambda: B.correlated_pair_probability(1 - D.CORRELATED_SPECIFICITY, D.CORRELATION_WEIGHT),
        "correlated_pair_probability",
    )
    correlated = attempt(lambda: B.posterior_general(D.PREVALENCE, corr_tp, corr_fp), "posterior_general")
    assert naive > correlated


# --------------------------------------------------------------------------
# Exercise 8 -- Naive Bayes with Laplace smoothing
# --------------------------------------------------------------------------


def _trained_model():
    return NB.train({"spam": D.SPAM_DOCS, "ham": D.HAM_DOCS})


def test_8_clear_documents_classify_correctly():
    model = attempt(_trained_model, "train")
    spam_winner, _ = attempt(lambda: NB.classify(model, D.HELD_OUT_CLEAR_SPAM, alpha=D.LAPLACE_ALPHA), "classify")
    ham_winner, _ = attempt(lambda: NB.classify(model, D.HELD_OUT_CLEAR_HAM, alpha=D.LAPLACE_ALPHA), "classify")
    assert (spam_winner, ham_winner) == ("spam", "ham")


def test_8_veto_case_classifies_ham_with_smoothing():
    model = attempt(_trained_model, "train")
    winner, _ = attempt(lambda: NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=D.LAPLACE_ALPHA), "classify")
    assert winner == "ham"


def test_8_veto_case_ham_score_is_exactly_zero_without_smoothing():
    model = attempt(_trained_model, "train")
    _, scores = attempt(lambda: NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=0), "classify")
    assert scores["ham"] == 0


def test_8_word_probability_with_smoothing_never_zero():
    model = attempt(_trained_model, "train")
    prob = attempt(lambda: NB.word_probability(model, "watches", "ham", alpha=1), "word_probability")
    assert prob > 0


# --------------------------------------------------------------------------
# Exercise 9 -- log space
# --------------------------------------------------------------------------


def test_9_500_factors_of_001_underflow_to_exactly_zero():
    result = attempt(lambda: NB.multiply_probabilities([0.01] * 500), "multiply_probabilities")
    assert result == 0.0


def test_9_sum_of_logs_stays_finite():
    import math

    result = attempt(lambda: NB.sum_of_logs([0.01] * 500), "sum_of_logs")
    assert math.isfinite(result)


def test_9_sum_of_logs_matches_the_measured_value():
    result = attempt(lambda: NB.sum_of_logs([0.01] * 500), "sum_of_logs")
    assert result == D.UNDERFLOW_LOG_SUM


# --------------------------------------------------------------------------
# Predictions -- fifteen numbers, checked against the real answers
# --------------------------------------------------------------------------

EXPECTED: dict[str, object] = {
    "opening_posterior": float(D.OPENING_POSTERIOR_EXACT),
    "opening_posterior_below_naive_guess": True,
    "natural_frequency_total_positives": 1098.0,
    "simulation_within_tolerance": True,
    "prevalence_half_posterior": 0.99,
    "prevalence_sweep_increasing": True,
    "likelihood_ratio_value": 99.0,
    "odds_form_matches_direct": True,
    "sequential_two_test_posterior": float(Fraction(1045, 1267)),
    "sequential_order_independent": True,
    "correlated_naive_posterior": 0.9075,
    "correlated_naive_overstates": True,
    "veto_case_smoothed_class": "ham",
    "veto_case_unsmoothed_ham_is_zero": True,
    "underflow_to_exactly_zero": True,
}

HINTS: dict[str, str] = {
    "opening_posterior": (
        "Bayes' theorem: (prior x sensitivity) / evidence, where evidence "
        "is the law of total probability over {condition, no condition}."
    ),
    "correlated_naive_posterior": (
        "This is the SAME formula as exercise 6, just with sensitivity "
        "squared and (1 - specificity) squared as the two likelihoods."
    ),
    "veto_case_smoothed_class": (
        "Three of the four words in this document are ham-associated; "
        "only 'watches' points toward spam."
    ),
}


@pytest.mark.parametrize("key", sorted(EXPECTED))
def test_predictions(key):
    got = need(answers.ANSWERS.get(key), f"answers.ANSWERS[{key!r}]")
    want = EXPECTED[key]
    hint = HINTS.get(key, "")
    if isinstance(want, (bool, str)):
        assert got == want, f"{key}: your answer {got!r}, expected {want!r}. {hint}"
    else:
        assert abs(float(got) - float(want)) < 1e-6, (
            f"{key}: your answer {got!r}, expected {want!r}. {hint}"
        )


def test_every_answer_key_is_still_present():
    missing = sorted(set(EXPECTED) - set(answers.ANSWERS))
    assert not missing, f"answers.py is missing these keys: {missing}"

#!/usr/bin/env bash
# Tests for the Day 115 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:
#
#   * the opening posterior -- 99% sensitive, 99% specific, 1-in-1000
#     prevalence -- is exactly 99/1098 (about 0.0902), and NOT 0.99;
#   * the 100,000-person natural-frequency table's TP/(TP+FP) matches the
#     formula exactly;
#   * a 2,000,000-person seeded simulation lands within 3 standard errors
#     of the exact posterior;
#   * the posterior is strictly increasing in prevalence, and at a
#     prevalence of 1/2 it is EXACTLY 0.99 -- the naive guess, correct for
#     a different question;
#   * posterior odds equal prior odds times the likelihood ratio, exactly;
#   * two different tests, updated sequentially, give an identical
#     posterior regardless of the order they are applied in;
#   * a correlated pair of same-sample tests has a naive (assumes-
#     independence) posterior that is strictly, and substantially, higher
#     than the correct correlation-aware one;
#   * a from-scratch naive Bayes classifier with Laplace smoothing
#     correctly classifies held-out documents, and without smoothing a
#     single absent word collapses a class's probability to exactly zero;
#   * multiplying 500 factors of 0.01 underflows to exactly 0.0 in
#     float64, while the corresponding sum of logs stays finite;
#   * 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 115 — Bayes You Can Trust"
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_opening_posterior 02_natural_frequencies 03_simulation \
              04_prevalence_sweep 05_odds_form 06_sequential_updating \
              07_correlated_tests 08_naive_bayes_smoothing 09_log_space; 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 60 ]; then
  check "the reference suite ran at least 60 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 60 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 `bayes`,
# `simulate`, `naive_bayes`, `dataset` 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 math
from fractions import Fraction

import numpy as np

import bayes as B
import dataset as D
import naive_bayes as NB
import simulate as S

result = B.posterior(D.PREVALENCE, D.SENSITIVITY, D.SPECIFICITY)
print("opening_posterior", result)
print("opening_posterior_rounded", round(float(result), 4))
print("opening_posterior_not_naive", result != Fraction(99, 100))

total_positive = D.NATURAL_FREQUENCY_TP + D.NATURAL_FREQUENCY_FP
natural_ratio = Fraction(D.NATURAL_FREQUENCY_TP, total_positive)
print("natural_total_positive", total_positive)
print("natural_ratio_matches_formula", natural_ratio == result)

rng = np.random.default_rng(D.SIMULATION_SEED)
counts = S.simulate_population(
    rng, D.SIMULATION_POPULATION, float(D.PREVALENCE), float(D.SENSITIVITY), float(D.SPECIFICITY)
)
exact = float(result)
tol = 3.0 * D.standard_error(exact, counts.positives)
print("simulation_within_tolerance", abs(counts.empirical_posterior - exact) < tol)

sweep = [B.posterior(p, D.SENSITIVITY, D.SPECIFICITY) for p in D.PREVALENCE_SWEEP]
print("sweep_increasing", all(sweep[i] < sweep[i + 1] for i in range(len(sweep) - 1)))
print("sweep_half_is_099", B.posterior(Fraction(1, 2), D.SENSITIVITY, D.SPECIFICITY) == Fraction(99, 100))

prior_odds = B.probability_to_odds(D.PREVALENCE)
ratio = B.likelihood_ratio(D.SENSITIVITY, D.SPECIFICITY)
posterior_odds = B.update_odds(prior_odds, ratio)
print("odds_form_matches", posterior_odds == prior_odds * ratio)
print("odds_form_probability_matches_direct", B.odds_to_probability(posterior_odds) == result)

test_a = (D.TEST_A_SENSITIVITY, D.TEST_A_SPECIFICITY)
test_b = (D.TEST_B_SENSITIVITY, D.TEST_B_SPECIFICITY)
seq_ab = B.sequential_posterior(D.PREVALENCE, [test_a, test_b])
seq_ba = B.sequential_posterior(D.PREVALENCE, [test_b, test_a])
print("sequential_posterior", seq_ab)
print("sequential_order_independent", seq_ab == seq_ba)

naive_tp = B.independent_pair_probability(D.CORRELATED_SENSITIVITY)
naive_fp = B.independent_pair_probability(1 - D.CORRELATED_SPECIFICITY)
naive_post = B.posterior_general(D.PREVALENCE, naive_tp, naive_fp)
corr_tp = B.correlated_pair_probability(D.CORRELATED_SENSITIVITY, D.CORRELATION_WEIGHT)
corr_fp = B.correlated_pair_probability(1 - D.CORRELATED_SPECIFICITY, D.CORRELATION_WEIGHT)
corr_post = B.posterior_general(D.PREVALENCE, corr_tp, corr_fp)
print("correlated_naive_posterior", naive_post)
print("correlated_correct_posterior", corr_post)
print("correlated_naive_higher", naive_post > corr_post)

model = NB.train({"spam": D.SPAM_DOCS, "ham": D.HAM_DOCS})
smoothed_winner, _ = NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=D.LAPLACE_ALPHA)
_, unsmoothed_scores = NB.classify(model, D.HELD_OUT_VETO_CASE, alpha=0)
print("veto_smoothed_winner", smoothed_winner)
print("veto_unsmoothed_ham_zero", unsmoothed_scores["ham"] == 0)

factors = [0.01] * 500
product = NB.multiply_probabilities(factors)
logsum = NB.sum_of_logs(factors)
print("underflow_is_zero", product == 0.0)
print("logsum_finite", math.isfinite(logsum))
print("logsum_matches_measured", logsum == D.UNDERFLOW_LOG_SUM)
PY
)"

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

check_eq "the opening posterior is exactly 99/1098, reduced to 11/122" "11/122" "$(get opening_posterior)"
check_eq "which rounds to 0.0902" "0.0902" "$(get opening_posterior_rounded)"
check_eq "and is NOT the naive 0.99 guess" "True" "$(get opening_posterior_not_naive)"
check_eq "the natural-frequency table has 1098 total positives" "1098" "$(get natural_total_positive)"
check_eq "TP/(TP+FP) matches the formula's answer exactly" "True" "$(get natural_ratio_matches_formula)"
check_eq "the 2,000,000-person simulation lands within 3 standard errors" "True" "$(get simulation_within_tolerance)"
check_eq "the prevalence sweep is strictly increasing" "True" "$(get sweep_increasing)"
check_eq "at prevalence 1/2 the posterior is exactly 0.99" "True" "$(get sweep_half_is_099)"
check_eq "posterior odds equal prior odds times the likelihood ratio, exactly" "True" "$(get odds_form_matches)"
check_eq "the odds-form probability matches the direct formula exactly" "True" "$(get odds_form_probability_matches_direct)"
check_eq "the two-test sequential posterior is exactly 1045/1267" "1045/1267" "$(get sequential_posterior)"
check_eq "updating in either order gives an identical result" "True" "$(get sequential_order_independent)"
check_eq "the naive correlated-test posterior is exactly 363/400" "363/400" "$(get correlated_naive_posterior)"
check_eq "the correct correlated posterior is exactly 2189/13400" "2189/13400" "$(get correlated_correct_posterior)"
check_eq "the naive posterior is strictly higher than the correct one" "True" "$(get correlated_naive_higher)"
check_eq "the veto-case document classifies ham with smoothing" "ham" "$(get veto_smoothed_winner)"
check_eq "and ham's score is exactly zero without smoothing" "True" "$(get veto_unsmoothed_ham_zero)"
check_eq "500 factors of 0.01 underflow to exactly 0.0" "True" "$(get underflow_is_zero)"
check_eq "the corresponding sum of logs is finite" "True" "$(get logsum_finite)"
check_eq "and matches the independently measured value" "True" "$(get logsum_matches_measured)"

# --------------------------------------------------------------------------
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 naive-correlated-posterior comparison flipped -- 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_naive_posterior_is_strictly_higher_than_the_correlated_posterior" "${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 naive > correlated\n\n\ndef test_correlated_posterior_still_exceeds_a_single_tests_posterior"
replacement = (
    "    assert naive < correlated  # DELIBERATELY WRONG: self-test only\n\n\n"
    "def test_correlated_posterior_still_exceeds_a_single_tests_posterior"
)
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_naive_posterior_is_strictly_higher_than_the_correlated_posterior"*"failed"*|*"FAILED"*"test_naive_posterior_is_strictly_higher_than_the_correlated_posterior"*)
      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 ]