| |
| """ |
| protocol_engine.py -- protocol design engine. Python standard library only. |
| |
| Input : one JSON config describing a registered construct test |
| (stages/arms, effect thresholds, alpha, target power, gate rules). |
| Output: a markdown design report (stdout) containing, per arm: |
| 1. required sample size, closed form; |
| 2. Monte Carlo power at that size (cross-check of the closed form); |
| 3. power at the currently registered size; |
| plus a falsifier-fire audit, an input-separation gate budget, |
| a scale table, and a registration skeleton filled from the config. |
| |
| Deterministic: one seed in the config; every simulation stream derives |
| from it. Re-running the same command reproduces every number. |
| |
| Run: python3 protocol_engine.py config_summarization_bias.json |
| """ |
|
|
| import json |
| import math |
| import random |
| import sys |
| from statistics import NormalDist |
|
|
| ND = NormalDist() |
|
|
|
|
| |
| |
| |
|
|
| def z_alpha(alpha, sides): |
| return ND.inv_cdf(1 - alpha / sides) |
|
|
|
|
| def z_power(power): |
| return ND.inv_cdf(power) |
|
|
|
|
| def cohens_h(p1, p2): |
| return 2 * math.asin(math.sqrt(p1)) - 2 * math.asin(math.sqrt(p2)) |
|
|
|
|
| def rate_from_h(p_base, h): |
| """Rate p2 > p_base such that cohens_h(p2, p_base) == h.""" |
| return math.sin(math.asin(math.sqrt(p_base)) + h / 2.0) ** 2 |
|
|
|
|
| def rate_from_or(p_base, odds_ratio): |
| odds = odds_ratio * p_base / (1.0 - p_base) |
| return odds / (1.0 + odds) |
|
|
|
|
| def binom_pmf(k, n, p): |
| return math.comb(n, k) * p ** k * (1.0 - p) ** (n - k) |
|
|
|
|
| def binom_cdf(k, n, p): |
| if k < 0: |
| return 0.0 |
| if k >= n: |
| return 1.0 |
| return sum(binom_pmf(i, n, p) for i in range(k + 1)) |
|
|
|
|
| def binom_draw(rng, n, p): |
| if hasattr(rng, "binomialvariate"): |
| return rng.binomialvariate(n, p) |
| return sum(1 for _ in range(n) if rng.random() < p) |
|
|
|
|
| def fmt(x, digits=3): |
| return f"{x:.{digits}f}" |
|
|
|
|
| |
| |
| |
|
|
| def n_two_proportions(h, alpha, sides, power): |
| """Per-group n for a two-proportion test, arcsine (Cohen's h) form.""" |
| z = z_alpha(alpha, sides) + z_power(power) |
| return math.ceil(2.0 * (z / h) ** 2) |
|
|
|
|
| def n_one_proportion(p0, p1, alpha, sides, power): |
| za = z_alpha(alpha, sides) |
| zb = z_power(power) |
| num = za * math.sqrt(p0 * (1 - p0)) + zb * math.sqrt(p1 * (1 - p1)) |
| return math.ceil((num / abs(p1 - p0)) ** 2) |
|
|
|
|
| def _beta_cf(a, b, x): |
| """Continued fraction for the incomplete beta function.""" |
| qab, qap, qam = a + b, a + 1.0, a - 1.0 |
| c = 1.0 |
| d = 1.0 - qab * x / qap |
| if abs(d) < 1e-300: |
| d = 1e-300 |
| d = 1.0 / d |
| h = d |
| for m in range(1, 201): |
| m2 = 2 * m |
| aa = m * (b - m) * x / ((qam + m2) * (a + m2)) |
| d = 1.0 + aa * d |
| if abs(d) < 1e-300: |
| d = 1e-300 |
| c = 1.0 + aa / c |
| if abs(c) < 1e-300: |
| c = 1e-300 |
| d = 1.0 / d |
| h *= d * c |
| aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) |
| d = 1.0 + aa * d |
| if abs(d) < 1e-300: |
| d = 1e-300 |
| c = 1.0 + aa / c |
| if abs(c) < 1e-300: |
| c = 1e-300 |
| d = 1.0 / d |
| delta = d * c |
| h *= delta |
| if abs(delta - 1.0) < 3e-16: |
| break |
| return h |
|
|
|
|
| def _reg_incomplete_beta(x, a, b): |
| if x <= 0.0: |
| return 0.0 |
| if x >= 1.0: |
| return 1.0 |
| bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) |
| + a * math.log(x) + b * math.log(1.0 - x)) |
| if x < (a + 1.0) / (a + b + 2.0): |
| return bt * _beta_cf(a, b, x) / a |
| return 1.0 - bt * _beta_cf(b, a, 1.0 - x) / b |
|
|
|
|
| def t_cdf(t, df): |
| """CDF of the central Student t distribution (incomplete beta).""" |
| x = df / (df + t * t) |
| ib = _reg_incomplete_beta(x, df / 2.0, 0.5) |
| return 1.0 - 0.5 * ib if t >= 0 else 0.5 * ib |
|
|
|
|
| def t_crit(alpha, sides, df): |
| target = 1.0 - alpha / sides |
| lo, hi = 0.0, 1000.0 |
| while hi - lo > 1e-12: |
| mid = (lo + hi) / 2.0 |
| if t_cdf(mid, df) < target: |
| lo = mid |
| else: |
| hi = mid |
| return (lo + hi) / 2.0 |
|
|
|
|
| def _chi_pdf(v, df): |
| if v <= 0.0: |
| return 0.0 |
| k = df / 2.0 |
| return math.exp((k - 1.0) * math.log(v) - v / 2.0 |
| - k * math.log(2.0) - math.lgamma(k)) |
|
|
|
|
| def _adaptive_simpson(f, a, b, eps, depth): |
| c = (a + b) / 2.0 |
| fa, fb, fc = f(a), f(b), f(c) |
| whole = (b - a) * (fa + 4.0 * fc + fb) / 6.0 |
|
|
| def step(a, b, eps, whole, fa, fb, fc, depth): |
| c = (a + b) / 2.0 |
| lm, rm = (a + c) / 2.0, (c + b) / 2.0 |
| flm, frm = f(lm), f(rm) |
| left = (c - a) * (fa + 4.0 * flm + fc) / 6.0 |
| right = (b - c) * (fc + 4.0 * frm + fb) / 6.0 |
| delta = left + right - whole |
| if depth <= 0 or abs(delta) <= 15.0 * eps: |
| return left + right + delta / 15.0 |
| return (step(a, c, eps / 2.0, left, fa, fc, flm, depth - 1) |
| + step(c, b, eps / 2.0, right, fc, fb, frm, depth - 1)) |
|
|
| return step(a, b, eps, whole, fa, fb, fc, depth) |
|
|
|
|
| def nct_cdf(t, df, ncp): |
| """CDF of the noncentral t distribution, from its definition |
| T = (Z + ncp) / sqrt(V/df): P(T <= t) = E_V[Phi(t sqrt(V/df) - ncp)] |
| over the chi-square density of V, by adaptive Simpson quadrature.""" |
| if t < 0: |
| return 1.0 - nct_cdf(-t, df, -ncp) |
| f = lambda v: ND.cdf(t * math.sqrt(v / df) - ncp) * _chi_pdf(v, df) |
| hi = df + 20.0 * math.sqrt(2.0 * df) + 60.0 |
| return _adaptive_simpson(f, 0.0, hi, 1e-12, 40) |
|
|
|
|
| def paired_t_power_exact(n, dz, alpha, sides): |
| """Exact power of the paired (one-sample) t test on differences: |
| noncentral t at df = n-1, noncentrality dz*sqrt(n).""" |
| df = float(n - 1) |
| ncp = dz * math.sqrt(n) |
| crit = t_crit(alpha, sides, df) |
| power = 1.0 - nct_cdf(crit, df, ncp) |
| if sides == 2: |
| power += nct_cdf(-crit, df, ncp) |
| return power |
|
|
|
|
| def n_paired_dz(dz, alpha, sides, power): |
| """Smallest n whose exact paired-t power reaches the target -- |
| an exact search over the noncentral t power, not the normal |
| closed form.""" |
| n = 2 |
| while paired_t_power_exact(n, dz, alpha, sides) < power: |
| n += 1 |
| return n |
|
|
|
|
| |
| |
| |
|
|
| def two_prop_z_pvalue(x1, n1, x2, n2): |
| pooled = (x1 + x2) / (n1 + n2) |
| se = math.sqrt(pooled * (1 - pooled) * (1.0 / n1 + 1.0 / n2)) |
| if se == 0.0: |
| return 1.0 |
| z = (x1 / n1 - x2 / n2) / se |
| return 2.0 * (1.0 - ND.cdf(abs(z))) |
|
|
|
|
| def mc_power_two_prop(n, p1, p2, alpha, reps, seed): |
| rng = random.Random(seed) |
| hits = 0 |
| for _ in range(reps): |
| x = binom_draw(rng, n, p1) |
| y = binom_draw(rng, n, p2) |
| if two_prop_z_pvalue(x, n, y, n) < alpha: |
| hits += 1 |
| p = hits / reps |
| return p, math.sqrt(p * (1 - p) / reps) |
|
|
|
|
| def exact_power_two_prop_z(n, p1, p2, alpha): |
| """Exact power of the pooled z test by enumerating all outcome pairs. |
| Practical for small n (used at the registered n = 20).""" |
| probs1 = [binom_pmf(x, n, p1) for x in range(n + 1)] |
| probs2 = [binom_pmf(y, n, p2) for y in range(n + 1)] |
| total = 0.0 |
| for x in range(n + 1): |
| for y in range(n + 1): |
| if two_prop_z_pvalue(x, n, y, n) < alpha: |
| total += probs1[x] * probs2[y] |
| return total |
|
|
|
|
| def fisher_pvalue(x, n, y, m): |
| """Two-sided Fisher exact p-value for [[x, n-x], [y, m-y]].""" |
| s = x + y |
| lo = max(0, s - m) |
| hi = min(n, s) |
| denom = math.comb(n + m, s) |
|
|
| def hyper(k): |
| return math.comb(n, k) * math.comb(m, s - k) / denom |
|
|
| p_obs = hyper(x) |
| return sum(hyper(k) for k in range(lo, hi + 1) |
| if hyper(k) <= p_obs * (1 + 1e-12)) |
|
|
|
|
| def exact_power_fisher(n, p1, p2, alpha): |
| probs1 = [binom_pmf(x, n, p1) for x in range(n + 1)] |
| probs2 = [binom_pmf(y, n, p2) for y in range(n + 1)] |
| total = 0.0 |
| for x in range(n + 1): |
| for y in range(n + 1): |
| if fisher_pvalue(x, n, y, n) < alpha: |
| total += probs1[x] * probs2[y] |
| return total |
|
|
|
|
| |
| |
| |
|
|
| def one_sample_critical(n, p0, alpha, sides): |
| """Equal-tail exact critical counts. Reject if X <= lo or X >= hi.""" |
| if sides == 2: |
| lo = -1 |
| for k in range(n + 1): |
| if binom_cdf(k, n, p0) <= alpha / 2: |
| lo = k |
| else: |
| break |
| hi = n + 1 |
| for k in range(n, -1, -1): |
| if 1.0 - binom_cdf(k - 1, n, p0) <= alpha / 2: |
| hi = k |
| else: |
| break |
| return lo, hi |
| hi = n + 1 |
| for k in range(n, -1, -1): |
| if 1.0 - binom_cdf(k - 1, n, p0) <= alpha: |
| hi = k |
| else: |
| break |
| return -1, hi |
|
|
|
|
| def exact_power_one_sample(n, p0, p1, alpha, sides): |
| lo, hi = one_sample_critical(n, p0, alpha, sides) |
| power = 0.0 |
| if lo >= 0: |
| power += binom_cdf(lo, n, p1) |
| if hi <= n: |
| power += 1.0 - binom_cdf(hi - 1, n, p1) |
| return power |
|
|
|
|
| |
| |
| |
|
|
| def mc_power_paired_dz(n, dz, alpha, sides, reps, seed): |
| """Independent check of the exact paired-t size: simulate the |
| paired t test itself. Differences D ~ Normal(dz, 1) are drawn, |
| the variance is ESTIMATED from each sample (sample sd, ddof=1), |
| and the statistic is compared with the t critical value -- the |
| t behaviour the closed form must answer to. Because the check |
| shares no formula with n_paired_dz, it can falsify it: run at |
| the old normal-approximation size n = 32 (dz = 0.5, alpha 0.05) |
| it returns ~0.78, below the 0.80 target the old size claimed.""" |
| rng = random.Random(seed) |
| crit = t_crit(alpha, sides, float(n - 1)) |
| hits = 0 |
| for _ in range(reps): |
| diffs = [rng.gauss(dz, 1.0) for _ in range(n)] |
| mean = sum(diffs) / n |
| var = sum((d - mean) ** 2 for d in diffs) / (n - 1) |
| if var <= 0.0: |
| continue |
| t = mean / math.sqrt(var / n) |
| if (sides == 2 and abs(t) > crit) or (sides == 1 and t > crit): |
| hits += 1 |
| p = hits / reps |
| return p, math.sqrt(p * (1 - p) / reps) |
|
|
|
|
| |
| |
| |
|
|
| def gate_build_count(n_required, fail_rate, assurance): |
| """Smallest build count B with P(Binomial(B, 1-fail) >= n_required) |
| >= assurance. Exact search over the binomial CDF.""" |
| p_pass = 1.0 - fail_rate |
| B = n_required |
| while True: |
| if 1.0 - binom_cdf(n_required - 1, B, p_pass) >= assurance: |
| return B |
| B += 1 |
|
|
|
|
| def mc_gate_check(B, n_required, fail_rate, reps, seed): |
| rng = random.Random(seed) |
| hits = sum(1 for _ in range(reps) |
| if binom_draw(rng, B, 1.0 - fail_rate) >= n_required) |
| p = hits / reps |
| return p, math.sqrt(p * (1 - p) / reps) |
|
|
|
|
| def exact_search_one_sample(p0, p1, alpha, sides, target, n_start): |
| """Smallest n >= n_start whose exact binomial power reaches target. |
| The exact test is conservative at small n, so this sits at or above |
| the closed-form value. Search is capped at n_start + 120.""" |
| n = n_start |
| while n <= n_start + 120: |
| if exact_power_one_sample(n, p0, p1, alpha, sides) >= target: |
| return n |
| n += 1 |
| return n_start |
|
|
|
|
| |
| |
| |
| |
| |
| |
|
|
| DECISIVE_ALPHA = 0.01 |
| DECISIVE_POWER = 0.95 |
|
|
|
|
| def run(cfg): |
| alpha = cfg["alpha"] |
| sides = cfg.get("sides", 2) |
| power_t = cfg["target_power"] |
| reps = cfg.get("mc_reps", 20000) |
| seed0 = cfg.get("seed", 1) |
| out = [] |
| out.append(f"# Design report -- {cfg['construct_name']}") |
| out.append("") |
| out.append(f"Source protocol: {cfg['source_protocol']}") |
| out.append("") |
| out.append(f"- Standard design point: alpha = {alpha} ({sides}-sided), " |
| f"power = {power_t}") |
| out.append(f"- Decisive design point: alpha = {DECISIVE_ALPHA} " |
| f"({sides}-sided), power = {DECISIVE_POWER}") |
| out.append(f"- Monte Carlo reps per cell = {reps}; base seed = {seed0}; " |
| f"one stream per row. Re-running reproduces every number.") |
| out.append("") |
|
|
| |
| out.append("## 1. Required sample size per arm") |
| out.append("") |
| out.append("| Arm | Test | Effect as specified | n at standard point | " |
| "MC power at that n | n at decisive point | Current n | " |
| "Power at current n |") |
| out.append("|---|---|---|---|---|---|---|---|") |
| res = {} |
|
|
| |
| s3 = cfg["stage3"] |
| p_base = s3["human_told_preference_rate"] |
| p_h = rate_from_h(p_base, s3["registered_thresholds"]["cohens_h"]) |
| p_or = rate_from_or(p_base, s3["registered_thresholds"]["odds_ratio"]) |
| h_or = abs(cohens_h(p_or, p_base)) |
| h_gov = min(s3["registered_thresholds"]["cohens_h"], h_or) |
| p_alt = p_h if s3["registered_thresholds"]["cohens_h"] <= h_or else p_or |
| n_std = n_two_proportions(h_gov, alpha, sides, power_t) |
| n_dec = n_two_proportions(h_gov, DECISIVE_ALPHA, sides, DECISIVE_POWER) |
| mc_std, se_std = mc_power_two_prop(n_std, p_alt, p_base, alpha, |
| reps, seed0 + 101) |
| n_cur = s3["current_pairs"] |
| mc_cur, se_cur = mc_power_two_prop(n_cur, p_alt, p_base, alpha, |
| reps, seed0 + 102) |
| res["stage3"] = dict(n_std=n_std, n_dec=n_dec, p_alt=p_alt, p_base=p_base, |
| h=h_gov) |
| out.append( |
| f"| Stage 3, evaluative (core) | two-proportion z, pooled | " |
| f"h = {s3['registered_thresholds']['cohens_h']} -> rate {p_base} vs " |
| f"{fmt(p_h)}; OR = {s3['registered_thresholds']['odds_ratio']} -> rate " |
| f"{fmt(p_or)} (h = {fmt(h_or)}); governing h = {fmt(h_gov)} | " |
| f"{n_std} per group | {fmt(mc_std)} +/- {fmt(se_std)} | {n_dec} per group | " |
| f"{n_cur} per group | {fmt(mc_cur)} +/- {fmt(se_cur)} |") |
|
|
| |
| s2 = cfg["stage2"] |
| dz = s2["planning_effect_dz"] |
| n2_std = n_paired_dz(dz, alpha, sides, power_t) |
| n2_dec = n_paired_dz(dz, DECISIVE_ALPHA, sides, DECISIVE_POWER) |
| mc2, se2 = mc_power_paired_dz(n2_std, dz, alpha, sides, reps, seed0 + 201) |
| mc2c, se2c = mc_power_paired_dz(s2["current_matrices"], dz, alpha, sides, |
| reps, seed0 + 202) |
| res["stage2"] = dict(n_std=n2_std, n_dec=n2_dec) |
| out.append( |
| f"| Stage 2, generative | paired t on SI differences | dz = {dz} " |
| f"(worked value -- the paper fixes no SI effect size) | {n2_std} pairs | " |
| f"{fmt(mc2)} +/- {fmt(se2)} | {n2_dec} pairs | " |
| f"{s2['current_matrices']} pairs | {fmt(mc2c)} +/- {fmt(se2c)} |") |
|
|
| |
| da = cfg["density_arm"] |
| p1d = rate_from_h(da["null_preference_rate"], |
| da["planning_effect_h_vs_null"]) |
| nd_cf = n_one_proportion(da["null_preference_rate"], p1d, alpha, sides, |
| power_t) |
| nd_std = exact_search_one_sample(da["null_preference_rate"], p1d, alpha, |
| sides, power_t, nd_cf) |
| nd_cf_dec = n_one_proportion(da["null_preference_rate"], p1d, |
| DECISIVE_ALPHA, sides, DECISIVE_POWER) |
| nd_dec = exact_search_one_sample(da["null_preference_rate"], p1d, |
| DECISIVE_ALPHA, sides, DECISIVE_POWER, |
| nd_cf_dec) |
| exd = exact_power_one_sample(nd_std, da["null_preference_rate"], p1d, |
| alpha, sides) |
| exdd = exact_power_one_sample(nd_dec, da["null_preference_rate"], p1d, |
| DECISIVE_ALPHA, sides) |
| exdc = exact_power_one_sample(20, da["null_preference_rate"], p1d, |
| alpha, sides) |
| res["density"] = dict(n_std=nd_std, n_dec=nd_dec) |
| out.append( |
| f"| Density calibration arm | one-sample exact binomial vs chance | " |
| f"h = {da['planning_effect_h_vs_null']} -> rate " |
| f"{da['null_preference_rate']} vs {fmt(p1d)} (exact search on " |
| f"the unrounded rate {fmt(p1d, 5)}) | {nd_std} pairs " |
| f"(exact power {fmt(exd)}; closed form {nd_cf}) | -- | {nd_dec} pairs " |
| f"(exact power {fmt(exdd)}; closed form {nd_cf_dec}) | 20 pairs " |
| f"(exact power {fmt(exdc)}) |") |
|
|
| |
| rv = cfg["reader_validation"] |
| nr_cf = n_one_proportion(rv["null_convergence_rate"], |
| rv["planning_target_rate"], alpha, sides, power_t) |
| nr_std = exact_search_one_sample(rv["null_convergence_rate"], |
| rv["planning_target_rate"], alpha, |
| sides, power_t, nr_cf) |
| nr_cf_dec = n_one_proportion(rv["null_convergence_rate"], |
| rv["planning_target_rate"], |
| DECISIVE_ALPHA, sides, DECISIVE_POWER) |
| nr_dec = exact_search_one_sample(rv["null_convergence_rate"], |
| rv["planning_target_rate"], |
| DECISIVE_ALPHA, sides, DECISIVE_POWER, |
| nr_cf_dec) |
| exr = exact_power_one_sample(nr_std, rv["null_convergence_rate"], |
| rv["planning_target_rate"], alpha, sides) |
| exrd = exact_power_one_sample(nr_dec, rv["null_convergence_rate"], |
| rv["planning_target_rate"], |
| DECISIVE_ALPHA, sides) |
| exrc = exact_power_one_sample(20, rv["null_convergence_rate"], |
| rv["planning_target_rate"], alpha, sides) |
| res["reader"] = dict(n_std=nr_std, n_dec=nr_dec) |
| out.append( |
| f"| SI reader validation | one-sample exact binomial vs chance | " |
| f"rate {rv['null_convergence_rate']} vs " |
| f"{rv['planning_target_rate']} (worked values -- chance depends on the " |
| f"coding scheme) | {nr_std} items (exact power {fmt(exr)}; closed " |
| f"form {nr_cf}) | -- | {nr_dec} items (exact power {fmt(exrd)}; " |
| f"closed form {nr_cf_dec}) | 20 items " |
| f"(exact power {fmt(exrc)}) |") |
|
|
| |
| ra = cfg["register_arm"] |
| p_alt_r = rate_from_h(ra["human_baseline_rate"], |
| ra["planning_effect_h_between_registers"]) |
| nreg_std = n_two_proportions(ra["planning_effect_h_between_registers"], |
| alpha, sides, power_t) |
| nreg_dec = n_two_proportions(ra["planning_effect_h_between_registers"], |
| DECISIVE_ALPHA, sides, DECISIVE_POWER) |
| mcr, ser = mc_power_two_prop(nreg_std, p_alt_r, ra["human_baseline_rate"], |
| alpha, reps, seed0 + 301) |
| mcrc, serc = mc_power_two_prop(20, p_alt_r, ra["human_baseline_rate"], |
| alpha, reps, seed0 + 302) |
| res["register"] = dict(n_std=nreg_std, n_dec=nreg_dec) |
| out.append( |
| f"| Register arm | two-proportion z per register stratum | " |
| f"h = {ra['planning_effect_h_between_registers']} -> rate " |
| f"{ra['human_baseline_rate']} vs {fmt(p_alt_r)} | {nreg_std} contents " |
| f"per stratum ({nreg_std * 4} versions) | {fmt(mcr)} +/- {fmt(ser)} | " |
| f"{nreg_dec} contents per stratum | 20 contents per stratum | " |
| f"{fmt(mcrc)} +/- {fmt(serc)} |") |
| out.append("") |
|
|
| |
| approx = ND.cdf(h_gov * math.sqrt(n_cur / 2) - z_alpha(alpha, sides)) |
| mc_f, mc_f_se = mc_power_two_prop(n_cur, p_alt, p_base, alpha, |
| reps, seed0 + 999) |
| ex_z = exact_power_two_prop_z(n_cur, p_alt, p_base, alpha) |
| fisher = exact_power_fisher(n_cur, p_alt, p_base, alpha) |
| one_s = exact_power_one_sample(n_cur, p_base, p_or, alpha, sides) |
| paired_mc, _ = mc_power_paired_dz(n_cur, |
| s3["registered_thresholds"]["cohens_h"], |
| alpha, sides, reps, seed0 + 998) |
| out.append("## 2. Falsifier-fire audit -- Stage 3 at the registered size") |
| out.append("") |
| out.append(f"Setup: {n_cur} pairs per group; true effect set exactly at the " |
| f"registered threshold (h = {fmt(h_gov)}, rates {p_base} vs " |
| f"{fmt(p_alt)}). 'Fires' = the test is not significant, so the " |
| f"registered falsifier triggers against a true effect.") |
| out.append("") |
| out.append("| Reading | Test | Power | P(falsifier fires) |") |
| out.append("|---|---|---|---|") |
| out.append(f"| Registered (model rate vs human rate) | normal approx | " |
| f"{fmt(approx)} | {fmt(1 - approx)} |") |
| out.append(f"| Registered | pooled z, Monte Carlo | {fmt(mc_f)} +/- " |
| f"{fmt(mc_f_se)} | {fmt(1 - mc_f)} |") |
| out.append(f"| Registered | pooled z, exact enumeration | {fmt(ex_z)} | " |
| f"{fmt(1 - ex_z)} |") |
| out.append(f"| Registered | Fisher exact, exact enumeration | " |
| f"{fmt(fisher)} | {fmt(1 - fisher)} |") |
| out.append(f"| Alternate (model rate vs chance, OR threshold) | exact " |
| f"binomial | {fmt(one_s)} | {fmt(1 - one_s)} |") |
| out.append(f"| Alternate (paired per-item contrast) | paired t, Monte " |
| f"Carlo | {fmt(paired_mc)} | {fmt(1 - paired_mc)} |") |
| out.append("") |
| fires = [1 - approx, 1 - mc_f, 1 - ex_z, 1 - fisher] |
| out.append(f"Registered reading, across test variants: P(fire) = " |
| f"{fmt(min(fires))} to {fmt(max(fires))}.") |
| out.append("") |
|
|
| |
| gate = cfg["gate"] |
| for label, target in (("standard", n_std), ("decisive", n_dec)): |
| out.append(f"## 3{label == 'decisive' and 'b' or 'a'}. " |
| f"{gate['name']} -- build budget ({label} point)") |
| out.append("") |
| out.append(f"Target: >= {target} pairs pass the gate. Assurance: " |
| f"P(passing >= {target}) >= {gate['assurance']}. " |
| f"Exact binomial search; last column re-checks by Monte Carlo.") |
| out.append("") |
| out.append("| Pair failure rate | Pairs to build (exact) | " |
| "Expected passing | MC assurance |") |
| out.append("|---|---|---|---|") |
| for j, f in enumerate(gate["failure_rate_scenarios"]): |
| B = gate_build_count(target, f, gate["assurance"]) |
| mc_a, mc_a_se = mc_gate_check(B, target, f, reps, |
| seed0 + 5000 + (0 if label == "standard" else 100) + j) |
| out.append(f"| {f} | {B} | {fmt(B * (1 - f), 1)} | " |
| f"{fmt(mc_a)} +/- {fmt(mc_a_se)} |") |
| out.append("") |
|
|
| |
| sc = cfg["scale"] |
| judges = sc["model_families"] + sc["human_raters_per_pair"] |
| out.append("## 4. Scale table") |
| out.append("") |
| out.append(f"Per pair: {sc['versions_per_pair']} stimulus versions " |
| f"({sc['versions_per_pair'] * 2} with the register arm), " |
| f"{judges} Stage 3 judge calls ({sc['model_families']} model " |
| f"families + {sc['human_raters_per_pair']} human rater), " |
| f"SI coding passes = versions x 2 raters. Stage 4 = two further " |
| f"Stage 3 passes (length-mismatched, framing-varied).") |
| out.append("") |
| out.append("| Pairs | Versions | Versions with register arm | " |
| "Stage 3 judge calls | Stage 4 judge calls | SI coding passes |") |
| out.append("|---|---|---|---|---|---|") |
| for P in sc["pair_counts"]: |
| v = P * sc["versions_per_pair"] |
| calls = P * judges |
| out.append(f"| {P} | {v} | {v * 2} | {calls} | {calls * 2} | {v * 2} |") |
| out.append("") |
|
|
| |
| out.append("## 5. Registration skeleton (filled from this config)") |
| out.append("") |
| out.append("1. Construct: summarization bias -- the systematic tendency of a " |
| "model, in generation or in evaluation, to replace or prefer " |
| "reconstructable shown-mode structure with the abstract summary " |
| "label that names its content (paper S4.3).") |
| out.append(f"2. Thresholds, fixed in advance: model told-preference rate " |
| f"exceeds human rate with Cohen's h >= " |
| f"{s3['registered_thresholds']['cohens_h']} or odds ratio >= " |
| f"{s3['registered_thresholds']['odds_ratio']} (paper S5.3). " |
| f"Alpha = {alpha} ({sides}-sided) at the standard point; " |
| f"alpha = {DECISIVE_ALPHA}, power = {DECISIVE_POWER} at the " |
| f"decisive point.") |
| out.append(f"3. Stimuli: build the gate-budget count from Table 3 for the " |
| f"expected failure rate (standard point) -- word count +/-5%, " |
| f"grade matched (paper S5.1).") |
| out.append("4. Coding gate: SI coded by two independent raters per version; " |
| "no category proceeds at Cohen's kappa < 0.60 (paper S5.1).") |
| out.append(f"5. Input-separation gate: a pair enters Stage 3 only if its two " |
| f"members separate on SI by the registered minimum; failed pairs " |
| f"are counted, reported, and replaced from the build budget " |
| f"(thread post #5).") |
| out.append(f"6. Stage 2 (generative): >= {sc['model_families']} model " |
| f"families, shown-mode instruction, SI compared to the human " |
| f"shown target; n = {n2_std} pairs at the standard point, " |
| f"{n2_dec} at the decisive point (Table 1).") |
| out.append(f"7. Stage 3 (evaluative, core): pairs judged by the same models " |
| f"and by a human panel, order counterbalanced; n = {n_std} pairs " |
| f"per group at the standard point, {n_dec} at the decisive " |
| f"point (Table 1).") |
| out.append("8. Stage 4 (discriminant): Stage 3 re-run (a) length-mismatched " |
| "and (b) framing-varied; a preference that vanishes under either " |
| "is reclassified as that known bias (paper S5.4).") |
| out.append(f"9. Added arms, sized in Table 1: density calibration " |
| f"({nd_std} pairs), SI reader validation ({nr_std} items), " |
| f"register arm ({nreg_std} contents per stratum). No arm is " |
| f"registered without its own size from this table.") |
| out.append("10. Decision rules (paper S5.5): Stage 3 supports and Stage 4 " |
| "survives -> confirmed at the scope tested. Stage 3 falsifies " |
| "-> construct withdrawn, published with the same prominence as " |
| "a confirmation. Post-hoc qualifiers are prohibited.") |
| out.append("") |
| return "\n".join(out) + "\n" |
|
|
|
|
| KNOWN_SCENARIOS = ("null", "registered_alternative") |
|
|
|
|
| def check_scenario(cfg): |
| """Refuse an unknown truth.scenario by name, on any config that |
| carries one -- an unknown string must never silently select the |
| world a verdict is about. Configs without a truth block (the |
| sizing configs this program normally reads) are unaffected.""" |
| truth = cfg.get("truth") |
| if isinstance(truth, dict) and "scenario" in truth: |
| if truth["scenario"] not in KNOWN_SCENARIOS: |
| sys.exit(f"config refused: truth.scenario " |
| f"{truth['scenario']!r} is not a known scenario " |
| f"(expected 'null' or 'registered_alternative')") |
|
|
|
|
| def main(): |
| if len(sys.argv) != 2: |
| sys.exit("usage: python3 protocol_engine.py <config.json>") |
| with open(sys.argv[1], "r", encoding="utf-8") as fh: |
| cfg = json.load(fh) |
| check_scenario(cfg) |
| if "versions_per_pair" not in cfg.get("scale", {}): |
| cfg["scale"]["versions_per_pair"] = 2 |
| if "mc_reps" not in cfg: |
| cfg["mc_reps"] = 20000 |
| sys.stdout.write(run(cfg)) |
|
|
|
|
| if __name__ == "__main__": |
| main() |
|
|
