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 exposes how one noisy realization and the mass-distance-eccentricity degeneracy affect recovery; a single 68% interval need not contain every planted parameter.
JPL comparison (6-parameter): fit the simplified Neptune model to JPL Horizons Uranus longitudes with a stationary AR(1) discrepancy model. Its free σ and ρ describe model-conditional residual scale and correlation; 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
Code
import loggingfrom pathlib import Pathimport matplotlib.pyplot as pltimport numpy as npfrom IPython.display import displayfrom discoverneptune.bayesian import ( AR1_MCMC_BOUNDS, AR1_MCMC_PARAM_NAMES, run_mcmc_fixed_sigma, run_prior_comparison,)from discoverneptune.data import PLANET_IDS, fetch_planet_vectorsfrom discoverneptune.objective import BOUNDSfrom discoverneptune.plot_style import ( apply_style, category_cycle, dual_render, fig_size, palette, truth_line,)from discoverneptune.plotting import ( plot_corner, plot_residuals, plot_traces,)from discoverneptune.plotting import ( plot_prior_comparison as plot_prior_cmp,)from discoverneptune.sensitivity import inject_noisefrom discoverneptune.simulation import ( MASS_NEPTUNE_APPROX, NeptuneCandidate, StateVector, simulate_uranus_longitudes_from_vectors,)from discoverneptune.solver import parameters_near_bounds, solvefrom discoverneptune.synthetic import build_synthetic_problemapply_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.parentBAYESIAN_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 Siliconprint(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.pireturnfloat(np.degrees(np.sqrt(np.mean(residuals**2))) *3600.0)def _print_result_summary(result, *, label: str) ->None:print(label)for i, name inenumerate(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.0print(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 >0elsefloat("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.05else"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}")
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)?
Synthetic observations: 66 annual epochs
Injected sigma : 2.0 arcsec
MLE : a = 26.167 AU e = 0.1399 l = 119.53 deg
Synthetic RMS : 1.67 arcsec
Code
# Residuals before and after Neptune — visual context for what the MCMC is fittingneptune_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",))
When the simplified 4-parameter Neptune model is confronted with JPL Horizons Uranus longitudes, its residuals include structured model discrepancy from simplifications such as coplanar orbits, fixed perihelion, and missing planets. The JPL fit therefore uses a stationary AR(1) residual model with free marginal scale σ and correlation ρ. These are model-conditional discrepancy parameters, not a measurement of an irreducible physical noise floor.
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.
The two regimes answer different questions by design:
Synthetic benchmark: with σ fixed at the known 2 arcsec injection level, the posterior is a controlled calibration example, not a guarantee that every planted value lies in a 68% interval. In this committed realization the mass and eccentricity truths fall just below their 16th percentiles, while the MLE eccentricity is high (0.1399 versus the 0.15 cap) but not within the diagnostic’s 1% boundary tolerance. The notebook reports both boundary checks and the interval misses explicitly; the mass–semi-major axis correlation reveals the degeneracy ridge that makes the noisy inverse problem weakly identified.
JPL comparison: the six-parameter fit uses a stationary AR(1) discrepancy model because the independent-Gaussian residuals are strongly autocorrelated. Its intervals remain model-conditional: convergence, boundary sensitivity, SBC, and posterior-predictive checks determine whether they are calibrated. The reported rho is near the upper bound, and modern Neptune lies outside the displayed orbital intervals. This is evidence of simplified-model discrepancy, not a physical recovery of Neptune’s orbit. The optimizer RMS is reported only within its stated search domain, not as an irreducible floor.
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.