""" Joint Lithwick-Wu + broad-Farey mixture model. Three components on the fractional-part distribution {r} ∈ [0, 1): p(x) = α₁ · g(x) (Gauss baseline) + α₂ · f_LW(x | ε_LW, σ_LW, γ_LW) (narrow asymmetric pairs at 4 first-order MMRs) + α₃ · f_Farey(x | σ_F) (broad bumps at 22 Farey-12 (p,q) centers) with α₁ + α₂ + α₃ = 1, all ≥ 0. LW parameters are held at the Test 9 MLE (ε=0.017, σ_LW=0.0036, γ=5.0). The broad-Farey component uses 22 per-(p,q) centers (not collapsed) with a single free bandwidth σ_F. Free parameters: α₁, α₂, σ_F (with α₃ = 1 − α₁ − α₂). This tests whether the data prefers BOTH components together over either alone, and quantifies the relative weight of the narrow asymmetric LW component vs the broad Farey shoulders. """ 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) # Test 9 MLE for the LW pair component LW_EPS = 0.0166 LW_SIGMA = 0.0036 LW_GAMMA = 5.0 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"] # ---------------- I/O ------------------------------------------------ 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 farey_pairs(max_sum): out = [] for q in range(1, max_sum): for p in range(q + 1, max_sum - q + 1): if math.gcd(p, q) == 1: out.append((p, q)) return out def farey12_per_pq_centers(): """Per-(p,q) centers and labels for Farey-12 (no collapsing).""" out = [] for p, q in farey_pairs(12): c = (p % q) / q if q > 1 else 0.0 out.append((p, q, c)) return out # ---------------- densities ------------------------------------------ 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, gamma): K = len(centers) total = np.zeros_like(x) norm = 1.0 / (1.0 + gamma) for c in centers: total += norm * wrapped_gaussian_vec(x, c - eps, sigma) total += gamma * norm * wrapped_gaussian_vec(x, c + eps, sigma) return total / K def farey_broad_density_vec(x, centers_arr, sigma): """Equal-weighted broad Farey-12 component.""" K = len(centers_arr) total = np.zeros_like(x) for c in centers_arr: total += wrapped_gaussian_vec(x, c, sigma) return total / K def normalize(density_fn, n_grid=2001): xs = np.linspace(0, 1, n_grid) ys = density_fn(xs) return float(np.trapezoid(ys, xs)) # ---------------- main routine --------------------------------------- def main(): rs = load_ratios("ps_dm.csv") fracs = rs - np.floor(rs) fracs = fracs[(fracs >= 0) & (fracs < 1)] n = len(fracs) print(f"n fractional parts: {n}") g_at_data = gauss_density(fracs) ll_pure_gauss = float(np.log(g_at_data).sum()) print(f"pure Gauss LL: {ll_pure_gauss:.2f}") # LW component (fixed) -- compute density at data and normalize once lw_at_data_raw = lw_density_vec(fracs, FIRST_ORDER_CENTERS, LW_EPS, LW_SIGMA, LW_GAMMA) lw_norm_fn = lambda x: lw_density_vec(x, FIRST_ORDER_CENTERS, LW_EPS, LW_SIGMA, LW_GAMMA) Z_lw = normalize(lw_norm_fn, n_grid=4001) lw_at_data = lw_at_data_raw / Z_lw print(f"LW norm Z_LW = {Z_lw:.4f} (≈ 1 expected)") # Broad-Farey component (per-(p,q) centers, free σ_F) centers_pq = farey12_per_pq_centers() centers_arr = np.array([c for _, _, c in centers_pq]) print(f"Farey-12 centers (per-(p,q)): {len(centers_arr)}") # ---------- compare component-only baselines ---------- print("\n--- baselines ---") # LW only (α₁ = 0, α₃ = 0, α₂ = 1): just LW component p_lw_only = lw_at_data if np.any(p_lw_only <= 0): ll_lw_only = -np.inf else: ll_lw_only = float(np.log(p_lw_only).sum()) print(f"LW-only (α₂ = 1): LL = {ll_lw_only:.2f} (probably -∞: LW pairs are sharp)") # Sweep σ_F to find the best Farey-broad-only baseline sigma_F_grid = np.exp(np.linspace(math.log(0.01), math.log(0.5), 31)) best_F_only = (-np.inf, None) for sF in sigma_F_grid: f_F = farey_broad_density_vec(fracs, centers_arr, sF) Z_F = normalize(lambda x: farey_broad_density_vec(x, centers_arr, sF)) if Z_F <= 0: continue f_F_norm = f_F / Z_F if np.any(f_F_norm <= 0): continue ll = float(np.log(f_F_norm).sum()) if ll > best_F_only[0]: best_F_only = (ll, sF) print(f"Farey-broad-only (α₃ = 1): best σ_F = {best_F_only[1]:.4f} LL = {best_F_only[0]:.2f} " f"ΔLL = {best_F_only[0] - ll_pure_gauss:+.2f}") # ---------- joint 3-component fit ---------- print("\n--- joint Gauss + LW + Farey-broad mixture fit ---") # 3-D grid: (α₁, α₂, σ_F). α₃ = 1 - α₁ - α₂. # Coarse alpha1_grid = np.linspace(0.0, 1.0, 31) alpha2_grid = np.linspace(0.0, 0.20, 21) # LW component should be small (~0.02 from Test 9) sigma_F_grid_c = np.exp(np.linspace(math.log(0.01), math.log(0.5), 21)) best = (-np.inf, None, None, None) print(f"coarse 3-D grid: {len(alpha1_grid) * len(alpha2_grid) * len(sigma_F_grid_c)} cells") for sF in sigma_F_grid_c: f_F = farey_broad_density_vec(fracs, centers_arr, sF) Z_F = normalize(lambda x: farey_broad_density_vec(x, centers_arr, sF)) if Z_F <= 0: continue f_F_norm = f_F / Z_F for a1 in alpha1_grid: for a2 in alpha2_grid: a3 = 1.0 - a1 - a2 if a3 < -1e-9: continue a3 = max(a3, 0.0) p = a1 * g_at_data + a2 * lw_at_data + a3 * f_F_norm if np.any(p <= 0): continue ll = float(np.log(p).sum()) if ll > best[0]: best = (ll, a1, a2, sF) ll_best, a1_b, a2_b, sF_b = best a3_b = 1.0 - a1_b - a2_b print(f"coarse best: α₁ = {a1_b:.3f} α₂ = {a2_b:.3f} α₃ = {a3_b:.3f} σ_F = {sF_b:.4f}") print(f" LL = {ll_best:.2f} ΔLL vs pure Gauss = {ll_best - ll_pure_gauss:+.2f}") # Refined a1_lo = max(0.0, a1_b - 0.1); a1_hi = min(1.0, a1_b + 0.1) a2_lo = max(0.0, a2_b - 0.05); a2_hi = min(0.3, a2_b + 0.05) sF_lo = max(0.005, sF_b * 0.7); sF_hi = min(0.6, sF_b * 1.4) alpha1_g2 = np.linspace(a1_lo, a1_hi, 21) alpha2_g2 = np.linspace(a2_lo, a2_hi, 21) sigma_F_g2 = np.linspace(sF_lo, sF_hi, 21) best2 = (-np.inf, None, None, None) print(f"refined 3-D grid: {len(alpha1_g2) * len(alpha2_g2) * len(sigma_F_g2)} cells") for sF in sigma_F_g2: f_F = farey_broad_density_vec(fracs, centers_arr, sF) Z_F = normalize(lambda x: farey_broad_density_vec(x, centers_arr, sF)) if Z_F <= 0: continue f_F_norm = f_F / Z_F for a1 in alpha1_g2: for a2 in alpha2_g2: a3 = 1.0 - a1 - a2 if a3 < -1e-9: continue a3 = max(a3, 0.0) p = a1 * g_at_data + a2 * lw_at_data + a3 * f_F_norm if np.any(p <= 0): continue ll = float(np.log(p).sum()) if ll > best2[0]: best2 = (ll, a1, a2, sF) ll_best, a1_b, a2_b, sF_b = best2 a3_b = 1.0 - a1_b - a2_b print(f"refined best: α₁ = {a1_b:.3f} α₂ = {a2_b:.3f} α₃ = {a3_b:.3f} σ_F = {sF_b:.4f}") print(f" LL = {ll_best:.2f} ΔLL vs pure Gauss = {ll_best - ll_pure_gauss:+.2f}") # ---------- profile-likelihood for α₂ (LW component weight) ---------- print("\n--- profile likelihood for α₂ (LW weight), other params fixed ---") f_F_at_best = farey_broad_density_vec(fracs, centers_arr, sF_b) Z_F_best = normalize(lambda x: farey_broad_density_vec(x, centers_arr, sF_b)) f_F_norm_best = f_F_at_best / Z_F_best a2_profile = np.linspace(0.0, 0.20, 81) lls_a2 = [] for a2 in a2_profile: a1 = a1_b # keep alpha_1 fixed a3 = 1.0 - a1 - a2 if a3 < 0: lls_a2.append(-np.inf) continue p = a1 * g_at_data + a2 * lw_at_data + a3 * f_F_norm_best if np.any(p <= 0): lls_a2.append(-np.inf) continue lls_a2.append(float(np.log(p).sum())) lls_a2 = np.array(lls_a2) ll_max = lls_a2.max() a2_mle_prof = a2_profile[np.argmax(lls_a2)] in_ci68 = a2_profile[lls_a2 > ll_max - 0.5] in_ci95 = a2_profile[lls_a2 > ll_max - 1.92] a2_ci = (in_ci68.min(), in_ci68.max()) if len(in_ci68) else (a2_mle_prof, a2_mle_prof) a2_95 = (in_ci95.min(), in_ci95.max()) if len(in_ci95) else (a2_mle_prof, a2_mle_prof) print(f"α₂ MLE (profile, α₁ fixed): {a2_mle_prof:.4f}") print(f"68% CI: [{a2_ci[0]:.4f}, {a2_ci[1]:.4f}]") print(f"95% CI: [{a2_95[0]:.4f}, {a2_95[1]:.4f}]") # Compare to α₂ = 0 (no LW component) a3_no_lw = 1 - a1_b p_no_lw = a1_b * g_at_data + a3_no_lw * f_F_norm_best ll_no_lw = float(np.log(p_no_lw).sum()) if np.all(p_no_lw > 0) else -np.inf print(f"\nLL at α₂ = 0 (no LW component): {ll_no_lw:.2f}") print(f"ΔLL(joint MLE − Gauss+Farey only) = {ll_best - ll_no_lw:+.2f}") print(f"2·ΔLL = {2 * (ll_best - ll_no_lw):.2f} (chi² on 1 dof; critical@0.001 ≈ 10.83)") print(f"z_eff for LW component presence: {math.sqrt(2 * max(0, ll_best - ll_no_lw)):.2f}") # ---------- visualizations ---------- xs = np.linspace(0.0001, 0.9999, 2000) ys_gauss = gauss_density(xs) ys_lw_raw = lw_density_vec(xs, FIRST_ORDER_CENTERS, LW_EPS, LW_SIGMA, LW_GAMMA) ys_lw = ys_lw_raw / Z_lw ys_farey_raw = farey_broad_density_vec(xs, centers_arr, sF_b) ys_farey = ys_farey_raw / Z_F_best ys_total = a1_b * ys_gauss + a2_b * ys_lw + a3_b * ys_farey fig, ax = plt.subplots(figsize=(10, 5.5)) ax.hist(fracs, bins=80, density=True, color="#7f8c8d", alpha=0.4, edgecolor="white", linewidth=0.5, label=f"data (n={n})") ax.plot(xs, a1_b * ys_gauss, color="#2c3e50", lw=1.5, alpha=0.7, label=fr"α₁·Gauss ({a1_b:.3f})") ax.plot(xs, a2_b * ys_lw, color="#c0392b", lw=1.4, alpha=0.85, label=fr"α₂·LW-pair ({a2_b:.3f})") ax.plot(xs, a3_b * ys_farey, color="#9b59b6", lw=1.4, alpha=0.7, label=fr"α₃·Farey-broad ({a3_b:.3f}, σ_F={sF_b:.3f})") ax.plot(xs, ys_total, color="#16a085", lw=2.5, label=fr"total mixture, ΔLL = {ll_best - ll_pure_gauss:+.1f}") # Mark first-order MMR centers for c, lbl in zip(FIRST_ORDER_CENTERS, FIRST_ORDER_LABELS): ax.axvline(c, color="#c0392b", lw=0.5, alpha=0.4) ax.text(c + 0.005, ax.get_ylim()[1] * 0.93, lbl, fontsize=8, color="#c0392b") ax.set_xlabel(r"$\{r\} = r - \lfloor r \rfloor$") ax.set_ylabel("density") ax.set_title(rf"Joint Gauss + LW-pair + Farey-broad mixture ($\alpha_1, \alpha_2, \alpha_3, \sigma_F$)") ax.legend(loc="upper left", fontsize=9) ax.set_xlim(0, 1) fig.tight_layout() fig.savefig("c1_joint_fit.png", dpi=140) # Profile likelihood for α₂ fig, ax = plt.subplots(figsize=(8, 4.5)) ax.plot(a2_profile, lls_a2, color="#2c3e50", lw=2) ax.axvline(a2_b, color="#c0392b", lw=1.5, ls="--", label=fr"MLE α₂ = {a2_b:.4f}") ax.axvline(0.0, color="#888", lw=1, ls=":", label="α₂ = 0 (no LW)") ax.axhline(ll_max - 0.5, color="#bbb", lw=0.8, ls=":", label="68%") ax.axhline(ll_max - 1.92, color="#ddd", lw=0.8, ls=":", label="95%") ax.axvspan(a2_ci[0], a2_ci[1], color="#16a085", alpha=0.20, label=f"68%: [{a2_ci[0]:.3f}, {a2_ci[1]:.3f}]") ax.axvspan(a2_95[0], a2_95[1], color="#16a085", alpha=0.08, label=f"95%: [{a2_95[0]:.3f}, {a2_95[1]:.3f}]") ax.set_xlabel(r"$\alpha_2$ (LW component weight)") ax.set_ylabel("log-likelihood") ax.set_title(rf"Profile log-likelihood for $\alpha_2$ (joint fit, $\alpha_1$, $\sigma_F$ at MLE)") ax.legend(fontsize=8) fig.tight_layout() fig.savefig("c1_joint_a2_profile.png", dpi=140) # Bar comparison: each model's LL fig, ax = plt.subplots(figsize=(8, 4.5)) models = ["pure Gauss", "Farey-broad\nonly", "LW only\n(Test 9)", "joint\n(this fit)"] lls_models = [ll_pure_gauss, best_F_only[0], -25.85, ll_best] # -25.85 from Test 9 MLE dlls = [ll - ll_pure_gauss for ll in lls_models] colors = ["#2c3e50", "#9b59b6", "#c0392b", "#16a085"] bars = ax.bar(models, dlls, color=colors, alpha=0.85) for b, dll in zip(bars, dlls): ax.annotate(f"{dll:+.1f}", (b.get_x() + b.get_width()/2, b.get_height()), ha="center", va="bottom" if dll > 0 else "top", fontsize=9) ax.axhline(0, color="#999", lw=0.7) ax.set_ylabel(r"$\Delta$LL vs pure Gauss") ax.set_title("Model comparison: each component vs the joint mixture") fig.tight_layout() fig.savefig("c1_joint_model_compare.png", dpi=140) print("\nplots: c1_joint_fit.png c1_joint_a2_profile.png c1_joint_model_compare.png") if __name__ == "__main__": main()