# 10,000 simulated experiments per scenario, 10 mice per group, same model as the worked example (model.py).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import simulate, DAYS
rng = np.random.default_rng(20261111); R = 10_000; n = 10

def mw(a, b):
    return 1.0 if len(set(a) | set(b)) == 1 else stats.mannwhitneyu(b, a, method="asymptotic").pvalue

for label, drug in (("No drug effect", (0.9, 3.0, 0)), ("Drug effect as in the worked example", (0.7, 2.2, 1))):
    hits = np.zeros(3)
    for _ in range(R):
        v, g = simulate(rng, n, 0.9, 3.0), simulate(rng, n, *drug)
        hits[0] += mw(v.sum(1), g.sum(1)) < 0.05
        hits[1] += mw(v.max(1), g.max(1)) < 0.05
        hits[2] += min(mw(v[:, j], g[:, j]) for j in range(len(DAYS))) < 0.05
    print(f"{label}: p<0.05 for cumulative score {hits[0] / R:.1%}, maximum score {hits[1] / R:.1%}, "
          f"at least one of {len(DAYS)} daily tests {hits[2] / R:.1%}")