1import json
  2import math
  3import numpy as np
  4
  5# IMM stale-feedback detector for the known scalar linear process
  6# z[t+1] = a*z[t] + w[t], y[t] = z[t-delay] + v[t].
  7# In this verification, the process state history is known to isolate the
  8# paper's exact likelihood/posterior mechanism from online system-ID error.
  9
 10def trial(delay_true, noise, a=0.85, q=0.04, D=4, n=1200, seed=0,
 11          persistence=0.985, threshold=0.01):
 12    rng = np.random.default_rng(seed)
 13    z = np.zeros(n + D + 1)
 14    # Stationary-ish initialization avoids an artificial transient.
 15    z[0] = rng.normal(0, math.sqrt(q / (1 - a*a)))
 16    for t in range(n + D):
 17        z[t + 1] = a * z[t] + rng.normal(0, math.sqrt(q))
 18    y = np.array([z[D + t - delay_true] + rng.normal(0, math.sqrt(noise))
 19                  for t in range(n)])
 20
 21    m = D + 1
 22    T = np.full((m, m), (1 - persistence) / (m - 1))
 23    np.fill_diagonal(T, persistence)
 24    pi = np.ones(m) / m
 25    posterior = []
 26    logodds = []
 27    kl = []
 28    llr = []
 29    norm_error = []
 30    alarm = None
 31    run = 0
 32    for t in range(n):
 33        # IMM mode mixing prior c_i = sum_j T_ji*pi_j.
 34        c = T.T @ pi
 35        means = z[D + t - np.arange(m)]
 36        ll = -0.5 * ((y[t] - means)**2 / noise + math.log(2 * math.pi * noise))
 37        un = c * np.exp(ll - np.max(ll))
 38        pi = un / un.sum()
 39        posterior.append(pi.copy())
 40        norm_error.append(abs(pi.sum() - 1.0))
 41        logodds.append(math.log(max(pi[delay_true], 1e-300) /
 42                              max(pi[0], 1e-300)))
 43        # Equal-variance Gaussian innovation KL for true vs no-delay mode.
 44        delta = means[delay_true] - means[0]
 45        kl.append(delta * delta / (2 * noise))
 46        # Raw cumulative innovation log-likelihood ratio avoids posterior saturation.
 47        llr.append(((y[t] - means[0])**2 - (y[t] - means[delay_true])**2) / (2 * noise))
 48        run = run + 1 if pi[0] <= threshold else 0
 49        if alarm is None and run >= 3:
 50            alarm = t + 1
 51
 52    posterior = np.asarray(posterior)
 53    logodds = np.asarray(logodds)
 54    kl = np.asarray(kl)
 55    llr = np.asarray(llr)
 56    start = 10
 57    # Use an early/mid window before numerical posterior saturation.
 58    stop = min(n, 250)
 59    slope = float(np.polyfit(np.arange(start, stop), logodds[start:stop], 1)[0])
 60    cumulative_llr = np.cumsum(llr)
 61    llr_slope = float(np.polyfit(np.arange(start, stop), cumulative_llr[start:stop], 1)[0])
 62    mean_kl = float(np.mean(kl[start:stop]))
 63    dom = next((i + 1 for i, p in enumerate(posterior)
 64                if p[delay_true] > .5 and p[delay_true] > p[0]), None)
 65    return dict(slope=slope, llr_slope=llr_slope, mean_kl=mean_kl, ratio=llr_slope / max(mean_kl, 1e-12),
 66                alarm=alarm, dominant=dom, final_true=float(posterior[-1, delay_true]),
 67                final_zero=float(posterior[-1, 0]), max_norm_error=max(norm_error))
 68
 69
 70def aggregate(delay, noise, reps=30, **kw):
 71    vals = [trial(delay, noise, seed=7000 + k, **kw) for k in range(reps)]
 72    def med(k):
 73        x = [v[k] for v in vals if v[k] is not None]
 74        return None if not x else float(np.median(x))
 75    return {"delay": delay, "noise": noise, "slope": med("slope"),
 76            "mean_KL": med("mean_kl"), "cumulative_llr_slope": med("llr_slope"), "slope_over_KL": med("ratio"),
 77            "alarm_median": med("alarm"), "dominant_median": med("dominant"),
 78            "final_true_posterior": med("final_true"),
 79            "final_zero_posterior": med("final_zero"),
 80            "max_normalization_error": max(v["max_norm_error"] for v in vals)}
 81
 82
 83def main():
 84    # Prediction 1: persistent delay separates from no-delay.
 85    delayed = aggregate(2, .03)
 86    clean = aggregate(0, .03)
 87    # Prediction 2: expected log posterior-odds slope is the innovation KL.
 88    kl_sweep = [aggregate(2, r) for r in (.015, .03, .06, .12)]
 89    # Prediction 3: lower KL takes longer to cross the fixed alarm threshold.
 90    detect_sweep = [aggregate(2, r, reps=20) for r in (.015, .03, .06, .12)]
 91    out = {
 92      "prediction_checks": {
 93        "P1_mode_separation": {
 94          "prediction": "persistent d=2 gives p(d=2)>0.5 while clean d=0 retains p(d=0)",
 95          "delayed": delayed, "clean": clean},
 96        "P2_KL_slope": {
 97          "prediction": "early cumulative innovation log-likelihood slope approximately equals mean innovation KL",
 98          "noise_sweep": kl_sweep},
 99        "P3_detection_scaling": {
100          "prediction": "alarm time increases as KL decreases (larger observation noise)",
101          "noise_sweep": detect_sweep},
102        "P4_normalization": {
103          "prediction": "posterior sums to one to numerical precision",
104          "max_abs_sum_error": max(x["max_normalization_error"] for x in [delayed, clean] + kl_sweep)}
105      },
106      "settings": {"a": .85, "q": .04, "D": 4, "persistence": .985,
107                   "steps": 1200, "threshold": .01, "replicates": 30,
108                   "likelihood": "exact Gaussian innovation with known state history"}
109    }
110    print(json.dumps(out, indent=2))
111
112if __name__ == '__main__':
113    main()