# Checks the Bonferroni-based sample sizes against Dunnett's test, the analysis actually planned.
# Means 100, 95, 90, 85 (vehicle, low, mid, high), SD 12; 10,000 simulated experiments per row.
# Python 3.13, NumPy, SciPy (scipy.stats.dunnett, with the same generator passed for its numerical integration). Seed fixed before running.
import numpy as np
from scipy import stats
rng = np.random.default_rng(20261111); R = 10_000; mu = (100, 95, 90, 85)
for n, which, lab in ((15, 2, "high dose (difference 15)"), (31, 1, "mid dose (difference 10)"), (6, 2, "high dose (difference 15)")):
    hits = 0
    for _ in range(R):
        g = [rng.normal(m, 12, n) for m in mu]
        hits += stats.dunnett(*g[1:], control=g[0], rng=rng).pvalue[which] < 0.05
    print(f"n = {n} per group: Dunnett detects the {lab} in {hits / R:.1%} of experiments")