1#!/usr/bin/env python32"""Reproducible exact-sufficient-statistics Monte Carlo for rate-parameter exponential MLE.3Each worker is partitioned into four quarters; sums are sampled from their exact Gamma law.4No data movement/communication/network benchmark is claimed.5"""6import csv, numpy as np7from pathlib import Path8import matplotlib; matplotlib.use('Agg')9import matplotlib.pyplot as plt10ROOT=Path(__file__).parent11SEED=2026102812N=2**2013THETA=1.014REPS=50015M_LIST=(256,64,16)16rng=np.random.default_rng(SEED)1718def newton(x, mle):19 return 2.0*x-x*x/mle2021rows=[]; summary=[]22for m in M_LIST:23 K=N//m24 obs={key: [] for key in ('raw','one_step_raw','two_step_raw','jack2','one_step_jack2','jack3','one_step_jack3')}25 for rep in range(REPS):26 # Gamma is a *sum* of m/4 exponential(rate=theta) observations.27 quarters=rng.gamma(shape=m//4, scale=1/THETA,size=(K,4))28 halves=quarters.reshape(K,2,2).sum(axis=2)29 totals=halves.sum(axis=1)30 full_mle=N/np.sum(totals)31 local_full=m/totals32 local_halves=(m/2)/halves33 local_quarts=(m/4)/quarters34 raw=float(local_full.mean())35 J2=float((2*local_full - local_halves.mean(axis=1)).mean())36 J3=float(((8/3)*local_full -2*local_halves.mean(axis=1) +(1/3)*local_quarts.mean(axis=1)).mean())37 values={'raw':raw,38 'one_step_raw':newton(raw,full_mle),39 'two_step_raw':newton(newton(raw,full_mle),full_mle),40 'jack2':J2,'one_step_jack2':newton(J2,full_mle),41 'jack3':J3,'one_step_jack3':newton(J3,full_mle)}42 for key,val in values.items():43 delta=np.sqrt(N)*(val-full_mle)44 obs[key].append(delta)45 rows.append((m,K,rep,key,delta))46 for name, data in obs.items():47 a=np.array(data)48 summary.append(dict(m=m,K=K,method=name,mean_scaled_error=float(a.mean()),49 mc_standard_error=float(a.std(ddof=1)/np.sqrt(REPS)),50 rmse_scaled_error=float(np.sqrt(np.mean(a*a))),51 p95_abs_scaled_error=float(np.quantile(np.abs(a),.95))))52 print('m',m,'K',K,[(k,round(np.mean(v),6)) for k,v in obs.items()],flush=True)53with open(ROOT/'simulation_replicates.csv','w',newline='') as f:54 w=csv.writer(f); w.writerow(['m','K','rep','method','scaled_error_relative_to_centralized_mle']); w.writerows(rows)55with open(ROOT/'simulation_summary.csv','w',newline='') as f:56 w=csv.DictWriter(f,fieldnames=list(summary[0]));w.writeheader();w.writerows(summary)5758palette={'one_step_raw':'#A55040','two_step_raw':'#9B7C5A','one_step_jack2':'#176F86','one_step_jack3':'#407B58'}59labels={'one_step_raw':'Raw averaging + 1 Newton','two_step_raw':'Raw averaging + 2 Newton','one_step_jack2':'Half-block jackknife + 1 Newton','one_step_jack3':'Quarter-block jackknife + 1 Newton'}60fig,axs=plt.subplots(1,2,figsize=(11.5,4.0),layout='constrained')61ks=np.array([N//m for m in M_LIST])62for method in palette:63 vals=[next(r for r in summary if r['K']==int(k) and r['method']==method) for k in ks]64 axs[0].plot(ks,[abs(x['mean_scaled_error']) for x in vals],'-o',lw=2,ms=5,65 color=palette[method],label=labels[method])66 axs[1].plot(ks,[x['rmse_scaled_error'] for x in vals],'-o',lw=2,ms=5,color=palette[method])67for ax in axs:68 ax.set_xscale('log',base=2)69 ax.set_yscale('log')70 ax.set_xticks(ks)71 ax.set_xticklabels([f'{k:,}' for k in ks])72 ax.grid(alpha=.2,which='both')73 ax.set_xlabel('Number of partitions K')74axs[0].set_ylabel('Absolute mean scaled gap')75axs[1].set_ylabel('RMSE of scaled gap')76axs[0].set_title('Bias relative to centralized MLE')77axs[1].set_title('Error relative to centralized MLE')78fig.legend(loc='outside lower center',bbox_to_anchor=(.5,-.09),ncol=2,frameon=False,fontsize=9)79fig.savefig(ROOT/'figures'/'mc_comparison.pdf',bbox_inches='tight')80fig.savefig(ROOT/'figures'/'mc_comparison.png',dpi=170,bbox_inches='tight')81plt.close(fig)82print('wrote outputs')83