Bayesian Inference for Neptune’s Orbit

This notebook maps the posterior distribution of Neptune’s orbital parameters using emcee’s affine-invariant ensemble sampler, in two complementary regimes:

• Synthetic benchmark (4-parameter): inject known Gaussian noise into synthetic Uranus longitudes and run MCMC with σ fixed at the injected value. This tests whether the posterior recovers the planted Neptune with correct uncertainties.
• JPL comparison (5-parameter): fit the simplified Neptune model to real JPL Horizons Uranus longitudes and treat σ as a free parameter. Here σ absorbs the mismatch between the 4-parameter model and the modern ephemeris, and the flat-vs-Bode prior comparison reveals how prior assumptions reshape the mass-distance ridge.

Each MCMC run targets the 50τ convergence heuristic. Convergence is assessed via ESS (effective sample size) and split-R̂ from ArviZ. Chains are checkpointed to data/cache/mcmc/*.h5 so interrupted runs can resume without resampling.

Setup

from IPython.display import display
from pathlib import Path
import logging

import matplotlib.pyplot as plt
import numpy as np

from discoverneptune.bayesian import run_mcmc_fixed_sigma, run_prior_comparison
from discoverneptune.data import PLANET_IDS, fetch_planet_vectors
from discoverneptune.sensitivity import inject_noise
from discoverneptune.simulation import (
    MASS_NEPTUNE_APPROX,
    NeptuneCandidate,
    StateVector,
    simulate_uranus_longitudes_from_vectors,
)
from discoverneptune.solver import solve
from discoverneptune.synthetic import build_synthetic_problem
from discoverneptune.plotting import (
    plot_corner,
    plot_prior_comparison as plot_prior_cmp,
    plot_residuals,
    plot_traces,
)

# plot_style imports added by transform script
from discoverneptune.plot_style import (
    apply_style, palette, category_cycle, dual_render,
    SEQUENTIAL_CMAP, truth_line, annotate_stat, fig_size,
)
apply_style()


logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s  %(levelname)-8s  %(name)s  %(message)s",
    datefmt="%H:%M:%S",
)
logger = logging.getLogger("bayesian_nb")

REPO_ROOT = Path.cwd()
if REPO_ROOT.name == "docs":
    REPO_ROOT = REPO_ROOT.parent

BAYESIAN_CACHE_DIR = REPO_ROOT / "data/cache/mcmc"
BAYESIAN_CACHE_DIR.mkdir(parents=True, exist_ok=True)
SYNTHETIC_POSTERIOR_CACHE = BAYESIAN_CACHE_DIR / "bayesian_synthetic_fixed_sigma.h5"
JPL_FLAT_POSTERIOR_CACHE = BAYESIAN_CACHE_DIR / "bayesian_jpl_flat.h5"
JPL_BODE_POSTERIOR_CACHE = BAYESIAN_CACHE_DIR / "bayesian_jpl_bode.h5"
MCMC_WORKER_OPTIONS = [1, 4, 6, 8]
MCMC_WORKERS = 4  # chosen from the local 1/4/6/8-worker benchmark on Apple Silicon

print(f"MCMC cache dir      : {BAYESIAN_CACHE_DIR}")
print(f"Synthetic cache     : {SYNTHETIC_POSTERIOR_CACHE}")
print(f"JPL flat cache      : {JPL_FLAT_POSTERIOR_CACHE}")
print(f"JPL Bode cache      : {JPL_BODE_POSTERIOR_CACHE}")
print(
    f"MCMC workers        : {MCMC_WORKERS} (process-based, candidate set {MCMC_WORKER_OPTIONS})"
)


def _sv_from_df(df) -> StateVector:
    return StateVector(
        x=float(df["x"].to_numpy()[0]),
        y=float(df["y"].to_numpy()[0]),
        z=float(df["z"].to_numpy()[0]),
        vx=float(df["vx"].to_numpy()[0]),
        vy=float(df["vy"].to_numpy()[0]),
        vz=float(df["vz"].to_numpy()[0]),
    )


def _rms_arcsec(truth: np.ndarray, sim: np.ndarray) -> float:
    residuals = ((truth - sim + np.pi) % (2 * np.pi)) - np.pi
    return float(np.degrees(np.sqrt(np.mean(residuals**2))) * 3600.0)


def _print_result_summary(result, *, label: str) -> None:
    print(label)
    for i, name in enumerate(result.labels):
        med = result.median_params[i]
        lo = result.ci_16[i]
        hi = result.ci_84[i]
        if name == "sigma_rad":
            med = np.degrees(med) * 3600.0
            lo = np.degrees(lo) * 3600.0
            hi = np.degrees(hi) * 3600.0
            print(f"  sigma_arcsec : {med:.3f}  [{lo:.3f}, {hi:.3f}]")
        else:
            print(f"  {name:12s}: {med:.6g}  [{lo:.6g}, {hi:.6g}]")
    print(f"  Acceptance   : {result.acceptance_fraction:.3f}")
    print(f"  Autocorr tau : {', '.join(f'{t:.0f}' for t in result.autocorr_times)}")
    print(f"  ESS          : {', '.join(f'{v:.0f}' for v in result.ess)}")
    print(f"  Split-R_hat  : {', '.join(f'{v:.3f}' for v in result.rhat)}")
    print(f"  Eff. samples : {len(result.flat_samples):,}")



def _print_convergence_summary(result, *, label: str) -> None:
    tau_max = float(np.max(result.autocorr_times))
    steps_per_tau = float(result.n_steps) / tau_max if tau_max > 0 else float("inf")
    max_rhat = float(np.max(result.rhat))
    min_ess = float(np.min(result.ess))
    passed_50tau = steps_per_tau > 50.0
    verdict = "PASS" if passed_50tau and max_rhat < 1.05 else "CHECK"

    print(label)
    print(f"  Max tau      : {tau_max:.1f}")
    print(f"  Steps / tau  : {steps_per_tau:.1f}x")
    print(f"  50τ target   : {'passed' if passed_50tau else 'not reached'}")
    print(f"  Max R_hat    : {max_rhat:.3f}")
    print(f"  Min ESS      : {min_ess:.0f}")
    print(f"  Verdict      : {verdict}")

MCMC cache dir      : <local-path>
Synthetic cache     : <local-path>
JPL flat cache      : <local-path>
JPL Bode cache      : <local-path>
MCMC workers        : 4 (process-based, candidate set [1, 4, 6, 8])

Synthetic benchmark posterior

Running MCMC on the noiseless synthetic benchmark would push σ toward the lower prior bound — the forward model can reproduce its own observations exactly, so there is no residual floor to infer. To make the Bayesian problem well-posed, we inject 2 arcsec Gaussian noise into the synthetic Uranus longitudes and hold σ fixed at that known value.

The posterior then answers a clean question: given known measurement noise and 66 annual observations spanning 1781–1846, how tightly can MCMC constrain the four Neptune parameters (mass, semi-major axis, eccentricity, mean longitude)?

TRUE_NEPTUNE = NeptuneCandidate(
    mass_solar=MASS_NEPTUNE_APPROX,
    a=30.07,
    e=0.009,
    l_rad=np.radians(132.0),
)
SIGMA_SYNTHETIC_ARCSEC = 2.0

logger.info("Building synthetic benchmark problem ...")
problem = build_synthetic_problem(TRUE_NEPTUNE)
synthetic_rng = np.random.default_rng(42)
synthetic_truth_clean = problem["truth_longs"]
synthetic_truth_longs = inject_noise(
    synthetic_truth_clean,
    sigma_arcsec=SIGMA_SYNTHETIC_ARCSEC,
    rng=synthetic_rng,
)

synthetic_mle = solve(
    jd_times=problem["jd_times"],
    truth_longs=synthetic_truth_longs,
    start_jd=problem["start_jd"],
    jupiter_sv=problem["jupiter_sv"],
    saturn_sv=problem["saturn_sv"],
    uranus_sv=problem["uranus_sv"],
    de_maxiter=30,
    de_popsize=5,
    seed=42,
)
synthetic_rms = np.degrees(
    np.sqrt(synthetic_mle.sse / len(problem["jd_times"]))
) * 3600.0

print(f"Synthetic observations: {len(problem['jd_times'])} annual epochs")
print(f"Injected sigma       : {SIGMA_SYNTHETIC_ARCSEC:.1f} arcsec")
print(
    f"MLE                  : a = {synthetic_mle.neptune.a:.3f} AU   "
    f"e = {synthetic_mle.neptune.e:.4f}   "
    f"l = {np.degrees(synthetic_mle.neptune.l_rad):.2f} deg"
)
print(f"Synthetic RMS        : {synthetic_rms:.2f} arcsec")

Synthetic observations: 66 annual epochs
Injected sigma       : 2.0 arcsec
MLE                  : a = 31.173 AU   e = 0.0126   l = 133.99 deg
Synthetic RMS        : 1.57 arcsec

# Residuals before and after Neptune — visual context for what the MCMC is fitting
neptune_absent = NeptuneCandidate(
    mass_solar=1e-12, a=30.07, e=0.009, l_rad=synthetic_mle.neptune.l_rad,
)
sim_baseline = simulate_uranus_longitudes_from_vectors(
    problem["jd_times"], problem["start_jd"],
    problem["jupiter_sv"], problem["saturn_sv"], problem["uranus_sv"],
    neptune_absent,
)
sim_best = simulate_uranus_longitudes_from_vectors(
    problem["jd_times"], problem["start_jd"],
    problem["jupiter_sv"], problem["saturn_sv"], problem["uranus_sv"],
    synthetic_mle.neptune,
)
display(dual_render(
    plot_residuals,
    problem["jd_times"], synthetic_truth_longs, sim_baseline, sim_best,
    alt="Uranus residuals before and after Neptune perturbation",
))

logger.info("Running fixed-sigma MCMC on noisy synthetic data ...")
synthetic_mcmc = run_mcmc_fixed_sigma(
    mle_result=synthetic_mle,
    jd_times=problem["jd_times"],
    truth_longs=synthetic_truth_longs,
    start_jd=problem["start_jd"],
    jupiter_sv=problem["jupiter_sv"],
    saturn_sv=problem["saturn_sv"],
    uranus_sv=problem["uranus_sv"],
    sigma_rad=np.radians(SIGMA_SYNTHETIC_ARCSEC / 3600.0),
    n_walkers=32,
    n_steps=40000,
    seed=42,
    backend_path=SYNTHETIC_POSTERIOR_CACHE,
    resume=True,
    n_workers=MCMC_WORKERS,
    progress=False,
)
_print_result_summary(synthetic_mcmc, label="Synthetic fixed-σ posterior")
_print_convergence_summary(synthetic_mcmc, label="Synthetic convergence summary")

synthetic_truths = [
    TRUE_NEPTUNE.mass_solar,
    TRUE_NEPTUNE.a,
    TRUE_NEPTUNE.e,
    TRUE_NEPTUNE.l_rad,
]
display(dual_render(plot_traces, synthetic_mcmc, alt="MCMC trace plots — synthetic fixed-sigma posterior"))
display(dual_render(plot_corner, synthetic_mcmc, truths=synthetic_truths, alt="Corner plot — synthetic fixed-sigma posterior"))

Synthetic fixed-σ posterior
  mass_solar  : 6.86549e-05  [5.1137e-05, 8.65407e-05]
  a           : 30.2095  [28.3883, 32.0713]
  e           : 0.0532538  [0.0152844, 0.102868]
  l_rad       : 2.28886  [2.20991, 2.3593]
  Acceptance   : 0.270
  Autocorr tau : 802, 650, 331, 568
  ESS          : 727, 1312, 3732, 1902
  Split-R_hat  : 1.032, 1.022, 1.010, 1.018
  Eff. samples : 7,424
Synthetic convergence summary
  Max tau      : 802.4
  Steps / tau  : 49.9x
  50τ target   : not reached
  Max R_hat    : 1.032
  Min ESS      : 727
  Verdict      : CHECK

JPL comparison posterior

When the simplified 4-parameter Neptune model is confronted with real JPL Horizons Uranus longitudes, the residual floor is no longer zero — it reflects genuine model-fidelity limitations (coplanar orbits, fixed perihelion, missing planets). Treating σ as a free fifth parameter lets the posterior absorb that mismatch floor, giving σ a physically meaningful interpretation rather than an arbitrary fixed value.

This is also the natural place to compare flat and Bode-style priors. The mass-distance degeneracy ridge is broad enough in the JPL regime that prior pressure visibly shifts the posterior — the Bode prior (centered near 38.8 AU) pulls the semi-major axis estimate higher, illustrating why Le Verrier’s distance was biased while his longitude prediction remained accurate.

logger.info("Loading cached JPL vectors ...")
df_uranus = fetch_planet_vectors(PLANET_IDS["uranus"], use_cache=True)
df_jupiter = fetch_planet_vectors(PLANET_IDS["jupiter"], use_cache=True)
df_saturn = fetch_planet_vectors(PLANET_IDS["saturn"], use_cache=True)

jpl_jd_times = df_uranus["datetime_jd"].to_numpy(dtype=float)
jpl_truth_longs = df_uranus["longitude_rad"].to_numpy(dtype=float)
jpl_start_jd = float(jpl_jd_times[0])
jpl_jupiter_sv = _sv_from_df(df_jupiter)
jpl_saturn_sv = _sv_from_df(df_saturn)
jpl_uranus_sv = _sv_from_df(df_uranus)

logger.info("Running JPL MLE baseline ...")
jpl_mle = solve(
    jd_times=jpl_jd_times,
    truth_longs=jpl_truth_longs,
    start_jd=jpl_start_jd,
    jupiter_sv=jpl_jupiter_sv,
    saturn_sv=jpl_saturn_sv,
    uranus_sv=jpl_uranus_sv,
    de_maxiter=30,
    de_popsize=5,
    seed=42,
)
jpl_mle_rms = np.degrees(np.sqrt(jpl_mle.sse / len(jpl_jd_times))) * 3600.0
print(
    f"JPL MLE              : a = {jpl_mle.neptune.a:.3f} AU   "
    f"e = {jpl_mle.neptune.e:.4f}   "
    f"l = {np.degrees(jpl_mle.neptune.l_rad):.2f} deg"
)
print(f"JPL RMS floor        : {jpl_mle_rms:.2f} arcsec")
print()

logger.info("Running JPL prior comparison MCMC ...")
jpl_comparison = run_prior_comparison(
    mle_result=jpl_mle,
    jd_times=jpl_jd_times,
    truth_longs=jpl_truth_longs,
    start_jd=jpl_start_jd,
    jupiter_sv=jpl_jupiter_sv,
    saturn_sv=jpl_saturn_sv,
    uranus_sv=jpl_uranus_sv,
    n_walkers=32,
    n_steps=26000,
    seed=42,
    flat_backend_path=JPL_FLAT_POSTERIOR_CACHE,
    bode_backend_path=JPL_BODE_POSTERIOR_CACHE,
    resume=True,
    n_workers=MCMC_WORKERS,
    progress=False,
)
jpl_flat = jpl_comparison.flat
jpl_bode = jpl_comparison.bode

_print_result_summary(jpl_flat, label="JPL posterior (flat prior, free σ)")
_print_convergence_summary(jpl_flat, label="JPL flat-prior convergence summary")
print()
_print_result_summary(jpl_bode, label="JPL posterior (Bode prior, free σ)")
_print_convergence_summary(jpl_bode, label="JPL Bode-prior convergence summary")

display(dual_render(plot_traces, jpl_flat, alt="MCMC trace plots — JPL flat prior"))
display(dual_render(plot_corner, jpl_flat, alt="Corner plot — JPL flat prior"))
display(dual_render(plot_prior_cmp, jpl_comparison, alt="Flat vs Bode prior comparison"))

JPL MLE              : a = 30.140 AU   e = 0.0645   l = 196.92 deg
JPL RMS floor        : 2.87 arcsec

JPL posterior (flat prior, free σ)
  mass_solar  : 0.000187401  [0.000173726, 0.000195895]
  a           : 48.3298  [46.5193, 49.4715]
  e           : 0.442375  [0.418453, 0.457724]
  l_rad       : 4.85359  [4.76827, 4.90676]
  sigma_arcsec : 1.919  [1.759, 2.107]
  Acceptance   : 0.444
  Autocorr tau : 467, 520, 515, 517, 147
  ESS          : 2554, 2132, 2181, 2156, 5798
  Split-R_hat  : 1.011, 1.012, 1.012, 1.012, 1.005
  Eff. samples : 10,912
JPL flat-prior convergence summary
  Max tau      : 519.8
  Steps / tau  : 50.0x
  50τ target   : passed
  Max R_hat    : 1.012
  Min ESS      : 2132
  Verdict      : PASS

JPL posterior (Bode prior, free σ)
  mass_solar  : 0.000181947  [0.00016395, 0.000193872]
  a           : 47.4092  [45.0809, 49.0072]
  e           : 0.430254  [0.39953, 0.45167]
  l_rad       : 4.81035  [4.69863, 4.88574]
  sigma_arcsec : 1.939  [1.775, 2.132]
  Acceptance   : 0.427
  Autocorr tau : 411, 475, 470, 471, 120
  ESS          : 2112, 1596, 1593, 1575, 4572
  Split-R_hat  : 1.014, 1.017, 1.017, 1.017, 1.006
  Eff. samples : 12,288
JPL Bode-prior convergence summary
  Max tau      : 474.6
  Steps / tau  : 50.6x
  50τ target   : passed
  Max R_hat    : 1.017
  Min ESS      : 1575
  Verdict      : PASS

# σ posterior marginal — the model mismatch floor
sigma_samples_arcsec = np.degrees(jpl_flat.flat_samples[:, -1]) * 3600.0
med_sigma = float(np.median(sigma_samples_arcsec))


def _build_sigma_hist(theme="light"):
    apply_style(theme)
    pal = palette(theme)
    cy = category_cycle(theme)
    fig, ax = plt.subplots(figsize=fig_size("wide"))
    ax.hist(sigma_samples_arcsec, bins=60, density=True, color=cy[2], alpha=0.7,
            edgecolor="white", linewidth=0.5)
    ax.axvline(jpl_mle_rms, color=cy[3], ls="--", lw=1.5,
               label=f'MLE RMS: {jpl_mle_rms:.2f}"')
    ax.axvline(med_sigma, color=cy[0], ls="-", lw=1.5,
               label=f'Posterior median: {med_sigma:.2f}"')
    ax.set_xlabel("σ (arcsec)")
    ax.set_ylabel("Density")
    ax.set_title("Inferred observation noise — JPL flat prior")
    ax.legend()
    fig.tight_layout()
    return fig


display(dual_render(_build_sigma_hist, alt="Sigma posterior marginal histogram — JPL flat prior"))

# Posterior predictive band — draw Neptune parameters from the posterior,
# simulate Uranus longitudes, and overlay the envelope on the observations.
n_draws = 100
rng = np.random.default_rng(0)
draw_indices = rng.choice(len(jpl_flat.flat_samples), size=n_draws, replace=False)

jpl_years = (jpl_jd_times - jpl_jd_times[0]) / 365.25 + 1781
residual_draws = np.empty((n_draws, len(jpl_jd_times)))

for i, idx in enumerate(draw_indices):
    sample = jpl_flat.flat_samples[idx]
    candidate = NeptuneCandidate(
        mass_solar=float(sample[0]),
        a=float(sample[1]),
        e=float(sample[2]),
        l_rad=float(sample[3]),
    )
    sim_longs = simulate_uranus_longitudes_from_vectors(
        jpl_jd_times, jpl_start_jd,
        jpl_jupiter_sv, jpl_saturn_sv, jpl_uranus_sv, candidate,
    )
    delta = ((jpl_truth_longs - sim_longs + np.pi) % (2 * np.pi)) - np.pi
    residual_draws[i] = np.degrees(delta) * 3600.0

q10 = np.percentile(residual_draws, 10, axis=0)
q50 = np.percentile(residual_draws, 50, axis=0)
q90 = np.percentile(residual_draws, 90, axis=0)

# Observed residuals at MLE for reference
mle_sim = simulate_uranus_longitudes_from_vectors(
    jpl_jd_times, jpl_start_jd,
    jpl_jupiter_sv, jpl_saturn_sv, jpl_uranus_sv, jpl_mle.neptune,
)
mle_delta = ((jpl_truth_longs - mle_sim + np.pi) % (2 * np.pi)) - np.pi
mle_residuals = np.degrees(mle_delta) * 3600.0


def _build_ppc_band(theme="light"):
    apply_style(theme)
    pal = palette(theme)
    cy = category_cycle(theme)
    fig, ax = plt.subplots(figsize=fig_size("double"))
    ax.fill_between(jpl_years, q10, q90, alpha=0.3, color=cy[0], label="80% posterior band")
    ax.plot(jpl_years, q50, lw=1.2, color=cy[0], label="Posterior median")
    ax.scatter(jpl_years, mle_residuals, s=12, color=cy[1], zorder=5, label="MLE residuals")
    truth_line(ax, 0, theme=theme)
    ax.set_xlabel("Year")
    ax.set_ylabel("Residual (arcsec)")
    ax.set_title("Posterior predictive residuals — JPL flat prior")
    ax.legend(loc="upper left")
    return fig


display(dual_render(_build_ppc_band, alt="Posterior predictive residuals with 80% band — JPL flat prior"))

Interpretation

The two regimes answer different questions by design:

• Synthetic benchmark: with σ fixed at the known 2 arcsec injection level, the posterior tests whether MCMC recovers the planted Neptune. The corner plot should show the true parameters inside the credible regions, and the mass–semi-major axis correlation reveals the degeneracy ridge even in this controlled setting.

• JPL comparison: with σ free, the posterior quantifies both Neptune’s orbital uncertainties and the model’s residual floor. The inferred σ (typically ~2.5 arcsec) matches the optimizer’s RMS residual, confirming that the likelihood is well-calibrated. The flat-vs-Bode prior comparison shows that prior pressure on semi-major axis shifts the mass estimate along the degeneracy ridge without meaningfully changing the longitude — the same qualitative finding as Le Verrier’s historical experience.

• Mass-distance degeneracy: the strong positive correlation between mass and semi-major axis appears in both regimes, but it is more scientifically interesting in the JPL case, where the ridge is broad enough for prior choice to matter. A more distant Neptune must be more massive to produce the same Uranus perturbations, so the two parameters trade off while longitude stays tightly constrained.