# Number of independent experiments for 80% power to detect KO 0.030 vs WT 0.035 per hour (model.py), using the
# paired t-test on per-experiment growth rates fitted over 0-72 h. The SD of the per-experiment difference is
# estimated from 20,000 simulated experiments; power then follows the noncentral t distribution.
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import experiment, TIMES
rng = np.random.default_rng(7)
def rate(x): m = np.log(x).mean(0); k = TIMES <= 72; return np.polyfit(TIMES[k], m[k], 1)[0]
d = np.array([rate(experiment(rng, 0.030)) - rate(experiment(rng, 0.035)) for _ in range(20_000)])
print(f"Mean difference in fitted rate {d.mean():.5f} per hour; SD across experiments {d.std(ddof=1):.5f}")
for n in range(3, 13):
    df = n - 1; ncp = d.mean() / (d.std(ddof=1) / np.sqrt(n)); tc = stats.t.ppf(0.975, df)
    print(f"{n} experiments: power {stats.nct.cdf(-tc, df, ncp) + stats.nct.sf(tc, df, ncp):.1%}")