scripts/claim3_cvar.py

"""Claim 3 -- Proposition 2.6: the family bounds and connects to CVaR, recovering CVaR
and LogSumExp as limiting cases:

  CVaR_rho(phi;mu) + lam(log rho - 1) <= lam F_rho(phi/lam; mu)
                                      <= CVaR_rho(phi;mu) + lam(log rho - 1 + 1/rho)

Independent tests:
  T1  CVaR cross-check: variational form (9) vs. tail-average/quantile definition.
  T2  the two-sided bound (10) over thousands of random (phi, mu, lam, rho).
  T3  tightness audit: shrink the upper additive constant lam/rho -> c*lam/rho and grow
      the lower one, and locate the constant at which each side starts to fail.
  T4  LIMIT 1 (CVaR): lam -> 0 with rho fixed; fitted convergence rate in lam.
  T5  LIMIT 2 (LogSumExp): rho -> 0 with lam fixed; lam*F_rho(phi/lam) -> lam*log E e^{phi/lam},
      i.e. the LogSumExp/log-partition functional; fitted rate in rho.
  T6  end-to-end interpolation table: one instance, grid over (lam, rho), reporting the
      distance to CVaR_rho and to the LogSumExp value.
  T7  boundary audit: rho -> 1 (CVaR_1 = mean; lam F_1 = mean - lam).
"""

import json
import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from common import F_rho, F_true, cvar

OUT = os.path.join(
    os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "outputs"
)
os.makedirs(OUT, exist_ok=True)
SEED = 20260725
rng = np.random.default_rng(SEED)
res = {"seed": SEED}


def rand_phi(rng):
    kind = rng.choice(["gauss", "heavy", "uniform", "bimodal", "spiky"])
    n = int(rng.integers(20, 500))
    if kind == "gauss":
        phi = rng.normal(0, rng.uniform(0.3, 3), n)
    elif kind == "heavy":
        phi = rng.standard_t(3, n) * rng.uniform(0.5, 2)
    elif kind == "uniform":
        phi = rng.uniform(-5, 5, n)
    elif kind == "bimodal":
        phi = np.concatenate(
            [rng.normal(-2, 0.4, n // 2), rng.normal(2, 0.4, n - n // 2)]
        )
    else:
        phi = rng.normal(0, 0.3, n)
        phi[0] += rng.uniform(2, 8)
    w = rng.dirichlet(np.ones(n))
    return phi, w, str(kind)


# ------------------------------------------------------------------ T1 CVaR sanity
t1 = {"max_abs_err_vs_tail_average": 0.0, "n": 0}
for _ in range(300):
    n = int(rng.integers(50, 400))
    phi = rng.normal(0, 2, n)
    w = np.full(n, 1.0 / n)
    for rho in [0.05, 0.1, 0.25, 0.5]:
        k = int(round(rho * n))
        if k < 1:
            continue
        if abs(k / n - rho) > 1e-12:
            continue  # only exact-quantile cases are directly comparable
        tail = np.mean(np.sort(phi)[-k:])
        v = cvar(phi, rho, w)
        t1["max_abs_err_vs_tail_average"] = max(
            t1["max_abs_err_vs_tail_average"], abs(tail - v)
        )
        t1["n"] += 1
res["T1_cvar_crosscheck"] = t1

# ------------------------------------------------------------------ T2 the bound
t2 = {
    "n": 0,
    "lower_violations": 0,
    "upper_violations": 0,
    "min_lower_slack": np.inf,
    "min_upper_slack": np.inf,
    "max_rel_position": 0.0,
    "min_rel_position": np.inf,
    "examples": [],
}
for _ in range(4000):
    phi, w, kind = rand_phi(rng)
    lam = float(np.exp(rng.uniform(np.log(1e-3), np.log(30))))
    rho = float(rng.choice([1e-3, 1e-2, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9]))
    C = cvar(phi, rho, w)
    lF = lam * F_rho(phi / lam, rho, w)[0]
    lo = C + lam * (np.log(rho) - 1)
    hi = C + lam * (np.log(rho) - 1 + 1 / rho)
    t2["n"] += 1
    if lF < lo - 1e-9 * max(1, abs(lo)):
        t2["lower_violations"] += 1
    if lF > hi + 1e-9 * max(1, abs(hi)):
        t2["upper_violations"] += 1
    t2["min_lower_slack"] = min(t2["min_lower_slack"], float(lF - lo))
    t2["min_upper_slack"] = min(t2["min_upper_slack"], float(hi - lF))
    pos = (lF - lo) / (hi - lo)
    t2["max_rel_position"] = max(t2["max_rel_position"], float(pos))
    t2["min_rel_position"] = min(t2["min_rel_position"], float(pos))
    if len(t2["examples"]) < 4:
        t2["examples"].append(
            {
                "kind": kind,
                "n": int(phi.size),
                "lam": lam,
                "rho": rho,
                "CVaR": C,
                "lam_F_rho": lF,
                "lower": float(lo),
                "upper": float(hi),
            }
        )
t2["min_lower_slack"] = float(t2["min_lower_slack"])
t2["min_upper_slack"] = float(t2["min_upper_slack"])
res["T2_proposition26_bounds"] = t2

# ------------------------------------------------------------------ T3 tightness audit
t3 = {}
for c in [1.0, 0.9, 0.7, 0.5, 0.25]:
    viol = 0
    n = 0
    rng2 = np.random.default_rng(7 + int(100 * c))
    for _ in range(1500):
        phi, w, _ = rand_phi(rng2)
        lam = float(np.exp(rng2.uniform(np.log(1e-2), np.log(20))))
        rho = float(rng2.choice([1e-2, 0.05, 0.1, 0.25, 0.5, 0.9]))
        C = cvar(phi, rho, w)
        lF = lam * F_rho(phi / lam, rho, w)[0]
        hi = C + lam * (np.log(rho) - 1 + c / rho)
        n += 1
        if lF > hi + 1e-9:
            viol += 1
    t3[f"upper_const_c_times_lam_over_rho__c={c}"] = {"violations": viol, "n": n}
for d in [0.0, 0.05, 0.2, 0.5]:
    viol = 0
    n = 0
    rng2 = np.random.default_rng(99 + int(100 * d))
    for _ in range(1500):
        phi, w, _ = rand_phi(rng2)
        lam = float(np.exp(rng2.uniform(np.log(1e-2), np.log(20))))
        rho = float(rng2.choice([1e-2, 0.05, 0.1, 0.25, 0.5, 0.9]))
        C = cvar(phi, rho, w)
        lF = lam * F_rho(phi / lam, rho, w)[0]
        lo = C + lam * (np.log(rho) - 1) + d * lam
        n += 1
        if lF < lo - 1e-9:
            viol += 1
    t3[f"lower_shifted_by_d_times_lam__d={d}"] = {"violations": viol, "n": n}
res["T3_tightness_audit"] = t3

# ------------------------------------------------------------------ T4 CVaR limit
lams = np.array([1.0, 0.3, 0.1, 0.03, 1e-2, 3e-3, 1e-3, 3e-4, 1e-4])
t4 = {"instances": [], "rate_exponent_mean": None}
exps = []
for inst in range(60):
    phi, w, kind = rand_phi(rng)
    rho = float(rng.choice([0.05, 0.1, 0.25, 0.5]))
    C = cvar(phi, rho, w)
    errs = []
    for lam in lams:
        val = lam * F_rho(phi / lam, rho, w)[0] - lam * (np.log(rho) - 1)
        errs.append(abs(val - C))
    errs = np.array(errs)
    p = np.polyfit(np.log(lams[-5:]), np.log(np.maximum(errs[-5:], 1e-300)), 1)
    exps.append(float(p[0]))
    if inst < 3:
        t4["instances"].append(
            {
                "kind": kind,
                "rho": rho,
                "CVaR": C,
                "lam_grid": lams.tolist(),
                "abs_err": errs.tolist(),
            }
        )
t4["rate_exponent_mean"] = float(np.mean(exps))
t4["rate_exponent_std"] = float(np.std(exps))
t4["rate_exponent_min"] = float(np.min(exps))
t4["note"] = (
    "lam*F_rho(phi/lam) - lam(log rho - 1) -> CVaR_rho ; exponent of err ~ lam^p"
)
res["T4_limit_CVaR"] = t4

# ------------------------------------------------------------------ T5 LogSumExp limit
rhos = np.array([0.5, 0.1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6])
t5 = {
    "instances": [],
}
exps = []
for inst in range(60):
    phi, w, kind = rand_phi(rng)
    lam = float(rng.choice([0.5, 1.0, 2.0]))
    target = lam * F_true(phi / lam, w)  # = lam log int e^{phi/lam} dmu  (LogSumExp)
    errs = []
    for rho in rhos:
        errs.append(abs(lam * F_rho(phi / lam, rho, w)[0] - target))
    errs = np.array(errs)
    p = np.polyfit(np.log(rhos[-4:]), np.log(np.maximum(errs[-4:], 1e-300)), 1)
    exps.append(float(p[0]))
    if inst < 3:
        t5["instances"].append(
            {
                "kind": kind,
                "lam": lam,
                "logsumexp_value": float(target),
                "rho_grid": rhos.tolist(),
                "abs_err": errs.tolist(),
            }
        )
t5["rate_exponent_mean"] = float(np.mean(exps))
t5["rate_exponent_std"] = float(np.std(exps))
t5["note"] = (
    "F_rho -> F (log-partition / LogSumExp) as rho -> 0 ; exponent of err ~ rho^p"
)
res["T5_limit_LogSumExp"] = t5

# ------------------------------------------------------------------ T6 interpolation table
rng3 = np.random.default_rng(2026)
phi = rng3.normal(0, 1.5, 300)
w = np.full(300, 1 / 300)
table = []
for lam in [3.0, 1.0, 0.3, 0.1, 0.03, 0.01]:
    for rho in [0.5, 0.1, 0.01, 1e-3]:
        C = cvar(phi, rho, w)
        lse = lam * F_true(phi / lam, w)
        lF = lam * F_rho(phi / lam, rho, w)[0]
        table.append(
            {
                "lam": lam,
                "rho": rho,
                "lam_F_rho": lF,
                "CVaR_rho_shifted": float(C + lam * (np.log(rho) - 1)),
                "logsumexp_lam": float(lse),
                "dist_to_CVaR_shifted": float(abs(lF - (C + lam * (np.log(rho) - 1)))),
                "dist_to_logsumexp": float(abs(lF - lse)),
            }
        )
res["T6_interpolation_table"] = {
    "mean_phi": float(np.sum(w * phi)),
    "n": 300,
    "rows": table,
}

# ------------------------------------------------------------------ T7 boundary rho->1
t7 = []
phi = rng3.normal(0, 1.0, 400)
w = np.full(400, 1 / 400)
for rho in [0.9, 0.99, 0.999]:
    for lam in [1.0, 0.1]:
        t7.append(
            {
                "rho": rho,
                "lam": lam,
                "lam_F_rho": float(lam * F_rho(phi / lam, rho, w)[0]),
                "mean_minus_lam": float(np.sum(w * phi) - lam),
                "CVaR_rho": float(cvar(phi, rho, w)),
                "mean": float(np.sum(w * phi)),
            }
        )
res["T7_boundary_rho_to_1"] = t7

with open(os.path.join(OUT, "claim3_cvar.json"), "w") as fh:
    json.dump(res, fh, indent=2, default=float)
print(
    json.dumps(
        {k: v for k, v in res.items() if k != "T6_interpolation_table"},
        indent=2,
        default=float,
    )[:5000]
)