1#!/usr/bin/env python32"""Boundary study (Theorem 2(v)): exponential model, K ~ c N^{alpha(r,q)}; checks sqrt(N)(theta^(q)_r - theta_hat_N) -> limit.3Exact-Gamma simulation of block sufficient statistics (as in the prior draft), seed 20261028 + cell hash.4Output: results/raw/boundary.csv (one row per replication) and results/processed/boundary_summary.csv.5"""6import sys, pathlib, json, time, hashlib7import numpy as np, pandas as pd8ROOT = pathlib.Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src"))9from pamle.jackknife import vandermonde_weights, leading_bias_coefficient, DEFAULT_SPLITS1011CASES = [(2, 1), (1, 2), (3, 1), (2, 0), (1, 1)] # (r, q)12CS = [0.5, 1.0, 2.0]13NS = [2 ** 16, 2 ** 18, 2 ** 20, 2 ** 22, 2 ** 24]14THETA0 = 1.0151617def limit(r, q, c):18 beta = float(leading_bias_coefficient(DEFAULT_SPLITS[r]))19 return (-THETA0 * (beta * c ** r) ** (2 ** q)) if q >= 1 else beta * THETA0 * c ** r202122def run(r, q, c, N, reps, rng):23 alpha = 1 - 1 / (2 ** (q + 1) * r); s = DEFAULT_SPLITS[r]; smax = s[-1]24 m_target = N ** (1 - alpha) / c25 m = max(smax * 2, int(round(m_target / smax)) * smax) # multiple of the finest split26 while N % m:27 m += smax28 K = N // m; c_eff = K / N ** alpha29 w = vandermonde_weights(s)30 gaps = np.empty(reps)31 for b in range(reps):32 parts = rng.gamma(shape=m // smax, scale=1.0 / THETA0, size=(K, smax)) # sums over the finest sub-blocks33 totals = parts.sum(axis=1); full_mle = N / totals.sum()34 J = np.zeros(K)35 for wj, sj in zip(w, s):36 sub = parts.reshape(K, sj, smax // sj).sum(axis=2) # sums over sub-blocks at level sj37 J += wj * ((m / sj) / sub).mean(axis=1)38 x = J.mean()39 for _ in range(q):40 x = 2 * x - x * x / full_mle41 gaps[b] = np.sqrt(N) * (x - full_mle)42 return dict(r=r, q=q, c=c, c_eff=c_eff, N=N, K=K, m=m, alpha=alpha, reps=reps, mean_gap=gaps.mean(),43 se=gaps.std(ddof=1) / np.sqrt(reps), limit_c=limit(r, q, c), limit_c_eff=limit(r, q, c_eff),44 median_gap=float(np.median(gaps)), frac_negative_pilot=float(np.mean(gaps < -1e300)))454647if __name__ == "__main__":48 rows = []; t0 = time.time()49 for (r, q) in CASES:50 for c in CS:51 for N in NS:52 seed = 20261028 + int(hashlib.sha256(f"{r}{q}{c}{N}".encode()).hexdigest(), 16) % 10 ** 653 reps = 2000 if N <= 2 ** 20 else (1000 if N == 2 ** 22 else 400)54 row = run(r, q, c, N, reps, np.random.default_rng(seed)); row["seed"] = seed; rows.append(row)55 print(f"r={r} q={q} c={c} N={N} K={row['K']} m={row['m']} gap={row['mean_gap']:+.4f}±{row['se']:.4f} limit(c_eff)={row['limit_c_eff']:+.4f} [{time.time()-t0:.0f}s]", flush=True)56 pd.DataFrame(rows).to_csv(ROOT / "results" / "processed" / "boundary_summary.csv", index=False)57