# 10,000 simulated experiments per scenario, 12 mice per group, same model as the worked example (model.py).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats, optimize
from model import simulate, LOD
rng = np.random.default_rng(20261111); R = 10_000; n = 12

def mle(y, ev, xg):
    def nll(q):
        z = (y - q[0] - q[1] * xg) / np.exp(q[2])
        return -(ev * (stats.norm.logpdf(z) - q[2]) + (1 - ev) * stats.norm.logcdf(z)).sum()
    f = optimize.minimize(nll, [np.mean(y), 0, np.log(0.8)], method="BFGS")
    E = np.eye(3) * 1e-4; p0 = f.x
    H = np.array([[(nll(p0 + a + b) - nll(p0 + a - b) - nll(p0 - a + b) + nll(p0 - a - b)) / 4e-8 for b in E] for a in E])
    return p0[1], 2 * stats.norm.sf(abs(p0[1] / np.sqrt(np.linalg.inv(H)[1, 1])))

for label, mc, mt in (("Moderate censoring, no effect (medians 8 and 8)", 8, 8), ("Moderate censoring, twofold effect (8 and 16)", 8, 16),
                      ("Heavy censoring, no effect (medians 5 and 5)", 5, 5), ("Heavy censoring, twofold effect (5 and 10)", 5, 10)):
    sig = np.zeros(3); rs = np.zeros(2); ratios = [[], []]; sdc = []; cens = 0
    for _ in range(R):
        a, ea = simulate(rng, n, mc); b, eb = simulate(rng, n, mt)
        la, lb = np.log(np.where(ea == 1, a, LOD / 2)), np.log(np.where(eb == 1, b, LOD / 2))
        sig[0] += stats.ttest_ind(lb, la, equal_var=False).pvalue < 0.05
        sig[1] += stats.mannwhitneyu(b, a, method="asymptotic").pvalue < 0.05
        y = np.log(np.r_[a, b]); ev = np.r_[ea, eb]; xg = np.r_[np.zeros(n), np.ones(n)]
        b1, p = mle(y, ev, xg); sig[2] += p < 0.05
        ratios[0].append(np.exp(lb.mean() - la.mean())); ratios[1].append(np.exp(b1)); sdc.append(la.std(ddof=1)); cens += (1 - ea).sum()
    print(f"{label}: control values below the LOD {cens / (R * n):.0%}; p<0.05 with LOD/2 substitution {sig[0] / R:.1%}, "
          f"Mann-Whitney {sig[1] / R:.1%}, left-censored model {sig[2] / R:.1%}; median ratio estimate: substitution "
          f"{np.median(ratios[0]):.2f}, model {np.median(ratios[1]):.2f}; median control log-SD after substitution {np.median(sdc):.2f} (true 0.80)")