1"""Seeded runs are bit-reproducible, and the library reproduces the prior draft's exact-Gamma Monte Carlo numbers."""2import csv, pathlib3import numpy as np4from pamle.jackknife import vandermonde_weights56ROOT = pathlib.Path(__file__).resolve().parents[1]789def _prior_summary():10 f = ROOT / "references" / "initial_pamle" / "unpacked" / "simulation_summary.csv"11 rows = list(csv.DictReader(open(f)))12 return {(int(r['m']), r['method']): float(r['mean_scaled_error']) for r in rows}131415def _run_prior_design(reps=500, seed=20261028, N=2 ** 20, m_list=(256, 64, 16)):16 rng = np.random.default_rng(seed); w2 = vandermonde_weights((1, 2)); w3 = vandermonde_weights((1, 2, 4))17 out = {}18 for m in m_list:19 K = N // m; acc = {k: [] for k in ('raw', 'one_step_raw', 'two_step_raw', 'jack2', 'one_step_jack2', 'jack3', 'one_step_jack3')}20 for _ in range(reps):21 quarters = rng.gamma(shape=m // 4, scale=1.0, size=(K, 4))22 halves = quarters.reshape(K, 2, 2).sum(axis=2); totals = halves.sum(axis=1)23 full_mle = N / totals.sum()24 lf, lh, lq = m / totals, (m / 2) / halves, (m / 4) / quarters25 raw = lf.mean()26 J2 = (w2[0] * lf + w2[1] * lh.mean(axis=1)).mean()27 J3 = (w3[0] * lf + w3[1] * lh.mean(axis=1) + w3[2] * lq.mean(axis=1)).mean()28 newton = lambda x: 2 * x - x * x / full_mle29 vals = {'raw': raw, 'one_step_raw': newton(raw), 'two_step_raw': newton(newton(raw)), 'jack2': J2,30 'one_step_jack2': newton(J2), 'jack3': J3, 'one_step_jack3': newton(J3)}31 for k, v in vals.items():32 acc[k].append(np.sqrt(N) * (v - full_mle))33 for k, v in acc.items():34 out[(m, k)] = float(np.mean(v))35 return out363738def test_prior_draft_numbers_are_reproduced():39 prior = _prior_summary(); mine = _run_prior_design()40 for key, val in prior.items():41 assert abs(mine[key] - val) <= 1e-9 * max(1.0, abs(val)), (key, mine[key], val)424344def test_seeded_determinism():45 a = _run_prior_design(reps=5, m_list=(64,)); b = _run_prior_design(reps=5, m_list=(64,))46 assert a == b47