"""Simulation for /guides/anova-assumptions-not-met/: skewed, positive (log-normal) data in three groups.

Question: for concentration-like data whose SD grows with the mean, which analysis keeps the
false-positive rate near 5% and which detects a real multiplicative difference most often?
Methods: classical one-way ANOVA on raw values, Welch's ANOVA on raw values, Welch's ANOVA on
log values, Kruskal-Wallis on raw values. log-scale SD 0.6 in every group (as in the worked example).
No difference: all groups share the same distribution. Difference: geometric means 1, 1, 1.8.
20,000 runs per scenario, seed 20261110 (fixed in advance). Python 3.13, SciPy 1.18.1.
Run: python simulation.py  (writes simulation_output.json)
"""
import json

import numpy as np
from scipy import stats

RUNS = 20_000
SIGMA = 0.6
rng = np.random.default_rng(20261110)

def one(groups):
    logs = [np.log(g) for g in groups]
    return np.array([
        stats.f_oneway(*groups).pvalue,
        stats.f_oneway(*groups, equal_var=False).pvalue,
        stats.f_oneway(*logs, equal_var=False).pvalue,
        stats.kruskal(*groups).pvalue,
    ]) < 0.05

scenarios = []
for sizes in [(8, 8, 8), (5, 8, 11), (11, 8, 5)]:
    for label, gms in [("no difference", (1.0, 1.0, 1.0)), ("third group 1.8 times higher", (1.0, 1.0, 1.8))]:
        hits = np.zeros(4)
        for _ in range(RUNS):
            groups = [np.exp(rng.normal(np.log(gm), SIGMA, n)) for gm, n in zip(gms, sizes)]
            hits += one(groups)
        rates = hits / RUNS
        scenarios.append({"sizes": list(sizes), "truth": label,
                          "classical_raw": rates[0], "welch_raw": rates[1],
                          "welch_log": rates[2], "kruskal_raw": rates[3]})
        print(sizes, label, np.round(rates * 100, 1))

json.dump({"runs": RUNS, "sigma_log": SIGMA, "seed": 20261110, "scenarios": scenarios},
          open("simulation_output.json", "w"), indent=2)