1import json, math
 2import numpy as np
 3from scipy.integrate import quad
 4from shape_private_noise import select_shape, hockey_divergence, sample_noise, moment
 5
 6
 7def monte_carlo_moments(p, b, seed=7, n=500000):
 8    z = sample_noise(np.random.default_rng(seed), p, b, n)
 9    return float(np.mean(np.abs(z))), float(np.mean(z*z)), float(np.std(z))
10
11
12def numerical_divergence_mc(p, b, d, eps, seed=9, n=1000000):
13    # Importance sampling under P: D_e(P||Q)=E_P[(1-exp(eps+log q-log p))_+].
14    rng = np.random.default_rng(seed)
15    x = sample_noise(rng, p, b, n)
16    lp = math.log(p) - math.log(2*b) - math.lgamma(1/p) - (np.abs(x)/b)**p
17    lq = math.log(p) - math.log(2*b) - math.lgamma(1/p) - (np.abs(x-d)/b)**p
18    val = np.maximum(0, 1-np.exp(np.minimum(700, eps+lq-lp)))
19    return float(np.mean(val)), float(np.std(val)/math.sqrt(n))
20
21
22def main():
23    grid=[1,1.25,1.5,1.75,2,2.5,3,4,6,8]
24    regimes=[(0.1,1e-5),(0.5,1e-5),(1,1e-5),(2,1e-5),(1,1e-3),(1,1e-2),(5,1e-5)]
25    sweep=[]
26    for eps,delta in regimes:
27        best, rows=select_shape(epsilon=eps,delta=delta,grid=grid)
28        lap=next(x for x in rows if x['p']==1)
29        gauss=next(x for x in rows if x['p']==2)
30        sweep.append({'epsilon':eps,'delta':delta,'best':best,
31                      'laplace_moment':lap['moment'],'gaussian_moment':gauss['moment'],
32                      'gain_vs_laplace':1-best['moment']/lap['moment'],
33                      'gain_vs_gaussian':1-best['moment']/gauss['moment']})
34    b,rows=select_shape(epsilon=1,delta=1e-5,grid=grid)
35    moment_check=[]
36    for r in rows[:3]:
37        mc1,mc2,sd=monte_carlo_moments(r['p'],r['b'])
38        moment_check.append({'p':r['p'],'formula_m1':moment(r['p'],r['b'],1),
39          'mc_m1':mc1,'formula_m2':moment(r['p'],r['b'],2),'mc_m2':mc2,'std':sd,
40          'quad_divergence':r['divergence'],'mc_divergence':numerical_divergence_mc(r['p'],r['b'],1,1)})
41    result={'sweep':sweep,'checks':moment_check}
42    open('sweep_results.json','w').write(json.dumps(result,indent=2))
43    print(json.dumps(result,indent=2))
44main()