import json, numpy as np, pandas as pd, affect_sim as s from scipy.stats import chi2, t as tdist OUT = "../DATA/" res = {} # 1. Main Monte Carlo comparison (nominal: noise x1, dropout 0.2, outliers 0.05) NMC = 200 rows = [] rng = np.random.default_rng(2026) for i in range(NMC): r = s.run_methods(rng, 1.0, 0.2, 0.05) for m, v in r.items(): rows.append(dict(run=i, method=m, V=v["rmse"][0], A=v["rmse"][1], D=v["rmse"][2], overall=v["overall"], quad=v["oct"])) df = pd.DataFrame(rows); df.to_csv(OUT + "mc_main.csv", index=False) def ci(x): x = np.asarray(x); h = tdist.ppf(0.975, len(x) - 1) * x.std(ddof=1) / np.sqrt(len(x)) return float(x.mean()), float(h) summ = {} for m in s.METHODS: d = df[df.method == m] summ[m] = {k: ci(d[k]) for k in ["V", "A", "D", "overall", "quad"]} res["main"] = summ # paired comparison Fused KF vs Late fusion+EMA and gated vs ungated from scipy.stats import ttest_rel, wilcoxon a = df[df.method == "Fused KF"].overall.values; b = df[df.method == "Late fusion+EMA"].overall.values g = df[df.method == "Gated fused KF"].overall.values res["paired"] = dict(kf_vs_ema_meandiff=float((b - a).mean()), kf_vs_ema_rel=float(((b - a) / b).mean()), kf_vs_ema_p=float(wilcoxon(a, b).pvalue), kf_better_frac=float((a < b).mean()), gated_vs_kf_meandiff=float((g - a).mean()), gated_vs_kf_p=float(wilcoxon(g, a).pvalue)) # 2. Monte Carlo convergence of mean overall RMSE (Fused KF) conv = [] for n in [25, 50, 100, 200]: m, h = ci(a[:n]); conv.append(dict(n=n, mean=m, ci95=h, rel=h / m)) res["convergence"] = conv # 3. Dropout sweep (noise x1, outliers 0.05) sweep = [] rng = np.random.default_rng(7) for p in [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]: for i in range(60): r = s.run_methods(rng, 1.0, p, 0.05, which=["Late fusion+EMA", "Fused KF", "Gated fused KF", "Face-only KF"]) for m, v in r.items(): sweep.append(dict(p=p, method=m, overall=v["overall"])) sw = pd.DataFrame(sweep); sw.to_csv(OUT + "dropout_sweep.csv", index=False) # 4. Verification: NEES with matched truth model (6-dim augmented state) rng = np.random.default_rng(11) nees_all = [] for i in range(50): x, mu = s.gen_truth(rng, matched=True) obs = s.gen_obs(rng, x, 1.0, 0.2, 0.0) o = s.kalman(obs, return_P=True) e = np.hstack([o["x"] - x, o["mu"] - mu]) nees = np.einsum("ki,kij,kj->k", e, np.linalg.inv(o["P"]), e) nees_all.append(nees) nees_all = np.array(nees_all) # [runs, N] avg_t = nees_all.mean(0) lo, hi = chi2.ppf([0.025, 0.975], 6 * 50) / 50 res["nees"] = dict(mean=float(avg_t[50:].mean()), lo=float(lo), hi=float(hi), frac_inside=float(np.mean((avg_t[50:] > lo) & (avg_t[50:] < hi)))) np.save(OUT + "nees_avg.npy", avg_t) # 5. 2^4 full-factorial DOE, 10 replicates rng = np.random.default_rng(99) levels = dict(noise=(0.5, 1.5), drop=(0.1, 0.4), out=(0.0, 0.10), gate=(0, 1)) doe = [] for rep in range(10): for A in (-1, 1): for B in (-1, 1): for C in (-1, 1): for Dg in (-1, 1): ns = levels["noise"][A > 0]; pdp = levels["drop"][B > 0]; po = levels["out"][C > 0] x, mu = s.gen_truth(rng); obs = s.gen_obs(rng, x, ns, pdp, po) e = s.kalman(obs, gate=bool(Dg > 0), noise_scale=ns)["x"] r = s.rmse(e, x) doe.append(dict(rep=rep, X1_noise=A, X2_dropout=B, X3_outlier=C, X4_gating=Dg, noise_scale=ns, p_drop=pdp, p_out=po, gating=int(Dg > 0), Y_rmse=float(np.sqrt(np.mean(r ** 2))))) dd = pd.DataFrame(doe); dd.to_csv(OUT + "doe_minitab.csv", index=False) import statsmodels.formula.api as smf from statsmodels.stats.anova import anova_lm fit = smf.ols("Y_rmse ~ X1_noise*X2_dropout*X3_outlier*X4_gating", data=dd).fit() an = anova_lm(fit, typ=2) coef = fit.params.drop("Intercept"); se = fit.bse.drop("Intercept") tstat = (coef / se) res["doe"] = dict(r2=float(fit.rsquared), r2adj=float(fit.rsquared_adj), F=float(fit.fvalue), Fp=float(fit.f_pvalue), dfres=float(fit.df_resid), effects={k: dict(effect=float(2 * coef[k]), t=float(tstat[k]), p=float(fit.pvalues[k])) for k in coef.index}, anova={k: dict(SS=float(an.loc[k, "sum_sq"]), df=float(an.loc[k, "df"]), F=(None if np.isnan(an.loc[k, "F"]) else float(an.loc[k, "F"])), p=(None if np.isnan(an.loc[k, "PR(>F)"]) else float(an.loc[k, "PR(>F)"]))) for k in an.index}) # reduced model for residual check fit2 = smf.ols("Y_rmse ~ X1_noise*X2_dropout*X3_outlier*X4_gating", data=dd).fit() from scipy.stats import shapiro res["doe"]["shapiro_p"] = float(shapiro(fit2.resid).pvalue) np.save(OUT + "doe_resid.npy", fit2.resid.values) cells = dd.groupby(["X3_outlier", "X4_gating"]).Y_rmse.mean() res["doe"]["cell_out_gate"] = {f"{k[0]},{k[1]}": float(v) for k, v in cells.items()} # 6. Behaviour generation trade-off (lambda sweep), 30 runs, nominal conditions rng = np.random.default_rng(5) lams = [0, 1, 2, 4, 8, 16, 32, 64] trials = [] for i in range(30): x, mu = s.gen_truth(rng); obs = s.gen_obs(rng, x, 1.0, 0.2, 0.05) o = s.kalman(obs, gate=True, return_P=True); trials.append((x, o)) beh = [] for lam in lams: for rate in (False, True): for tau in (None, 0.15): mets = [s.behaviour(o["x"], o["P"], x, lam, tau, rate)[1] for x, o in trials] beh.append(dict(lam=lam, rate=rate, tau=tau, **{k: float(np.mean([m[k] for m in mets])) for k in mets[0]})) bh = pd.DataFrame(beh); bh.to_csv(OUT + "behaviour_sweep.csv", index=False) # baseline: late-fusion+EMA estimate, lam=0, no rate, no conf bmets = [] for x, o in trials: pass res["behaviour"] = beh json.dump(res, open(OUT + "results.json", "w"), indent=1) print(json.dumps({k: res[k] for k in ["main", "paired", "convergence", "nees"]}, indent=1)) print("DOE", res["doe"]["r2"], res["doe"]["r2adj"], res["doe"]["F"], res["doe"]["shapiro_p"]) for k, v in res["doe"]["effects"].items(): print(k, round(v["effect"], 4), round(v["t"], 2), f"{v['p']:.2e}") print(res["doe"]["cell_out_gate"]) print(bh.to_string()) print(sw.groupby(["p", "method"]).overall.mean().unstack().round(3))