# Sample size. Primary outcome CE = |H - 0.80| over 20 items (H = share of 80 %
# intervals that hit). ANCOVA on pre-test CE (correlation rho); alpha 0.05 two-sided,
# power 0.90. SD of CE from a simulated control arm: hit chance ~ Beta(mean 0.50, SD
# 0.15), an assumption.
import math, random
from statistics import NormalDist, mean, stdev

rng = random.Random(15)
N = NormalDist()

def beta_ab(m, sd):
    k = m * (1 - m) / sd ** 2 - 1
    return m * k, (1 - m) * k

def ce_sample(n, m, sd=0.15, items=20):
    a, b = beta_ab(m, sd)
    out = []
    for _ in range(n):
        h = rng.betavariate(a, b)
        hits = sum(rng.random() < h for _ in range(items))
        out.append(abs(hits / items - 0.80))
    return out

def per_arm(delta, sigma, rho, alpha=0.05, power=0.90):
    z = N.inv_cdf(1 - alpha / 2) + N.inv_cdf(power)
    return math.ceil(2 * z * z * sigma ** 2 * (1 - rho ** 2) / delta ** 2) + 1

if __name__ == "__main__":
    ctrl, trt = ce_sample(20000, 0.50), ce_sample(20000, 0.57)
    sigma = stdev(ctrl)
    print(f"simulated control: mean CE {mean(ctrl):.3f}, SD {sigma:.3f};"
          f" hit chance 0.57 gives mean CE {mean(trt):.3f}")
    print("analysable participants per arm | rho = 0.0 | 0.3 | 0.5")
    for d in (0.05, 0.06, 0.08):
        ns = [per_arm(d, sigma, r) for r in (0.0, 0.3, 0.5)]
        print(f"   CE difference {d:.2f} (d = {d / sigma:.2f})  |"
              f" {ns[0]:9d} | {ns[1]:3d} | {ns[2]:3d}")
    n = per_arm(0.06, sigma, 0.5)
    rec = math.ceil(n / 0.90)
    print(f"chosen: 0.06, rho = 0.5 -> {n} per arm")
    print(f"   allowing 10 % loss (incomplete sessions, failed attention checks):")
    print(f"   recruit {rec} per arm, {2 * rec} in all")
    n0, reps, hits = per_arm(0.06, sigma, 0.0), 1000, 0
    for _ in range(reps):
        a, b = ce_sample(n0, 0.50), ce_sample(n0, 0.57)
        se = math.sqrt(stdev(a) ** 2 / n0 + stdev(b) ** 2 / n0)
        hits += abs(mean(a) - mean(b)) / se > N.inv_cdf(0.975)
    print(f"check without covariate: {n0} per arm, simulated power {hits / reps:.2f}")
    pk = [math.comb(20, k) * 0.8 ** k * 0.2 ** (20 - k) for k in range(21)]
    floor = sum(p * abs(k / 20 - 0.8) for k, p in enumerate(pk))
    print(f"floor: a perfectly calibrated person (h = 0.80) has expected CE {floor:.3f}"
          " on 20 items")
# simulated control: mean CE 0.304, SD 0.173; hit chance 0.57 gives mean CE 0.243
# analysable participants per arm | rho = 0.0 | 0.3 | 0.5
#    CE difference 0.05 (d = 0.29)  |       253 | 231 | 190
#    CE difference 0.06 (d = 0.35)  |       176 | 161 | 133
#    CE difference 0.08 (d = 0.46)  |       100 |  91 |  75
# chosen: 0.06, rho = 0.5 -> 133 per arm
#    allowing 10 % loss (incomplete sessions, failed attention checks):
#    recruit 148 per arm, 296 in all
# check without covariate: 176 per arm, simulated power 0.93
# floor: a perfectly calibrated person (h = 0.80) has expected CE 0.070 on 20 items

# Decision rule fixed in advance. D = arm A minus arm B in post-test CE (ANCOVA-adjusted);
# negative = better calibrated in A; SESOI 0.06. Manipulation check: A has lower base-rate
# error than B (one-sided p < 0.05). Integrity: at most 15 % excluded.
SESOI = 0.06

def verdict(d, lo, hi, manip_ok, excluded):
    if excluded > 0.15:
        return f"D = {d:+.3f}: NOT INTERPRETABLE - {excluded:.0%} excluded (limit 15 %)"
    if not manip_ok:
        return f"D = {d:+.3f}: NOT INTERPRETABLE - manipulation check failed"
    if hi <= -SESOI:
        v = "improves calibration by at least 0.06"
    elif hi < 0:
        v = "improves calibration; not shown to reach 0.06"
    elif lo > 0:
        v = "training WORSENS calibration"
    elif lo >= -SESOI and hi <= SESOI:
        v = "no difference worth acting on"
    else:
        v = "unclear: the interval is too wide to decide"
    return f"D = {d:+.3f} ({lo:+.3f} to {hi:+.3f}): {v}"

for case in ((-0.095, -0.130, -0.060, True, 0.08), (-0.040, -0.075, -0.005, True, 0.07),
             (-0.010, -0.045, 0.025, True, 0.09), (0.045, 0.010, 0.080, True, 0.06),
             (-0.020, -0.080, 0.040, True, 0.10), (-0.110, -0.145, -0.075, False, 0.05),
             (-0.090, -0.125, -0.055, True, 0.21)):
    print(verdict(*case))
print("secondary (signed gap 0.80 - H, base-rate error, interval width, 2-week CE):")
print("   reported with intervals; none decides the verdict (the base-rate error is used")
print("   only as the manipulation check, which decides interpretability)")
# D = -0.095 (-0.130 to -0.060): improves calibration by at least 0.06
# D = -0.040 (-0.075 to -0.005): improves calibration; not shown to reach 0.06
# D = -0.010 (-0.045 to +0.025): no difference worth acting on
# D = +0.045 (+0.010 to +0.080): training WORSENS calibration
# D = -0.020 (-0.080 to +0.040): unclear: the interval is too wide to decide
# D = -0.110: NOT INTERPRETABLE - manipulation check failed
# D = -0.090: NOT INTERPRETABLE - 21% excluded (limit 15 %)
# secondary (signed gap 0.80 - H, base-rate error, interval width, 2-week CE):
#    reported with intervals; none decides the verdict (the base-rate error is used
#    only as the manipulation check, which decides interpretability)