# How often each approach concludes "the knockout changes anxiety" (p < 0.05): 5,000 simulated studies of 10 + 10 mice,
# with no effect, with less anxiety (latent level +0.8 or +1.2 SD) and with more activity only (+30%, anxiety unchanged).
# Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np
from scipy import stats
from model import study, COLS
rng = np.random.default_rng(20261112); R = 5_000
ANX = [COLS.index(c) for c in ["pct_open_time", "pct_open_entries", "open_entries", "pct_centre_time", "centre_entries"]]
Z = [COLS.index(c) for c in ["pct_open_time", "pct_open_entries", "pct_centre_time"]]
COUNTS = [COLS.index(c) for c in ["open_entries", "centre_entries"]]

def one(anx, act):
    s = study(rng, anx, act); x = np.array([r[2:] for r in s], float); g = np.array([r[0] == "knockout" for r in s])
    wt, ko = x[~g], x[g]
    p = stats.ttest_ind(ko, wt, axis=0, equal_var=False).pvalue
    pa = p[ANX]; holm = min(1.0, pa.min() * len(pa))         # smallest Holm-adjusted p = smallest p x number of tests
    z = ((x[:, Z] - wt[:, Z].mean(0)) / wt[:, Z].std(0, ddof=1)).mean(1)
    pd_ = stats.ttest_ind(z[g], z[~g], equal_var=False).pvalue
    return np.array([(pa < 0.05).any(), holm < 0.05, p[COLS.index("pct_open_time")] < 0.05, pd_ < 0.05, (p[COUNTS] < 0.05).any()])
for lab, anx, act in (("No effect", 0.0, 1.0), ("Less anxious by 0.8 SD", 0.8, 1.0), ("Less anxious by 1.2 SD", 1.2, 1.0), ("30% more active, anxiety unchanged", 0.0, 1.3)):
    r = np.array([one(anx, act) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A any of 5 anxiety measures {r[0]:.1%}; B Holm over 5 {r[1]:.1%}; C primary % open-arm time {r[2]:.1%}; "
          f"D composite of 3 percentages {r[3]:.1%}; count measures (open-arm or centre entries) {r[4]:.1%}")