import json import math import numpy as np # Fisher-floor-corrected DSM: toy verification and mini validation experiment. # All quantities are evaluated for a 1D Gaussian-mixture data prior, whose # posterior and marginal score are available exactly by finite mixture sums. def mixture_params(): # Unequal, separated components make posterior uncertainty nontrivial. 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]) def sample_y(n, rng): means, probs, _ = mixture_params() z = rng.choice(len(means), size=n, p=probs) return means[z] + rng.normal(size=n) * np.array([0.35, 0.5, 0.4])[z] def posterior_quantities(x, alpha, sigma): """Exact posterior mean/variance and marginal score for mixture prior.""" means, probs, scales = mixture_params() # x | component k is N(alpha*mean_k, sigma^2 + alpha^2*scale_k^2) varx = sigma * sigma + (alpha * scales) ** 2 logp = np.log(probs) - 0.5 * (np.log(2*np.pi*varx) + (x[:, None]-alpha*means)**2/varx) logp -= logp.max(axis=1, keepdims=True) w = np.exp(logp); w /= w.sum(axis=1, keepdims=True) post_var_k = (scales**2 * sigma**2) / varx post_mean_k = (means * sigma**2 + alpha * scales**2 * x[:, None]) / varx mu = (w * post_mean_k).sum(axis=1) second = (w * (post_var_k + post_mean_k**2)).sum(axis=1) vy = second - mu**2 # Conditional-score mean equals marginal score (Tweedie identity). score = (alpha * mu - x) / (sigma*sigma) floor = (alpha*alpha / sigma**4) * vy return mu, vy, score, floor def bank_floor(x, alpha, sigma, bank): logits = -(x[:, None] - alpha*bank[None, :])**2 / (2*sigma*sigma) logits -= logits.max(axis=1, keepdims=True) p = np.exp(logits); p /= p.sum(axis=1, keepdims=True) mu = p @ bank v = (p * (bank[None, :] - mu[:, None])**2).sum(axis=1) return alpha*alpha / sigma**4 * v def decomposition_check(rng): n = 180000; alpha, sigma = 0.8, 0.9 y = sample_y(n, rng); x = alpha*y + sigma*rng.normal(size=n) target = (alpha*y-x)/sigma**2 _, _, marginal, floor = posterior_quantities(x, alpha, sigma) # Deliberately imperfect score, so equality is not only checked at optimum. model = marginal + 0.35*np.sin(1.7*x) + 0.12*rng.normal(size=n) raw = np.mean((model-target)**2) ideal = np.mean((model-marginal)**2) f = np.mean(floor) return {"raw": float(raw), "ideal": float(ideal), "floor": float(f), "raw_minus_ideal": float(raw-ideal), "abs_identity_error": float(abs(raw-ideal-f)), "relative_identity_error": float(abs(raw-ideal-f)/raw)} def alpha_scaling(rng): # At high noise, posterior is close to prior; F ~= Var(Y)*alpha^2/sigma^4. sigma = 5.0; n = 220000; y = sample_y(n, rng); eps = rng.normal(size=n) vals = [] for a in np.array([0.10, 0.20, 0.35, 0.50, 0.70]): x = a*y + sigma*eps vals.append(np.mean(posterior_quantities(x, a, sigma)[3])) aa = np.array([.10,.20,.35,.50,.70]) slope = float(np.sum(aa**2*np.array(vals))/np.sum(aa**4)) predicted = slope*aa**2 rel = np.max(np.abs(np.array(vals)-predicted)/np.maximum(np.array(vals),1e-15)) return {"alphas": aa.tolist(), "floors": np.array(vals).tolist(), "quadratic_fit_relative_max_error": float(rel), "observed_ratio_F(a=.7)/F(a=.1)": float(vals[-1]/vals[0]), "predicted_ratio_(.7/.1)^2": 49.0} def alpha_zero_and_noise(rng): n=120000; y=sample_y(n,rng) rows=[] for a,s in [(0.0,0.8),(0.8,0.8),(0.8,5.0),(0.8,10.0)]: x=a*y+s*rng.normal(size=n) rows.append((a,s,float(np.mean(posterior_quantities(x,a,s)[3])))) return {"cases":[{"alpha":a,"sigma":s,"floor":f} for a,s,f in rows], "alpha_zero_floor": rows[0][2], "zero_prediction": 0.0} def bank_convergence(rng): # Prediction: iid bank Monte Carlo error decreases approximately B^-1/2. n=9000; alpha=.8; sigma=.9 y=sample_y(n,rng); x=alpha*y+sigma*rng.normal(size=n) exact=posterior_quantities(x,alpha,sigma)[3] out=[] for b in [16,32,64,128,256,512,1024]: errs=[] for _ in range(5): bank=sample_y(b,rng) errs.append(np.mean(np.abs(bank_floor(x,alpha,sigma,bank)-exact))) out.append((b,float(np.mean(errs)))) bs=np.array([z[0] for z in out],float); es=np.array([z[1] for z in out]) slope=float(np.polyfit(np.log(bs),np.log(es),1)[0]) return {"mean_absolute_errors":[{"bank":b,"mae":e} for b,e in out], "loglog_error_slope":slope,"predicted_slope":-0.5} def ranking_experiment(rng): # Same checkpoints (score perturbations), two schedules with different # additive DSM floors. Corrected estimates target marginal-score error. n=100000; y=sample_y(n,rng) checkpoints=[] xref=0.8*y+0.9*rng.normal(size=n) _,_,sref,_=posterior_quantities(xref,.8,.9) for j,noise in enumerate([.04,.10,.18,.28,.40,.60]): # perturbation grows with x and independent score noise pred=sref + noise*(0.7*np.tanh(xref)+rng.normal(size=n)) checkpoints.append(np.mean((pred-sref)**2)) schedules=[(.35,1.5),(.9,.65)] # (alpha,sigma), deliberately different floors rows=[] for a,s in schedules: x=a*y+s*rng.normal(size=n) _,_,m,f=posterior_quantities(x,a,s) for j,noise in enumerate([.04,.10,.18,.28,.40,.60]): pred=m + noise*(0.7*np.tanh(x)+rng.normal(size=n)) raw=float(np.mean((pred-(a*y-x)/s**2)**2)) corr=raw-float(np.mean(f)) ideal=float(np.mean((pred-m)**2)) rows.append({"schedule":[a,s],"checkpoint":j,"raw":raw,"corrected":corr,"ideal":ideal,"floor":float(np.mean(f))}) def rank_corr(vals, ideal): return float(np.corrcoef(np.argsort(np.argsort(vals)),np.argsort(np.argsort(ideal)))[0,1]) summary=[] for sch in schedules: rr=[r for r in rows if r['schedule']==list(sch)] summary.append({"schedule":sch,"floor":rr[0]['floor'], "raw_ideal_rank_corr":rank_corr([r['raw'] for r in rr],[r['ideal'] for r in rr]), "corrected_ideal_rank_corr":rank_corr([r['corrected'] for r in rr],[r['ideal'] for r in rr])}) return {"schedule_summary":summary,"rows":rows} def main(): rng=np.random.default_rng(2713) result={"decomposition":decomposition_check(rng),"alpha_scaling":alpha_scaling(rng), "alpha_zero_and_noise":alpha_zero_and_noise(rng),"bank_convergence":bank_convergence(rng), "ranking":ranking_experiment(rng)} with open("results.json","w") as f: json.dump(result,f,indent=2) print(json.dumps(result,indent=2)) if __name__ == '__main__': main()