1"""Monte Carlo summaries, confidence intervals and coverage."""2from __future__ import annotations3import numpy as np4from scipy import stats567def summarize(errors: np.ndarray):8 """errors: (B,) or (B,p) array -> dict with bias, sd, rmse, MC standard errors."""9 e = np.asarray(errors, float)10 B = e.shape[0]11 mean = e.mean(axis=0); sd = e.std(axis=0, ddof=1)12 rmse = np.sqrt(np.mean(e ** 2, axis=0))13 return {"n": B, "bias": mean, "sd": sd, "rmse": rmse, "bias_se": sd / np.sqrt(B),14 "rmse_se": np.sqrt(np.var(e ** 2, axis=0, ddof=1) / B) / (2 * np.maximum(rmse, 1e-300)),15 "p95_abs": np.quantile(np.abs(e), 0.95, axis=0)}161718def wald_ci(theta, cov, N, level=0.95):19 z = stats.norm.ppf(0.5 + level / 2)20 se = np.sqrt(np.diag(cov) / N)21 return theta - z * se, theta + z * se222324def coverage(theta_hat, cov_hat, N, theta0, levels=(0.90, 0.95, 0.99)):25 """theta_hat (B,p), cov_hat (B,p,p) [asymptotic covariance of sqrt(N)(theta-theta0)] -> coverage per level and param."""26 out = {}27 for lv in levels:28 z = stats.norm.ppf(0.5 + lv / 2)29 se = np.sqrt(np.einsum('bpp->bp', cov_hat) / N)30 hit = np.abs(theta_hat - theta0) <= z * se31 out[lv] = {"coverage": hit.mean(axis=0), "mean_length": (2 * z * se).mean(axis=0)}32 return out33