"""Reproduce the Reichenau retrospective binomial analysis with Python 3.

Save this script and reichenau_67_audit.json in the same directory, then run:
    python3 statistical_study.py
Outputs statistical_study.json beside the script. Standard library only.
Each included row is one Bernoulli trial under the assumed model. Status
filters retain the author's original side; they do not validate that side.
"""
import hashlib
import json
from decimal import Decimal, localcontext
from fractions import Fraction
from math import comb
from pathlib import Path

ROOT = Path(__file__).resolve().parent
SOURCE = ROOT / "reichenau_67_audit.json"

def tail(n, k, q):
    with localcontext() as ctx:
        ctx.prec = 70
        q = Decimal(str(q))
        return sum(Decimal(comb(n, j)) * q**j * (1-q)**(n-j)
                   for j in range(k, n+1))

def analyze():
    rows = json.loads(SOURCE.read_text())
    assert len(rows) == len({r["id"] for r in rows}) == 67
    included = [r for r in rows if r["included_in_original_tally"]]
    assert all(r["author_side_original"] in {"late", "classical"} for r in included)
    assert len(included) == 66
    assert sum(r["author_side_original"] == "late" for r in included) == 49
    groups = [
        ("all_counted", included),
        ("supported_or_qualified", [r for r in included if r["audit_status"] in
                                    {"supported_lexeme", "qualified_relation"}]),
        ("supported_lexeme", [r for r in included if r["audit_status"] == "supported_lexeme"]),
    ]
    results = []
    for name, group in groups:
        n = len(group)
        k = sum(r["author_side_original"] == "late" for r in group)
        exact = Fraction(sum(comb(n, j) for j in range(k, n+1)), 2**n)
        computed = tail(n, k, "0.5")
        assert abs(float(exact) - float(computed)) < 1e-16
        results.append(dict(subset=name, n=n, k=k, proportion=k/n,
                            p_one_sided_q_half=str(computed),
                            p_exact_fraction=str(exact), row_ids=[r["id"] for r in group]))
    scenarios = [dict(q=q, p_one_sided=str(tail(66, 49, q)))
                 for q in ["0.50", "0.60", "0.65", "0.75"]]
    lo, hi = Decimal("0.5"), Decimal("0.75")
    for _ in range(80):
        mid = (lo + hi)/2
        if tail(66, 49, mid) < Decimal("0.05"):
            lo = mid
        else:
            hi = mid
    return dict(
        analysis_date="2026-09-09", source_file=SOURCE.name,
        source_sha256=hashlib.sha256(SOURCE.read_bytes()).hexdigest(),
        design="Retrospective analysis of a selected published list; no new lexical coding.",
        null="Independent Bernoulli trials, constant q=0.5; alternative q>0.5.",
        alpha="0.05", empirical_null_available=False,
        interpretation="Conditional binomial probabilities; no date or historical-cause estimate.",
        subset_rule="Filter existing dictionary statuses, retaining author_side_original.",
        subset_tests="Exploratory sensitivity analyses on nested, non-independent subsets.",
        subsets=results, null_sensitivity=scenarios,
        q_at_one_sided_p_005=str((lo+hi)/2),
        method="Exact binomial upper tail, Decimal precision 70; q=0.5 checked with rational arithmetic.",
        references=["https://www.jhss.ro/downloads/30/articles/vol%2015%20no%202%20(30)%202024-105-117.pdf",
                    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binomtest.html"])

if __name__ == "__main__":
    result = analyze()
    (ROOT / "statistical_study.json").write_text(json.dumps(result, ensure_ascii=False, indent=2)+"\n")
    print(json.dumps({"subsets": [{k:v for k,v in s.items() if k != "row_ids"} for s in result["subsets"]],
                      "null_sensitivity": result["null_sensitivity"],
                      "q_at_one_sided_p_005": result["q_at_one_sided_p_005"]}, indent=2))