Validation: Detectability, Hindcast, SBC, and Posterior Predictive Checks
This notebook asks four questions about the Neptune inverse problem:
Detectability — before we validate any particular Neptune fit, is the Uranus anomaly itself large enough and structured enough to rule out noise? This is the precondition that motivates solving an inverse problem at all.
Hindcast validation — fit on 1781..cutoff-1, predict the cutoff..1846 holdout. This is the most important scientific check: does an early Neptune fit actually generalize out of sample?
Simulation-based calibration (SBC) — draw synthetic problems from a narrow prior, fit each one, and check whether posterior intervals behave sensibly. The checked-in run is intentionally lightweight and demonstrative, not statistically rigorous.
Posterior predictive checks — draw posterior samples, replicate the Uranus longitude series, and compare the implied residual structure against what was observed.
Runtime notes for a fresh run: - detectability: under a second - hindcast: about 30 seconds - SBC: a few minutes - posterior predictive (synthetic): under a minute - posterior predictive (JPL): optional, cache-guarded
Setup
The notebook works from either the repository root or the docs/ directory. We keep cache artifacts under data/cache/validation/ and rebuild the lightweight trial objects from cached parquet files when present.
Code
import loggingfrom pathlib import Pathimport emceeimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdfrom IPython.display import displayfrom discoverneptune.bayesian import run_mcmc_fixed_sigmafrom discoverneptune.data import PLANET_IDS, fetch_planet_vectorsfrom discoverneptune.plot_style import ( apply_style, dual_render, fig_size, palette, truth_line,)from discoverneptune.plotting import ( plot_hindcast_1846_error, plot_hindcast_rms, plot_posterior_predictive_bands, plot_posterior_predictive_summaries, plot_sbc_coverage,)from discoverneptune.simulation import ( MASS_NEPTUNE_APPROX, NeptuneCandidate, StateVector,)from discoverneptune.solver import solvefrom discoverneptune.synthetic import build_synthetic_problemfrom discoverneptune.validation import ( HindcastTrial, SBCTrial, hindcast_trials_to_frame, run_hindcast_experiment, run_posterior_predictive_check, run_sbc_experiment, sbc_trials_to_frame, summarize_sbc_trials,)apply_style()logging.basicConfig( level=logging.INFO,format="%(asctime)s%(levelname)-8s%(name)s%(message)s", datefmt="%H:%M:%S",)logger = logging.getLogger("validation_nb")cwd = Path.cwd()candidate_repo_roots = [cwd, cwd.parent]repo_root =next( (candidate for candidate in candidate_repo_roots if (candidate /"pyproject.toml").exists()),None,)if repo_root isNone:raiseFileNotFoundError("Could not find the repository root from the current working directory." )CACHE_DIR = repo_root /"data"/"cache"/"validation"CACHE_DIR.mkdir(parents=True, exist_ok=True)def _sv_from_df(df: pd.DataFrame) -> 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 _hindcast_trial_from_row(row: pd.Series) -> HindcastTrial:return HindcastTrial( mode=str(row["mode"]), cutoff_year=int(row["cutoff_year"]), n_fit=int(row["n_fit"]), n_holdout=int(row["n_holdout"]), in_sample_rms_arcsec=float(row["in_sample_rms_arcsec"]), holdout_rms_arcsec=float(row["holdout_rms_arcsec"]), holdout_max_abs_arcsec=float(row["holdout_max_abs_arcsec"]), longitude_1846_error_deg=float(row["longitude_1846_error_deg"]), best_params=np.array( [row["best_mass"], row["best_a"], row["best_e"], row["best_l_rad"]], dtype=float, ), best_sse=float(row["best_sse"]), converged=bool(row["converged"]), )def _sbc_trial_from_row(row: pd.Series) -> SBCTrial: residual_model =str(row.get("residual_model", "iid")) params = ["mass", "a", "e", "l_rad"]if residual_model =="ar1": params += ["sigma_rad", "rho"]return SBCTrial( noise_arcsec=float(row["noise_arcsec"]), truth_params=np.array([row[f"truth_{name}"] for name in params], dtype=float), mle_params=np.array( [row[f"mle_{name}"] for name in params[:4]], dtype=float ), converged=bool(row["converged"]), rank_by_param={name: int(row[f"rank_{name}"]) for name in params}, coverage_by_level={ name: {0.5: bool(row[f"cov_{name}_50"]),0.8: bool(row[f"cov_{name}_80"]),0.95: bool(row[f"cov_{name}_95"]), }for name in params }, interval_width_by_level={ name: {0.5: float(row[f"width_{name}_50"]),0.8: float(row[f"width_{name}_80"]),0.95: float(row[f"width_{name}_95"]), }for name in params }, max_rhat=float(row["max_rhat"]), min_ess=float(row["min_ess"]), residual_model=residual_model, )
Detectability: Is the Anomaly Real?
Before validating any particular Neptune fit, we should verify the Uranus anomaly is statistically distinguishable from noise. If it were not, our whole inverse problem would be fitting sky-position noise.
We simulate Uranus’s longitude under a no-Neptune 4-body model (Sun + Jupiter + Saturn + Uranus, all pinned at JPL truth at 1781-01-01), compare to the bundled annual JPL-derived proxy longitudes, and ask two questions of the residual sequence:
Is the RMS larger than plausible measurement noise? We score \(\chi^2_\mathrm{red} = \sum r_i^2 / (\sigma^2 (N - k))\) with \(\sigma = 1\,\mathrm{arcsec}\) as an illustrative reference scale and \(k = 0\) free parameters. This is not a likelihood for the historical measurement process. A value near 1 would mean the no-Neptune model is consistent with noise; values far above 1 rule it out.
Are the residuals structured, not random? Lag-1 autocorrelation is near 0 for white noise (std \(\approx 1/\sqrt{N} \approx 0.12\) for \(N=66\)). A value near \(+1\) means a smoothly drifting signal — the fingerprint of a missing perturbing body.
No-Neptune residuals over 66 JPL-derived proxy points (1781–1846):
RMS = 31.5 arcsec
Peak-to-peak = 131.4 arcsec
chi2_red (sigma=1 arcsec) = 994 (expected ~1 under noise-only)
Lag-1 autocorrelation = +0.900 (expected ~0 +/- 0.12 under white noise)
Hindcast Validation
We fit the simplified Neptune model on an early prefix of the synthetic observation window and then score its predictions on the held-out suffix. The key question is how quickly the holdout RMS and 1846-position error improve as more years accumulate.
SBC repeatedly simulates synthetic problems from a narrow prior near the true Neptune regime, reruns the usual solve + MCMC pipeline, and checks whether posterior intervals cover the truth at roughly the advertised rates. The notebook run below is lightweight: it uses one fixed IID noise level plus a two-replicate AR(1) path check. These surfaces exercise both calibration paths but are not statistically definitive.
Posterior predictive checks ask what the fitted model says Uranus residuals should look like under repeated draws from the posterior. In synthetic mode the residual band should stay narrow and centered near zero; in JPL mode the same machinery exposes the residual structure that the simplified 4-parameter model cannot capture.
This cell reuses the flat-prior AR(1) HDF posterior produced by bayesian.ipynb when that cache is available. It uses the usual JPL Uranus/Jupiter/Saturn state-vector path, then feeds the thinned six-parameter posterior draws through the same predictive machinery as the synthetic case. The output states explicitly when the cache is absent.
Detectability proxy. Relative to the illustrative 1-arcsecond scale, the no-Neptune JPL-derived proxy residuals have a large chi2_red and lag-1 autocorrelation near +1. This demonstrates that omitting Neptune creates a strong, structured signal in the computational benchmark. It is not evidence about the original historical measurement errors; those require the separately sourced historical residual series.
Hindcast. The interesting quantity is the holdout error: if a fit learned a real Neptune-like signal, holdout RMS and the 1846-position error should improve as later cutoff years include more of the synthetic benchmark interval.
SBC. With only a handful of IID replicates and two AR(1) path-check replicates, empirical coverage is noisy. The point is to execute both calibration mechanisms and produce reusable result surfaces, not to claim definitive frequentist coverage.
Posterior predictive. The upper-tail probability reports P(T_rep >= T_obs). The committed synthetic RMS value of 1.0 means every replicated RMS exceeded its paired observed RMS. Here the fitted residual RMS is 1.67 arcsec against a fixed 2 arcsec replication scale, so the extreme value can reflect conditioning on an over-fit point estimate; it is not confirmation of calibration or proof of model over-dispersion. The synthetic AR(1) result is explicitly an implementation path check with a fixed rho, not a posterior inference. The cached JPL run did execute here: its lag-1 upper-tail probability is 0.00 and its absolute-drift probability is 0.07, so even the fitted AR(1) discrepancy fails to reproduce important residual structure. Those low tail probabilities reinforce, rather than resolve, the boundary-sensitive model-mismatch warning. JPL results should be interpreted only when the matching cached Bayesian posterior was loaded and the cell ran.