"""Worked example for /guides/when-can-you-exclude-an-outlier/ (constructed data, built by hand, no seed).

IL-6 (pg/mL), 6 control and 6 treated mice; treated mouse T6 reads 118 pg/mL.
Python 3.13, NumPy 2.5, SciPy 1.18.1. Run from this folder: python example.py
"""
import json

import numpy as np
from scipy import stats

control = np.array([40, 44, 37, 46, 42, 39.0])
treated = np.array([52, 58, 49, 61, 55, 118.0])
out = {}


def welch(a, b):
    r = stats.ttest_ind(a, b, equal_var=False)
    ci = r.confidence_interval(0.95)
    return {"diff": float(a.mean() - b.mean()), "t": float(r.statistic), "df": float(r.df),
            "p": float(r.pvalue), "ci": [float(ci.low), float(ci.high)]}


out["means"] = {"control": float(control.mean()), "treated": float(treated.mean()),
                "treated_without_T6": float(treated[:-1].mean())}
out["sds"] = {"control": float(control.std(ddof=1)), "treated": float(treated.std(ddof=1)),
              "treated_without_T6": float(treated[:-1].std(ddof=1))}
out["welch_all"] = welch(treated, control)
out["welch_without_T6"] = welch(treated[:-1], control)
mw = stats.mannwhitneyu(treated, control, alternative="two-sided", method="exact")
out["mann_whitney_all"] = {"U": float(mw.statistic), "p": float(mw.pvalue)}
yu = stats.ttest_ind(treated, control, equal_var=False, trim=0.2)
out["yuen_20pct_all"] = {"t": float(yu.statistic), "df": float(yu.df), "p": float(yu.pvalue),
                         "trimmed_diff": float(stats.trim_mean(treated, 0.2) - stats.trim_mean(control, 0.2))}
lg = stats.ttest_ind(np.log(treated), np.log(control), equal_var=False)
ci = lg.confidence_interval(0.95)
out["welch_log_all"] = {"t": float(lg.statistic), "df": float(lg.df), "p": float(lg.pvalue),
                        "ratio": float(np.exp(np.log(treated).mean() - np.log(control).mean())),
                        "ratio_ci": [float(np.exp(ci.low)), float(np.exp(ci.high))]}

lw = stats.ttest_ind(np.log(treated[:-1]), np.log(control), equal_var=False)
ciw = lw.confidence_interval(0.95)
out["welch_log_without_T6"] = {"t": float(lw.statistic), "df": float(lw.df), "p": float(lw.pvalue),
                               "ratio": float(np.exp(np.log(treated[:-1]).mean() - np.log(control).mean())),
                               "ratio_ci": [float(np.exp(ciw.low)), float(np.exp(ciw.high))]}

# Screens on the treated group
z = (treated - treated.mean()) / treated.std(ddof=1)
out["z_treated"] = [float(v) for v in z]
n = len(treated)
out["max_possible_abs_z_n6"] = float((n - 1) / np.sqrt(n))
# Grubbs (two-sided, alpha 0.05), NIST e-Handbook 1.3.5.17.1
G = float(np.max(np.abs(z)))
tcrit = stats.t.ppf(1 - 0.05 / (2 * n), n - 2)
G_crit = float((n - 1) / np.sqrt(n) * np.sqrt(tcrit**2 / (n - 2 + tcrit**2)))
out["grubbs"] = {"G": G, "G_crit": G_crit}
med = np.median(treated)
mad = 1.4826 * np.median(np.abs(treated - med))
out["mad_scores_treated"] = [float(v) for v in (treated - med) / mad]

print(json.dumps(out, indent=2))
with open("example_output.json", "w") as f:
    json.dump(out, f, indent=2)