""" Quadrature calculations for the noisy-sensor calibration analysis. Computes numerically evaluated expectations for the physical sensor model (interpolator record-averaged fresh risk, cos^2(M0,a), and the ray oracle's record-averaged risk for the physical target Y, with two labelled comparisons to the clean-score ray optimum d0/(d0+gamma^2) of the deterministic stylized design: cross_target_cross_experiment_ratio = ray_oracle_fresh / clean_score_exact -- numerator: the ray oracle on the RANDOM physical calibration ray M0, target Y; denominator: the stylized ray optimum on its DETERMINISTIC design, target y. Different targets AND different designs; NOT an efficiency or oracle-gap ratio within one problem; tends to (1+s2)/s2. same_target_cross_experiment_ratio = ray_oracle_fresh / clean_score_physical, clean_score_physical = (d0+gamma)/(d0+gamma^2), the stylized ray optimum converted to a predictor of Y via Y = sqrt(1-1/gamma) y + eps with eps ~ N(0, 1/gamma) independent of the readings and of the calibration data. Same target, still different designs; tends to 1. On ONE ray the conversion is an increasing affine map of the risk, so a genuine same-ray ratio would be exactly 1 record by record.), by 2-D quadrature over b = M0[1] ~ N(sqrt(k), s2) (Gauss-Hermite) and Q = ||M0_perp||^2 / s2 ~ chi2_{k-1} (Gauss-Legendre in probability space). Writes a "quadrature" block to sensor_calibration_results.json beside this script (updating that receipt when it exists, otherwise creating it), and an "exact_cdf" block holding the three noncentral chi-square lower-tail probabilities recorded at k = 16, sigma^2 = 4, where P/sigma^2 = ||M0||^2/sigma^2 ~ chi^2_k(noncentrality k/sigma^2): nonpositive_floor = F(k) = Pr(P - k sigma^2 <= 0) ~ 0.307 nonpositive_tuned = F(12) = Pr(P <= 48), 48 = tuned magnitude ~ 0.109 stress_window = F(k) - F(k/1.05) = Pr(-0.05 < (P - k sigma^2)/P < 0) ~ 0.044 (the near-pole band; used only by the numerical-guard audit) negative_denominator_regular = F_k(k/1.05) = Pr((P - k sigma^2)/P <= -0.05) at k = 16 (~0.263) and k = 100 (~0.034): the marginal probability of a negative corrected denominator with five-percent margin, i.e. the denominator part of the blog's sign-reversed realization gate (which adds J*Y > 0, TLS regularity, a slope-magnitude window, and legibility on top) with rounding assertions on the three quoted decimals (an "exact_scalars" block adds the remaining recorded closed-form scalars -- correlator pairs, window edges, k_min, null quantiles, anti-alignment -- each with a rounding assertion); the same block records the flip probabilities F(k) at (k, sigma^2) = (16, 4), (100, 4), (1000, 4) and the smallest k with F(k) < 0.01 at sigma^2 = 2, 4, 9 (the 1% thresholds 71, 232, 1013). """ import json import numpy as np from pathlib import Path from scipy.special import roots_hermitenorm, roots_legendre from scipy.stats import chi2 as chi2_dist, ncx2 NH, NQ = 120, 400 def nodes(k, s2): xh, wh = roots_hermitenorm(NH) wh = wh / np.sqrt(2.0 * np.pi) b = np.sqrt(k) + np.sqrt(s2) * xh xg, wg = roots_legendre(NQ) u = 0.5 * (xg + 1.0) wq = 0.5 * wg Q = chi2_dist.ppf(u, k - 1) return b[:, None], wh[:, None], Q[None, :], wq[None, :] def block(k, s2): b, wb, Q, wq = nodes(k, s2) W = wb * wq theta = np.sqrt(k) * b # M0^T a P = b * b + s2 * Q # ||M0||^2 interp = (theta / P - 1.0) ** 2 + s2 / P # conditional fresh risk oracle = s2 * P / (theta * theta + s2 * P) # ray-oracle conditional risk cos2 = theta * theta / (P * k) g = 1.0 + k / s2 d0 = k - 1.0 oracle_fresh = float((W * oracle).sum()) clean_score_exact = d0 / (d0 + g * g) # clean-score target y clean_score_physical = (d0 + g) / (d0 + g * g) # converted to target Y return dict(k=k, sigma2=s2, interpolator_fresh=float((W * interp).sum()), ray_oracle_fresh=oracle_fresh, cos2=float((W * cos2).sum()), clean_score_exact=clean_score_exact, clean_score_physical=clean_score_physical, cross_target_cross_experiment_ratio=oracle_fresh / clean_score_exact, same_target_cross_experiment_ratio=oracle_fresh / clean_score_physical, nodes=dict(hermite=NH, legendre=NQ)) def exact_cdf_block(k=16, s2=4.0): """Exact lower-tail probabilities of P/s2 ~ chi^2_k(k/s2), P = ||M0||^2. Thresholds on the P/s2 scale: the full floor k*s2 maps to k; the tuned magnitude k*s2*(k-1)/(k+s2) maps to t_tuned = k*(k-1)/(k+s2) (12 at k=16, s2=4); the stress window -0.05 < (P-k*s2)/P < 0 maps to (k/1.05, k). The decimal assertions apply to the recorded case (k, s2) = (16, 4.0) only. """ nc = k / s2 F = lambda x: float(ncx2.cdf(x, k, nc)) t_tuned = k * (k - 1) / (k + s2) out = dict(k=k, sigma2=s2, df=k, noncentrality=nc, floor_threshold=float(k), tuned_threshold=float(t_tuned), stress_upper=float(k), stress_lower=float(k / 1.05), nonpositive_floor=F(k), nonpositive_tuned=F(t_tuned), stress_window=F(k) - F(k / 1.05)) # near-pole band (audit only) # Pr(D_EIV/P <= -0.05): negative corrected denominator with five-percent margin. # This is a statement about the denominator only -- not about slope orientation, # TLS regularity, or legibility, which the blog's selection gate adds on top. out["negative_denominator_regular"] = [ dict(k=kk, sigma2=s2, value=float(ncx2.cdf(kk / 1.05, kk, kk / s2))) for kk in (16, 100)] if (k, s2) == (16, 4.0): assert out["tuned_threshold"] == 12.0, out assert round(out["nonpositive_floor"], 3) == 0.307, out assert round(out["nonpositive_tuned"], 3) == 0.109, out assert round(out["stress_window"], 3) == 0.044, out assert [round(r["value"], 3) for r in out["negative_denominator_regular"]] == [0.263, 0.034] # flip probabilities F_{chi^2_k(k/s2)}(k) at other k, and 1% thresholds Fk = lambda kk, ss: float(ncx2.cdf(kk, kk, kk / ss)) out["flip_probability"] = [dict(k=kk, sigma2=ss, value=Fk(kk, ss)) for (kk, ss) in ((16, 4.0), (100, 4.0), (1000, 4.0))] assert [round(r["value"], 3) for r in out["flip_probability"]] == [0.307, 0.066, 0.0] assert 5.145e-7 <= out["flip_probability"][2]["value"] < 5.155e-7 # quoted as 5.15e-7 thresholds = [] for ss in (2.0, 4.0, 9.0): kk = 1 while Fk(kk, ss) >= 0.01: kk += 1 thresholds.append(dict(sigma2=ss, k_min=kk, value=Fk(kk, ss), value_below=Fk(kk - 1, ss))) assert [t["k_min"] for t in thresholds] == [71, 232, 1013], thresholds assert all(t["value"] < 0.01 <= t["value_below"] for t in thresholds) out["one_percent_threshold"] = thresholds return out def exact_scalars_block(): """Recorded closed-form scalars that neither simulation block nor the quadrature grid emits elsewhere: the known-strength sensor correlator's exact expected fresh/replay pair at sigma^2 = 4 for k in {16, 25, 100, 1000, 10000}; the two window edges in sigma^2 (misfit lower root of s4(1+2/k)+4 s2/k = 1 and the prediction upper edge sqrt(k+1)-1) for k in {16, 100, 10000}; the smallest k with k > s2(2+s2) at s2 = 4; the null rule's exact chi^2_1 loss quantiles; and the anti-alignment probability Phi(-sqrt(k/s2)) at (16, 4). Each carries a rounding assertion on its recorded decimals.""" from scipy.stats import chi2 as chi2_law, norm s2 = 4.0 pairs = [dict(k=k, sigma2=s2, fresh=s2 * (2 + s2) / k, replay=s2 ** 2 + (4 * s2 + 2 * s2 ** 2) / k) for k in (16, 25, 100, 1000, 10_000)] # (fresh to 4 decimals, replay to 2 recorded decimals; 10^4 to 4) quoted = {16: (1.5000, 19.00, 2), 25: (0.9600, 17.92, 2), 100: (0.2400, 16.48, 2), 1000: (0.0240, 16.05, 2), 10_000: (0.0024, 16.0048, 4)} for r in pairs: f, rp, nd = quoted[r["k"]] assert round(r["fresh"], 4) == f and round(r["replay"], nd) == rp, r edges = [] for k in (16, 100, 10_000): a, b, c = 1.0 + 2.0 / k, 4.0 / k, -1.0 lower = (-b + np.sqrt(b * b - 4 * a * c)) / (2 * a) edges.append(dict(k=k, misfit_lower_root=float(lower), prediction_upper_edge=float(np.sqrt(k + 1.0) - 1.0))) assert [round(e["misfit_lower_root"], 3) for e in edges[:2]] == [0.838, 0.971] assert round(edges[2]["misfit_lower_root"], 4) == 0.9997 assert [round(e["prediction_upper_edge"], 3) for e in edges] == [3.123, 9.050, 99.005] k_min = int(np.floor(s2 * (2 + s2))) + 1 assert k_min == 25 null_q = dict(median=float(chi2_law.ppf(0.5, 1)), q90=float(chi2_law.ppf(0.9, 1))) assert round(null_q["median"], 4) == 0.4549 and round(null_q["q90"], 4) == 2.7055 anti = float(norm.cdf(-np.sqrt(16 / s2))) assert round(anti, 5) == 0.02275 return dict(correlator_exact_pairs=pairs, window_edges=edges, fresh_window_k_min=dict(sigma2=s2, k_min=k_min), null_rule_loss_quantiles=null_q, anti_alignment_probability=dict(k=16, sigma2=s2, value=anti)) def has_simulation_results(data): """Recognize computed Monte Carlo records, not merely an existing receipt.""" blocks = data.get("blocks") if isinstance(blocks, list) and any( isinstance(block, dict) and isinstance(block.get("draws"), (int, float)) and block["draws"] > 0 and isinstance(block.get("rules"), dict) and block["rules"] for block in blocks): return True multi = data.get("multi_record") return isinstance(multi, list) and any( isinstance(block, dict) and isinstance(block.get("draws"), (int, float)) and block["draws"] > 0 and all(key in block for key in ("fresh_sim", "replay_sim", "oracle_sim")) for block in multi) def main(): cases = [(16, 4.0), (100, 4.0), (1000, 4.0), (10_000, 4.0), (1_000_000, 4.0), (10_000, 2.0)] out = [block(k, s2) for (k, s2) in cases] # The published site ships this script without the simulation receipt. # Update the receipt when it is present; otherwise start a new one. path = Path(__file__).with_name("sensor_calibration_results.json") data = json.loads(path.read_text()) if path.exists() else {} simulation_receipt = has_simulation_results(data) data["quadrature"] = out data["exact_cdf"] = exact_cdf_block() data["exact_scalars"] = exact_scalars_block() import platform, scipy data["regeneration"] = dict( python=platform.python_version(), numpy=np.__version__, scipy=scipy.__version__, note="Existing Monte Carlo simulation results are retained. This run updates " "only the quadrature, exact-CDF and closed-form blocks. The exact-CDF " "and closed-form quoted values carry rounding or ordering assertions; " "quadrature values are numerical evaluations on the declared grid. " "Byte identity is not a cross-environment guarantee." if simulation_receipt else "Written by sensor_calibration_quadrature.py alone, without the simulation " "receipt: it holds only the quadrature, exact-CDF and closed-form blocks. " "The exact-CDF and closed-form quoted values carry rounding or ordering " "assertions; quadrature values are numerical evaluations on the declared grid. " "Byte identity is not a cross-environment guarantee.") path.write_text(json.dumps(data, indent=2, allow_nan=False)) e = data["exact_cdf"] print(f"exact_cdf k={e['k']} s2={e['sigma2']:.1f} F(16)={e['nonpositive_floor']:.10f}" f" F(12)={e['nonpositive_tuned']:.10f} stress={e['stress_window']:.10f}") for r in out: print(f"k={r['k']:>8} s2={r['sigma2']:.1f} interp={r['interpolator_fresh']:.7f}" f" oracle={r['ray_oracle_fresh']:.7f} cos2={r['cos2']:.4f}" f" clean={r['clean_score_exact']:.7f}" f" cross_target_cross_experiment_ratio={r['cross_target_cross_experiment_ratio']:.4f}" f" same_target_cross_experiment_ratio={r['same_target_cross_experiment_ratio']:.4f}") if __name__ == "__main__": main()