# How often each analysis gives p < 0.05: 5,000 simulated studies of 8 + 8 mice with 9 scored fields each, with no
# treatment effect and with the treatment lowering the underlying inflammation by 1 SD of the between-mouse variation.
# Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np, warnings
from scipy import stats, optimize
from model import study, MICE, FIELDS
warnings.simplefilter("ignore")
rng = np.random.default_rng(20261112); R = 5_000

def polr_p(score, x):
    lev = np.sort(np.unique(score)); y = np.searchsorted(lev, score); K = len(lev)
    if K < 2: return 1.0
    def nll(p, use_x):
        zeta = np.cumsum(np.r_[p[0], np.exp(p[1:K - 1])]); eta = p[K - 1] * x if use_x else 0 * x
        cum = np.column_stack([np.zeros(len(y)), 1 / (1 + np.exp(-(zeta[None, :] - eta[:, None]))), np.ones(len(y))])
        return -np.sum(np.log(np.clip(cum[np.arange(len(y)), y + 1] - cum[np.arange(len(y)), y], 1e-300, None)))
    p0 = np.r_[-1.0, np.zeros(K - 2)]
    f1 = optimize.minimize(nll, np.r_[p0, 0.0], args=(True,), method="BFGS"); f0 = optimize.minimize(lambda p: nll(np.r_[p, 0.0], False), p0, method="BFGS")
    return stats.chi2.sf(max(2 * (f0.fun - f1.fun), 0), 1)

def one(shift):
    s = np.array([r[3] for r in study(rng, shift)], float).reshape(2, MICE, FIELDS)   # [group, mouse, field]
    pa = stats.mannwhitneyu(s[1].ravel(), s[0].ravel(), alternative="two-sided", method="asymptotic").pvalue
    med = np.median(s, axis=2); mean = s.mean(axis=2)
    pb = stats.ttest_ind(med[1], med[0], equal_var=False).pvalue
    pc = stats.mannwhitneyu(med[1], med[0], alternative="two-sided", method="asymptotic").pvalue
    pd_ = polr_p(med.ravel(), np.r_[np.zeros(MICE), np.ones(MICE)])
    pe = stats.ttest_ind(mean[1], mean[0], equal_var=False).pvalue
    return np.array([pa, pb, pc, pe, pd_]) < 0.05
for lab, shift in (("No effect", 0.0), ("Treatment lowers inflammation by 1 SD", -1.0)):
    r = np.array([one(shift) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: p<0.05 with A fields as n {r[0]:.1%}; B mouse medians, t-test {r[1]:.1%}; C mouse medians, Mann-Whitney {r[2]:.1%}; "
          f"D mean of fields per mouse, t-test {r[3]:.1%}; E proportional odds on mouse medians {r[4]:.1%}")