1import json
  2import math
  3import numpy as np
  4
  5# Fisher-floor-corrected DSM: toy verification and mini validation experiment.
  6# All quantities are evaluated for a 1D Gaussian-mixture data prior, whose
  7# posterior and marginal score are available exactly by finite mixture sums.
  8
  9
 10def mixture_params():
 11    # Unequal, separated components make posterior uncertainty nontrivial.
 12    return np.array([-2.0, 0.5, 2.5]), np.array([0.25, 0.45, 0.30]), np.array([0.38, 0.55, 0.42])
 13
 14
 15def sample_y(n, rng):
 16    means, probs, _ = mixture_params()
 17    z = rng.choice(len(means), size=n, p=probs)
 18    return means[z] + rng.normal(size=n) * np.array([0.35, 0.5, 0.4])[z]
 19
 20
 21def posterior_quantities(x, alpha, sigma):
 22    """Exact posterior mean/variance and marginal score for mixture prior."""
 23    means, probs, scales = mixture_params()
 24    # x | component k is N(alpha*mean_k, sigma^2 + alpha^2*scale_k^2)
 25    varx = sigma * sigma + (alpha * scales) ** 2
 26    logp = np.log(probs) - 0.5 * (np.log(2*np.pi*varx) + (x[:, None]-alpha*means)**2/varx)
 27    logp -= logp.max(axis=1, keepdims=True)
 28    w = np.exp(logp); w /= w.sum(axis=1, keepdims=True)
 29    post_var_k = (scales**2 * sigma**2) / varx
 30    post_mean_k = (means * sigma**2 + alpha * scales**2 * x[:, None]) / varx
 31    mu = (w * post_mean_k).sum(axis=1)
 32    second = (w * (post_var_k + post_mean_k**2)).sum(axis=1)
 33    vy = second - mu**2
 34    # Conditional-score mean equals marginal score (Tweedie identity).
 35    score = (alpha * mu - x) / (sigma*sigma)
 36    floor = (alpha*alpha / sigma**4) * vy
 37    return mu, vy, score, floor
 38
 39
 40def bank_floor(x, alpha, sigma, bank):
 41    logits = -(x[:, None] - alpha*bank[None, :])**2 / (2*sigma*sigma)
 42    logits -= logits.max(axis=1, keepdims=True)
 43    p = np.exp(logits); p /= p.sum(axis=1, keepdims=True)
 44    mu = p @ bank
 45    v = (p * (bank[None, :] - mu[:, None])**2).sum(axis=1)
 46    return alpha*alpha / sigma**4 * v
 47
 48
 49def decomposition_check(rng):
 50    n = 180000; alpha, sigma = 0.8, 0.9
 51    y = sample_y(n, rng); x = alpha*y + sigma*rng.normal(size=n)
 52    target = (alpha*y-x)/sigma**2
 53    _, _, marginal, floor = posterior_quantities(x, alpha, sigma)
 54    # Deliberately imperfect score, so equality is not only checked at optimum.
 55    model = marginal + 0.35*np.sin(1.7*x) + 0.12*rng.normal(size=n)
 56    raw = np.mean((model-target)**2)
 57    ideal = np.mean((model-marginal)**2)
 58    f = np.mean(floor)
 59    return {"raw": float(raw), "ideal": float(ideal), "floor": float(f),
 60            "raw_minus_ideal": float(raw-ideal), "abs_identity_error": float(abs(raw-ideal-f)),
 61            "relative_identity_error": float(abs(raw-ideal-f)/raw)}
 62
 63
 64def alpha_scaling(rng):
 65    # At high noise, posterior is close to prior; F ~= Var(Y)*alpha^2/sigma^4.
 66    sigma = 5.0; n = 220000; y = sample_y(n, rng); eps = rng.normal(size=n)
 67    vals = []
 68    for a in np.array([0.10, 0.20, 0.35, 0.50, 0.70]):
 69        x = a*y + sigma*eps
 70        vals.append(np.mean(posterior_quantities(x, a, sigma)[3]))
 71    aa = np.array([.10,.20,.35,.50,.70])
 72    slope = float(np.sum(aa**2*np.array(vals))/np.sum(aa**4))
 73    predicted = slope*aa**2
 74    rel = np.max(np.abs(np.array(vals)-predicted)/np.maximum(np.array(vals),1e-15))
 75    return {"alphas": aa.tolist(), "floors": np.array(vals).tolist(),
 76            "quadratic_fit_relative_max_error": float(rel),
 77            "observed_ratio_F(a=.7)/F(a=.1)": float(vals[-1]/vals[0]),
 78            "predicted_ratio_(.7/.1)^2": 49.0}
 79
 80
 81def alpha_zero_and_noise(rng):
 82    n=120000; y=sample_y(n,rng)
 83    rows=[]
 84    for a,s in [(0.0,0.8),(0.8,0.8),(0.8,5.0),(0.8,10.0)]:
 85        x=a*y+s*rng.normal(size=n)
 86        rows.append((a,s,float(np.mean(posterior_quantities(x,a,s)[3]))))
 87    return {"cases":[{"alpha":a,"sigma":s,"floor":f} for a,s,f in rows],
 88            "alpha_zero_floor": rows[0][2],
 89            "zero_prediction": 0.0}
 90
 91
 92def bank_convergence(rng):
 93    # Prediction: iid bank Monte Carlo error decreases approximately B^-1/2.
 94    n=9000; alpha=.8; sigma=.9
 95    y=sample_y(n,rng); x=alpha*y+sigma*rng.normal(size=n)
 96    exact=posterior_quantities(x,alpha,sigma)[3]
 97    out=[]
 98    for b in [16,32,64,128,256,512,1024]:
 99        errs=[]
100        for _ in range(5):
101            bank=sample_y(b,rng)
102            errs.append(np.mean(np.abs(bank_floor(x,alpha,sigma,bank)-exact)))
103        out.append((b,float(np.mean(errs))))
104    bs=np.array([z[0] for z in out],float); es=np.array([z[1] for z in out])
105    slope=float(np.polyfit(np.log(bs),np.log(es),1)[0])
106    return {"mean_absolute_errors":[{"bank":b,"mae":e} for b,e in out],
107            "loglog_error_slope":slope,"predicted_slope":-0.5}
108
109
110def ranking_experiment(rng):
111    # Same checkpoints (score perturbations), two schedules with different
112    # additive DSM floors. Corrected estimates target marginal-score error.
113    n=100000; y=sample_y(n,rng)
114    checkpoints=[]
115    xref=0.8*y+0.9*rng.normal(size=n)
116    _,_,sref,_=posterior_quantities(xref,.8,.9)
117    for j,noise in enumerate([.04,.10,.18,.28,.40,.60]):
118        # perturbation grows with x and independent score noise
119        pred=sref + noise*(0.7*np.tanh(xref)+rng.normal(size=n))
120        checkpoints.append(np.mean((pred-sref)**2))
121    schedules=[(.35,1.5),(.9,.65)] # (alpha,sigma), deliberately different floors
122    rows=[]
123    for a,s in schedules:
124        x=a*y+s*rng.normal(size=n)
125        _,_,m,f=posterior_quantities(x,a,s)
126        for j,noise in enumerate([.04,.10,.18,.28,.40,.60]):
127            pred=m + noise*(0.7*np.tanh(x)+rng.normal(size=n))
128            raw=float(np.mean((pred-(a*y-x)/s**2)**2))
129            corr=raw-float(np.mean(f))
130            ideal=float(np.mean((pred-m)**2))
131            rows.append({"schedule":[a,s],"checkpoint":j,"raw":raw,"corrected":corr,"ideal":ideal,"floor":float(np.mean(f))})
132    def rank_corr(vals, ideal):
133        return float(np.corrcoef(np.argsort(np.argsort(vals)),np.argsort(np.argsort(ideal)))[0,1])
134    summary=[]
135    for sch in schedules:
136        rr=[r for r in rows if r['schedule']==list(sch)]
137        summary.append({"schedule":sch,"floor":rr[0]['floor'],
138                        "raw_ideal_rank_corr":rank_corr([r['raw'] for r in rr],[r['ideal'] for r in rr]),
139                        "corrected_ideal_rank_corr":rank_corr([r['corrected'] for r in rr],[r['ideal'] for r in rr])})
140    return {"schedule_summary":summary,"rows":rows}
141
142
143def main():
144    rng=np.random.default_rng(2713)
145    result={"decomposition":decomposition_check(rng),"alpha_scaling":alpha_scaling(rng),
146            "alpha_zero_and_noise":alpha_zero_and_noise(rng),"bank_convergence":bank_convergence(rng),
147            "ranking":ranking_experiment(rng)}
148    with open("results.json","w") as f: json.dump(result,f,indent=2)
149    print(json.dumps(result,indent=2))
150
151if __name__ == '__main__': main()