import json, time from pathlib import Path import numpy as np SEED = 2399 rng = np.random.default_rng(SEED) def phi(x, K=8): x = np.asarray(x).reshape(-1) # Frequencies are deterministic, as required by the proposal. out = [np.ones_like(x)] for k in range(1, K + 1): out += [np.cos(np.pi * k * x), np.sin(np.pi * k * x)] return np.stack(out, axis=1) def true_fn(x): x = np.asarray(x) return 0.4 + 0.8*np.sin(np.pi*x) - 0.35*np.cos(2*np.pi*x) + 0.25*np.sin(4*np.pi*x) def update_stats(X, y, lam=1e-2, K=8): P = phi(X, K) V = lam*np.eye(P.shape[1]) + P.T @ P b = P.T @ y return V, b def predict(V, b, X, sigma, beta=2.5, K=8): P = phi(X, K) theta = np.linalg.solve(V, b) # solve once, then form diagonal quadratic forms VinvP = np.linalg.solve(V, P.T).T q = np.maximum(0.0, np.sum(P * VinvP, axis=1)) mu = P @ theta sd = sigma*np.sqrt(q) return mu, sd def online_stats(X, y, lam=1e-2, K=8): m = 2*K+1 V = lam*np.eye(m) b = np.zeros(m) for x, yy in zip(X, y): p = phi([x], K)[0] V += np.outer(p, p) b += p*yy return V, b def contraction_sweep(): # Repeated identical observations give q_n=q_0/(1+n*q_0/lambda-like scale), # hence posterior sd is asymptotically proportional to n^-1/2. sigma, lam, K = 0.12, 1e-2, 8 x0 = 0.15 ns = np.array([1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]) sds = [] for n in ns: X = np.full(n, x0) y = true_fn(X) # y values do not affect posterior variance V, b = update_stats(X, y, lam, K) sds.append(predict(V, b, [x0], sigma, K=K)[1][0]) # Fit slope only in the data-dominated regime. slope = np.polyfit(np.log(ns[2:]), np.log(sds[2:]), 1)[0] # Uncovered point after observations at x0. X = np.full(1024, x0); y = true_fn(X) V, b = update_stats(X, y, lam, K) covered = predict(V, b, [x0], sigma, K=K)[1][0] uncovered = predict(V, b, [-0.93], sigma, K=K)[1][0] return {"n": ns.tolist(), "sd": np.array(sds).tolist(), "loglog_slope": float(slope), "predicted_slope": -0.5, "covered_sd_n1024": float(covered), "uncovered_sd_n1024": float(uncovered), "uncovered_to_covered": float(uncovered/covered)} def coverage_sweep(): # Repeated trials assess simultaneous grid coverage of beta*s. K, lam, sigma, beta = 8, 1e-2, 0.12, 2.5 grid = np.linspace(-1, 1, 201) rows = [] for n, region in [(40, "concentrated"), (160, "concentrated"), (40, "broad"), (160, "broad")]: covers = [] widths = [] rmses = [] for t in range(100): rr = np.random.default_rng(SEED + 10000 + t + n + (0 if region == "concentrated" else 1000)) X = rr.uniform(-0.25, 0.25, n) if region == "concentrated" else rr.uniform(-1, 1, n) y = true_fn(X) + rr.normal(0, sigma, n) V, b = update_stats(X, y, lam, K) mu, sd = predict(V, b, grid, sigma, beta, K) err = np.abs(true_fn(grid)-mu) covers.append(float(np.all(err <= beta*sd))) widths.append(float(np.mean(2*beta*sd))) rmses.append(float(np.sqrt(np.mean((true_fn(grid)-mu)**2)))) rows.append({"n": n, "region": region, "simultaneous_coverage": float(np.mean(covers)), "mean_interval_width": float(np.mean(widths)), "grid_rmse": float(np.mean(rmses)), "target_coverage": 0.95, "beta": beta}) return rows def baseline_comparison(): # Same feature family and same data: baseline gives only a point prediction; # proposed head additionally supplies a calibrated, location-dependent interval. rr = np.random.default_rng(SEED) n = 160 X = rr.uniform(-0.25, 0.25, n) sigma = 0.12 y = true_fn(X) + rr.normal(0, sigma, n) grid = np.linspace(-1, 1, 401) V, b = update_stats(X, y, 1e-2, 8) mu, sd = predict(V, b, grid, sigma, beta=2.5, K=8) baseline_rmse = float(np.sqrt(np.mean((true_fn(grid)-mu)**2))) # deterministic ridge has identical mean but no uncertainty/rejection signal feature_dim = 17 baseline_interval_width = 0.0 proposed_width = float(np.mean(2*2.5*sd)) unsafe_fraction = float(np.mean((mu - 2.5*sd > true_fn(grid)) | (mu + 2.5*sd < true_fn(grid)))) return {"deterministic_ridge_grid_rmse": baseline_rmse, "certificate_grid_rmse": baseline_rmse, "deterministic_interval_width": baseline_interval_width, "certificate_mean_interval_width": proposed_width, "certificate_grid_miss_fraction": unsafe_fraction, "feature_dim": feature_dim} def main(): # Algebraic/numerical identity: rank-one updates and batch sufficient statistics. rr = np.random.default_rng(SEED) X = rr.uniform(-1, 1, 300); y = true_fn(X) + rr.normal(0, .12, len(X)) Vb, bb = update_stats(X, y) Vo, bo = online_stats(X, y) identity = {"max_abs_V_difference": float(np.max(np.abs(Vb-Vo))), "max_abs_b_difference": float(np.max(np.abs(bb-bo))), "max_abs_theta_difference": float(np.max(np.abs(np.linalg.solve(Vb,bb)-np.linalg.solve(Vo,bo))))} result = {"seed": SEED, "identity_check": identity, "contraction_sweep": contraction_sweep(), "coverage_sweep": coverage_sweep(), "baseline_comparison": baseline_comparison(), "parameter_scaling_sweeps": parameter_scaling_sweeps()} Path("results.json").write_text(json.dumps(result, indent=2)) print(json.dumps(result, indent=2)) def parameter_scaling_sweeps(): rr = np.random.default_rng(SEED + 77) X = rr.uniform(-1, 1, 240) y0 = true_fn(X) + rr.normal(0, 1, len(X)) # Prediction: posterior s is exactly linear in the supplied noise scale sigma. V, b = update_stats(X, y0) sigmas = np.array([0.03, 0.06, 0.12, 0.24]) mean_s = [] for sig in sigmas: mean_s.append(float(np.mean(predict(V, b, X, sig, K=8)[1]))) sigma_ratio = np.array(mean_s) / mean_s[2] expected_sigma_ratio = sigmas / sigmas[2] # Prediction: increasing lambda increases uncertainty, while reducing data fit strength. lambdas = np.array([1e-3, 1e-2, 1e-1, 1.0]) grid = np.linspace(-1, 1, 201) lambda_s = [] lambda_rmse = [] for lam in lambdas: Vl, bl = update_stats(X, true_fn(X) + rr.normal(0, .12, len(X)), lam) mu, sd = predict(Vl, bl, grid, .12, K=8) lambda_s.append(float(np.mean(sd))) lambda_rmse.append(float(np.sqrt(np.mean((mu-true_fn(grid))**2)))) return {"sigma": sigmas.tolist(), "mean_sd": mean_s, "observed_sigma_ratios": sigma_ratio.tolist(), "predicted_sigma_ratios": expected_sigma_ratio.tolist(), "lambda": lambdas.tolist(), "mean_sd_by_lambda": lambda_s, "grid_rmse_by_lambda": lambda_rmse} if __name__ == "__main__": main()