1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6SEED = 2747
  7
  8
  9def empirical_cvar_lower(z, gamma):
 10    """Max over eta of eta - (gamma*K)^-1 sum(max(eta-z,0))."""
 11    z = np.asarray(z, dtype=float)
 12    K = len(z)
 13    if not (0 < gamma <= 1) or K == 0:
 14        raise ValueError("gamma and samples")
 15    # The objective is piecewise linear; a maximizer is an order statistic.
 16    candidates = np.r_[z, z.min() - 1.0, z.max() + 1.0]
 17    vals = np.array([eta - np.maximum(eta - z, 0).sum() / (gamma * K)
 18                     for eta in candidates])
 19    # Return the optimized objective, not eta itself.
 20    return float(vals.max())
 21
 22
 23def empirical_bottom_cvar(z, gamma):
 24    """Equivalent weighted mean of the lower empirical tail (fraction gamma)."""
 25    z = np.sort(np.asarray(z, dtype=float))
 26    K = len(z)
 27    h = gamma * K
 28    m = int(math.floor(h))
 29    frac = h - m
 30    # Integral of the empirical quantile over [0,gamma], divided by gamma.
 31    total = z[:m].sum()
 32    if frac > 1e-12 and m < K:
 33        total += frac * z[m]
 34    return float(total / h)
 35
 36
 37def conformal_q(residuals, alpha):
 38    r = np.sort(np.asarray(residuals, dtype=float))
 39    k = int(math.ceil((len(r) + 1) * (1 - alpha)))
 40    if k > len(r):
 41        return float("inf")
 42    return float(r[k - 1])
 43
 44
 45def episodes(rng, n, K=16, gamma=0.25, shift=0.0):
 46    """Selected-trajectory episodes: latent true clearance and noisy biased samples."""
 47    # Heteroskedastic scenes make the calibration nontrivial while retaining exchangeability.
 48    y = rng.normal(0.72 + shift, 0.20, n)
 49    y = np.clip(y, -0.8, 1.5)
 50    bias = 0.12 + 0.12 * np.maximum(0, 0.5 - y)  # predictor overstates danger clearance
 51    # Candidate/stochastic predictions for the already-selected trajectory.
 52    z = y[:, None] + bias[:, None] + rng.normal(0, 0.13, (n, K))
 53    # Vectorized empirical lower-tail mean, equal to the eta maximization.
 54    zs = np.sort(z, axis=1)
 55    h = gamma * K
 56    mfull = int(math.floor(h))
 57    frac = h - mfull
 58    if mfull == 0:
 59        m = zs[:, 0].copy()
 60    else:
 61        m = zs[:, :mfull].sum(axis=1)
 62        if frac > 1e-12 and mfull < K:
 63            m = m + frac * zs[:, mfull]
 64        m = m / h
 65    nominal = z.mean(axis=1)
 66    return y, m, nominal
 67
 68
 69def one_trial(rng, ncal, ntest, alpha, gamma, shift_test=0.0):
 70    yc, mc, nc = episodes(rng, ncal, gamma=gamma)
 71    q = conformal_q(mc - yc, alpha)
 72    qt = conformal_q(nc - yc, alpha)
 73    yt, mt, nt = episodes(rng, ntest, gamma=gamma, shift=shift_test)
 74    cert = mt - q
 75    cert_nom = nt - qt
 76    # Uncalibrated is C=M (the common baseline).
 77    return {
 78        "q_cvar": q, "q_nominal": qt,
 79        "coverage_cvar": float(np.mean(yt >= cert)),
 80        "coverage_nominal_cal": float(np.mean(yt >= cert_nom)),
 81        "coverage_uncalibrated": float(np.mean(yt >= mt)),
 82        "nonnegative_cvar": float(np.mean(cert >= 0)),
 83        "nonnegative_nominal": float(np.mean(cert_nom >= 0)),
 84    }
 85
 86
 87def mean_trials(ncal, alpha, gamma, trials=30, ntest=700, shift=0.0):
 88    out = []
 89    for t in range(trials):
 90        out.append(one_trial(np.random.default_rng(SEED + 10000*t + ncal), ncal,
 91                             ntest, alpha, gamma, shift))
 92    keys = out[0].keys()
 93    return {k: float(np.mean([x[k] for x in out])) for k in keys}
 94
 95
 96def main():
 97    rng = np.random.default_rng(SEED)
 98    # Math sanity: CVaR optimizer versus direct lower-tail integral.
 99    cvar_rows = []
100    for K in [8, 16, 64, 127]:
101        z = rng.normal(size=K)
102        for gamma in [0.1, 0.25, 0.5, 1.0]:
103            a = empirical_cvar_lower(z, gamma)
104            b = empirical_bottom_cvar(z, gamma)
105            cvar_rows.append({"K": K, "gamma": gamma, "optimizer": a,
106                              "lower_tail_mean": b, "abs_error": abs(a-b)})
107
108    # Prediction 1: exact finite-sample conformal rank coverage is approximately
109    # 1-alpha (and never systematically below it) under exchangeability.
110    coverage_sweep = {}
111    for alpha in [0.05, 0.10, 0.20, 0.30]:
112        coverage_sweep[str(alpha)] = mean_trials(200, alpha, 0.25, trials=80)
113
114    # Prediction 2: increasing calibration n makes the quantile/certificate more
115    # stable; the expected rank granularity is 1/(n+1).
116    n_sweep = {}
117    for n in [20, 50, 100, 200, 500]:
118        n_sweep[str(n)] = mean_trials(n, 0.10, 0.25, trials=80)
119
120    # Prediction 3: the lower-tail mechanism is conservative relative to mean
121    # prediction when prediction noise is nonzero, increasing the required q.
122    gamma_sweep = {}
123    for gamma in [0.10, 0.25, 0.50, 1.0]:
124        gamma_sweep[str(gamma)] = mean_trials(200, 0.10, gamma, trials=80)
125
126    # Deliberate distribution shift: exchangeability assumption fails.
127    shifted = mean_trials(200, 0.10, 0.25, trials=30, shift=-0.45)
128    result = {
129        "seed": SEED,
130        "cvar_math_max_abs_error": max(x["abs_error"] for x in cvar_rows),
131        "cvar_math_rows": cvar_rows,
132        "coverage_sweep": coverage_sweep,
133        "calibration_size_sweep": n_sweep,
134        "gamma_sweep": gamma_sweep,
135        "shifted_test": shifted,
136        "predictions": {
137            "exchangeable_coverage": "coverage should be >= 1-alpha in aggregate; finite n gives rank granularity about 1/(n+1)",
138            "alpha_monotonicity": "larger alpha lowers q and target coverage toward 1-alpha",
139            "lower_tail_conservatism": "smaller gamma lowers M_CVaR and should require a smaller/equal residual correction than gamma=1"
140        }
141    }
142    Path("results.json").write_text(json.dumps(result, indent=2))
143    print(json.dumps({"cvar_max_error": result["cvar_math_max_abs_error"],
144                      "coverage": result["coverage_sweep"],
145                      "n": result["calibration_size_sweep"],
146                      "gamma": result["gamma_sweep"],
147                      "shifted": shifted}, indent=2))
148
149
150if __name__ == "__main__":
151    main()