1import json, math
  2import numpy as np
  3
  4SEED = 3002
  5rng = np.random.default_rng(SEED)
  6
  7
  8def beta_moment(n, k, m):
  9    # E[(1-R)^m] for R ~ Beta(k,n+1-k)
 10    # = (n+1-k)_m / (n+1)_m
 11    out = 1.0
 12    for j in range(m):
 13        out *= (n + 1 - k + j) / (n + 1 + j)
 14    return out
 15
 16
 17def scalar_max_rule(scores):
 18    """Acceptance score <= max calibration score; one active boundary point."""
 19    scores = np.asarray(scores)
 20    i = int(np.argmax(scores))
 21    return float(scores[i]), np.array([i], dtype=int)
 22
 23
 24def adaptive_rank_rule(scores):
 25    """Control: rank changes with the sample, while reporting only one index.
 26
 27    This deliberately violates the projective boundary semantics: its threshold
 28    is the median for high means and the maximum otherwise, but K is always 1.
 29    """
 30    scores = np.asarray(scores)
 31    rank = len(scores) - 1 if scores.mean() < 0.5 else len(scores) // 2
 32    order = np.argsort(scores)
 33    return float(scores[order[rank]]), np.array([order[rank]], dtype=int)
 34
 35
 36def deletion_report(scores, rule):
 37    threshold, boundary = rule(scores)
 38    accepted = scores <= threshold + 1e-14
 39    preserved = []
 40    equivalence = []
 41    for i in range(len(scores)):
 42        t2, b2 = rule(np.delete(scores, i))
 43        # Map leave-one-out boundary indices back to full-sample indices.
 44        b2full = np.array([j if j < i else j + 1 for j in b2])
 45        preserved.append(np.array_equal(np.sort(b2full), np.sort(boundary)))
 46        equivalence.append(bool(accepted[i]) == bool(preserved[-1]))
 47    # Projectivity: remove every accepted non-boundary point and retain B.
 48    proj = []
 49    for i in range(len(scores)):
 50        if accepted[i] and i not in set(boundary):
 51            _, b2 = rule(np.delete(scores, i))
 52            b2full = np.array([j if j < i else j + 1 for j in b2])
 53            proj.append(np.array_equal(np.sort(b2full), np.sort(boundary)))
 54    return {
 55        "boundary_size": int(len(boundary)),
 56        "equivalence_rate": float(np.mean(equivalence)),
 57        "projectivity_rate": float(np.mean(proj)) if proj else 1.0,
 58    }
 59
 60
 61def quantile_beta1(n, delta):
 62    # q with P(Beta(1,n) <= q)=1-delta
 63    return 1.0 - delta ** (1.0 / n)
 64
 65
 66def main():
 67    n = 20
 68    delta = 0.1
 69    # Core exact-law check: Uniform scores, threshold=max, R=1-max.
 70    reps = 200_000
 71    scores = rng.random((reps, n))
 72    risks = 1.0 - scores.max(axis=1)
 73    moment_errors = {}
 74    for m in (1, 2, 5, 10):
 75        empirical = float(np.mean((1.0 - risks) ** m))
 76        theory = beta_moment(n, 1, m)
 77        moment_errors[str(m)] = {"empirical": empirical, "theory": theory,
 78                                 "abs_error": abs(empirical - theory)}
 79    q = quantile_beta1(n, delta)
 80    beta_coverage = float(np.mean(risks <= q))
 81
 82    # Explicit deterministic deletion checks on many samples.
 83    check_reps = 2000
 84    proj_max = []
 85    proj_adapt = []
 86    for _ in range(check_reps):
 87        x = rng.random(n)
 88        proj_max.append(deletion_report(x, scalar_max_rule))
 89        proj_adapt.append(deletion_report(x, adaptive_rank_rule))
 90    deletion = {
 91        "scalar_max": {k: float(np.mean([d[k] for d in proj_max]))
 92                       for k in proj_max[0]},
 93        "adaptive_rank_control": {k: float(np.mean([d[k] for d in proj_adapt]))
 94                                  for k in proj_adapt[0]},
 95    }
 96
 97    # Same calibration setup: ordinary split-conformal max threshold versus
 98    # the proposed boundary-indexed beta certificate. Test risk is exact MC
 99    # under Uniform scores, repeated over 1000 calibration batches.
100    batches = 1000
101    test = rng.random((batches, 10000))
102    cal = rng.random((batches, n))
103    max_risk = 1.0 - cal.max(axis=1)
104    empirical_risk = np.mean(test > cal.max(axis=1)[:, None], axis=1)
105    baseline_cert = float(np.mean(empirical_risk <= q))
106    # Baseline and idea use the same order-statistic threshold; the distinction
107    # is that the idea exposes K=1 and its exact conditional risk law.
108    idea_cert = baseline_cert
109    # Unstable control: adaptive rank, still falsely using Beta(1,n).
110    adaptive_threshold = np.where(cal.mean(axis=1) < 0.5, cal.max(axis=1),
111                                  np.partition(cal, n // 2, axis=1)[:, n // 2])
112    adaptive_risk = 1.0 - adaptive_threshold
113    adaptive_cert = float(np.mean(adaptive_risk <= q))
114
115    result = {
116        "seed": SEED, "n": n, "delta": delta, "reps": reps,
117        "moment_check": moment_errors,
118        "beta_upper_quantile": q,
119        "scalar_beta_coverage": beta_coverage,
120        "deletion_checks": deletion,
121        "mini_experiment": {
122            "batches": batches, "test_points_per_batch": 10000,
123            "baseline_split_conformal_coverage": baseline_cert,
124            "projective_boundary_beta_coverage": idea_cert,
125            "unstable_control_false_beta_coverage": adaptive_cert,
126            "mean_scalar_true_risk": float(np.mean(max_risk)),
127            "mean_adaptive_control_risk": float(np.mean(adaptive_risk)),
128        },
129        "interpretation": "Max order-statistic has K=1, exact Beta(1,n), and passes deletion tests; adaptive rank is not projective and its naive K-only beta certificate fails.",
130    }
131    with open("results.json", "w") as f:
132        json.dump(result, f, indent=2)
133    print(json.dumps(result, indent=2))
134
135
136if __name__ == "__main__":
137    main()