1import json
 2import math
 3import numpy as np
 4
 5
 6def skew_balance_check(seed=0, n=10000):
 7    rng = np.random.default_rng(seed)
 8    errors = []
 9    for _ in range(n):
10        U = rng.normal(0, 5)
11        du = rng.normal(0, 4)
12        p = abs(rng.normal()) + 1e-12
13        qf = p * math.exp(-du / 2.0)
14        qr = p * math.exp(du / 2.0)
15        errors.append(abs((math.log(qf) - U) -
16                          (math.log(qr) - (U + du))))
17    return float(max(errors)), float(np.mean(errors))
18
19
20def simulate(rho, events=160000, burn=5000, seed=1, nring=128):
21    rng = np.random.default_rng(seed)
22    p = rng.normal()
23    x = 0
24    signs, waits, positions = [], [], []
25    for k in range(events + burn):
26        rate = abs(p)
27        dt = rng.exponential(1.0 / rate) if rate > 1e-14 else 0.0
28        direction = 1 if p >= 0 else -1
29        x = (x + direction) % nring
30        p = rho * p + math.sqrt(max(0.0, 1.0 - rho * rho)) * rng.normal()
31        if k >= burn:
32            signs.append(direction)
33            waits.append(dt)
34            positions.append(x)
35    signs, waits, positions = map(np.asarray, (signs, waits, positions))
36    corr = float(np.mean(signs[:-1] * signs[1:]))
37    run_length = float(1.0 / (1.0 - corr))
38    centered = signs - signs.mean()
39    var = np.mean(centered * centered)
40    ac = [float(np.mean(centered[:-lag] * centered[lag:]) / var)
41          for lag in range(1, 100)]
42    sign_tau = float(1.0 + 2.0 * sum(a for a in ac if a > 0))
43    angle = 2 * math.pi * positions / nring
44    obs = np.cos(angle)
45    c = obs - obs.mean()
46    v = np.mean(c*c)
47    acx = [float(np.mean(c[:-lag] * c[lag:]) / v)
48           for lag in range(1, 100)]
49    x_tau = float(1.0 + 2.0 * sum(a for a in acx if a > 0))
50    return {
51        "rho": rho, "sign_corr": corr, "run_length": run_length,
52        "sign_tau": sign_tau, "position_tau": x_tau,
53        "mean_wait": float(waits.mean()),
54        "event_rate_per_time": float(len(waits) / waits.sum())
55    }
56
57
58def main():
59    max_err, mean_err = skew_balance_check()
60    # For Gaussian OU-refresh momenta, corr(sign_t, sign_{t+1})
61    # = 2 asin(rho)/pi. Hence mean consecutive same-sign run length is
62    # 1/(1-corr). These are direct mechanism predictions.
63    rows = []
64    for rho in [0.0, 0.25, 0.5, 0.75, 0.9]:
65        out = simulate(rho, seed=100 + int(100 * rho))
66        out["pred_corr"] = 2.0 / math.pi * math.asin(rho)
67        out["pred_run_length"] = 1.0 / (1.0 - out["pred_corr"])
68        out["corr_abs_error"] = abs(out["sign_corr"] - out["pred_corr"])
69        rows.append(out)
70    result = {
71        "skew_balance_max_abs_log_error": max_err,
72        "skew_balance_mean_abs_log_error": mean_err,
73        "prediction_rows": rows,
74        "baseline_rho0_position_tau": rows[0]["position_tau"],
75        "idea_rho075_position_tau": rows[3]["position_tau"],
76        "note": "Event-rate is reported empirically; per-event Gaussian refresh has heavy-tailed waiting times, so E|p| is not the event rate."
77    }
78    with open("results.json", "w") as f:
79        json.dump(result, f, indent=2)
80    print(json.dumps(result, indent=2))
81
82
83if __name__ == "__main__":
84    main()