""" Asymmetric Lithwick-Wu pair model. For each first-order mean-motion resonance r = (k+1)/k (k = 1, 2, 3, 4), model the local density structure as an asymmetric Gaussian pair: f_LW(x | c) = (1/(1+γ)) · N(x; c−ε, σ) + (γ/(1+γ)) · N(x; c+ε, σ) where x = {r} is the fractional part, c is the resonance's fractional center, ε is the offset of the pile-ups from the exact resonance, σ is the pile-up width, and γ is the inside-vs-outside amplitude ratio. γ = 1 : symmetric (no LW asymmetry) γ > 1 : enhanced "above" — the LW prediction γ < 1 : enhanced "below" Total model: p(x) = α · g(x) + (1−α) · (1/K) Σ_res f_LW(x | c_res) with g(x) = Gauss measure on [0,1), K = 4 first-order resonances. Free parameters: α, ε, σ, γ (four). We compare: (i) MLE of the four-parameter LW model (ii) Symmetric version (γ = 1, three parameters) (iii) Per-resonance γ (four extra parameters) """ import math import csv from collections import defaultdict import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt LN2 = math.log(2.0) FIRST_ORDER_RES = [(2, 1), (3, 2), (4, 3), (5, 4)] FIRST_ORDER_CENTERS = np.array([0.0, 0.5, 1.0 / 3.0, 0.25]) FIRST_ORDER_LABELS = ["2:1", "3:2", "4:3", "5:4"] def load_ratios(path): systems = defaultdict(list) with open(path) as f: for row in csv.DictReader(f): try: p = float(row["pl_orbper"]) except (ValueError, TypeError): continue if math.isfinite(p) and p > 0: systems[row["hostname"]].append(p) out = [] for periods in systems.values(): if len(periods) < 2: continue periods = sorted(periods) for i in range(len(periods) - 1): r = periods[i + 1] / periods[i] if r > 1 and math.isfinite(r): out.append(r) return np.array(out) def gauss_density(x): return 1.0 / ((1.0 + x) * LN2) def wrapped_gaussian_vec(x, mu, sigma, n_wraps=3): ks = np.arange(-n_wraps, n_wraps + 1) delta = x[:, None] - mu + ks[None, :] coef = 1.0 / (sigma * math.sqrt(2 * math.pi)) return (coef * np.exp(-0.5 * (delta / sigma) ** 2)).sum(axis=1) def lw_density_vec(x, centers, eps, sigma, gammas): """Density at points x given (centers, eps, sigma) and per-center gammas (array same length as centers).""" K = len(centers) total = np.zeros_like(x) for c, gamma in zip(centers, gammas): norm = 1.0 / (1.0 + gamma) total = total + norm * wrapped_gaussian_vec(x, c - eps, sigma) total = total + gamma * norm * wrapped_gaussian_vec(x, c + eps, sigma) return total / K def lw_normalizer(centers, eps, sigma, gammas, n_grid=2001): xs = np.linspace(0, 1, n_grid) ys = lw_density_vec(xs, centers, eps, sigma, gammas) return float(np.trapezoid(ys, xs)) def log_likelihood(fracs, alpha, eps, sigma, gammas, centers, gauss_base): f = lw_density_vec(fracs, centers, eps, sigma, gammas) Z = lw_normalizer(centers, eps, sigma, gammas) if Z <= 0: return -np.inf f_norm = f / Z p = alpha * gauss_base + (1 - alpha) * f_norm if np.any(p <= 0): return -np.inf return float(np.log(p).sum()) def fit_shared_gamma(fracs, centers, alpha_grid, eps_grid, sigma_grid, gamma_grid, gauss_base): """4-D grid search with shared γ across centers. Z and f depend only on (eps, sigma, gamma) — recompute them in the outer triple loop, then sweep α cheaply inside. """ best = (-np.inf, None, None, None, None) for eps in eps_grid: for sigma in sigma_grid: for gamma in gamma_grid: gammas = np.full(len(centers), gamma) f = lw_density_vec(fracs, centers, eps, sigma, gammas) Z = lw_normalizer(centers, eps, sigma, gammas) if Z <= 0: continue f_norm = f / Z for alpha in alpha_grid: p = alpha * gauss_base + (1 - alpha) * f_norm if np.any(p <= 0): continue ll = float(np.log(p).sum()) if ll > best[0]: best = (ll, alpha, eps, sigma, gamma) return best def main(): rs = load_ratios("ps_dm.csv") fracs = rs - np.floor(rs) fracs = fracs[(fracs >= 0) & (fracs < 1)] print(f"n fractional parts: {len(fracs)}") gauss_base = gauss_density(fracs) ll_pure_gauss = float(np.log(gauss_base).sum()) print(f"pure Gauss LL: {ll_pure_gauss:.2f}") centers = FIRST_ORDER_CENTERS # ----------- Coarse grid search (LW-scale parameters) ------------ # Lithwick–Wu (2012) found the pile-up at ε ≈ 0.01–0.04 (1–4% offset) # with width σ ≈ 0.005–0.015 (sharp). We constrain the grid to these # ranges so the fit tests the LW prediction at the right physical # scale, not the broader shoulder structure picked up by Test 8. alpha_grid = np.linspace(0.50, 0.99, 25) eps_grid = np.linspace(0.005, 0.04, 18) sigma_grid = np.exp(np.linspace(math.log(0.003), math.log(0.020), 16)) gamma_grid = np.exp(np.linspace(math.log(0.2), math.log(5.0), 21)) n = len(alpha_grid) * len(eps_grid) * len(sigma_grid) * len(gamma_grid) print(f"\ncoarse 4-D grid: {n} cells") ll, alpha, eps, sigma, gamma = fit_shared_gamma( fracs, centers, alpha_grid, eps_grid, sigma_grid, gamma_grid, gauss_base ) print(f"best: α={alpha:.3f} ε={eps:.4f} σ={sigma:.4f} γ={gamma:.3f}") print(f"LL = {ll:.2f} ΔLL vs pure Gauss = {ll - ll_pure_gauss:+.2f}") # ----------- Refined grid ------------ a_lo, a_hi = max(0.5, alpha - 0.05), min(0.99, alpha + 0.05) eps_lo, eps_hi = max(0.002, eps * 0.7), min(0.05, eps * 1.3) sig_lo, sig_hi = max(0.002, sigma * 0.7), min(0.03, sigma * 1.3) gam_lo, gam_hi = max(0.1, gamma * 0.5), min(8.0, gamma * 2.0) alpha_grid2 = np.linspace(a_lo, a_hi, 15) eps_grid2 = np.linspace(eps_lo, eps_hi, 15) sigma_grid2 = np.linspace(sig_lo, sig_hi, 15) gamma_grid2 = np.linspace(gam_lo, gam_hi, 15) print(f"\nrefined grid: {15**4} cells") ll, alpha, eps, sigma, gamma = fit_shared_gamma( fracs, centers, alpha_grid2, eps_grid2, sigma_grid2, gamma_grid2, gauss_base ) print(f"refined best: α={alpha:.3f} ε={eps:.4f} σ={sigma:.4f} γ={gamma:.3f}") print(f"LL = {ll:.2f} ΔLL vs pure Gauss = {ll - ll_pure_gauss:+.2f}") # ----------- Compare to symmetric (γ=1) at same (α, ε, σ) ------------ gammas_sym = np.full(len(centers), 1.0) ll_sym = log_likelihood(fracs, alpha, eps, sigma, gammas_sym, centers, gauss_base) dll_asym = ll - ll_sym print(f"\nSymmetric (γ=1) at same (α, ε, σ): LL = {ll_sym:.2f}") print(f"ΔLL_asym (asymmetry vs symmetry) = {dll_asym:+.2f}") print(f"2·ΔLL_asym = {2*dll_asym:.2f} (chi² on 1 dof; critical@0.001 = 10.83)") # ----------- Best symmetric fit (γ=1), profiling over (α, ε, σ) ------------ print("\n--- finding best SYMMETRIC fit (γ=1 fixed) ---") best_sym = (-np.inf, None, None, None) for e in eps_grid2: for s in sigma_grid2: f_s = lw_density_vec(fracs, centers, e, s, gammas_sym) Z_s = lw_normalizer(centers, e, s, gammas_sym) if Z_s <= 0: continue f_n_s = f_s / Z_s for a in alpha_grid2: p_s = a * gauss_base + (1 - a) * f_n_s if np.any(p_s <= 0): continue ll_s = float(np.log(p_s).sum()) if ll_s > best_sym[0]: best_sym = (ll_s, a, e, s) ll_best_sym, a_sym, e_sym, s_sym = best_sym print(f"best symmetric: α={a_sym:.3f} ε={e_sym:.4f} σ={s_sym:.4f} LL = {ll_best_sym:.2f}") print(f"ΔLL(asym − best sym) = {ll - ll_best_sym:+.2f} (4 vs 3 params, 1 dof)") # ----------- Profile likelihood for γ ------------ print("\n--- profile likelihood for γ (other params fixed at MLE) ---") gamma_profile = np.exp(np.linspace(math.log(0.3), math.log(5.0), 81)) lls_gamma = [] for g in gamma_profile: gammas = np.full(len(centers), g) ll_g = log_likelihood(fracs, alpha, eps, sigma, gammas, centers, gauss_base) lls_gamma.append(ll_g) lls_gamma = np.array(lls_gamma) ll_max = lls_gamma.max() g_mle = gamma_profile[np.argmax(lls_gamma)] in_ci68 = gamma_profile[lls_gamma > ll_max - 0.5] in_ci95 = gamma_profile[lls_gamma > ll_max - 1.92] print(f"γ MLE: {g_mle:.3f}") print(f"68% CI: [{in_ci68.min():.3f}, {in_ci68.max():.3f}]") print(f"95% CI: [{in_ci95.min():.3f}, {in_ci95.max():.3f}]") # Test against γ=1 ll_at_1 = log_likelihood(fracs, alpha, eps, sigma, np.full(len(centers), 1.0), centers, gauss_base) print(f"LL at γ=1: {ll_at_1:.2f}, ΔLL(MLE − γ=1) = {ll_max - ll_at_1:+.2f}") print(f"Significance of asymmetry: z_eff = √(2·ΔLL) = {math.sqrt(2 * max(0, ll_max - ll_at_1)):.2f}") # ----------- Per-resonance γ fit (one γ per center) ------------ print("\n--- per-resonance γ fit (keep α, ε, σ at shared MLE) ---") gamma_grid_pr = np.exp(np.linspace(math.log(0.2), math.log(5.0), 25)) per_res_gammas = [] per_res_lls = [] for i, center in enumerate(centers): best_g = None best_ll_at = -np.inf for g in gamma_grid_pr: gammas = np.ones(len(centers)) # others at γ=1 gammas[i] = g ll_g = log_likelihood(fracs, alpha, eps, sigma, gammas, centers, gauss_base) if ll_g > best_ll_at: best_ll_at = ll_g best_g = g per_res_gammas.append(best_g) per_res_lls.append(best_ll_at) print(f" {FIRST_ORDER_LABELS[i]:>4} best γ (with others=1) = {best_g:.3f} " f"LL = {best_ll_at:.2f}") # Joint per-resonance γ fit (4-D grid in γ alone, others fixed) print("\n joint 4-γ fit (15^4 = 50625 cells)...") g_grid = np.exp(np.linspace(math.log(0.3), math.log(4.0), 15)) best_joint = (-np.inf, None) for g1 in g_grid: for g2 in g_grid: for g3 in g_grid: for g4 in g_grid: gammas = np.array([g1, g2, g3, g4]) ll_j = log_likelihood(fracs, alpha, eps, sigma, gammas, centers, gauss_base) if ll_j > best_joint[0]: best_joint = (ll_j, gammas.copy()) print(f" joint best: γ = {best_joint[1]} LL = {best_joint[0]:.2f}") print(f" ΔLL(per-resonance − shared γ) = {best_joint[0] - ll_max:+.2f} (3 extra dof)") print(f" 2·ΔLL = {2*(best_joint[0] - ll_max):.2f} (critical chi²(3)@0.05 ≈ 7.81)") # ----------- Plots ------------ # Density plot xs = np.linspace(0.0001, 0.9999, 2000) gammas_mle = np.full(len(centers), gamma) ys_gauss = gauss_density(xs) ys_lw = lw_density_vec(xs, centers, eps, sigma, gammas_mle) Z = lw_normalizer(centers, eps, sigma, gammas_mle) ys_lw /= Z ys_mix = alpha * ys_gauss + (1 - alpha) * ys_lw # also LW with γ=1 ys_lw_sym = lw_density_vec(xs, centers, eps, sigma, np.full(len(centers), 1.0)) Zsym = lw_normalizer(centers, eps, sigma, np.full(len(centers), 1.0)) ys_lw_sym /= Zsym fig, ax = plt.subplots(figsize=(9.5, 5)) ax.hist(fracs, bins=80, density=True, color="#7f8c8d", alpha=0.4, edgecolor="white", linewidth=0.5, label=f"data (n={len(fracs)})") ax.plot(xs, ys_gauss, color="#2c3e50", lw=1.5, alpha=0.6, label="pure Gauss") ax.plot(xs, ys_lw_sym, color="#8e44ad", lw=1.3, ls="--", alpha=0.7, label=fr"symmetric LW pairs (γ=1, ε={eps:.3f}, σ={sigma:.3f})") ax.plot(xs, ys_lw, color="#c0392b", lw=1.5, alpha=0.85, label=fr"asymmetric LW pairs (γ={gamma:.2f})") ax.plot(xs, ys_mix, color="#16a085", lw=2.5, label=fr"best mixture: α={alpha:.3f}, ΔLL = {ll - ll_pure_gauss:.1f}") for c, lbl in zip(centers, FIRST_ORDER_LABELS): ax.axvline(c, color="#c0392b", lw=0.6, alpha=0.4) ax.text(c + 0.005, ax.get_ylim()[1] * 0.93, lbl, fontsize=8, color="#c0392b") ax.set_xlabel(r"$\{r\}$") ax.set_ylabel("density") ax.set_title(rf"Lithwick–Wu pair model: $\alpha = {alpha:.3f}$, $\varepsilon = {eps:.4f}$, " rf"$\sigma = {sigma:.4f}$, $\gamma = {gamma:.2f}$") ax.legend(loc="upper left", fontsize=8) ax.set_xlim(0, 1) fig.tight_layout() fig.savefig("c1_lw_pair_fit.png", dpi=140) # γ profile plot fig, ax = plt.subplots(figsize=(8, 4.5)) ax.plot(gamma_profile, lls_gamma, color="#2c3e50", lw=2) ax.axvline(1.0, color="#888", lw=1, ls=":", label="γ = 1 (symmetric)") ax.axvline(gamma, color="#c0392b", lw=1.5, ls="--", label=fr"MLE γ = {gamma:.3f}") ax.axhline(ll_max - 0.5, color="#bbb", lw=0.8, ls=":") ax.axhline(ll_max - 1.92, color="#ddd", lw=0.8, ls=":") ax.axvspan(in_ci68.min(), in_ci68.max(), color="#16a085", alpha=0.20, label=f"68%: [{in_ci68.min():.2f}, {in_ci68.max():.2f}]") ax.axvspan(in_ci95.min(), in_ci95.max(), color="#16a085", alpha=0.08, label=f"95%: [{in_ci95.min():.2f}, {in_ci95.max():.2f}]") ax.set_xscale("log") ax.set_xlabel(r"$\gamma$ (above/below amplitude ratio)") ax.set_ylabel("log-likelihood") ax.set_title(rf"Profile log-likelihood for $\gamma$ (α, ε, σ at MLE)") ax.legend(fontsize=8) fig.tight_layout() fig.savefig("c1_lw_gamma_profile.png", dpi=140) # Per-resonance γ bar fig, ax = plt.subplots(figsize=(7, 4.5)) pr_gammas = np.array(per_res_gammas) joint_gammas = best_joint[1] x_pos = np.arange(len(FIRST_ORDER_LABELS)) w = 0.35 ax.bar(x_pos - w/2, pr_gammas, w, color="#c0392b", alpha=0.8, label="single-resonance fit (others γ=1)") ax.bar(x_pos + w/2, joint_gammas, w, color="#16a085", alpha=0.8, label="joint 4-γ fit") ax.axhline(1.0, color="#888", lw=1, ls=":", label="γ = 1 (symmetric)") ax.axhline(gamma, color="#2c3e50", lw=1, ls="--", label=f"shared MLE γ = {gamma:.2f}") ax.set_xticks(x_pos) ax.set_xticklabels(FIRST_ORDER_LABELS) ax.set_xlabel("first-order resonance") ax.set_ylabel(r"$\gamma$") ax.set_title("Per-resonance LW asymmetry") ax.set_yscale("log") ax.legend(fontsize=8) fig.tight_layout() fig.savefig("c1_lw_per_res.png", dpi=140) print("\nplots: c1_lw_pair_fit.png c1_lw_gamma_profile.png c1_lw_per_res.png") if __name__ == "__main__": main()