1import json
  2import numpy as np
  3from scipy.special import ndtr
  4
  5SEED = 2816
  6BETA = 0.99
  7
  8
  9def cvar_rockafellar(y, beta=BETA):
 10    y = np.asarray(y)
 11    eta = np.quantile(y, beta, method='linear')
 12    return float(eta + np.mean(np.maximum(y - eta, 0.0)) / (1-beta))
 13
 14
 15def ridge_fit(P, z, lam):
 16    P = np.asarray(P); z = np.asarray(z)
 17    return np.linalg.solve(P.T @ P + lam*np.eye(P.shape[1]), P.T @ z)
 18
 19
 20def phi(x):
 21    x = np.asarray(x)
 22    bump = np.exp(-0.5*((x-0.90)/0.025)**2)
 23    return np.column_stack([np.ones_like(x), x, bump, np.sin(6*np.pi*x)])
 24
 25
 26def low_members(x):
 27    # Five cheap models: close globally, but uncertain near the tail boundary.
 28    x = np.asarray(x)
 29    base = 0.15 + 0.35*x + 0.08*np.sin(4*np.pi*x)
 30    bump = np.exp(-0.5*((x-0.90)/0.025)**2)
 31    offsets = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
 32    return np.array([base + 0.025*o*(0.25 + 2.0*bump) for o in offsets])
 33
 34
 35def high_oracle(x):
 36    x = np.asarray(x)
 37    # Narrow, high-impact failure region missing from the cheap model.
 38    return 0.15 + 0.35*x + 0.08*np.sin(4*np.pi*x) + 1.25*np.exp(-0.5*((x-0.90)/0.025)**2) + 0.04*np.sin(18*np.pi*x)
 39
 40
 41def residual_uncertainty(x, Xh, zh, lam):
 42    if len(Xh) < 3: return np.ones(len(x))
 43    P = phi(Xh); w = ridge_fit(P, zh, lam)
 44    pred = P @ w
 45    rng = np.random.default_rng(123)
 46    ws = []
 47    for _ in range(20):
 48        idx = rng.integers(0, len(Xh), len(Xh))
 49        ws.append(ridge_fit(P[idx], zh[idx], lam))
 50    return np.std(np.asarray(ws) @ phi(x).T, axis=0)
 51
 52
 53def acquire(policy, cand, Xh, zh, n, eps=1.0, rho=0.5, lam=1e-3):
 54    members = low_members(cand); mu = members.mean(0); s = members.std(0, ddof=1)
 55    q = np.quantile(mu, BETA)
 56    region = np.abs(mu-q) <= eps*np.maximum(s, 1e-5)
 57    if policy == 'random':
 58        score = np.random.default_rng(999+len(Xh)).random(len(cand))
 59        return cand[np.argsort(score)[-n:]]
 60    if policy == 'global_uncertainty':
 61        score = s
 62    else:
 63        # stage 1: half from tail boundary by uncertainty; stage 2: residual/correction impact
 64        ids = np.where(region)[0]
 65        if len(ids) == 0: ids = np.arange(len(cand))
 66        n1 = n//2
 67        take1 = ids[np.argsort(s[ids])[-min(n1,len(ids)):]]
 68        sr = residual_uncertainty(cand, Xh, zh, lam)
 69        if len(Xh): w = ridge_fit(phi(Xh), zh, lam); impact = np.abs(phi(cand) @ w)
 70        else: impact = np.zeros(len(cand))
 71        score = (sr/(sr.std()+1e-9) + rho*impact/(impact.std()+1e-9))
 72        available = np.ones(len(cand), dtype=bool); available[take1] = False
 73        ids2 = np.where(available)[0]
 74        take2 = ids2[np.argsort(score[ids2])[-(n-n1):]]
 75        return np.concatenate([cand[take1], cand[take2]])
 76    return cand[np.argsort(score)[-n:]]
 77
 78
 79def mechanism_checks():
 80    # Prediction 1: Rockafellar minimizer is the beta quantile.
 81    y = np.sort(np.random.default_rng(SEED).normal(size=10001))
 82    etas = np.linspace(y[int(.985*len(y))], y[int(.995*len(y))], 401)
 83    obj = [e + np.mean(np.maximum(y-e,0))/(1-BETA) for e in etas]
 84    eta_obs = etas[np.argmin(obj)]; q = np.quantile(y, BETA)
 85    cvar_direct = y[int(np.ceil(BETA*len(y))):].mean()
 86    # Prediction 2: in a single orthogonal direction, ridge correction ratio is 1/(1+lambda).
 87    z = np.array([2.0]); P = np.array([[1.0]])
 88    lams = np.array([0., .01, .1, 1., 10.])
 89    ratios = np.array([ridge_fit(P,z,l)[0]/2.0 for l in lams])
 90    pred_ratios = 1/(1+lams)
 91    # Prediction 3: tail-boundary region probability for standard normal standardized distance is 2Phi(eps)-1.
 92    # Generate standardized distances directly, matching the definition |mu-q|/s.
 93    d = np.random.default_rng(SEED+1).normal(size=200000)
 94    eps = np.array([0.5, 1., 2.])
 95    observed = np.array([(np.abs(d)<=e).mean() for e in eps])
 96    predicted = 2*ndtr(eps)-1
 97    return {
 98        'cvar_quantile_check': {'predicted_eta_quantile':float(q), 'observed_eta':float(eta_obs), 'abs_error':float(abs(q-eta_obs)), 'direct_cvar':float(cvar_direct), 'rockafellar_at_eta':float(min(obj))},
 99        'ridge_shrinkage_check': {'lambda':lams.tolist(), 'predicted_ratio':pred_ratios.tolist(), 'observed_ratio':ratios.tolist(), 'max_abs_error':float(np.max(abs(ratios-pred_ratios)))},
100        'tail_region_check': {'epsilon':eps.tolist(), 'predicted_fraction':predicted.tolist(), 'observed_fraction':observed.tolist(), 'max_abs_error':float(np.max(abs(observed-predicted)))}
101    }
102
103
104def mini_experiment():
105    rng = np.random.default_rng(SEED)
106    pool = np.linspace(0,1,5000,endpoint=False)+0.5/5000
107    truth = high_oracle(pool); true_cvar = cvar_rockafellar(truth)
108    results = {}
109    for policy in ['random','global_uncertainty','tail_residual']:
110        # same initial design across policies
111        init = rng.choice(len(pool), 24, replace=False)
112        Xh = pool[init].copy(); zh = high_oracle(Xh)-low_members(Xh).mean(0)
113        rows=[]
114        for rnd in range(5):
115            # corrected predictor uses low mean plus fitted residual
116            w = ridge_fit(phi(Xh), zh, 1e-3)
117            pred = low_members(pool).mean(0) + phi(pool)@w
118            rows.append({'calls':int(len(Xh)), 'cvar_abs_error':float(abs(cvar_rockafellar(pred)-true_cvar)), 'global_mse':float(np.mean((pred-truth)**2)), 'tail_mean_error':float(abs(np.mean(np.sort(pred)[-50:])-np.mean(np.sort(truth)[-50:])))})
119            cand = pool[rng.choice(len(pool), 2500, replace=False)]
120            new = acquire(policy, cand, Xh, zh, 16)
121            Xh = np.concatenate([Xh,new]); zh = np.concatenate([zh, high_oracle(new)-low_members(new).mean(0)])
122        results[policy] = rows
123    return {'true_cvar':true_cvar, 'runs':results}
124
125
126if __name__ == '__main__':
127    out = {'seed':SEED, 'beta':BETA, 'mechanism_checks':mechanism_checks(), 'mini_experiment':mini_experiment()}
128    with open('results.json','w') as f: json.dump(out,f,indent=2)
129    print(json.dumps(out, indent=2))