# How often each analysis gives p < 0.05, and the average estimated difference in log10 CFU: 5,000 simulated studies
# of 8 + 8 mice, with no effect (burden 3.0 log10), no effect at a low burden (1.8 log10, more mice below the limit
# of detection), and a treatment that lowers burden by 1.5 log10 (3.0 to 1.5). Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np, warnings
from scipy import stats, optimize
from model import study
warnings.simplefilter("ignore")
rng = np.random.default_rng(20261112); R = 5_000; LOD = 10.0

def tobit(y, cens, x):
    def nll(p, use_x):
        mu = p[0] + (p[1] * x if use_x else 0); s = np.exp(p[-1])
        return -(np.sum(stats.norm.logpdf(y[~cens], (mu[~cens] if use_x else mu), s)) + np.sum(stats.norm.logcdf(((y[cens] - (mu[cens] if use_x else mu)) / s))))
    f1 = optimize.minimize(nll, [y.mean(), 0, np.log(y.std() + 0.1)], args=(True,), method="BFGS")
    f0 = optimize.minimize(nll, [y.mean(), np.log(y.std() + 0.1)], args=(False,), method="BFGS")
    return f1.x[1], stats.chi2.sf(max(2 * (f0.fun - f1.fun), 0), 1)

def one(mc, mt):
    rows = study(rng, mc, mt); g = np.array([r[0] == "treated" for r in rows]); dil = np.array([r[2] for r in rows]); col = np.array([r[3] for r in rows], float)
    cfu = col * 10.0 ** dil / 0.1; below = col == 0
    pa = stats.ttest_ind(cfu[g], cfu[~g], equal_var=False).pvalue
    lh = np.log10(np.where(below, LOD / 2, cfu)); pb = stats.ttest_ind(lh[g], lh[~g], equal_var=False).pvalue; db = lh[g].mean() - lh[~g].mean()
    pc = stats.mannwhitneyu(cfu[g], cfu[~g], alternative="two-sided", method="asymptotic").pvalue
    dd, pd_ = tobit(np.log10(np.where(below, LOD, cfu)), below, g.astype(float)) if below.sum() < len(below) - 1 else (np.nan, 1.0)
    return np.array([pa < 0.05, pb < 0.05, pc < 0.05, pd_ < 0.05, db, dd, below.mean()])
for lab, mc, mt in (("No effect, burden 3.0 log10", 3.0, 3.0), ("No effect, burden 1.8 log10", 1.8, 1.8), ("Treatment lowers burden by 1.5 log10", 3.0, 1.5)):
    r = np.nanmean(np.array([one(mc, mt) for _ in range(R)]), axis=0)
    print(f"{lab} ({r[6]:.0%} of mice below the LOD): p<0.05 with A raw CFU {r[0]:.1%}; B log10, LOD/2 {r[1]:.1%}; C Mann-Whitney {r[2]:.1%}; "
          f"D censored model {r[3]:.1%}; average estimated difference B {r[4]:+.2f}, D {r[5]:+.2f} log10")