"""False-positive rates for /guides/superplots/: cell-level Welch t-test vs paired t-test on
experiment means, when the drug has NO average effect. Design as in the worked example: 4 independent
experiments, 50 cells per condition per experiment, cell SD 2, day-to-day baseline SD 1.
Scenario A: the drug effect varies between experiments (SD 0.5) around a true mean of 0.
Scenario B: no variation in the effect (only the baseline differs between days).
Python 3.13.18, NumPy 2.5, SciPy 1.18.1. Run: python simulation.py -> simulation_output.json"""
import json

import numpy as np
from scipy import stats

N_SIM, N_EXP, N_CELLS, SD_CELL, SD_DAY = 20_000, 4, 50, 2.0, 1.0


def run(sd_effect, seed):
    rng = np.random.default_rng(seed)
    hits_cell = hits_mean = 0
    for _ in range(N_SIM):
        day = rng.normal(0, SD_DAY, N_EXP)
        eff = rng.normal(0, sd_effect, N_EXP) if sd_effect > 0 else np.zeros(N_EXP)
        ctrl = day[:, None] + rng.normal(0, SD_CELL, (N_EXP, N_CELLS))
        drug = (day + eff)[:, None] + rng.normal(0, SD_CELL, (N_EXP, N_CELLS))
        if stats.ttest_ind(drug.ravel(), ctrl.ravel(), equal_var=False).pvalue < 0.05:
            hits_cell += 1
        if stats.ttest_rel(drug.mean(axis=1), ctrl.mean(axis=1)).pvalue < 0.05:
            hits_mean += 1
    return {"cell_level_welch": hits_cell / N_SIM, "paired_t_on_experiment_means": hits_mean / N_SIM}


out = {"n_sim": N_SIM, "design": {"experiments": N_EXP, "cells_per_condition_per_experiment": N_CELLS,
                                  "sd_cell": SD_CELL, "sd_day": SD_DAY},
       "A_effect_varies_sd_0.5": run(0.5, 1),
       "B_effect_constant": run(0.0, 2)}
with open("simulation_output.json", "w") as f:
    json.dump(out, f, indent=2)
print(json.dumps(out, indent=2))