#!/usr/bin/env python3
"""Analytic Student-four tail and independent Gaussian reference experiment."""
import argparse
import json
import math
from pathlib import Path
import numpy as np


def verify(reference):
    report = json.loads(reference.read_text())
    threshold = report["threshold"]
    u = threshold / math.sqrt(threshold * threshold + 4)
    # Integrating 3/8 * (1+t²/4)^(-5/2), with t=2tan(theta), gives this tail.
    exact_tail = .5 - .75 * u + .25 * u ** 3
    rng = np.random.default_rng(report["monte_carlo_seed"])
    values = rng.normal(size=(report["monte_carlo_samples"], 6))
    target, comparisons = values[:, 0], values[:, 1:]
    statistic = (target - comparisons.mean(axis=1)) / (comparisons.std(axis=1, ddof=1) * math.sqrt(1 + 1 / 5))
    exceedances = int(np.sum(statistic > threshold))
    assert abs(exact_tail - report["exact_student4_one_sided_tail"]) < 1e-14
    assert exceedances == report["monte_carlo_exceedances"]
    print("PASS Student-four reference", exact_tail, "and", exceedances, "of", report["monte_carlo_samples"], "independent Gaussian experiments")


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--reference", type=Path, default=Path("small-sample-reference.json"))
    args = parser.parse_args()
    verify(args.reference)