"""Repeated-sampling experiments; results include Monte Carlo uncertainty."""
from pathlib import Path
import json
import numpy as np
from scipy.stats import norm, binom

ROOT = Path(__file__).resolve().parents[1]


def interval_experiment(n: int, p: float, repeats: int, rng) -> dict:
    counts = rng.binomial(n, p, repeats)
    phat = counts / n
    z = norm.ppf(0.975)
    se = np.sqrt(phat * (1 - phat) / n)
    denominator = 1 + z * z / n
    center = (phat + z * z / (2 * n)) / denominator
    half = z * np.sqrt(phat * (1 - phat) / n + z * z / (4 * n * n)) / denominator
    wald = (phat - z * se <= p) & (p <= phat + z * se)
    wilson = (center - half <= p) & (p <= center + half)
    def exact_coverage(method):
        k = np.arange(n + 1)
        sample_p = k / n
        if method == "wald":
            h = z * np.sqrt(sample_p * (1 - sample_p) / n)
            lo, hi = sample_p - h, sample_p + h
        else:
            mid = (sample_p + z * z / (2 * n)) / denominator
            h = z * np.sqrt(sample_p * (1 - sample_p) / n + z * z / (4 * n * n)) / denominator
            lo, hi = mid - h, mid + h
        return float(binom.pmf(k[(lo <= p) & (p <= hi)], n, p).sum())
    result = {"n": n, "p": p, "repeats": repeats}
    for method, event in [("wald", wald), ("wilson", wilson)]:
        coverage = float(event.mean())
        result[method] = {"empirical_coverage": coverage, "exact_binomial_coverage": exact_coverage(method),
                          "monte_carlo_se": float(np.sqrt(coverage * (1 - coverage) / repeats))}
    return result


if __name__ == "__main__":
    rng = np.random.default_rng(20261012)
    repeats = 20000
    experiments = [interval_experiment(n, p, repeats, rng) for n, p in [(20, 0.05), (100, 0.30)]]
    null_pvalues = rng.uniform(size=(repeats, 20))
    familywise = {}
    for name, alpha in [("unadjusted", 0.05), ("bonferroni", 0.05 / 20)]:
        rate = float((null_pvalues < alpha).any(axis=1).mean())
        familywise[name] = {"empirical_fwer": rate, "theoretical_fwer_under_independence": 1 - (1 - alpha) ** 20,
                            "monte_carlo_se": float(np.sqrt(rate * (1 - rate) / repeats))}
    result = {"seed": 20261012, "repeats": repeats, "ci_coverage": experiments,
              "all_null_independent_uniform_pvalues_20_tests": familywise,
              "note": "Uniform independent null p-values are an idealized model, not a universal description of tests."}
    (ROOT / "results").mkdir(exist_ok=True)
    output = ROOT / "results/simulation.json"
    output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n", encoding="utf-8")
    print(f"Simulation complete: {output}")