files / problems / union-closed / compute / best_on_ray.py
"""Exact envelope of the {b,1} family for iid + Example-4 (and + Example-5).
On support {b,1} with P(S=1)=a,
ratio = (1-a) * Hmix(b,β) / h(b)
where Hmix = (1-β) h(2b-b²) + β h_or_proto(b,b) (independent C3).
Equality forces a = 1 - h(b)/Hmix and
mean = 1 - (1-b) h(b) / Hmix(b,β).
We maximise that mean over (b,β), then check the mean-preserving
derivative (Liu's stationarity) and a fine verification grid.
"""
from __future__ import annotations
import json
import math
from pathlib import Path
import numpy as np
from entropy import a_example4, h, h_or_example4, h_or_example5, h_or_indep, pi_example4
def Hmix_ex4(b, beta):
return (1.0 - beta) * h_or_indep(b, b) + beta * h_or_example4(b, b)
def Hmix_ex5(b, beta):
return (1.0 - beta) * h_or_indep(b, b) + beta * h_or_example5(b, b)
def envelope(b, beta, proto="ex4"):
Hm = Hmix_ex4(b, beta) if proto == "ex4" else Hmix_ex5(b, beta)
hb = h(b)
if Hm <= hb * 0.5 or Hm <= 1e-18 or hb <= 1e-18:
return None
one_minus_a = hb / Hm
a = 1.0 - one_minus_a
if a < -1e-12 or a > 0.5 + 1e-12:
return None
a = min(max(a, 0.0), 0.5)
mean = 1.0 - (1.0 - b) * one_minus_a
return dict(b=b, beta=beta, a=a, mean=mean, Hmix=Hm, hb=hb)
def grid_envelope(proto="ex4", nb=2000, nbeta=800):
best = None
# also record best for each beta
by_beta = []
for j in range(nbeta + 1):
beta = j / nbeta
loc = None
for i in range(nb + 1):
b = 0.05 + 0.44 * i / nb
e = envelope(b, beta, proto)
if e is None:
continue
if loc is None or e["mean"] > loc["mean"]:
loc = e
if best is None or e["mean"] > best["mean"]:
best = e
if loc:
by_beta.append(loc)
return best, by_beta
def verify_box(proto, beta, c, nb=2500, na=2000):
"""Min ratio on {b,1} among mean ≤ c."""
fn = h_or_example4 if proto == "ex4" else h_or_example5
min_r = 10.0
arg = None
n_checked = 0
for i in range(nb + 1):
b = 0.02 + 0.48 * i / nb
hb = h(b)
hid = h_or_indep(b, b)
hp = fn(b, b)
Hm = (1.0 - beta) * hid + beta * hp
if hb <= 1e-18:
continue
for j in range(na + 1):
a = 0.5 * j / na
mean = a + (1.0 - a) * b
if mean > c + 1e-15:
continue
eh = (1.0 - a) * hb
if eh <= 1e-18:
continue
# independent C3 for proto; also correlated
e_ind = (1.0 - a) ** 2 * hp
e_cor = (1.0 - a) * hp
ep = min(e_ind, e_cor)
eiid = (1.0 - a) ** 2 * hid
r = ((1.0 - beta) * eiid + beta * ep) / eh
n_checked += 1
if r < min_r:
min_r = r
arg = dict(a=a, b=b, mean=mean, ratio=r)
return dict(min_ratio=min_r, arg=arg, n_checked=n_checked)
def stationarity_beta(b, proto="ex4"):
"""Liu-style β that kills the mean-preserving derivative at the
equality 2-point for this b (a determined by iid-equality or mix-equality).
We use the envelope a(b,β) and pick β so d mean / db = 0 along equality,
equivalently maximise mean(b,β) in b for that β — already done by the grid.
Here we just compute a finite-difference d(ratio)/db along d(mean)=0
at the envelope point, to confirm it's ~0 at the reported maximiser.
"""
return None
def main():
report = {}
for proto in ("ex5", "ex4"):
print("envelope", proto, flush=True)
best, by_beta = grid_envelope(proto, nb=2500, nbeta=1000)
print(" best", best, flush=True)
# verify a few c's just below best.mean
ver = {}
c_best = best["mean"]
beta = best["beta"]
for dc, label in [(0.0, "at_best"), (2e-6, "minus_2e-6"), (1e-5, "minus_1e-5"),
(c_best - 0.38270908792, "liu_c")]:
c = c_best - dc
print(" verify", label, c, "beta", beta, flush=True)
ver[label] = {"c": c, **verify_box(proto, beta, c, nb=1800, na=1400)}
print(" min_r", ver[label]["min_ratio"], ver[label]["arg"], flush=True)
# also verify Liu's published β for ex5 at Liu's c
report[proto] = {
"best_envelope": best,
"verify": ver,
"beta_curve_sample": by_beta[::50],
}
# Liu published point replay
print("liu published replay", flush=True)
report["liu_published_replay"] = verify_box(
"ex5", 0.100052559862974, 0.382709087918735, nb=2000, na=1600
)
print(report["liu_published_replay"], flush=True)
path = Path(__file__).resolve().parent / "best_on_ray.json"
path.write_text(json.dumps(report, indent=2) + "\n")
print("wrote", path)
if __name__ == "__main__":
main()