# 10,000 simulated experiments with NO treatment effect: 20 mice, outcome correlated with body weight (r = 0.5).
# (a) Random allocation: 10 mice to each group at random.
# (b) "Haphazard" allocation: mice are taken from the cage one by one and the first 10 caught go to treatment;
#     heavier, less active mice tend to be caught first (catch order correlated with weight, rank r about 0.5).
# Share of experiments with Welch's t-test p < 0.05.
import numpy as np
from scipy import stats
rng = np.random.default_rng(4)
runs, hit_r, hit_h = 10_000, 0, 0
for _ in range(runs):
    w = rng.normal(0, 1, 20)
    y = 0.5 * w + rng.normal(0, np.sqrt(1 - 0.25), 20)
    idx = rng.permutation(20)
    hit_r += stats.ttest_ind(y[idx[:10]], y[idx[10:]], equal_var=False).pvalue < 0.05
    order = np.argsort(-(0.58 * w + rng.normal(0, 1, 20)))          # catch order: heavier first, with noise
    hit_h += stats.ttest_ind(y[order[:10]], y[order[10:]], equal_var=False).pvalue < 0.05
print(f"runs = {runs}; false positives with no treatment effect: random allocation {hit_r / runs:.1%}, first-caught to treatment {hit_h / runs:.1%}")