1import json
 2import numpy as np
 3from scipy.integrate import quad
 4from nonlinear_drift import scalar_pairwise_rate
 5
 6LAM, D = 1.0, 0.5
 7
 8
 9def exact_variance(beta):
10    def unnorm(x):
11        return np.exp(-(LAM*x*x/2.0 + beta*x**4/4.0)/D)
12    z = quad(unnorm, -np.inf, np.inf, epsabs=1e-12, epsrel=1e-12)[0]
13    m2 = quad(lambda x: x*x*unnorm(x), -np.inf, np.inf,
14              epsabs=1e-12, epsrel=1e-12)[0] / z
15    return m2
16
17
18def simulate(beta, seed=123, n=4000, burn=3000, sample_steps=7000, dt=.005):
19    rng = np.random.default_rng(seed)
20    x = np.zeros(n)
21    for _ in range(burn):
22        x += (-LAM*x - beta*x**3)*dt + np.sqrt(2*D*dt)*rng.normal(size=n)
23    vals = []
24    for k in range(sample_steps):
25        x += (-LAM*x - beta*x**3)*dt + np.sqrt(2*D*dt)*rng.normal(size=n)
26        if k % 10 == 0:
27            vals.append(np.mean(x*x))
28    return float(np.mean(vals)), float(np.std(vals) / np.sqrt(len(vals)))
29
30
31def main():
32    wide_betas = np.array([0., .02, .05, .1, .2, .5, 1.0])
33    exact_wide = np.array([exact_variance(b) for b in wide_betas])
34    baseline = D / LAM
35    # First-order perturbation: Var/baseline = 1 - 3 beta D/lambda^2 + O(beta^2).
36    tiny_betas = np.array([0., .0001, .0005, .001, .002, .005])
37    exact_tiny = np.array([exact_variance(b) for b in tiny_betas])
38    relative = exact_tiny / baseline - 1.0
39    slope = float(np.polyfit(tiny_betas, relative, 1)[0])
40    predicted_slope = -3 * D / LAM**2
41
42    # Prediction 1: pairwise rate is exactly >= lambda, equality only at x=y=0.
43    rng = np.random.default_rng(7)
44    x, y = rng.normal(size=100000), rng.normal(size=100000)
45    rates = scalar_pairwise_rate(x, y, LAM, .5)
46    contraction = {
47        "predicted_min_rate": LAM,
48        "observed_min_rate": float(rates.min()),
49        "observed_fraction_below_lambda": float(np.mean(rates < LAM - 1e-12)),
50        "rate_at_equal_radius_1": float(scalar_pairwise_rate(1., 1., LAM, .5)),
51        "predicted_rate_at_equal_radius_1": LAM + 3*.5
52    }
53    sims = {str(b): simulate(float(b), seed=100 + i) for i, b in enumerate([0., .5])}
54    out = {
55        "parameters": {"lambda": LAM, "D": D, "baseline_variance": baseline},
56        "prediction_1_contraction": contraction,
57        "prediction_2_small_beta_slope": {
58            "prediction": "relative variance slope = -3D/lambda^2 + O(beta)",
59            "predicted_relative_slope": predicted_slope,
60            "observed_relative_slope": slope,
61            "relative_error": abs(slope - predicted_slope) / abs(predicted_slope),
62            "betas": tiny_betas.tolist(),
63            "exact_variances": exact_tiny.tolist()
64        },
65        "prediction_3_monotonic_variance": {
66            "prediction": "variance strictly decreases as beta increases",
67            "betas": wide_betas.tolist(), "exact_variances": exact_wide.tolist(),
68            "strictly_decreasing": bool(np.all(np.diff(exact_wide) < 0)),
69            "relative_gap_at_beta_0.5": float(1 - exact_wide[5] / baseline)
70        },
71        "matched_euler_maruyama": {
72            "linear_beta_0": sims["0.0"], "cubic_beta_0.5": sims["0.5"],
73            "simulated_relative_gap": float(1 - sims["0.5"][0] / sims["0.0"][0]),
74            "dt": .005, "burn_steps": 3000, "sample_steps": 7000, "n_parallel": 4000
75        }
76    }
77    with open("results.json", "w") as f:
78        json.dump(out, f, indent=2)
79    print(json.dumps(out, indent=2))
80
81if __name__ == "__main__":
82    main()