# 10,000 simulated experiments: 12 mice per group, same model as the worked example, mice with < 20 s total
# exploration excluded. Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import simulate
rng = np.random.default_rng(20261111); R = 10_000
def one(mw, mk):
    out = []
    for m in (mw, mk):
        nov, fam = simulate(rng, 12, m); tot = nov + fam; keep = tot >= 20
        out.append(((nov - fam) / tot)[keep])
    w, k = out
    above = [stats.ttest_1samp(x, 0).pvalue < 0.05 and x.mean() > 0 for x in (w, k)]   # 'above chance': p < 0.05 and mean DI > 0
    return above[0] != above[1], stats.ttest_ind(k, w, equal_var=False).pvalue < 0.05
for lab, mw, mk in (("Same true DI in both groups (0.20)", 0.20, 0.20), ("No memory in either group (true DI 0)", 0.0, 0.0), ("KO lower (0.25 vs 0.05)", 0.25, 0.05)):
    r = np.array([one(mw, mk) for _ in range(R)])
    print(f"{lab}: one group 'above chance' and the other not {r[:, 0].mean():.1%}; Welch's t-test between groups p<0.05 {r[:, 1].mean():.1%}")