# 10,000 simulated experiments per row: vehicle and three doses, 6 mice per group, SD 4, alpha 0.05.
# "Any dose vs vehicle" = at least one Dunnett-adjusted p < 0.05; "trend" = linear contrast p < 0.05.
import numpy as np
from scipy import stats
rng = np.random.default_rng(7)
c = np.array([-3, -1, 1, 3])
def run(step, runs=10_000):
    hit_d = hit_t = 0
    for _ in range(runs):
        g = [100 - step * k + rng.normal(0, 4, 6) for k in range(4)]
        hit_d += stats.dunnett(*g[1:], control=g[0], rng=rng).pvalue.min() < 0.05
        mse = sum(((x - x.mean()) ** 2).sum() for x in g) / 20
        t = (c * [x.mean() for x in g]).sum() / np.sqrt(mse * (c ** 2 / 6).sum())
        hit_t += 2 * stats.t.sf(abs(t), 20) < 0.05
    return hit_d / runs, hit_t / runs
for step, label in [(0.0, "no effect (false positives)"), (2.2, "falls 2.2 per dose step (power)")]:
    d_, t_ = run(step)
    print(f"{label}: any dose vs vehicle (Dunnett) {d_:.1%}, linear trend {t_:.1%}")