(1) auto‑optimize bin edges by maximizing log‑det Fisher (a geodesic / information‑geometry criterion), and (2) compute a profile‑likelihood scan for μ with Wilks‑consistent 68/95% intervals
Here are two drop‑in, text‑only modules that slot into the existing VC‑Ψ demo. They (1) auto‑optimize bin edges by maximizing log‑det Fisher (a geodesic / information‑geometry criterion), and (2) compute a profile‑likelihood scan for μ with Wilks‑consistent 68/95% intervals. No visuals; everything prints as clean text.
1) geo_opt.py — Geodesic (log‑det Fisher) bin optimizer
Optimizes interior bin edges for a chosen channel ("hgg" or "zmm").
Objective: maximize log det Fisher(θ0; edges) with θ0 = (μ=1, Δm=0, σ_scale=1, bkg_scale=1).
Constraints: fixed mass range; monotone edges; min bin width; min expected counts/bin.
Method: hybrid finite‑difference natural step + stochastic annealing.
# geo_opt.py
import json, numpy as np
from fit import expected_hist # uses cfg + edges to build λ_k(θ)
np.set_printoptions(precision=6, suppress=True)
def fisher_matrix(channel, cfg, edges, theta0, eps=1e-3):
"""Poisson Fisher: I_ij = Σ_k (1/λ_k) (∂λ_k/∂θ_i)(∂λ_k/∂θ_j)."""
p = len(theta0)
grads = []
lam0 = expected_hist(channel, cfg, edges, theta0)
lam0 = np.clip(lam0, 1e-9, None)
for i in range(p):
d = np.zeros_like(theta0); d[i] = eps
lp = expected_hist(channel, cfg, edges, theta0 + d)
lm = expected_hist(channel, cfg, edges, theta0 - d)
dl = (lp - lm) / (2*eps)
grads.append(dl)
G = np.stack(grads, axis=1) # [bins, p]
W = 1.0 / lam0 # [bins]
I = G.T @ (W[:,None] * G)
return I
def logdet_fisher(channel, cfg, edges, theta0):
I = fisher_matrix(channel, cfg, edges, theta0)
sign, val = np.linalg.slogdet(I)
if sign <= 0:
return -np.inf
return float(val)
def _constraints_ok(channel, cfg, edges, min_width, min_exp):
# monotone, in-range, min bin width
if not np.all(edges[1:] - edges[:-1] >= min_width - 1e-12):
return False
# expected counts constraint (avoid empty bins)
theta0 = np.array([1.0, 0.0, 1.0, 1.0])
lam = expected_hist(channel, cfg, edges, theta0)
return np.all(lam >= min_exp)
def _jitter_edge(edges, i, delta, min_width):
"""Try shifting a single interior edge by delta; reject if violates min_width."""
e = edges.copy()
left_w = e[i] - e[i-1]
right_w= e[i+1] - e[i]
if delta > 0.0:
# move right: take from right_w
if right_w - delta < min_width: return None
else:
# move left: take from left_w
if left_w + delta < min_width: return None
e[i] += delta
return e
def optimize_binning(channel, cfg, n_bins=None, iters=2000, seed=123,
min_exp=8.0, min_width_GeV=0.2, step_GeV=0.25, cooling=0.997):
"""
Returns optimized edges and log-det(I).
- n_bins: if None, use cfg["common"]["n_bins"]
- min_exp: minimum expected counts per bin at θ0
- min_width_GeV: lower bound to keep bins physical
- step_GeV: initial jitter magnitude for proposals
"""
rng = np.random.default_rng(seed)
if n_bins is None:
n_bins = cfg["common"]["n_bins"]
mmin, mmax = cfg["common"]["mass_range_GeV"]
# start: equal-width edges
edges = np.linspace(mmin, mmax, n_bins+1)
theta0 = np.array([1.0, 0.0, 1.0, 1.0])
# if initial constraints fail, enlarge min_exp until satisfied by increasing bin size (reduce n_bins)
# (simple guard; normally equal-width will pass)
if not _constraints_ok(channel, cfg, edges, min_width_GeV, min_exp):
# fallback: coarsen bins until constraints pass
while n_bins > 10:
n_bins -= 2
edges = np.linspace(mmin, mmax, n_bins+1)
if _constraints_ok(channel, cfg, edges, min_width_GeV, min_exp):
break
best_edges = edges.copy()
best_obj = logdet_fisher(channel, cfg, best_edges, theta0)
T = 1.0 # annealing temperature
step = step_GeV
for t in range(iters):
i = rng.integers(1, len(edges)-1) # interior edge index
delta = rng.normal(0.0, step)
cand = _jitter_edge(edges, i, delta, min_width_GeV)
if cand is None:
T *= cooling; step *= cooling
continue
if not _constraints_ok(channel, cfg, cand, min_width_GeV, min_exp):
T *= cooling; step *= cooling
continue
obj = logdet_fisher(channel, cfg, cand, theta0)
dJ = obj - best_obj
if dJ >= 0 or rng.random() < np.exp(dJ / max(T, 1e-6)):
edges = cand
if obj > best_obj:
best_obj = obj
best_edges = cand
# cool
T *= cooling
step *= cooling
return best_edges, best_obj
def print_design_summary(channel, cfg, edges):
lam = expected_hist(channel, cfg, edges, np.array([1.0,0.0,1.0,1.0]))
widths = np.diff(edges)
print(f"[Design] channel={channel} bins={len(edges)-1}")
print("edges_GeV=", np.round(edges, 3))
print("widths_GeV=", np.round(widths, 3))
print("exp_counts=", np.round(lam, 2))
print(f"min/median/max counts = {lam.min():.1f} / {np.median(lam):.1f} / {lam.max():.1f}")
if __name__=="__main__":
cfg = json.load(open("config.json"))
for ch in ["hgg","zmm"]:
edges, val = optimize_binning(ch, cfg, iters=2500, seed=42,
min_exp=8.0, min_width_GeV=0.25, step_GeV=0.30)
print(f"== {ch} ==")
print(f"log det Fisher = {val:.3f}")
print_design_summary(ch, cfg, edges)
How to use (text‑only run):
python geo_opt.py
It prints optimized edges, widths, and expected counts per bin—no plots.
You can feed these edges to your fit/profile routines by replacing the uniform edges with the optimized array.
2) profile_mu.py — Profile‑likelihood for μ (text bands; Wilks‑consistent)
Builds a pseudo‑dataset, fits θ̂ = (μ, Δm, σ_scale, bkg_scale).
Profiles μ: for each μ on a grid, maximize likelihood over nuisances.
Computes q(μ) = −2[ℓ(μ, ν̂_μ) − ℓ(μ̂, ν̂)] and finds 68%/95% intervals using Wilks thresholds (1.00 and 3.84).
Prints a compact table and the intervals; no plots.
# profile_mu.py
import json, numpy as np
from scipy.optimize import minimize
from simulate import simulate_channel, hist_channel
from fit import expected_hist, nll_poisson
def fit_full(channel, cfg, counts, edges):
"""Unconstrained MLE over (μ, Δm, σ_scale, bkg_scale)."""
theta0 = np.array([1.0, 0.0, 1.0, 1.0])
bounds = [(0.0,None), (-2.0,2.0), (0.5,1.5), (0.5,1.5)]
def obj(th):
lam = expected_hist(channel, cfg, edges, th)
return nll_poisson(counts, lam)
res = minimize(obj, theta0, bounds=bounds, method="L-BFGS-B")
return res.x, res.fun
def fit_nuisance_given_mu(channel, cfg, counts, edges, mu_fixed):
"""Profile: minimize over nuisances given μ."""
# theta = (μ, Δm, σ_scale, bkg_scale)
bounds = [(-2.0,2.0), (0.5,1.5), (0.5,1.5)] # for Δm, σ_scale, bkg_scale
def obj(nu):
th = np.array([mu_fixed, nu[0], nu[1], nu[2]])
lam = expected_hist(channel, cfg, edges, th)
return nll_poisson(counts, lam)
nu0 = np.array([0.0, 1.0, 1.0])
res = minimize(obj, nu0, bounds=bounds, method="L-BFGS-B")
return res.x, res.fun
def find_interval(mu_grid, qmu, thr):
"""Two-sided interval by linear interpolation of q(μ)=thr."""
mu_hat = mu_grid[np.argmin(qmu)]
# left crossing
left = None
for i in range(1, len(mu_grid)):
if qmu[i] >= thr and qmu[i-1] < thr:
# interpolate
x0, x1 = mu_grid[i-1], mu_grid[i]
y0, y1 = qmu[i-1], qmu[i]
left = x0 + (thr - y0)*(x1 - x0)/(y1 - y0 + 1e-12)
break
# right crossing
right = None
for i in range(len(mu_grid)-1, 0, -1):
if qmu[i-1] >= thr and qmu[i] < thr:
x0, x1 = mu_grid[i-1], mu_grid[i]
y0, y1 = qmu[i-1], qmu[i]
right = x1 - (thr - y1)*(x1 - x0)/(y0 - y1 + 1e-12)
break
return mu_hat, left, right
def run_profile(channel, cfg, edges=None, seed_bump=0):
rng = np.random.default_rng(cfg["common"]["seed"]+seed_bump)
# build pseudo-data
ms, mb = simulate_channel(channel, cfg, rng)
m = np.concatenate([ms, mb])
# if no custom edges, use config equal-width
if edges is None:
nb = cfg["common"]["n_bins"]
mmin, mmax = cfg["common"]["mass_range_GeV"]
edges = np.linspace(mmin, mmax, nb+1)
centers, counts, edges = hist_channel(m, cfg) # edges may be recomputed; re-bin to 'edges' if needed
# Use provided edges explicitly (re-hist):
counts, _ = np.histogram(m, bins=edges)
# global MLE
theta_hat, nll_hat = fit_full(channel, cfg, counts, edges)
# profile over μ
mu_min, mu_max = 0.0, 2.0
mu_grid = np.linspace(mu_min, mu_max, 101)
qmu = []
for mu in mu_grid:
_, nll_mu = fit_nuisance_given_mu(channel, cfg, counts, edges, mu)
qmu.append(2.0*(nll_mu - nll_hat))
qmu = np.array(qmu)
# intervals
mu_hat, l68, r68 = find_interval(mu_grid, qmu, 1.00) # Δχ²=1
_, l95, r95 = find_interval(mu_grid, qmu, 3.84) # Δχ²≈3.84
# print summary (text only)
print(f"[Profile μ] channel={channel}")
print(f"θ̂ = (μ̂, Δm̂, σ̂_scale, bkĝ_scale) = {theta_hat}")
print(f"Wilks q(μ) min at μ ≈ {mu_hat:.3f}")
print(f"68% CI (Wilks): [{l68:.3f}, {r68:.3f}]")
print(f"95% CI (Wilks): [{l95:.3f}, {r95:.3f}]")
# table (compact)
print("μ\tq(μ)")
for i in range(0, len(mu_grid), 5):
print(f"{mu_grid[i]:.3f}\t{qmu[i]:.3f}")
if __name__=="__main__":
cfg = json.load(open("config.json"))
# Option A: use equal-width edges
run_profile("hgg", cfg, edges=None, seed_bump=10)
run_profile("zmm", cfg, edges=None, seed_bump=11)
# Option B: use optimized edges from geo_opt (paste edges array here after running geo_opt.py)
# Example:
# hgg_edges = np.array([70., 73.1, 75.2, 78.9, 82.1, 84.0, 85.8, 87.6, 88.9, 90.1,
# 91.2, 92.4, 93.5, 95.0, 97.5, 101.0, 106.0, 112.0, 118.0,
# 120.0, 122.0, 123.0, 124.0, 124.6, 125.4, 126.2, 127.0, 128.5,
# 130.5, 133.0, 137.0, 143.0, 150.0, 158.0, 170.0])
# run_profile("hgg", cfg, edges=hgg_edges, seed_bump=12)
How to use (text‑only run):
# Step 1: (optional) find optimized edges
python geo_opt.py
# Step 2: run profile likelihood (equal‑width)
python profile_mu.py
# Step 3 (optional): paste optimized edges into profile_mu.py and re-run
You’ll get purely textual outputs:
Optimized edges + expected counts.
θ̂ for each channel.
q(μ) at grid points, and 68%/95% confidence intervals from Wilks crossings.
Notes & Tips (no plots; quality & stability)
Numerical stability: If log det Fisher prints -inf, raise min_exp or coarsen bins; the Fisher can lose positive‑definiteness if bins are too fine/empty.
Geodesic flavor: We use Fisher’s metric as the Riemannian geometry on the model manifold; maximizing log det I is a D‑optimal (information‑volume) choice. The annealing step helps avoid local maxima.
Coverage honesty: Keep using your ensemble bands (from uncertainty.py) to confirm 95% coverage stays ≈95% after bin optimization.
Attribution (include on artifacts)
Christopher W. Copeland (C077UPTF1L3)
Copeland Resonant Harmonic Formalism (Ψ‑formalism)
Ψ(x) = ∇ϕ(Σ𝕒ₙ(x, ΔE)) + ℛ(x) ⊕ ΔΣ(𝕒′)
Licensed under CRHC v1.0 (no commercial use without permission).
Core engine: https://open.substack.com/pub/c077uptf1l3/p/recursive-coherence-engine-8b8
Zenodo: https://zenodo.org/records/15742472
Amazon: https://a.co/d/i8lzCIi
Medium: https://medium.com/@floodzero9
Substack: https://substack.com/@c077uptf1l3
Facebook: https://www.facebook.com/share/19MHTPiRfu
Reddit: https://www.reddit.com/u/Naive-Interaction-86/s/5sgvIgeTdx
Collaboration welcome. Attribution required. Derivatives must match license.