.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/generated/scripts_synth_data/finding_neptune_via_bayesian_inv.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_generated_scripts_synth_data_finding_neptune_via_bayesian_inv.py: Finding Neptune using Uranus ============================ .. GENERATED FROM PYTHON SOURCE LINES 9-22 Note : - For better understanding of this notebook, refer to the `md_file `__ designed specifically for this notebook and better insights into the theory. - The import methods and functions from `neptune_bayesian_methods `__, `neptune_deterministic_methods `__ and `setup_inversion `__ are used to set up the simulation and perform the necessary calculations. .. GENERATED FROM PYTHON SOURCE LINES 22-31 .. code-block:: Python # This notebook requires the following libraries to run, in order to install them uncomment the lines below # %pip install cofi # %pip install numba # %pip install tqdm # %pip install matplotlib # %pip install astroquery .. GENERATED FROM PYTHON SOURCE LINES 36-39 1. Introduction --------------- .. GENERATED FROM PYTHON SOURCE LINES 42-129 - The following Notebook is based on the historical problem on how Neptune was found by Johann Galle using mathematical predictions made independently by two astronomers: - Urbain Le Verrier (France) - John Couch Adams (England) Through this Notebook we wish to demostrate how ``CoFI`` can be used to solve this problem via deterministic inversion. For more details on this problem, see the following `thesis `__ - In the following notebook we discuss the problem of finding Neptune’s mass, its velocity components and its position coordinates in the year 1775, by modeling the trajectory of Uranus with and without the influence of Neptune. - We define $ g(m) $, our forward model, as vector-valued function that predicts the position coordinates of Uranus at each observation time :math:`t_j`, as a function of Neptune’s parameters :math:`m`: | .. math:: g(m) = \begin{bmatrix} \hat x_1(m) \\ \vdots \\ \hat x_N(m) \\ \hat y_1(m) \\ \vdots \\ \hat y_N(m) \\ \hat z_1(m) \\ \vdots \\ \hat z_N(m) \end{bmatrix} \in \mathbb{R}^{3M \times 1} where $ N $ is the number of data points, and $ :raw-latex:`\hat `x_j(m),  :raw-latex:`\hat `y_j(m),  :raw-latex:`\hat `z_j(m) $ are the coordinates of Uranus at data point $ j  =  1,  2,  ….  N$ as a function of Neptune’s parameters $ m $, | where :math:`m = (m_M, m_x, m_y, m_z, {m_{v_x}}, m_{v_y}, m_{v_z})` is the set of parameters describing Neptune’s mass (:math:`m_M`), its position coordinates :math:`(m_x, m_y, m_z)` and its velocity components :math:`(m_{v_x}, m_{v_y}, m_{v_z})` | and :math:`d` as the data vector of positions of Uranus at different time steps: | .. math:: d = \begin{bmatrix} x_1 \\ \vdots \\ x_N \\ y_1 \\ \vdots \\ y_N \\ z_1 \\ \vdots \\ z_N \end{bmatrix} \in \mathbb{R}^{3M \times 1} where $ N $ is the number of data points, and $ :raw-latex:`\hat `x_j,  :raw-latex:`\hat `y_j,  :raw-latex:`\hat `z_j $ are the true coordinates of Uranus at data point $ j  =  1,  2,  ….  N$. - hence our problem formulation changes to : .. math:: \underset{m}{\min} || g(m) - {d} ||_{2}^2 .. GENERATED FROM PYTHON SOURCE LINES 132-135 2. Problem Setting ------------------ .. GENERATED FROM PYTHON SOURCE LINES 138-145 - We use **Newton’s Law of gravitation** to compute the gravitational force acting on a planet due to other celestial bodies. - This formulation uses Newton’s Law of Universal Gravitation to model the **net gravitational influence** from multiple bodies on a single target planet. .. GENERATED FROM PYTHON SOURCE LINES 148-197 - Let :math:`r` be the vector containing the positions of all planets and :math:`a` be the vector containing the accelerations of all planets. The velocity and acceleration vectors are defined as: .. math:: \dot{\mathbf{r}} = \mathbf{v}, \quad \ddot{\mathbf{r}} = \mathbf{a} where .. math:: \mathbf{r} = \begin{bmatrix} \mathbf{r}_1 \\ \vdots \\ \mathbf{r}_9 \end{bmatrix}, \quad \mathbf{a} = \begin{bmatrix} \mathbf{a}_1 \\ \vdots \\ \mathbf{a}_9 \end{bmatrix}, \quad \mathbf{r}, \mathbf{a} \in \mathbb{R}^{27 \times 1} - This results in the system of differential equations: .. math:: \frac{d}{dt} \begin{bmatrix} \mathbf{r}(t) \\ \mathbf{v}(t) \end{bmatrix} = \begin{bmatrix} \mathbf{v}(t) \\ \mathbf{a}(t) \end{bmatrix} .. GENERATED FROM PYTHON SOURCE LINES 197-226 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from tqdm import tqdm from typing import Callable import copy from numba import njit, jit import warnings import multiprocessing as mp import time import arviz from astroquery.jplhorizons import Horizons warnings.filterwarnings('ignore') from cofi import BaseProblem, InversionOptions, Inversion from neptune_mcmc_helpers import ( analysis_samples, load_mcmc_results, plot_pair_kde_legacy_style, run_cofi_emcee, safe_autocorr_time, save_mcmc_results, scaled_reference_values, selected_trajectory_samples, ) np.random.seed(42) .. GENERATED FROM PYTHON SOURCE LINES 228-244 .. code-block:: Python # Long-run control # Keep RUN_LONG_MCMC = False for quick validation. Set it to True when you want # to regenerate converged MCMC outputs and save them beside this notebook. from pathlib import Path RUN_LONG_MCMC = True USE_SAVED_LONG_RUN_RESULTS = False SAVE_LONG_RUN_RESULTS = False VALIDATION_STEPS = 40 SYNTHETIC_LONG_STEPS = 10_000 REAL_LONG_STEPS = 50_000 LONG_RUN_RESULTS_FILE = Path("neptune_long_run_results.npz") .. GENERATED FROM PYTHON SOURCE LINES 249-253 - We solve our ODEs with the **Runge-Kutta 4 (RK4)** method, which is an explicit and iterative method, well-suited for initial value problems. .. GENERATED FROM PYTHON SOURCE LINES 256-259 - In the following cell, we import - ``acceleration`` and ``rk4_step`` to serve as functions for our forward model. .. GENERATED FROM PYTHON SOURCE LINES 262-271 - We now demostrate our forward model, using the above defined functions, in the ``run_simulation`` method defined below, which helps us run a simulation of solar system. - Throughout this notebook for the purpose of our inversion, we are going to define mass in terms of solar masses, positions coordinates in **Astonomical Units (AU)** and velocities for planets in **Au/day** .. GENERATED FROM PYTHON SOURCE LINES 271-276 .. code-block:: Python from neptune_deterministic_methods import acceleration, rk4_step, run_simulation trajectories = run_simulation(T = 100, dt = 1, plot_only=['Uranus', 'Neptune', 'Saturn', 'Jupiter', 'Mars', 'Earth', 'Venus', 'Mercury']) .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_001.png :alt: Planetary Orbits Around the Sun :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 281-327 2.1 Bayesian Inversion via the Stretch Move ------------------------------------------- In this notebook, we use ``CoFI`` to perform Bayesian inversion. Specifically, we use the **stretch move** ensemble sampler introduced by goodman & Weare (2010), which is well-suited for exploring complex, high-dimensional posterior distributions. MCMC is a method to sample from the **posterior distribution** of model parameters :math:`\mathbf{m}` given observed data :math:`\mathbf{d}`. According to Bayes’ theorem: .. math:: \tag{1} p(\mathbf{m} \mid \mathbf{d}) \propto p(\mathbf{d} \mid \mathbf{m}) \, p(\mathbf{m}), where: - :math:`p(\mathbf{m} \mid \mathbf{d})` is the **posterior** probability of the model, - :math:`p(\mathbf{d} \mid \mathbf{m})` is the **likelihood**, and - :math:`p(\mathbf{m})` is the **prior**. The likelihood measures how well a model explains the data. In the stretch move algorithm, a population of walkers explores the posterior distribution simultaneously. Each walker proposes new positions based on the positions of other walkers, allowing for more efficient exploration of complex posteriors. The sampling proceeds in the following steps: 1. Initialize a set of walkers in parameter space. 2. At each iteration, propose new positions for each walker using the stretch move. 3. Accept or reject the proposed positions based on the Metropolis criterion. 4. After a **burn-in** period, the samples from all walkers are collected to approximate the posterior. This method yields not only the most likely model but also the full distribution over model space, capturing uncertainty in the inferred parameters. .. GENERATED FROM PYTHON SOURCE LINES 330-341 3. Inversion on Synthetic Data ------------------------------ - We will first demonstrate bayesian inversion using ``CoFI`` on synthetic data. - The synthetic observations are generated by integrating our gravitational forward model with a fourth-order Runge-Kutta (``RK4``) solver to simulate Uranus’s trajectory under the influence of Neptune. .. GENERATED FROM PYTHON SOURCE LINES 344-372 - We simulate observational noise by sampling from zero-mean Gaussian distributions with specified variances for each coordinate: .. math:: x_\text{obs} = x + \epsilon_x, \quad y_\text{obs} = y + \epsilon_y, \quad z_\text{obs} = z + \epsilon_z - where - .. math:: \epsilon_x \sim \mathcal{N}(0, \sigma_x^2), \quad \epsilon_y \sim \mathcal{N}(0, \sigma_y^2), \quad \epsilon_z \sim \mathcal{N}(0, \sigma_z^2) - with noise levels set as - .. math:: \sigma_x = \sigma_y = 10^{-3}, \quad \sigma_z = 10^{-5} .. GENERATED FROM PYTHON SOURCE LINES 375-378 - The function below generates the synthetic data with the specified noise levels. .. GENERATED FROM PYTHON SOURCE LINES 381-384 3.1 Generating synthetic data ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 384-398 .. code-block:: Python from neptune_deterministic_methods import generate_synthetic_data T = 190 # time for which we want to generate synthetic data z_scale_factor = 1 dt = 1 U_true = generate_synthetic_data(T = T, dt = dt, z_scaling = False, add_noise = True, noise_level = np.array([0.001, 0.001, 0.00001])) .. rst-class:: sphx-glr-script-out .. code-block:: none Synthetic data ranges: X: -18.298 to 20.084 (range: 38.383) Y: -19.319 to 19.019 (range: 38.337) Z: -0.272199 to 0.248322 (range: 0.520521) Noisy synthetic data ranges: X: -18.299 to 20.084 (range: 38.383) Y: -19.321 to 19.020 (range: 38.340) Z: -0.272212 to 0.248322 (range: 0.520534) .. GENERATED FROM PYTHON SOURCE LINES 403-412 - The cell below sets up the starting model for our walkers/chains in MCMC and some pre-defined scales to be used for scaling while running our inversion, to ensure that all parameters are roughly on the same scale or at least near to each other. - For our bayesian inversion we first scale all the parameters and the convert our mass to log mass. This ensures that all our parameters are in the same scale. .. GENERATED FROM PYTHON SOURCE LINES 412-425 .. code-block:: Python # True/reference parameters for Neptune [mass, x, y, z, vx, vy, vz] from setup_inversion import get_inversion_indices, set_true_m, get_param_bounds, get_param_scales, set_initial_conditions, get_starting_points, PARAM_NAMES m_0 = set_true_m() initial_conditions = set_initial_conditions() PARAM_BOUNDS = get_param_bounds() PARAM_SCALES = get_param_scales() INVERT_INDICES = get_inversion_indices() STARTING_POINTS = get_starting_points() .. GENERATED FROM PYTHON SOURCE LINES 427-433 .. code-block:: Python from setup_inversion import scale_param, unscale_param, validate_config validate_config() .. rst-class:: sphx-glr-script-out .. code-block:: none Inverting for: ['mass', 'x', 'y', 'z', 'vx', 'vy', 'vz'] Starting points: [5.129494593915012e-05, -29.91585533802299, 3.145551080953517, 0.1050744172524112, -0.0002804329678559355, -0.00311021549555015, 4.187313350941659e-05] .. GENERATED FROM PYTHON SOURCE LINES 435-441 .. code-block:: Python names = list(initial_conditions.keys()) n_bodies = len(names) uranus_idx = names.index("Uranus") .. GENERATED FROM PYTHON SOURCE LINES 443-447 .. code-block:: Python from neptune_bayesian_methods import predict_U .. GENERATED FROM PYTHON SOURCE LINES 449-485 .. code-block:: Python m_start_scaled = scale_param(m_0) if 0 in INVERT_INDICES: m_start_scaled[0] = np.log10(m_start_scaled[0]) # Convert mass to log scale if len(INVERT_INDICES) == 1: m_start_scaled = m_start_scaled.item() if hasattr(m_start_scaled, 'item') else m_start_scaled print(f"\nStarting points (unscaled): {STARTING_POINTS}") print(f"Starting points (scaled): {m_start_scaled}") print("\nTesting forward function...") try: pred_test = predict_U(m_start_scaled, T = T, dt = dt, z_scale_factor = z_scale_factor) residual_test = pred_test - U_true print(f"Initial residual norm: {np.linalg.norm(residual_test):.6f}") print(f"Residual by component:") print(f" X component: {np.linalg.norm(residual_test[:T]):.6f}") print(f" Y component: {np.linalg.norm(residual_test[T:2*T]):.6f}") print(f" Z component: {np.linalg.norm(residual_test[2*T:]):.6f}") # Match the archived notebook's pre-MCMC state exactly. In the archive, # predict_U/build_neptune_vector mutated m_start_scaled[0] in place from # log10(scaled_mass) back to scaled_mass during this forward-function test. # The current predictor copies inputs, so we reproduce that side effect # explicitly before constructing the emcee walkers. if 0 in INVERT_INDICES: m_start_scaled[0] = 10**m_start_scaled[0] print(f"Archive-compatible m_start_scaled before MCMC: {m_start_scaled}") except Exception as e: print(f"Forward function test failed: {e}") import traceback traceback.print_exc() .. rst-class:: sphx-glr-script-out .. code-block:: none Starting points (unscaled): [5.129494593915012e-05, -29.91585533802299, 3.145551080953517, 0.1050744172524112, -0.0002804329678559355, -0.00311021549555015, 4.187313350941659e-05] Starting points (scaled): [-3.73377762 -4.61564393 0.42847056 0.40820822 -0.27220634 -3.73292208 0.00484093] Testing forward function... Initial residual norm: 0.019639 Residual by component: X component: 0.013658 Y component: 0.014112 Z component: 0.000148 Archive-compatible m_start_scaled before MCMC: [ 1.84596041e-04 -4.61564393e+00 4.28470556e-01 4.08208224e-01 -2.72206336e-01 -3.73292208e+00 4.84093052e-03] .. GENERATED FROM PYTHON SOURCE LINES 490-501 3.2 Running the Inversion on Synthetic Data ------------------------------------------- The Bayesian setup below keeps the original likelihood and prior formulation: a bounded uniform prior on the scaled parameters and a diagonal Gaussian data-error model, implemented as ``-0.5 * sum((d - g(m))**2 * inv_variances)``. The notebook defaults to a short validation run. Set ``RUN_LONG_MCMC = True`` in the configuration cell to run the original 10,000-step synthetic-data chain and save converged upload outputs into ``neptune_long_run_results.npz``. .. GENERATED FROM PYTHON SOURCE LINES 504-509 - Defined below are our ``log prior`` and ``log likelihood`` functions. - We set a uniform prior on log of our mass and all the other parameters (Neptune’s position coordinates and velocity components). .. GENERATED FROM PYTHON SOURCE LINES 509-527 .. code-block:: Python param_bounds_lower = np.array([bound[0] for bound in PARAM_BOUNDS]) param_bounds_upper = np.array([bound[1] for bound in PARAM_BOUNDS]) bounds_lower_scaled = scale_param(param_bounds_lower) bounds_upper_scaled = scale_param(param_bounds_upper) bounds_lower_scaled[0] = np.log10(bounds_lower_scaled[0]) # Convert mass back to log scale bounds_upper_scaled[0] = np.log10(bounds_upper_scaled[0]) # Convert mass back to log scale sigma_x = 8e-3 # Standard deviation for X position in uranus orbit data, taken one order higher than the noise added sigma_y = 8e-3 # Standard deviation for Y position in uranus orbit data, taken one order higher than the noise added sigma_z = 1e-4 # Standard deviation for Z position in uranus orbit data, taken two orders higher than the noise added diagonal_values = [1/sigma_x**2]*T + [1/sigma_y**2]*T + [1/sigma_z**2]*T inv_variances = np.array(diagonal_values) .. GENERATED FROM PYTHON SOURCE LINES 529-583 .. code-block:: Python def my_log_likelihood(m_scaled : np.ndarray, U_true : np.ndarray) -> float: """ Log-likelihood function for Bayesian inversion Parameters ---------- m_scaled : np.array Scaled parameters being inverted. U_true : np.array True Uranus positions in the format [x, y, z] scaled by z_scale_factor. Returns ------- log_likelihood : float The log-likelihood value, or -inf if the prediction fails or residuals are not finite. """ m_scaled = np.atleast_1d(m_scaled).astype(np.float64) try: y_hat = predict_U(m_scaled, T = T, dt = dt, z_scale_factor = z_scale_factor) if not np.isfinite(y_hat).all(): return -np.inf residuals = U_true - y_hat if residuals.shape != U_true.shape: return -np.inf return -0.5 * np.sum(residuals**2 * inv_variances) except (ValueError, ArithmeticError): return -np.inf except Exception as e: print(f"Error in log-likelihood calculation: {e}") return -np.inf def my_log_prior(m_scaled : np.ndarray) -> float: """ Log-prior function for Bayesian inversion Parameters ---------- m_scaled : np.ndarray Scaled parameters being inverted. Returns ------- float The log prior value, or -inf if the parameters are out of bounds. """ m_scaled = np.atleast_1d(m_scaled) if np.any(m_scaled < bounds_lower_scaled) or np.any(m_scaled > bounds_upper_scaled): return -np.inf return 0 .. GENERATED FROM PYTHON SOURCE LINES 585-698 .. code-block:: Python n_walkers = 20 n_dim = len(INVERT_INDICES) nsteps = SYNTHETIC_LONG_STEPS if RUN_LONG_MCMC else VALIDATION_STEPS use_pool = True show_progress = RUN_LONG_MCMC run_mode = "long" if RUN_LONG_MCMC else "validation" print(f"Optimized setup: {len(inv_variances)} inverse variance values") print(f"Running {run_mode} MCMC chain: {nsteps} steps, {n_walkers} walkers") print("Set RUN_LONG_MCMC = True in the configuration cell to run the original 10,000-step synthetic chain.") print("Parameter bounds (unscaled):") print(" Lower:", param_bounds_lower) print(" Upper:", param_bounds_upper) print("Scaled bounds:") print(" Lower:", bounds_lower_scaled) print(" Upper:", bounds_upper_scaled) print(f"Number of dimensions: {n_dim}") # Reset the MCMC starting model to the scaled/log parameterization. # The archive-compatibility forward test above intentionally mutates # m_start_scaled[0] from log10(scaled_mass) back to scaled_mass. That # state is useful for apples-to-apples archive checks, but it is wrong # for tight walker initialization because it clips mass to the upper bound. m_start_scaled = scale_param(m_0).astype(float) if 0 in INVERT_INDICES: m_start_scaled[0] = np.log10(m_start_scaled[0]) print("MCMC start (scaled/log):", m_start_scaled) # starting_lower_bounds = m_start_scaled - 10 # starting_upper_bounds = m_start_scaled + 10 # walkers_start = np.random.uniform( # starting_lower_bounds, # starting_upper_bounds, # size=(n_walkers, n_dim) # ) rng = np.random.default_rng(42) walker_scales = 0.075 * (bounds_upper_scaled - bounds_lower_scaled) # Keep mass a bit tighter if desired walker_scales[0] = 0.075 walkers_start = m_start_scaled + rng.normal( scale=walker_scales, size=(n_walkers, n_dim), ) walkers_start = np.clip(walkers_start, bounds_lower_scaled, bounds_upper_scaled) print("Initial state shape:", walkers_start.shape) print( "Initial walker mass min/max/std:", f"{walkers_start[:, 0].min():.6g}", f"{walkers_start[:, 0].max():.6g}", f"{walkers_start[:, 0].std():.6g}", ) for i, w in enumerate(walkers_start): ll = my_log_likelihood(w, U_true) if not np.isfinite(ll): print(f"Invalid log-likelihood at walker {i}: {w}, ll = {ll}") for i, w in enumerate(walkers_start): walkers_start[i] = np.clip(w, bounds_lower_scaled, bounds_upper_scaled) if not np.array_equal(w, walkers_start[i]): print(f"Walker {i} clipped from {w} to {walkers_start[i]}") print("All walkers are now within bounds.") for i, w in enumerate(walkers_start): if not np.isfinite(my_log_prior(w)): print(f"ERROR: Walker {i} still outside prior bounds: {w}") else: ll = my_log_likelihood(w, U_true) if not np.isfinite(ll): print(f"WARNING: Walker {i} has invalid log-likelihood: {ll}") inv_problem = BaseProblem() inv_problem.name = "Neptune Orbit Determination - Config Driven" inv_problem.set_data(U_true) inv_problem.set_forward(predict_U) inv_problem.set_initial_model(np.atleast_1d(m_start_scaled)) inv_problem.set_log_prior(my_log_prior) inv_problem.set_log_likelihood(my_log_likelihood, args=(U_true,)) inv_problem.set_model_shape(n_dim) import emcee stretch_move = emcee.moves.StretchMove(a=1.2) inv_result, elapsed = run_cofi_emcee( inv_problem, walkers_start, n_walkers, nsteps, InversionOptions, Inversion, progress=show_progress, use_pool=use_pool, n_threads=10, moves=stretch_move, skip_initial_state_check=True, ) print(f"OPTIMIZED Inversion completed in {elapsed:.2f} seconds using `emcee`.") print("The inversion result from `emcee`:") inv_result.summary() .. rst-class:: sphx-glr-script-out .. code-block:: none Optimized setup: 570 inverse variance values Running long MCMC chain: 10000 steps, 20 walkers Set RUN_LONG_MCMC = True in the configuration cell to run the original 10,000-step synthetic chain. Parameter bounds (unscaled): Lower: [ 1.e-08 -5.e+01 -5.e+01 -1.e+01 -1.e-03 -1.e-02 -1.e-03] Upper: [1.e-03 5.e+01 5.e+01 1.e+01 1.e-03 1.e-02 1.e-03] Scaled bounds: Lower: [ -7.44566917 -7.67844126 -6.8918876 -6.49935872 -0.79949471 -12.03003616 -0.06734446] Upper: [-2.44566917 7.67844126 6.8918876 6.49935872 0.79949471 12.03003616 0.06734446] Number of dimensions: 7 MCMC start (scaled/log): [-3.73377762 -4.61564393 0.42847056 0.40820822 -0.27220634 -3.73292208 0.00484093] Initial state shape: (20, 7) Initial walker mass min/max/std: -3.82945 -3.62403 0.0504557 All walkers are now within bounds. 0%| | 0/10000 [00:00 blob_names: ['log_likelihood', 'log_prior'] .. GENERATED FROM PYTHON SOURCE LINES 703-706 3.3 Getting the acceptance rates and auto-correlation times ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 706-784 .. code-block:: Python print("\n" + "="*50) print("RESULTS SUMMARY") print("="*50) param_names = PARAM_NAMES param_labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] sampler = inv_result.sampler reference_values_scaled = scaled_reference_values(scale_param, m_0) n_trajectory_samples = 10_000 discard_burn = 0 thin_interval = 10 saved_result = None if (not RUN_LONG_MCMC) and USE_SAVED_LONG_RUN_RESULTS: saved_result = load_mcmc_results(LONG_RUN_RESULTS_FILE, "synthetic") if saved_result is not None: analysis_chain = saved_result["chain"] flat_samples = analysis_chain[discard_burn::thin_interval].reshape(-1, analysis_chain.shape[-1]) acceptance_rates = saved_result["acceptance_fraction"] autocorr_times = saved_result["autocorr_times"] sample_source = f"saved long-run result ({LONG_RUN_RESULTS_FILE})" print(f"Loaded synthetic long-run results from {LONG_RUN_RESULTS_FILE}") else: analysis_chain, flat_samples, sample_source = analysis_samples( sampler, reference_values_scaled, bounds_lower_scaled, bounds_upper_scaled, thin=thin_interval, burn=discard_burn, seed=42, ) acceptance_rates = sampler.acceptance_fraction autocorr_times = safe_autocorr_time(sampler, fallback_chain=analysis_chain, discard=300) selected_samples = selected_trajectory_samples(flat_samples, n_trajectory_samples, seed=42) mean_acceptance = np.mean(acceptance_rates) max_autocorr = np.max(autocorr_times) mean_autocorr = np.mean(autocorr_times) if RUN_LONG_MCMC and SAVE_LONG_RUN_RESULTS: save_mcmc_results( LONG_RUN_RESULTS_FILE, "synthetic", analysis_chain, flat_samples, acceptance_rates, autocorr_times, nsteps, n_walkers, param_labels, ) print(f"Saved synthetic long-run results to {LONG_RUN_RESULTS_FILE}") print(f"Posterior samples for plots: {sample_source} ({len(flat_samples):,} samples)") print(f"Mean acceptance rate: {mean_acceptance:.3f}") print("Acceptance rate per walker:") for i, acc in enumerate(acceptance_rates): print(f" Walker {i:2d}: {acc:.3f}") print("\nAcceptance rate statistics:") print(f" Range: {np.min(acceptance_rates):.3f} - {np.max(acceptance_rates):.3f}") print(f" Std dev: {np.std(acceptance_rates):.3f}") print("Autocorrelation times by parameter:") for i, (label, tau) in enumerate(zip(param_labels, autocorr_times)): print(f" {label:>3s}: {tau:6.1f} steps") print("\nAutocorrelation statistics:") print(f" Mean τ: {mean_autocorr:.1f} steps") print(f" Max τ: {max_autocorr:.1f} steps") .. rst-class:: sphx-glr-script-out .. code-block:: none ================================================== RESULTS SUMMARY ================================================== Posterior samples for plots: emcee chain (20,000 samples) Mean acceptance rate: 0.265 Acceptance rate per walker: Walker 0: 0.461 Walker 1: 0.456 Walker 2: 0.455 Walker 3: 0.001 Walker 4: 0.460 Walker 5: 0.048 Walker 6: 0.015 Walker 7: 0.457 Walker 8: 0.003 Walker 9: 0.002 Walker 10: 0.454 Walker 11: 0.457 Walker 12: 0.002 Walker 13: 0.462 Walker 14: 0.106 Walker 15: 0.450 Walker 16: 0.467 Walker 17: 0.044 Walker 18: 0.449 Walker 19: 0.047 Acceptance rate statistics: Range: 0.001 - 0.467 Std dev: 0.214 Autocorrelation times by parameter: m: 718.9 steps x: 839.4 steps y: 794.0 steps z: 743.8 steps vx: 844.9 steps vy: 757.6 steps vz: 799.6 steps Autocorrelation statistics: Mean τ: 785.5 steps Max τ: 844.9 steps .. GENERATED FROM PYTHON SOURCE LINES 789-792 3.4 Plotting the Posteriors --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 795-800 - We now plot the posteriors and along with that, we alsomark the true parameter values of Neptune using a black ‘X’ mark. - Note that all our parameters are scaled. .. GENERATED FROM PYTHON SOURCE LINES 800-825 .. code-block:: Python reference_values_scaled = scaled_reference_values(scale_param, m_0) param_labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] ref_vals = [reference_values_scaled[i] for i in range(n_dim)] print("Creating corner plot with KDE contours and reference values...") fig = plot_pair_kde_legacy_style( flat_samples, param_labels, ref_vals, figsize=(16, 14), textsize=10, hdi_probs=(0.3, 0.6, 0.9), ) plt.suptitle(f'MCMC Parameter Posterior Distributions\n' f'(N_eff = {len(flat_samples):,}, burn-in = {discard_burn}, thin = {thin_interval})', fontsize=14, y=0.98) plt.tight_layout() plt.subplots_adjust(top=0.94) plt.show() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_002.png :alt: MCMC Parameter Posterior Distributions (N_eff = 20,000, burn-in = 0, thin = 10) :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Creating corner plot with KDE contours and reference values... .. GENERATED FROM PYTHON SOURCE LINES 830-833 3.5 Plotting the Chains ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 833-864 .. code-block:: Python fig, axes = plt.subplots(3, 3, figsize=(15, 12)) axes = axes.flatten() chain_full = analysis_chain labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] for i in range(n_dim): ax = axes[i] for walker in range(min(14, chain_full.shape[1])): ax.plot(chain_full[:, walker, i], alpha=0.6, linewidth=0.5) ax.axhline(reference_values_scaled[i], color='black', linestyle='--', alpha=0.8, linewidth=2, label='Reference Value' if i == 0 else '') ax.set_xlabel('Step') ax.set_ylabel(f'{labels[i]}') ax.set_title(f'Trace: {labels[i]}') ax.grid(True, alpha=0.3) if i == 0: ax.legend() for i in range(n_dim, len(axes)): axes[i].remove() plt.suptitle('MCMC Trace Plots (All Walkers)', fontsize=16) plt.tight_layout() plt.show() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_003.png :alt: MCMC Trace Plots (All Walkers), Trace: m, Trace: x, Trace: y, Trace: z, Trace: vx, Trace: vy, Trace: vz :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 869-894 4. Real Data Inversion ---------------------- - We then apply our deterministic inversion via CoFI on actual observational data obtained from `NASA JPL Horizons `__. - The data consists of geometric Cartesian position and velocity vectors of **Uranus**, relative to the **Solar System Barycenter**, with the following settings: - **Target body**: Uranus (799) - **Center body**: Solar System Barycenter (0) - **Reference frame**: Ecliptic of J2000.0 - **Time span**: A.D. 1775-Jan-01 to 2125-Jan-02 - **Step size**: 1 calendar year - **Output format**: Cartesian position and velocity (AU, AU/day) - **Output type**: GEOMETRIC states - **Calendar mode**: Mixed Julian/Gregorian - **Ephemeris source**: ``ura183_merged`` (Uranus), ``DE441`` (Solar System) This dataset provides real-world observations to test the robustness of our inversion pipeline. .. GENERATED FROM PYTHON SOURCE LINES 894-899 .. code-block:: Python #uncomment to install the astroquery package # !pip install astroquery .. GENERATED FROM PYTHON SOURCE LINES 904-907 4.1 Getting the actual data from NASA ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 907-912 .. code-block:: Python from neptune_deterministic_methods import get_actual_data U_true = get_actual_data(z_scaling = False, T = 190) .. rst-class:: sphx-glr-script-out .. code-block:: none Original data ranges: X: -18.297 to 20.085 (range: 38.382) Y: -19.319 to 19.013 (range: 38.332) Z: -0.272231 to 0.248352 (range: 0.520583) Z scaling is disabled, using original Z values. .. GENERATED FROM PYTHON SOURCE LINES 917-926 4.2 Setting the Inversion and Running it ---------------------------------------- The real-data section uses the same Bayesian objective as the synthetic section, but calibrates the likelihood with the effective per-coordinate residual scales obtained from the deterministic Morozov analysis. The MCMC walkers are initialized near that deterministic real-data solution rather than across the whole prior box. .. GENERATED FROM PYTHON SOURCE LINES 926-977 .. code-block:: Python param_bounds_lower = np.array([bound[0] for bound in PARAM_BOUNDS]) param_bounds_upper = np.array([bound[1] for bound in PARAM_BOUNDS]) bounds_lower_scaled = scale_param(param_bounds_lower) bounds_upper_scaled = scale_param(param_bounds_upper) bounds_lower_scaled[0] = np.log10(bounds_lower_scaled[0]) # Convert mass back to log scale bounds_upper_scaled[0] = np.log10(bounds_upper_scaled[0]) # Convert mass back to log scale # Effective real-data residual scales from the deterministic Morozov analysis. # These include observational noise plus model mismatch for this forward model. REAL_DETERMINISTIC_SIGMAS = np.array([7.8277e-02, 7.4224e-02, 9.7833e-04]) # sigma_x, sigma_y, sigma_z = REAL_DETERMINISTIC_SIGMAS sigma_x, sigma_y, sigma_z = np.array([5e-1, 5e-1, 5e-3]) diagonal_values = [1/sigma_x**2]*T + [1/sigma_y**2]*T + [1/sigma_z**2]*T inv_variances = np.array(diagonal_values) # Deterministic real-data solution from the Morozov run; used as the MCMC centre. REAL_DETERMINISTIC_M_UNSCALED = np.array([ 5.129472e-05, -2.991824e+01, 3.145458e+00, 1.050744e-01, -2.804390e-04, -3.109318e-03, 4.187314e-05, ]) real_mcmc_center_scaled = scale_param(REAL_DETERMINISTIC_M_UNSCALED).astype(float) real_mcmc_center_scaled[0] = np.log10(real_mcmc_center_scaled[0]) # real_mcmc_center_scaled = scale_param(m_0).astype(float) # real_mcmc_center_scaled[0] = np.log10(real_mcmc_center_scaled[0]) real_center_pred = predict_U(real_mcmc_center_scaled, T=T, dt=dt, z_scale_factor=z_scale_factor) real_center_residual = U_true - real_center_pred real_center_reduced_chi2 = np.sum(real_center_residual**2 * inv_variances) / len(inv_variances) print("Real-data effective sigmas:") print(f" sigma_x = {sigma_x:.4e}") print(f" sigma_y = {sigma_y:.4e}") print(f" sigma_z = {sigma_z:.4e}") print("Deterministic real-data MCMC centre (scaled/log):", real_mcmc_center_scaled) print("Residual RMS at deterministic centre:") print(f" X: {np.sqrt(np.mean(real_center_residual[:T]**2)):.4e}") print(f" Y: {np.sqrt(np.mean(real_center_residual[T:2*T]**2)):.4e}") print(f" Z: {np.sqrt(np.mean(real_center_residual[2*T:]**2)):.4e}") print(f"Reduced chi^2 at deterministic centre: {real_center_reduced_chi2:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Real-data effective sigmas: sigma_x = 5.0000e-01 sigma_y = 5.0000e-01 sigma_z = 5.0000e-03 Deterministic real-data MCMC centre (scaled/log): [-3.73559650e+00 -4.59450897e+00 4.33562860e-01 6.82916218e-02 -2.24209496e-01 -3.74052080e+00 2.81992383e-03] Residual RMS at deterministic centre: X: 7.9632e-02 Y: 7.4244e-02 Z: 9.9564e-04 Reduced chi^2 at deterministic centre: 0.0290 .. GENERATED FROM PYTHON SOURCE LINES 979-1100 .. code-block:: Python n_walkers = 28 n_dim = len(INVERT_INDICES) nsteps = REAL_LONG_STEPS if RUN_LONG_MCMC else VALIDATION_STEPS use_pool = True show_progress = RUN_LONG_MCMC run_mode = "long" if RUN_LONG_MCMC else "validation" print(f"Optimized setup: {len(inv_variances)} inverse variance values") print(f"Running {run_mode} MCMC chain: {nsteps} steps, {n_walkers} walkers") print("Set RUN_LONG_MCMC = True in the configuration cell to run the calibrated 80,000-step real-data chain.") print("Parameter bounds (unscaled):") print(" Lower:", param_bounds_lower) print(" Upper:", param_bounds_upper) print("Scaled bounds:") print(" Lower:", bounds_lower_scaled) print(" Upper:", bounds_upper_scaled) print(f"Number of dimensions: {n_dim}") rng = np.random.default_rng(84) walker_scales = 0.03 * (bounds_upper_scaled - bounds_lower_scaled) walker_scales[0] = 0.05 walkers_start = real_mcmc_center_scaled + rng.normal( scale=walker_scales, size=(n_walkers, n_dim), ) walkers_start = np.clip(walkers_start, bounds_lower_scaled, bounds_upper_scaled) print("Initial state shape:", walkers_start.shape) print( "Initial walker mass min/max/std:", f"{walkers_start[:, 0].min():.6g}", f"{walkers_start[:, 0].max():.6g}", f"{walkers_start[:, 0].std():.6g}", ) for i, w in enumerate(walkers_start): ll = my_log_likelihood(w, U_true) if not np.isfinite(ll): print(f"Invalid log-likelihood at walker {i}: {w}, ll = {ll}") for i, w in enumerate(walkers_start): walkers_start[i] = np.clip(w, bounds_lower_scaled, bounds_upper_scaled) if not np.array_equal(w, walkers_start[i]): print(f"Walker {i} clipped from {w} to {walkers_start[i]}") print("All walkers are now within bounds.") def my_log_prior(m_scaled): # The synthetic example can use a broad uniform prior because the data are # generated by this same simplified forward model. The real Horizons data # include full-Solar-System and ephemeris effects that this reduced model # cannot represent exactly. With only a uniform box prior, the sampler can # find unphysical Neptune parameters that overfit those model discrepancies. # This Gaussian prior is the Bayesian analogue of the deterministic # regularization term: it encodes external knowledge that Neptune's state is # near the MAP value determined in the deterministic inversion notebook # (finding_neptune_via_deterministic_inv) while still allowing uncertainty around it. if np.any(m_scaled < bounds_lower_scaled) or np.any(m_scaled > bounds_upper_scaled): return -np.inf prior_sigma_full = np.array([ 0.2, # log mass, atry 0.1 0.5, # x 0.5, # y 0.20, # z 0.1, # vx 0.1, # vy 0.20, # vz ]) prior_sigma = prior_sigma_full[INVERT_INDICES] dm = m_scaled - real_mcmc_center_scaled return -0.5 * np.sum((dm / prior_sigma)**2) for i, w in enumerate(walkers_start): if not np.isfinite(my_log_prior(w)): print(f"ERROR: Walker {i} still outside prior bounds: {w}") else: ll = my_log_likelihood(w, U_true) if not np.isfinite(ll): print(f"WARNING: Walker {i} has invalid log-likelihood: {ll}") inv_problem = BaseProblem() inv_problem.name = "Neptune Orbit Determination - Config Driven" inv_problem.set_data(U_true) inv_problem.set_forward(predict_U) inv_problem.set_initial_model(np.atleast_1d(real_mcmc_center_scaled)) inv_problem.set_log_prior(my_log_prior) inv_problem.set_log_likelihood(my_log_likelihood, args=(U_true,)) inv_problem.set_model_shape(n_dim) import emcee stretch_move = emcee.moves.StretchMove(a=1.2) inv_result, elapsed = run_cofi_emcee( inv_problem, walkers_start, n_walkers, nsteps, InversionOptions, Inversion, progress=show_progress, use_pool=use_pool, n_threads=14, moves=stretch_move, skip_initial_state_check=True, ) print(f"OPTIMIZED Inversion completed in {elapsed:.2f} seconds using `emcee`.") print("The inversion result from `emcee`:") inv_result.summary() .. rst-class:: sphx-glr-script-out .. code-block:: none Optimized setup: 570 inverse variance values Running long MCMC chain: 50000 steps, 28 walkers Set RUN_LONG_MCMC = True in the configuration cell to run the calibrated 80,000-step real-data chain. Parameter bounds (unscaled): Lower: [ 1.e-08 -5.e+01 -5.e+01 -1.e+01 -1.e-03 -1.e-02 -1.e-03] Upper: [1.e-03 5.e+01 5.e+01 1.e+01 1.e-03 1.e-02 1.e-03] Scaled bounds: Lower: [ -7.44566917 -7.67844126 -6.8918876 -6.49935872 -0.79949471 -12.03003616 -0.06734446] Upper: [-2.44566917 7.67844126 6.8918876 6.49935872 0.79949471 12.03003616 0.06734446] Number of dimensions: 7 Initial state shape: (28, 7) Initial walker mass min/max/std: -3.83146 -3.6507 0.0483306 All walkers are now within bounds. 0%| | 0/50000 [00:00 blob_names: ['log_likelihood', 'log_prior'] .. GENERATED FROM PYTHON SOURCE LINES 1105-1108 4.3 Getting the acceptance rates and auto-correlation times ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 1108-1186 .. code-block:: Python print("\n" + "="*50) print("RESULTS SUMMARY") print("="*50) param_names = PARAM_NAMES param_labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] sampler = inv_result.sampler reference_values_scaled = scaled_reference_values(scale_param, m_0) n_trajectory_samples = 10_000 discard_burn = 0 thin_interval = 1 saved_result = None if (not RUN_LONG_MCMC) and USE_SAVED_LONG_RUN_RESULTS: saved_result = load_mcmc_results(LONG_RUN_RESULTS_FILE, "real") if saved_result is not None: analysis_chain = saved_result["chain"] flat_samples = analysis_chain[discard_burn::thin_interval].reshape(-1, analysis_chain.shape[-1]) acceptance_rates = saved_result["acceptance_fraction"] autocorr_times = saved_result["autocorr_times"] sample_source = f"saved long-run result ({LONG_RUN_RESULTS_FILE})" print(f"Loaded real-data long-run results from {LONG_RUN_RESULTS_FILE}") else: analysis_chain, flat_samples, sample_source = analysis_samples( sampler, reference_values_scaled, bounds_lower_scaled, bounds_upper_scaled, thin=thin_interval, burn=discard_burn, seed=84, ) acceptance_rates = sampler.acceptance_fraction autocorr_times = safe_autocorr_time(sampler, fallback_chain=analysis_chain, discard=300) selected_samples = selected_trajectory_samples(flat_samples, n_trajectory_samples, seed=84) mean_acceptance = np.mean(acceptance_rates) max_autocorr = np.max(autocorr_times) mean_autocorr = np.mean(autocorr_times) if RUN_LONG_MCMC and SAVE_LONG_RUN_RESULTS: save_mcmc_results( LONG_RUN_RESULTS_FILE, "real", analysis_chain, flat_samples, acceptance_rates, autocorr_times, nsteps, n_walkers, param_labels, ) print(f"Saved real-data long-run results to {LONG_RUN_RESULTS_FILE}") print(f"Posterior samples for plots: {sample_source} ({len(flat_samples):,} samples)") print(f"Mean acceptance rate: {mean_acceptance:.3f}") print("Acceptance rate per walker:") for i, acc in enumerate(acceptance_rates): print(f" Walker {i:2d}: {acc:.3f}") print("\nAcceptance rate statistics:") print(f" Range: {np.min(acceptance_rates):.3f} - {np.max(acceptance_rates):.3f}") print(f" Std dev: {np.std(acceptance_rates):.3f}") print("Autocorrelation times by parameter:") for i, (label, tau) in enumerate(zip(param_labels, autocorr_times)): print(f" {label:>3s}: {tau:6.1f} steps") print("\nAutocorrelation statistics:") print(f" Mean τ: {mean_autocorr:.1f} steps") print(f" Max τ: {max_autocorr:.1f} steps") .. rst-class:: sphx-glr-script-out .. code-block:: none ================================================== RESULTS SUMMARY ================================================== Posterior samples for plots: emcee chain (1,400,000 samples) Mean acceptance rate: 0.637 Acceptance rate per walker: Walker 0: 0.615 Walker 1: 0.642 Walker 2: 0.669 Walker 3: 0.631 Walker 4: 0.662 Walker 5: 0.651 Walker 6: 0.623 Walker 7: 0.619 Walker 8: 0.659 Walker 9: 0.602 Walker 10: 0.640 Walker 11: 0.651 Walker 12: 0.637 Walker 13: 0.599 Walker 14: 0.639 Walker 15: 0.598 Walker 16: 0.633 Walker 17: 0.627 Walker 18: 0.648 Walker 19: 0.651 Walker 20: 0.650 Walker 21: 0.647 Walker 22: 0.622 Walker 23: 0.652 Walker 24: 0.649 Walker 25: 0.658 Walker 26: 0.597 Walker 27: 0.654 Acceptance rate statistics: Range: 0.597 - 0.669 Std dev: 0.020 Autocorrelation times by parameter: m: 578.7 steps x: 726.2 steps y: 785.4 steps z: 616.0 steps vx: 581.9 steps vy: 622.1 steps vz: 826.7 steps Autocorrelation statistics: Mean τ: 676.7 steps Max τ: 826.7 steps .. GENERATED FROM PYTHON SOURCE LINES 1191-1194 4.4 Plotting the Posteriors --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 1197-1202 - We now plot the posteriors and along with that, we alsomark the true parameter values of Neptune using a black ‘X’ mark. - Note that all our parameters are scaled. .. GENERATED FROM PYTHON SOURCE LINES 1202-1227 .. code-block:: Python reference_values_scaled = scaled_reference_values(scale_param, m_0) param_labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] ref_vals = [reference_values_scaled[i] for i in range(n_dim)] print("Creating corner plot with KDE contours and reference values...") fig = plot_pair_kde_legacy_style( flat_samples, param_labels, ref_vals, figsize=(16, 14), textsize=10, hdi_probs=(0.3, 0.6, 0.9), ) plt.suptitle(f'MCMC Parameter Posterior Distributions\n' f'(N_eff = {len(flat_samples):,}, burn-in = {discard_burn}, thin = {thin_interval})', fontsize=14, y=0.98) plt.tight_layout() plt.subplots_adjust(top=0.94) plt.show() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_004.png :alt: MCMC Parameter Posterior Distributions (N_eff = 1,400,000, burn-in = 0, thin = 1) :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Creating corner plot with KDE contours and reference values... .. GENERATED FROM PYTHON SOURCE LINES 1232-1235 4.5 Plotting the Chains ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 1235-1266 .. code-block:: Python fig, axes = plt.subplots(3, 3, figsize=(15, 12)) axes = axes.flatten() chain_full = analysis_chain labels = ["m" if idx == 0 else PARAM_NAMES[idx] for idx in INVERT_INDICES] for i in range(n_dim): ax = axes[i] for walker in range(min(14, chain_full.shape[1])): ax.plot(chain_full[:, walker, i], alpha=0.6, linewidth=0.5) ax.axhline(reference_values_scaled[i], color='black', linestyle='--', alpha=0.8, linewidth=2, label='Reference Value' if i == 0 else '') ax.set_xlabel('Step') ax.set_ylabel(f'{labels[i]}') ax.set_title(f'Trace: {labels[i]}') ax.grid(True, alpha=0.3) if i == 0: ax.legend() for i in range(n_dim, len(axes)): axes[i].remove() plt.suptitle('MCMC Trace Plots (All Walkers)', fontsize=16) plt.tight_layout() plt.show() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_005.png :alt: MCMC Trace Plots (All Walkers), Trace: m, Trace: x, Trace: y, Trace: z, Trace: vx, Trace: vy, Trace: vz :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_finding_neptune_via_bayesian_inv_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1268-1288 .. code-block:: Python true_scaled = scale_param(m_0).astype(float) true_scaled[0] = np.log10(true_scaled[0]) det_scaled = real_mcmc_center_scaled flat = sampler.get_chain(discard=10000, thin=100, flat=True) logp = sampler.get_log_prob(discard=10000, thin=100, flat=True) map_scaled = flat[np.argmax(logp)] print("log posterior at true:", my_log_prior(true_scaled) + my_log_likelihood(true_scaled, U_true)) print("log posterior at deterministic centre:", my_log_prior(det_scaled) + my_log_likelihood(det_scaled, U_true)) print("log posterior at chain MAP:", my_log_prior(map_scaled) + my_log_likelihood(map_scaled, U_true)) print("true:", true_scaled) print("det :", det_scaled) print("MAP :", map_scaled) .. rst-class:: sphx-glr-script-out .. code-block:: none log posterior at true: -9.975970292104922 log posterior at deterministic centre: -8.271254785759565 log posterior at chain MAP: -2.398142782038685 true: [-3.73377762 -4.61564393 0.42847056 0.40820822 -0.27220634 -3.73292208 0.00484093] det : [-3.73559650e+00 -4.59450897e+00 4.33562860e-01 6.82916218e-02 -2.24209496e-01 -3.74052080e+00 2.81992383e-03] MAP : [-3.57053498e+00 -3.70900076e+00 6.27588934e-01 9.68876835e-02 -2.24811733e-01 -3.69971220e+00 3.81295991e-04] .. GENERATED FROM PYTHON SOURCE LINES 1293-1298 5. Watermark ============ - For version of libraries used. .. GENERATED FROM PYTHON SOURCE LINES 1298-1320 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import tqdm import numba import arviz import astroquery import cofi def print_versions(): print("Library Versions:") print(f"{'numpy':<15}: {np.__version__}") print(f"{'matplotlib':<15}: {plt.matplotlib.__version__}") print(f"{'tqdm':<15}: {tqdm.__version__}") print(f"{'numba':<15}: {numba.__version__}") print(f"{'arviz':<15}: {arviz.__version__}") print(f"{'astroquery':<15}: {astroquery.__version__}") print(f"{'cofi':<15}: {cofi.__version__}") print_versions() .. rst-class:: sphx-glr-script-out .. code-block:: none Library Versions: numpy : 2.2.6 matplotlib : 3.10.9 tqdm : 4.67.3 numba : 0.65.1 arviz : 1.1.0 astroquery : 0.4.11 cofi : 0.2.11+68.gc319bc9 .. GENERATED FROM PYTHON SOURCE LINES 1321-1321 sphinx_gallery_thumbnail_number = -1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (244 minutes 28.758 seconds) .. _sphx_glr_download_examples_generated_scripts_synth_data_finding_neptune_via_bayesian_inv.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: finding_neptune_via_bayesian_inv.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: finding_neptune_via_bayesian_inv.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: finding_neptune_via_bayesian_inv.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_