Lab 1 — Simulate it: how tests behave

Author

Hicham Janati — FSML

This lab is about mastering (1) simulation with numpy and (2) statistical testing.

How this lab works. Four kinds of boxes:

🔮 Predict: write your guess before running the next cell. Being wrong is fine; guessing first is the point. Never use AI here.

🧠 Explain: a result just appeared; say in your own words what it means.

✍️ Your turn: a short coding task;

🤖 With AI: a prompt that sends you to find something the lab did not hand you, a concept or a piece of Python. AI is welcome for syntax and debugging throughout;

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

BLUE, RED, GREEN, GRAY = "#1f6fb2", "#c92a1a", "#1c7a43", "#8a8a8a"
plt.rcParams["figure.figsize"] = (6.5, 3.6)
plt.rcParams["axes.spines.top"] = plt.rcParams["axes.spines.right"] = False

Part 1 — Fake a thousand seasons

The lecture’s duel: Karim beat Yassine 15 games out of 20. Skill, or luck? H0: equals (θ = ½), H1: Karim is better (θ > ½). The lecture computed P(S ≥ 15) = 0.021 with the Binomial formula. Here we watch it happen.

1.1 Making randomness

A random-number generator is an object you create once, with a seed:

seed = 42
rng = np.random.default_rng(seed)
rng.random(5)

Five numbers, uniform between 0 and 1. A game Karim wins with probability ½ is a number below 0.5:

rng.random(5) < 0.5

🧠 Explain. Run the cell rng.random(5) above a second time, alone. Are the five numbers identical to the first run? Why?

✍️ Your answer:

…

Now put the seed line and the draw in the same cell:

seed = 42
rng = np.random.default_rng(seed)
rng.random(5)

🧠 Explain. Run this cell several times. Are the five numbers identical each time? Why?

✍️ Your answer:

…

🧠 Explain. Change seed = 42 to another value in the cell above and run it a few times. What happened? Say in one sentence what the seed is for.

✍️ Your answer:

…

✍️ Your turn. One season is 20 games. Draw 20 games with rng.random(20) < 0.5 and count Karim’s wins with .sum(). Run it a few times: does 15 ever show up?

# TODO: one season -> number of wins S

1.2 Ten thousand seasons at once

Instead of a loop, generate a whole table of games: 10 000 rows (seasons) × 20 columns (games), and sum each row.

✍️ Your turn. Compute S, Karim’s number of wins in each of the 10 000 seasons, by summing games along the right axis.

rng = np.random.default_rng(0)
games = rng.random((10_000, 20)) < 0.5      # True = Karim wins
# TODO: S = ...  (wins per season)
print(S.shape, S[:10])

🔮 Predict. These are 10 000 seasons between two equal players. What will the histogram of S look like: symmetric? centered around which value? roughly how wide?

✍️ Your answer:

…

plt.hist(S, bins=np.arange(-0.5, 21.5, 1), color=GRAY)
plt.axvline(14.5, color=RED, lw=2, label="S ≥ 15")
plt.xlabel("Karim's wins in a season of 20 (equal players)")
plt.ylabel("number of seasons"); plt.legend()
plt.show()

The exact counterpart of a simulated fraction is the cumulative distribution function (CDF, fonction de répartition): for the Binomial, stats.binom.cdf(k, n, θ) returns P(S ≤ k) for n games won with probability θ.

✍️ Your turn. Compute the fraction of seasons with S ≥ 15, and the same for S ≥ 17. Then verify both against the CDF.

# TODO: empirical_15 = ...    exact_15 = ...
# TODO: empirical_17 = ...    exact_17 = ...

print(f"S >= 15: simulated {empirical_15:.4f}   exact {exact_15:.4f}")
print(f"S >= 17: simulated {empirical_17:.4f}   exact {exact_17:.4f}")

🧠 Explain, step 1. Do simulation and formula agree? Rerun the table with 100 000 seasons instead of 10 000 and say why the agreement changes the way it does.

✍️ Your answer:

…

🧠 Explain, step 2. Assume both players really are the same level and you agreed beforehand on a significance threshold of 5 %. Computing P(S ≥ 14) as well, which score S makes you declare Karim the better player, and how often would you wrongfully declare it?

✍️ Your answer:

…

That frequency, declaring a difference when there is none, is the false-alarm rate, the test’s type I error. It is the thing a test promises to control: whatever the true θ ≤ ½, the rule “reject when S ≥ 15” fires in at most 5 % of seasons (here 2.1 %, because S is discrete and 14 would already be too lenient).

1.3 Misreading: “not significant, so they are equals”

The lecture’s other case: suppose the season had gone 12–8. Then P(S ≥ 12) = 0.252, nothing to reject. Yassine concludes: “see, we are equals”. Test his claim with the same tool, step by step.

🔮 Step 1. How many games out of 20 does Karim need to win to reach 5 % significance? (You found it in step 2 above.)

✍️ Your answer:

…

🔮 Step 2. Now suppose Karim really is better: he wins each game with probability 0.55. What are the chances of getting such a season (reaching the score of step 1) in 20 games?

✍️ Your answer:

…

✍️ Your turn. Rewrite section 1.2’s three lines from memory with 0.55 instead of 0.5, and compute the fraction of seasons with S ≥ 15. Verify with the CDF.

# TODO: S_better to be defined as the wins per season over 10 000 seasons with theta = 0.55, then the fraction with S_better >= 15
# TODO: S_better = ...
# TODO: print(...)

🧠 Explain. A genuinely better player clears the bar in barely one season out of twenty, so a 12–8 result is his normal outcome. What did “we fail to reject H0” actually prove about the two players, and was Yassine’s conclusion wrong?

✍️ Your answer:

…

Part 2 — The formula, at will

Simulation confirmed the formula; from here on we use the formula alone: 1 - stats.binom.cdf(k - 1, n, θ) gives P(S ≥ k) for any n and θ instantly.

2.1 Same 75 %, ten times more games

🔮 Predict. Karim wins 15 of 20: p = 0.021. Suppose the same 75 % win rate held over 200 games (150 of 200). Will the new p-value be smaller than 0.021, and by a little or by a lot?

✍️ Your answer:

…

print(f"15 of 20  : p = {1 - stats.binom.cdf(14, 20, 0.5):.4f}")
print(f"150 of 200: p = {1 - stats.binom.cdf(149, 200, 0.5):.2e}")

2.2 The rule, under both hypotheses

Fix α = 5 %. The rule “reject when S ≥ k” needs the smallest k whose false-alarm rate is at most 5 %:

def critical_k(n, alpha=0.05):
    """Smallest k such that P(S >= k) <= alpha when theta = 1/2."""
    k = np.arange(n + 1)
    return int(k[1 - stats.binom.cdf(k - 1, n, 0.5) <= alpha][0])

for n in [20, 200]:
    print(f"n = {n:3d}: reject H0 when S >= {critical_k(n)}")

Now suppose Karim’s true win probability were actually 70 %. How often would this same rule correctly declare him better?

print(f"P(S >= 15 | theta = 0.7) = {1 - stats.binom.cdf(14, 20, 0.7):.3f}")

Even when Karim really wins 70 % of games, a 20-game season detects it only about 42 % of the time. That probability, rejecting H0 when H1 is true, is the test’s power. Type I error is what happens under H0; power is what happens under H1; a good rule keeps the first below α and pushes the second toward 1.

✍️ Your turn. Same two numbers for n = 200 (its critical k is printed above): P(reject | θ = 0.5) and P(reject | θ = 0.7), with the CDF. Before running: which one must stay below 0.05, and how close to 1 do you expect the other?

# TODO: level and power to be defined as P(reject | theta = 0.5) and P(reject | theta = 0.7) at n = 200
# TODO: power = ...
print(f"n = 200, reject when S >= 113: level {level:.4f}   power {power:.4f}")

🧠 Explain. Compare the two seasons: what does going from 20 to 200 games improve, and what does it not change?

✍️ Your answer:

…

2.3 How many games to catch a 55 % player?

Karim is better, but only slightly: θ = 0.55. Part 1.3 showed 20 games have almost no chance of showing it.

🔮 Predict. Will the 5 % rule need more or fewer than 200 games to detect a 55 % player with probability at least 80 %? Why?

✍️ Your answer:

…

✍️ Your turn. Fill the loop: for each n, get the critical k, compute the power at θ = 0.55, and stop at the first n where it reaches 0.8.

for n in range(20, 3000, 10):
    # TODO: k = ...
    # TODO: power = ...
    # TODO: if power >= 0.8: print(n, power); break
    pass

🧠 Explain. Hundreds of games to see a 5-point edge. Write down the consequence for Part 1.3: with 20 games, what does a “not significant” verdict say about small differences?

✍️ Your answer:

…

🧠 Explain. Now the mirror image: with thousands of games, what kind of edge becomes detectable, and so what does “significant” tell you about the size of the edge?

✍️ Your answer:

…

Part 3 — Many tests at once

Yassine is now the club champion. One hundred challengers each play him a season of 100 games (enough for a test to have real power: run critical_k(100) to see the 5 % bar), and each one claims to be better. One test per challenger, at 5 %. What could go wrong?

3.1 One hundred equal challengers

🔮 Predict. Suppose every one of the 100 challengers is exactly Yassine’s level. How many of them will the 5 % rule declare “better”?

✍️ Your answer:

…

n, m = 100, 100                                   # games per duel, number of challengers
rng = np.random.default_rng(7)
S_all = (rng.random((m, n)) < 0.5).sum(axis=1)    # one score per challenger
pvals = 1 - stats.binom.cdf(S_all - 1, n, 0.5)    # one p-value per challenger
print(pvals[:10].round(3))

✍️ Your turn. rejected should be a True/False array, one entry per challenger: True when the p-value is at most 0.05. Count how many were rejected (.sum()) and print their scores S_all[rejected].

# TODO: rejected = ...
# TODO: print(rejected.sum(), S_all[rejected])

🧠 Explain. A hundred players of equal level, and several were declared better. Did any single test break its promise, and if not, where does the problem come from?

✍️ Your answer:

…

This is the multiple testing problem. Each test controls its own false-alarm rate α; nobody controls what happens to the family of 100 of them. The lecture’s exercise computed that family’s risk, the probability that at least one of m independent tests fires by mistake: 1 − (1 − α)^m, which for α = 5 % and m = 100 is 0.994. Near-certainty.

3.2 Bonferroni: divide α by m

The lecture’s fix: run each test at level α/m instead of α. The union bound then caps the probability of any false rejection at α.

🔮 Predict. Apply it to this tournament: 100 tests at 5 %, so what threshold must each p-value beat, and what score out of 100 does that demand?

✍️ Your answer:

…

✍️ Your turn. Redo 3.1’s count with the corrected threshold, on the same pvals. critical_k from Part 2 gives the new score bar.

alpha_bonf = 0.05 / m
print(f"Bonferroni bar: S >= {critical_k(n, alpha_bonf)} of {n}")
# TODO: rejected_bonf = ...  (on pvals, then .sum())
print(f"rejected in this tournament: {rejected_bonf}")

The price. Now let 20 of the 100 challengers really be better: they win each game with probability 0.65. The report tells the truth we never see in practice: who was declared better, how many are false positives, and how many of the 20 genuinely better players were found.

🔮 Predict. Of the 20 better players, which rule finds more, the 5 % rule or Bonferroni? Which one wrongly declares more equal players better?

✍️ Your answer:

…

theta = np.full(m, 0.5)
theta[:20] = 0.65                                  # challengers 0-19 really are better
better = theta > 0.5
rng = np.random.default_rng(9)
S_mix = (rng.random((m, n)) < theta[:, None]).sum(axis=1)
p_mix = 1 - stats.binom.cdf(S_mix - 1, n, 0.5)

def report(name, rejected):
    print(f"{name:12s} rejected {rejected.sum():2d}: {(rejected & ~better).sum()} false positives, "
          f"{(rejected & better).sum():2d} of the 20 better players found")

report("5 % rule", p_mix <= 0.05)
report("Bonferroni", p_mix <= 0.05 / m)

🧠 Explain. Bonferroni kept its promise and paid for it. In the words of Part 2: what happened to the power of each individual test when α became α/m?

✍️ Your answer:

…

3.3 A softer promise: the false discovery rate

Bonferroni promises “almost never a single false rejection”, and in a tournament of 100 nobody asked for that. A more useful promise: among the players we declare better, only a small fraction should in fact be equal. That fraction, false rejections divided by all rejections, is the false discovery rate (FDR), and controlling it instead of the family-wise risk is what restores the power.

The Benjamini–Hochberg (BH) procedure controls the FDR at level α for independent tests (no proof here):

  1. sort the p-values: \(p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}\);
  2. find the largest \(i\) such that \(p_{(i)} \le \dfrac{i}{m}\,\alpha\);
  3. reject the \(i\) tests with the smallest p-values.

The smallest p-value still has to beat α/m, Bonferroni’s bar, but the second only α · 2/m, the tenth α · 10/m, and so on: the bar rises with the rank. Many small p-values means many real effects, and BH spends that evidence. Here is the picture, the sorted p-values of the mixed tournament against the three bars:

def bh_plot(pvals, alpha=0.05, show=30):
    i = np.arange(1, show + 1)
    plt.plot(i, np.sort(pvals)[:show], "o", color=GRAY, label="sorted p-values")
    plt.plot(i, i * alpha / len(pvals), color=GREEN, lw=2, label="BH line  i·α/m")
    plt.axhline(alpha / len(pvals), color=RED, ls="--", label="Bonferroni  α/m")
    plt.axhline(alpha, color=BLUE, ls=":", label="5 % rule  α")
    plt.yscale("log"); plt.xlabel("rank i"); plt.ylabel("p-value (log scale)"); plt.legend()
    plt.show()

bh_plot(p_mix)

🧠 Explain. Read the picture: how many dots sit under the rising green BH line, and how many under the flat red Bonferroni line?

✍️ Your answer:

…

🤖 With AI. Implement the BH procedure in Python and apply it to p_mix below, then read the result against the picture and against Bonferroni’s results.

# TODO: import the function the assistant named

# TODO: rejected, qvalues, _, _ = ...

report("BH", rejected)
print("smallest adjusted p-values:", np.sort(qvalues)[:5].round(4))

✍️ Your answer:

…