1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def conformal_quantile(scores, alpha):
  7    scores = np.sort(np.asarray(scores, dtype=float))
  8    m = len(scores)
  9    k = int(math.ceil((m + 1) * (1 - alpha)))
 10    return float("inf") if k > m else float(scores[k - 1])
 11
 12
 13def kalman_run(y, q=0.04, r=0.25, gate=None, r_inflate=1.0):
 14    xhat, p = 0.0, 1.0
 15    estimates, innovations, triggered = [], [], []
 16    for obs in y:
 17        p_pred = p + q
 18        innov = float(obs - xhat)
 19        s = p_pred + r
 20        hit = gate is not None and abs(innov) / math.sqrt(s) > gate
 21        # Correct robust intervention: distrust the current suspicious observation.
 22        r_used = r * r_inflate if hit else r
 23        gain = p_pred / (p_pred + r_used)
 24        xhat = xhat + gain * innov
 25        p = (1.0 - gain) * p_pred
 26        estimates.append(xhat)
 27        innovations.append(innov / math.sqrt(s))
 28        triggered.append(bool(hit))
 29    return np.asarray(estimates), np.asarray(innovations), np.asarray(triggered)
 30
 31
 32def make_sequence(rng, n, q=0.04, r=0.25, outliers=False, outlier_prob=.08, outlier_scale=4.):
 33    x = np.zeros(n)
 34    for t in range(1, n):
 35        x[t] = x[t-1] + rng.normal(0, math.sqrt(q))
 36    y = x + rng.normal(0, math.sqrt(r), n)
 37    if outliers:
 38        mask = rng.random(n) < outlier_prob
 39        y[mask] += rng.normal(0, outlier_scale * math.sqrt(r), mask.sum())
 40    return x, y
 41
 42
 43def collect_scores(rng, n_scores, score_kind="state"):
 44    vals = []
 45    while len(vals) < n_scores:
 46        x, y = make_sequence(rng, 50)
 47        est, inn, _ = kalman_run(y)
 48        vals.extend(np.abs(x - est) if score_kind == "state" else np.abs(inn))
 49    return np.asarray(vals[:n_scores])
 50
 51
 52def coverage_for_q(rng, q, n_sequences=300):
 53    hits = []
 54    for _ in range(n_sequences):
 55        x, y = make_sequence(rng, 50)
 56        est, _, _ = kalman_run(y)
 57        hits.extend(np.abs(x - est) <= q)
 58    return float(np.mean(hits))
 59
 60
 61def run(seed=902):
 62    rng = np.random.default_rng(seed)
 63    alpha_values = (.05, .10, .20, .30)
 64    cal = collect_scores(rng, 5000)
 65    sigma = float(np.std(cal))
 66    # For |N(0,sigma^2)|, the exact reference is sigma*Phi^{-1}(1-alpha/2).
 67    zhalf = {0.05: 1.959964, 0.10: 1.644854, 0.20: 1.281552, 0.30: 1.036433}
 68    qvals = {str(a): conformal_quantile(cal, a) for a in alpha_values}
 69    coverage = {str(a): coverage_for_q(rng, qvals[str(a)]) for a in alpha_values}
 70    halfnormal_reference = {str(a): sigma * zhalf[a] for a in alpha_values}
 71
 72    # Calibration-size prediction: rank quantile converges with approximately O(m^-1/2)
 73    # sampling variation, while expected coverage remains near 1-alpha.
 74    size_sweep = {}
 75    for m in (100, 500, 2000, 5000):
 76        qs = []
 77        covs = []
 78        for rep in range(20):
 79            local = np.random.default_rng(seed + 10000 + m * 10 + rep)
 80            scores = collect_scores(local, m)
 81            q = conformal_quantile(scores, .10)
 82            qs.append(q)
 83            covs.append(coverage_for_q(local, q, 100))
 84        size_sweep[str(m)] = {
 85            "q_mean": float(np.mean(qs)), "q_std": float(np.std(qs)),
 86            "coverage_mean": float(np.mean(covs)), "coverage_std": float(np.std(covs)),
 87            "predicted_q_std_scaling": 1.0 / math.sqrt(m)
 88        }
 89
 90    # Deployment intervention: clean innovation gate, then contaminated observations.
 91    gate_scores = collect_scores(rng, 5000, "innovation")
 92    gate = conformal_quantile(gate_scores, .10)
 93    baseline, gated, triggers = [], [], []
 94    for _ in range(300):
 95        x, y = make_sequence(rng, 50, outliers=True, outlier_prob=.08, outlier_scale=4.)
 96        eb, _, _ = kalman_run(y)
 97        eg, _, tr = kalman_run(y, gate=gate, r_inflate=8.0)
 98        baseline.append(np.sqrt(np.mean((x-eb) ** 2)))
 99        gated.append(np.sqrt(np.mean((x-eg) ** 2)))
100        triggers.append(np.mean(tr))
101    result = {
102        "seed": seed, "calibration_scores": len(cal), "alpha": list(alpha_values),
103        "state_quantiles": qvals, "state_coverage": coverage,
104        "residual_std": sigma, "halfnormal_reference": halfnormal_reference,
105        "calibration_size_sweep": size_sweep,
106        "innovation_threshold": gate,
107        "contaminated_rmse_baseline": float(np.mean(baseline)),
108        "contaminated_rmse_gated": float(np.mean(gated)),
109        "contaminated_rmse_improvement_fraction": float(1 - np.mean(gated) / np.mean(baseline)),
110        "trigger_rate": float(np.mean(triggers)),
111        "predictions": {
112            "coverage": "exchangeable marginal coverage is at least approximately 1-alpha",
113            "alpha_scaling": "q decreases with alpha and follows the half-normal reference",
114            "calibration_scaling": "quantile variability decreases approximately as m^(-1/2)",
115            "robust_intervention": "R inflation on innovation-triggered observations reduces error under outliers"
116        }
117    }
118    Path("results.json").write_text(json.dumps(result, indent=2))
119    print(json.dumps(result, indent=2))
120
121
122if __name__ == "__main__":
123    run()