# Reviewer's independent check: exact coverage of the nominal 95% Wald interval
# phat +/- z sqrt(phat(1-phat)/n), z = 1.959963984540054, closed interval,
# n = 5..100, p = 0.001..0.999 step 0.001. Membership decided exactly with Fractions;
# coverage summed in floating point (scipy) and, near the thresholds, re-checked exactly.
from fractions import Fraction as F
import numpy as np
from scipy.stats import binom
Z = F("1.959963984540054")
below93 = below95 = pairs = 0
near = 0
for n in range(5, 101):
    k = np.arange(n + 1)
    for j in range(1, 1000):
        p = F(j, 1000)
        acc = []
        for kk in range(n + 1):
            ph = F(kk, n)
            d = ph - p
            # |d| <= z sqrt(ph(1-ph)/n)  <=>  d^2 * n <= z^2 ph(1-ph)
            acc.append(d * d * n <= Z * Z * ph * (1 - ph))
        acc = np.array(acc)
        cov = binom.pmf(k, n, j / 1000)[acc].sum()
        if min(abs(cov - 0.93), abs(cov - 0.95)) < 1e-9:
            near += 1
            num = sum(F(int(__import__('math').comb(n, kk))) * p**kk * (1 - p)**(n - kk) for kk in range(n + 1) if acc[kk])
            cov = float(num) if True else cov
            c93, c95 = num < F(93, 100), num < F(95, 100)
        else:
            c93, c95 = cov < 0.93, cov < 0.95
        pairs += 1; below93 += c93; below95 += c95
print({"pairs": pairs, "below_93": below93, "below_95": below95,
       "share_93": round(below93 / pairs, 6), "share_95": round(below95 / pairs, 6), "near_threshold_exact_rechecks": near})