.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/generated/scripts_synth_data/receiver_function_neighpy.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_generated_scripts_synth_data_receiver_function_neighpy.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_generated_scripts_synth_data_receiver_function_neighpy.py:


Receiver Function Inversion with the Neighbourhood Algorithm
============================================================

.. GENERATED FROM PYTHON SOURCE LINES 9-14

|Open In Colab|

.. |Open In Colab| image:: https://img.shields.io/badge/open%20in-Colab-b5e2fa?logo=googlecolab&style=flat-square&color=ffd670
   :target: https://colab.research.google.com/github/inlab-geo/cofi-examples/blob/main/examples/receiver_function/receiver_function_neighpy.ipynb


.. GENERATED FROM PYTHON SOURCE LINES 17-24

If you are running this notebook locally, make sure you’ve followed
`steps
here <https://github.com/inlab-geo/cofi-examples#run-the-examples-with-cofi-locally>`__
to set up the environment. (This
`environment.yml <https://github.com/inlab-geo/cofi-examples/blob/main/envs/environment.yml>`__
file specifies a list of packages required to run the notebooks)


.. GENERATED FROM PYTHON SOURCE LINES 27-29

--------------


.. GENERATED FROM PYTHON SOURCE LINES 32-38

1. Introduction 
----------------

1.1 Receiver functions 
~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 38-47

.. code-block:: Python


    # display theory on receiver functions
    from IPython.display import display, Markdown

    with open("../../theory/geo_receiver_function.md", "r") as f:
        content = f.read()

    display(Markdown(content))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    <IPython.core.display.Markdown object>




.. GENERATED FROM PYTHON SOURCE LINES 52-55

1.2 Neighbourhood Algorithm 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 55-64

.. code-block:: Python


    # display theory on the Neighbourhood Algorithm
    from IPython.display import display, Markdown

    with open("../../theory/neighbourhood_algorithm.md", "r") as f:
        content = f.read()

    display(Markdown(content))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    <IPython.core.display.Markdown object>




.. GENERATED FROM PYTHON SOURCE LINES 69-111

1.3 CoFI and the NA 
~~~~~~~~~~~~~~~~~~~~

The implementation of the NA that ``cofi`` wraps is called
```neighpy`` <https://github.com/inlab-geo/neighpy>`__. This
implementation implements both phases of the NA as described in the
original papers.

In this notebook we run the two stages of the NA **separately** using
the ``neighpyI`` and ``neighpyII`` tools in CoFI:

-  **Stage I** (``neighpyI``) — Direct search. An ensemble optimiser
   that minimises a user-defined **misfit** function to explore the
   parameter space.
-  **Stage II** (``neighpyII``) — Appraisal. An MCMC resampler that
   resamples the Stage I ensemble according to a user-supplied **log
   posterior probability density** (log-PPD). No additional forward
   model evaluations are required.

The key insight is that the log-PPD used in Stage II does not have to
equal the negative of the misfit used in Stage I. The relationship
between the two is a modelling choice that depends on the assumed
likelihood and prior. In this notebook we use
:math:`\log p(\mathbf{m}|\mathbf{d}) = -\Phi(\mathbf{m},\mathbf{d})`,
which corresponds to a Gaussian likelihood with a uniform prior, but we
make this choice explicit so the user can modify it.

We apply the NA to two receiver function inversion problems of
increasing complexity:

1. **Problem 1** — Invert for S-wave velocities only (4 parameters,
   layer depths fixed)
2. **Problem 2** — Invert for both layer depths and S-wave velocities (8
   parameters)

..

   **Note:** The hyperparameters in this notebook are intentionally set
   low for a quick demonstration run. Results may not be fully
   converged. To improve them, increase ``n_initial_samples``,
   ``n_iterations``, and ``n_resample`` in the cells below.


.. GENERATED FROM PYTHON SOURCE LINES 114-116

--------------


.. GENERATED FROM PYTHON SOURCE LINES 119-122

2. Import modules 
------------------


.. GENERATED FROM PYTHON SOURCE LINES 122-131

.. code-block:: Python


    # -------------------------------------------------------- #
    #                                                          #
    #     Uncomment below to set up environment on "colab"     #
    #                                                          #
    # -------------------------------------------------------- #

    # !pip install -U cofi pyrf96








.. GENERATED FROM PYTHON SOURCE LINES 133-143

.. code-block:: Python


    import math
    import numpy as np
    import matplotlib.pyplot as plt

    from cofi import BaseProblem, InversionOptions, Inversion
    import pyrf96

    np.random.seed(42)








.. GENERATED FROM PYTHON SOURCE LINES 148-150

--------------


.. GENERATED FROM PYTHON SOURCE LINES 153-164

3. Define the problem 
----------------------

We use a 4-layer Earth model in pyrf96’s Voronoi nuclei format
(``mtype=0``). Each row is ``[depth_km, Vs_km/s, Vp/Vs]``. The Vp/Vs
ratios are held fixed throughout the inversion — only depths and/or
velocities are inverted.

We generate synthetic observed data by computing the receiver function
for the true model and adding correlated noise.


.. GENERATED FROM PYTHON SOURCE LINES 164-184

.. code-block:: Python


    # True Earth model: 4 layers [depth_km, Vs_km/s, Vp/Vs]
    good_model = np.array([
        [ 1.0, 3.0, 1.7],
        [ 8.0, 3.2, 2.0],
        [20.0, 4.0, 1.7],
        [45.0, 4.2, 1.7],
    ])

    # Fixed Vp/Vs ratios
    vpvs = good_model[:, 2]

    # Generate synthetic data
    t, rfunc = pyrf96.rfcalc(good_model)                                  # noise-free
    t2, rfunc_noise = pyrf96.rfcalc(good_model, sn=0.5, seed=12345678)   # noisy observed
    observed_data = rfunc_noise

    # Inverse data covariance matrix for correlated noise
    Cdinv = pyrf96.InvDataCov(2.5, 0.01, len(rfunc))








.. GENERATED FROM PYTHON SOURCE LINES 186-197

.. code-block:: Python


    def plot_model_staircase(model, ax, max_depth=60.0, **kwargs):
        """Plot a pyrf96 Voronoi nuclei model as a Vs-depth staircase."""
        order = np.argsort(model[:, 0])
        depths = model[order, 0]
        vs = model[order, 1]
        interfaces = 0.5 * (depths[:-1] + depths[1:])
        plot_depths = np.concatenate([[0.0], np.repeat(interfaces, 2), [max_depth]])
        plot_vs = np.repeat(vs, 2)
        ax.plot(plot_vs, plot_depths, **kwargs)








.. GENERATED FROM PYTHON SOURCE LINES 199-222

.. code-block:: Python


    # Plot the true model and observed data
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5))

    # Velocity-depth profile
    plot_model_staircase(good_model, ax1, color="red", lw=2, label="True model")
    ax1.set_xlim(2.5, 5.0)
    ax1.set_ylim(60, 0)
    ax1.set_xlabel("Vs (km/s)")
    ax1.set_ylabel("Depth (km)")
    ax1.set_title("True Earth model")
    ax1.legend()

    # Receiver functions
    ax2.plot(t2, rfunc_noise, color="k", label="Observed (noisy)")
    ax2.plot(t, rfunc, color="r", ls="--", label="True (noise-free)")
    ax2.set_xlabel("Time (s)")
    ax2.set_ylabel("Amplitude")
    ax2.set_title("Synthetic receiver function")
    ax2.legend()

    plt.tight_layout()




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_001.png
   :alt: True Earth model, Synthetic receiver function
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 227-229

--------------


.. GENERATED FROM PYTHON SOURCE LINES 232-239

4. Problem 1: Invert for velocities only 
-----------------------------------------

In this first problem we fix the layer depths at their true values and
invert only for the 4 S-wave velocities. This is a simpler 4-dimensional
problem.


.. GENERATED FROM PYTHON SOURCE LINES 239-272

.. code-block:: Python


    # Conversion functions for velocity-only inversion
    fixed_depths = good_model[:, 0]  # [1, 8, 20, 45] — held fixed

    def get_vel_params(fullmodel):
        """Extract velocity parameters from full model."""
        return fullmodel[:, 1].copy()

    def make_model_vel(vel_params):
        """Reconstruct full model from velocity parameters (fixed depths and Vp/Vs)."""
        return np.column_stack([fixed_depths, vel_params, vpvs])

    # Misfit function for velocity-only inversion (used by Stage I)
    def misfit_vel(vel_params):
        model = make_model_vel(vel_params)
        t_pred, predicted_data = pyrf96.rfcalc(model)
        res = observed_data - predicted_data
        misfit_val = np.dot(res, np.transpose(np.dot(Cdinv, res))) / 2.0
        if math.isnan(misfit_val):
            return float("inf")
        return misfit_val

    # Log posterior function (used by Stage II)
    def log_posterior_vel(vel_params):
        """Log posterior = -misfit (uniform prior, Gaussian likelihood)."""
        return -misfit_vel(vel_params)

    # CoFI BaseProblem
    inv_problem_vel = BaseProblem()
    inv_problem_vel.name = "Receiver Function - Velocity only"
    inv_problem_vel.set_objective(misfit_vel)
    inv_problem_vel.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    =====================================================================
    Summary for inversion problem: Receiver Function - Velocity only
    =====================================================================
    Model shape: Unknown
    ---------------------------------------------------------------------
    List of functions/properties set by you:
    ['objective']
    ---------------------------------------------------------------------
    List of functions/properties created based on what you have provided:
    -- none --
    ---------------------------------------------------------------------
    List of functions/properties that can be further set for the problem:
    ( not all of these may be relevant to your inversion workflow )
    ['log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian', 'jacobian_times_vector', 'data_misfit', 'regularization', 'regularization_matrix', 'forward', 'data', 'data_covariance', 'data_covariance_inv', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']




.. GENERATED FROM PYTHON SOURCE LINES 277-286

4.1 Stage I: Direct Search 
~~~~~~~~~~~~~~~~~~~~~~~~~~~

Stage I of the Neighbourhood Algorithm is a derivative-free ensemble
optimiser. It minimises the misfit function defined above by iteratively
sampling new models in the neighbourhood of the best-fitting models
found so far. The ``neighpyI`` tool in CoFI runs only this direct search
phase.


.. GENERATED FROM PYTHON SOURCE LINES 286-304

.. code-block:: Python


    # NA-I hyperparameters for 4-D problem
    bounds_vel = [(2.0, 4.5)] * 4

    inv_options_vel = InversionOptions()
    inv_options_vel.set_tool("neighpyI")
    inv_options_vel.set_params(
        bounds=bounds_vel,
        n_initial_samples=2000,
        n_samples_per_iteration=50,
        n_cells_to_resample=10,
        n_iterations=100,
    )

    inv_vel = Inversion(inv_problem_vel, inv_options_vel)
    inv_result_vel = inv_vel.run()
    inv_result_vel.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    NAI - Initial Random Search

    NAI - Optimisation Loop:   0%|          | 0/100 [00:00<?, ?it/s]
    NAI - Optimisation Loop:   1%|          | 1/100 [00:01<02:04,  1.26s/it]
    NAI - Optimisation Loop:   3%|▎         | 3/100 [00:01<00:36,  2.62it/s]
    NAI - Optimisation Loop:   5%|▌         | 5/100 [00:01<00:21,  4.49it/s]
    NAI - Optimisation Loop:   7%|▋         | 7/100 [00:01<00:14,  6.25it/s]
    NAI - Optimisation Loop:   9%|▉         | 9/100 [00:01<00:11,  7.84it/s]
    NAI - Optimisation Loop:  11%|█         | 11/100 [00:02<00:09,  9.15it/s]
    NAI - Optimisation Loop:  13%|█▎        | 13/100 [00:02<00:08, 10.27it/s]
    NAI - Optimisation Loop:  15%|█▌        | 15/100 [00:02<00:07, 11.14it/s]
    NAI - Optimisation Loop:  17%|█▋        | 17/100 [00:02<00:07, 11.84it/s]
    NAI - Optimisation Loop:  19%|█▉        | 19/100 [00:02<00:06, 12.35it/s]
    NAI - Optimisation Loop:  21%|██        | 21/100 [00:02<00:06, 12.71it/s]
    NAI - Optimisation Loop:  23%|██▎       | 23/100 [00:02<00:05, 12.90it/s]
    NAI - Optimisation Loop:  25%|██▌       | 25/100 [00:03<00:05, 13.10it/s]
    NAI - Optimisation Loop:  27%|██▋       | 27/100 [00:03<00:05, 13.21it/s]
    NAI - Optimisation Loop:  29%|██▉       | 29/100 [00:03<00:05, 13.35it/s]
    NAI - Optimisation Loop:  31%|███       | 31/100 [00:03<00:05, 13.41it/s]
    NAI - Optimisation Loop:  33%|███▎      | 33/100 [00:03<00:04, 13.48it/s]
    NAI - Optimisation Loop:  35%|███▌      | 35/100 [00:03<00:04, 13.41it/s]
    NAI - Optimisation Loop:  37%|███▋      | 37/100 [00:03<00:04, 13.41it/s]
    NAI - Optimisation Loop:  39%|███▉      | 39/100 [00:04<00:04, 13.45it/s]
    NAI - Optimisation Loop:  41%|████      | 41/100 [00:04<00:04, 13.44it/s]
    NAI - Optimisation Loop:  43%|████▎     | 43/100 [00:04<00:04, 13.43it/s]
    NAI - Optimisation Loop:  45%|████▌     | 45/100 [00:04<00:04, 13.40it/s]
    NAI - Optimisation Loop:  47%|████▋     | 47/100 [00:04<00:03, 13.43it/s]
    NAI - Optimisation Loop:  49%|████▉     | 49/100 [00:04<00:03, 13.42it/s]
    NAI - Optimisation Loop:  51%|█████     | 51/100 [00:04<00:03, 13.43it/s]
    NAI - Optimisation Loop:  53%|█████▎    | 53/100 [00:05<00:03, 13.45it/s]
    NAI - Optimisation Loop:  55%|█████▌    | 55/100 [00:05<00:03, 13.53it/s]
    NAI - Optimisation Loop:  57%|█████▋    | 57/100 [00:05<00:03, 13.55it/s]
    NAI - Optimisation Loop:  59%|█████▉    | 59/100 [00:05<00:03, 13.55it/s]
    NAI - Optimisation Loop:  61%|██████    | 61/100 [00:05<00:02, 13.52it/s]
    NAI - Optimisation Loop:  63%|██████▎   | 63/100 [00:05<00:02, 13.48it/s]
    NAI - Optimisation Loop:  65%|██████▌   | 65/100 [00:06<00:02, 13.40it/s]
    NAI - Optimisation Loop:  67%|██████▋   | 67/100 [00:06<00:02, 13.40it/s]
    NAI - Optimisation Loop:  69%|██████▉   | 69/100 [00:06<00:02, 13.39it/s]
    NAI - Optimisation Loop:  71%|███████   | 71/100 [00:06<00:02, 13.39it/s]
    NAI - Optimisation Loop:  73%|███████▎  | 73/100 [00:06<00:02, 13.39it/s]
    NAI - Optimisation Loop:  75%|███████▌  | 75/100 [00:06<00:01, 13.43it/s]
    NAI - Optimisation Loop:  77%|███████▋  | 77/100 [00:06<00:01, 13.43it/s]
    NAI - Optimisation Loop:  79%|███████▉  | 79/100 [00:07<00:01, 13.45it/s]
    NAI - Optimisation Loop:  81%|████████  | 81/100 [00:07<00:01, 13.40it/s]
    NAI - Optimisation Loop:  83%|████████▎ | 83/100 [00:07<00:01, 13.39it/s]
    NAI - Optimisation Loop:  85%|████████▌ | 85/100 [00:07<00:01, 13.35it/s]
    NAI - Optimisation Loop:  87%|████████▋ | 87/100 [00:07<00:00, 13.32it/s]
    NAI - Optimisation Loop:  89%|████████▉ | 89/100 [00:07<00:00, 13.33it/s]
    NAI - Optimisation Loop:  91%|█████████ | 91/100 [00:07<00:00, 13.34it/s]
    NAI - Optimisation Loop:  93%|█████████▎| 93/100 [00:08<00:00, 13.39it/s]
    NAI - Optimisation Loop:  95%|█████████▌| 95/100 [00:08<00:00, 13.38it/s]
    NAI - Optimisation Loop:  97%|█████████▋| 97/100 [00:08<00:00, 13.38it/s]
    NAI - Optimisation Loop:  99%|█████████▉| 99/100 [00:08<00:00, 13.39it/s]
    NAI - Optimisation Loop: 100%|██████████| 100/100 [00:08<00:00, 11.60it/s]
    ============================
    Summary for inversion result
    ============================
    SUCCESS
    ----------------------------
    model: [3.08576275 3.15413977 4.05676828 4.47140278]
    samples: [[4.31768078 2.17658461 2.38897106 3.12093439]
     [3.98645265 3.85217428 4.22154291 3.43285316]
     [4.3300403  4.33965414 3.66323452 2.95553822]
     ...
     [3.08576269 3.1541397  4.05676856 4.47140271]
     [3.08576269 3.15413973 4.05676854 4.47140269]
     [3.08576271 3.15413972 4.05676856 4.47140269]]
    objectives: [3013.70661133 2142.70739802 3886.83915189 ...  445.86818384  445.86818384
      445.86818384]




.. GENERATED FROM PYTHON SOURCE LINES 306-316

.. code-block:: Python


    best_vel = inv_result_vel.model
    ds_samples_vel = inv_result_vel.samples
    ds_objectives_vel = inv_result_vel.objectives

    print("Best model (velocities):")
    print(make_model_vel(best_vel))
    print("\nTrue model:")
    print(good_model)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Best model (velocities):
    [[ 1.          3.08576275  1.7       ]
     [ 8.          3.15413977  2.        ]
     [20.          4.05676828  1.7       ]
     [45.          4.47140278  1.7       ]]

    True model:
    [[ 1.   3.   1.7]
     [ 8.   3.2  2. ]
     [20.   4.   1.7]
     [45.   4.2  1.7]]




.. GENERATED FROM PYTHON SOURCE LINES 321-324

Convergence
^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 324-342

.. code-block:: Python


    best_i_vel = np.argmin(ds_objectives_vel)

    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(ds_objectives_vel, marker=".", linestyle="", markersize=2, color="black")
    ax.scatter(best_i_vel, ds_objectives_vel[best_i_vel], c="g", s=30, zorder=10,
               label="Best model")
    ax.axvline(2000, c="grey", ls="--", lw=0.8)
    ax.set_yscale("log")
    ax.set_xlabel("Sample number")
    ax.set_ylabel("Misfit")
    ax.set_title("Problem 1: NA-I convergence (velocity-only)")
    ax.text(0.15, 0.95, "Initial random search", transform=ax.transAxes,
            ha="center", va="top")
    ax.text(0.65, 0.95, "Neighbourhood search", transform=ax.transAxes,
            ha="center", va="top")
    ax.legend()




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_002.png
   :alt: Problem 1: NA-I convergence (velocity-only)
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcb1c98cd0>



.. GENERATED FROM PYTHON SOURCE LINES 347-350

NA-I model ensemble
^^^^^^^^^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 350-368

.. code-block:: Python


    fig, ax = plt.subplots(figsize=(5, 8))

    # Plot initial random samples as faint ensemble
    for i in range(min(2000, len(ds_samples_vel))):
        m = make_model_vel(ds_samples_vel[i])
        plot_model_staircase(m, ax, color="k", alpha=0.01, lw=0.5)

    plot_model_staircase(good_model, ax, color="b", lw=2, label="True model")
    plot_model_staircase(make_model_vel(best_vel), ax, color="g", lw=2, ls="--",
                         label="Best model (NA-I)")
    ax.set_xlim(2.5, 5.0)
    ax.set_ylim(60, 0)
    ax.set_xlabel("Vs (km/s)")
    ax.set_ylabel("Depth (km)")
    ax.set_title("Problem 1: NA-I model ensemble")
    ax.legend(loc="lower left")




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_003.png
   :alt: Problem 1: NA-I model ensemble
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcb908aad0>



.. GENERATED FROM PYTHON SOURCE LINES 373-376

NA-I predicted data
^^^^^^^^^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 376-397

.. code-block:: Python


    fig, ax = plt.subplots(figsize=(10, 4))

    # Plot predicted RFs from initial samples
    for i in range(min(2000, len(ds_samples_vel))):
        m = make_model_vel(ds_samples_vel[i])
        try:
            t_pred, rf_pred = pyrf96.rfcalc(m)
            ax.plot(t_pred, rf_pred, color="k", alpha=0.01, lw=0.5)
        except Exception:
            pass

    ax.plot(t2, observed_data, color="k", lw=1.5, label="Observed", zorder=10)
    t_best, rf_best = pyrf96.rfcalc(make_model_vel(best_vel))
    ax.plot(t_best, rf_best, color="g", ls="--", lw=1.5, label="Best model (NA-I)",
            zorder=11)
    ax.set_xlabel("Time (s)")
    ax.set_ylabel("Amplitude")
    ax.set_title("Problem 1: Predicted receiver functions")
    ax.legend()




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_004.png
   :alt: Problem 1: Predicted receiver functions
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_004.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcca1c56d0>



.. GENERATED FROM PYTHON SOURCE LINES 402-422

4.2 Stage II: Appraisal 
~~~~~~~~~~~~~~~~~~~~~~~~

Stage II of the NA resamples the Voronoi cells constructed by Stage I
according to a **log posterior probability density** (log-PPD).
Crucially, this stage requires no additional forward model evaluations —
it only uses the geometry of the Stage I ensemble and the log-PPD values
assigned to each sample.

The user must provide ``log_ppd`` values for each sample in the Stage I
ensemble. The relationship between the log-PPD and the misfit is a
**modelling choice**:

.. math:: \log p(\mathbf{m}|\mathbf{d}) = -\Phi(\mathbf{m},\mathbf{d})

corresponds to a Gaussian likelihood with a uniform prior. More
generally, one could incorporate a non-uniform prior or a different
likelihood function, which would change the ``log_ppd`` values without
requiring Stage I to be re-run.


.. GENERATED FROM PYTHON SOURCE LINES 422-429

.. code-block:: Python


    # Compute log-PPD for each Stage I sample.
    # Since log_posterior_vel(s) = -misfit_vel(s), this is equivalent to:
    #   log_ppd_vel = np.array([log_posterior_vel(s) for s in ds_samples_vel])
    # but using the already-computed misfit values avoids redundant forward evaluations.
    log_ppd_vel = -ds_objectives_vel








.. GENERATED FROM PYTHON SOURCE LINES 431-447

.. code-block:: Python


    # NA-II hyperparameters
    inv_options_vel_app = InversionOptions()
    inv_options_vel_app.set_tool("neighpyII")
    inv_options_vel_app.set_params(
        bounds=bounds_vel,
        initial_ensemble=ds_samples_vel,
        log_ppd=log_ppd_vel,
        n_resample=20000,
        n_walkers=10,
    )

    inv_vel_app = Inversion(inv_problem_vel, inv_options_vel_app)
    inv_result_vel_app = inv_vel_app.run()
    inv_result_vel_app.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ============================
    Summary for inversion result
    ============================
    SUCCESS
    ----------------------------
    new_samples: [[3.08936785 3.14424749 4.04079552 4.44101329]
     [3.11039626 3.09111185 4.03866865 4.43234406]
     [3.10455271 3.17238427 4.02898751 4.49612852]
     ...
     [3.11025929 3.1584892  3.99314115 4.33662172]
     [3.09754188 3.15536398 4.03084395 4.48797973]
     [3.12008656 3.08642429 4.02306098 4.48868035]]




.. GENERATED FROM PYTHON SOURCE LINES 449-470

.. code-block:: Python


    appraisal_samples_vel = inv_result_vel_app.new_samples
    appraisal_mean_vel = appraisal_samples_vel.mean(axis=0)

    fig, ax = plt.subplots(figsize=(5, 8))
    for i in range(0, len(appraisal_samples_vel), 10):
        m = make_model_vel(appraisal_samples_vel[i])
        plot_model_staircase(m, ax, color="k", alpha=0.01, lw=0.5)

    plot_model_staircase(good_model, ax, color="b", lw=2, label="True model")
    plot_model_staircase(make_model_vel(best_vel), ax, color="g", lw=2, ls="--",
                         label="Best model (NA-I)")
    plot_model_staircase(make_model_vel(appraisal_mean_vel), ax, color="r", lw=2,
                         ls="--", label="Mean model (NA-II)")
    ax.set_xlim(2.5, 5.0)
    ax.set_ylim(60, 0)
    ax.set_xlabel("Vs (km/s)")
    ax.set_ylabel("Depth (km)")
    ax.set_title("Problem 1: NA-II appraisal ensemble")
    ax.legend(loc="lower left")




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_005.png
   :alt: Problem 1: NA-II appraisal ensemble
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_005.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcbc64c910>



.. GENERATED FROM PYTHON SOURCE LINES 472-511

.. code-block:: Python


    # 4x4 corner plot for velocity-only problem
    true_vel = get_vel_params(good_model)
    n_params_vel = 4
    var_names_vel = [f"Vs{i+1}" for i in range(n_params_vel)]

    fig, axes = plt.subplots(n_params_vel, n_params_vel, figsize=(8, 8),
                             tight_layout=True,
                             gridspec_kw={"hspace": 0, "wspace": 0})
    for i in range(n_params_vel):
        for j in range(n_params_vel):
            axes[i, j].set_xlim(bounds_vel[j])
            axes[i, j].set_xticks([])
            axes[i, j].set_yticks([])
            if i == j:
                axes[i, j].hist(appraisal_samples_vel[:, j], bins=100,
                                histtype="step", color="k")
                axes[i, j].axvline(best_vel[j], c="g", ls="--")
                axes[i, j].axvline(true_vel[j], c="b", ls="--")
                axes[i, j].axvline(appraisal_mean_vel[j], c="r", ls="--")
            elif j < i:
                axes[i, j].scatter(appraisal_samples_vel[:, j],
                                   appraisal_samples_vel[:, i],
                                   s=1, c="k", alpha=0.01)
                axes[i, j].scatter(best_vel[j], best_vel[i], c="g", s=10, zorder=10)
                axes[i, j].scatter(true_vel[j], true_vel[i], c="b", s=10, zorder=10)
                axes[i, j].scatter(appraisal_mean_vel[j], appraisal_mean_vel[i],
                                   c="r", s=10, zorder=10)
                axes[i, j].set_ylim(bounds_vel[i])
            else:
                axes[i, j].axis("off")
            if i == n_params_vel - 1 and j <= i:
                axes[i, j].set_xlabel(var_names_vel[j])
                axes[i, j].set_xticks(np.linspace(*bounds_vel[j], 3))
            if j == 0 and i > 0:
                axes[i, j].set_ylabel(var_names_vel[i])

    fig.suptitle("Problem 1: Posterior corner plot (velocity-only)")




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_006.png
   :alt: Problem 1: Posterior corner plot (velocity-only)
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_006.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 0.98, 'Problem 1: Posterior corner plot (velocity-only)')



.. GENERATED FROM PYTHON SOURCE LINES 516-518

--------------


.. GENERATED FROM PYTHON SOURCE LINES 521-528

5. Problem 2: Invert for depths and velocities 
-----------------------------------------------

Now we invert for both layer depths and S-wave velocities
simultaneously. This is a harder 8-dimensional problem with interleaved
parameters: ``[d1, vs1, d2, vs2, d3, vs3, d4, vs4]``.


.. GENERATED FROM PYTHON SOURCE LINES 528-559

.. code-block:: Python


    # Conversion functions for depth+velocity inversion
    def get_inversion_parameters(fullmodel):
        """Extract [depth, Vs] pairs as flat 8-vector."""
        return fullmodel[:, :2].flatten()

    def get_model_parameters(invmodel):
        """Reconstruct full [4,3] model from 8-vector, restoring fixed Vp/Vs."""
        return np.append(invmodel.reshape(len(vpvs), -1), vpvs[:, None], axis=1)

    # Misfit function for depth+velocity inversion (used by Stage I)
    def misfit_dv(imodel):
        model = get_model_parameters(imodel)
        t_pred, predicted_data = pyrf96.rfcalc(model)
        res = observed_data - predicted_data
        misfit_val = np.dot(res, np.transpose(np.dot(Cdinv, res))) / 2.0
        if math.isnan(misfit_val):
            return float("inf")
        return misfit_val

    # Log posterior function (used by Stage II)
    def log_posterior_dv(imodel):
        """Log posterior = -misfit (uniform prior, Gaussian likelihood)."""
        return -misfit_dv(imodel)

    # CoFI BaseProblem
    inv_problem_dv = BaseProblem()
    inv_problem_dv.name = "Receiver Function - Depth + Velocity"
    inv_problem_dv.set_objective(misfit_dv)
    inv_problem_dv.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    =====================================================================
    Summary for inversion problem: Receiver Function - Depth + Velocity
    =====================================================================
    Model shape: Unknown
    ---------------------------------------------------------------------
    List of functions/properties set by you:
    ['objective']
    ---------------------------------------------------------------------
    List of functions/properties created based on what you have provided:
    -- none --
    ---------------------------------------------------------------------
    List of functions/properties that can be further set for the problem:
    ( not all of these may be relevant to your inversion workflow )
    ['log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian', 'jacobian_times_vector', 'data_misfit', 'regularization', 'regularization_matrix', 'forward', 'data', 'data_covariance', 'data_covariance_inv', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']




.. GENERATED FROM PYTHON SOURCE LINES 564-567

5.1 Stage I: Direct Search 
~~~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 567-585

.. code-block:: Python


    # NA-I hyperparameters for 8-D problem (larger than Problem 1)
    bounds_dv = [(0.0, 60.0), (2.0, 4.5)] * 4

    inv_options_dv = InversionOptions()
    inv_options_dv.set_tool("neighpyI")
    inv_options_dv.set_params(
        bounds=bounds_dv,
        n_initial_samples=2000,
        n_samples_per_iteration=100,
        n_cells_to_resample=20,
        n_iterations=50,
    )

    inv_dv = Inversion(inv_problem_dv, inv_options_dv)
    inv_result_dv = inv_dv.run()
    inv_result_dv.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    NAI - Initial Random Search

    NAI - Optimisation Loop:   0%|          | 0/50 [00:00<?, ?it/s]
    NAI - Optimisation Loop:   2%|▏         | 1/50 [00:01<00:55,  1.13s/it]
    NAI - Optimisation Loop:   4%|▍         | 2/50 [00:01<00:26,  1.79it/s]
    NAI - Optimisation Loop:   6%|▌         | 3/50 [00:01<00:17,  2.67it/s]
    NAI - Optimisation Loop:   8%|▊         | 4/50 [00:01<00:13,  3.49it/s]
    NAI - Optimisation Loop:  10%|█         | 5/50 [00:01<00:10,  4.17it/s]
    NAI - Optimisation Loop:  12%|█▏        | 6/50 [00:01<00:09,  4.74it/s]
    NAI - Optimisation Loop:  14%|█▍        | 7/50 [00:02<00:08,  5.17it/s]
    NAI - Optimisation Loop:  16%|█▌        | 8/50 [00:02<00:07,  5.49it/s]
    NAI - Optimisation Loop:  18%|█▊        | 9/50 [00:02<00:07,  5.72it/s]
    NAI - Optimisation Loop:  20%|██        | 10/50 [00:02<00:06,  5.92it/s]
    NAI - Optimisation Loop:  22%|██▏       | 11/50 [00:02<00:06,  6.02it/s]
    NAI - Optimisation Loop:  24%|██▍       | 12/50 [00:02<00:06,  6.12it/s]
    NAI - Optimisation Loop:  26%|██▌       | 13/50 [00:03<00:05,  6.29it/s]
    NAI - Optimisation Loop:  28%|██▊       | 14/50 [00:03<00:05,  6.39it/s]
    NAI - Optimisation Loop:  30%|███       | 15/50 [00:03<00:05,  6.49it/s]
    NAI - Optimisation Loop:  32%|███▏      | 16/50 [00:03<00:05,  6.47it/s]
    NAI - Optimisation Loop:  34%|███▍      | 17/50 [00:03<00:05,  6.46it/s]
    NAI - Optimisation Loop:  36%|███▌      | 18/50 [00:03<00:04,  6.54it/s]
    NAI - Optimisation Loop:  38%|███▊      | 19/50 [00:03<00:04,  6.58it/s]
    NAI - Optimisation Loop:  40%|████      | 20/50 [00:04<00:04,  6.50it/s]
    NAI - Optimisation Loop:  42%|████▏     | 21/50 [00:04<00:04,  6.56it/s]
    NAI - Optimisation Loop:  44%|████▍     | 22/50 [00:04<00:04,  6.50it/s]
    NAI - Optimisation Loop:  46%|████▌     | 23/50 [00:04<00:04,  6.47it/s]
    NAI - Optimisation Loop:  48%|████▊     | 24/50 [00:04<00:04,  6.45it/s]
    NAI - Optimisation Loop:  50%|█████     | 25/50 [00:04<00:03,  6.51it/s]
    NAI - Optimisation Loop:  52%|█████▏    | 26/50 [00:04<00:03,  6.59it/s]
    NAI - Optimisation Loop:  54%|█████▍    | 27/50 [00:05<00:03,  6.65it/s]
    NAI - Optimisation Loop:  56%|█████▌    | 28/50 [00:05<00:03,  6.68it/s]
    NAI - Optimisation Loop:  58%|█████▊    | 29/50 [00:05<00:03,  6.59it/s]
    NAI - Optimisation Loop:  60%|██████    | 30/50 [00:05<00:03,  6.64it/s]
    NAI - Optimisation Loop:  62%|██████▏   | 31/50 [00:05<00:02,  6.64it/s]
    NAI - Optimisation Loop:  64%|██████▍   | 32/50 [00:05<00:02,  6.57it/s]
    NAI - Optimisation Loop:  66%|██████▌   | 33/50 [00:06<00:02,  6.59it/s]
    NAI - Optimisation Loop:  68%|██████▊   | 34/50 [00:06<00:02,  6.56it/s]
    NAI - Optimisation Loop:  70%|███████   | 35/50 [00:06<00:02,  6.53it/s]
    NAI - Optimisation Loop:  72%|███████▏  | 36/50 [00:06<00:02,  6.56it/s]
    NAI - Optimisation Loop:  74%|███████▍  | 37/50 [00:06<00:01,  6.53it/s]
    NAI - Optimisation Loop:  76%|███████▌  | 38/50 [00:06<00:01,  6.60it/s]
    NAI - Optimisation Loop:  78%|███████▊  | 39/50 [00:06<00:01,  6.52it/s]
    NAI - Optimisation Loop:  80%|████████  | 40/50 [00:07<00:01,  6.48it/s]
    NAI - Optimisation Loop:  82%|████████▏ | 41/50 [00:07<00:01,  6.46it/s]
    NAI - Optimisation Loop:  84%|████████▍ | 42/50 [00:07<00:01,  6.52it/s]
    NAI - Optimisation Loop:  86%|████████▌ | 43/50 [00:07<00:01,  6.49it/s]
    NAI - Optimisation Loop:  88%|████████▊ | 44/50 [00:07<00:00,  6.56it/s]
    NAI - Optimisation Loop:  90%|█████████ | 45/50 [00:07<00:00,  6.62it/s]
    NAI - Optimisation Loop:  92%|█████████▏| 46/50 [00:08<00:00,  6.56it/s]
    NAI - Optimisation Loop:  94%|█████████▍| 47/50 [00:08<00:00,  6.62it/s]
    NAI - Optimisation Loop:  96%|█████████▌| 48/50 [00:08<00:00,  6.68it/s]
    NAI - Optimisation Loop:  98%|█████████▊| 49/50 [00:08<00:00,  6.72it/s]
    NAI - Optimisation Loop: 100%|██████████| 50/50 [00:08<00:00,  6.66it/s]
    NAI - Optimisation Loop: 100%|██████████| 50/50 [00:08<00:00,  5.79it/s]
    ============================
    Summary for inversion result
    ============================
    SUCCESS
    ----------------------------
    model: [ 3.65446741  3.22530365  5.33884335  3.43241498 25.54124344  4.16511517
     41.47990624  4.49602258]
    samples: [[12.35888073  4.17986845 14.53017432 ...  3.29953543 27.92913793
       2.36780655]
     [54.77833827  3.25919094 54.25336229 ...  2.56894737 22.78114425
       3.13115686]
     [45.3030486   3.66195155 11.13520233 ...  3.19308944 55.49065324
       3.20351945]
     ...
     [ 3.65377737  3.22624301  5.32845836 ...  4.15493385 41.49363447
       4.49990418]
     [ 3.65387412  3.22655375  5.32909702 ...  4.15635468 41.49323695
       4.49834524]
     [ 3.65340715  3.22533198  5.32865265 ...  4.15686466 41.48756262
       4.49832081]]
    objectives: [4494.49996292 1122.94609476 1054.46175399 ...  494.62492545  494.87238816
      494.59962133]




.. GENERATED FROM PYTHON SOURCE LINES 587-597

.. code-block:: Python


    best_dv = inv_result_dv.model
    ds_samples_dv = inv_result_dv.samples
    ds_objectives_dv = inv_result_dv.objectives

    print("Best model (depth + velocity):")
    print(get_model_parameters(best_dv))
    print("\nTrue model:")
    print(good_model)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Best model (depth + velocity):
    [[ 3.65446741  3.22530365  1.7       ]
     [ 5.33884335  3.43241498  2.        ]
     [25.54124344  4.16511517  1.7       ]
     [41.47990624  4.49602258  1.7       ]]

    True model:
    [[ 1.   3.   1.7]
     [ 8.   3.2  2. ]
     [20.   4.   1.7]
     [45.   4.2  1.7]]




.. GENERATED FROM PYTHON SOURCE LINES 602-605

Convergence
^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 605-623

.. code-block:: Python


    best_i_dv = np.argmin(ds_objectives_dv)

    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(ds_objectives_dv, marker=".", linestyle="", markersize=2, color="black")
    ax.scatter(best_i_dv, ds_objectives_dv[best_i_dv], c="g", s=30, zorder=10,
               label="Best model")
    ax.axvline(2000, c="grey", ls="--", lw=0.8)
    ax.set_yscale("log")
    ax.set_xlabel("Sample number")
    ax.set_ylabel("Misfit")
    ax.set_title("Problem 2: NA-I convergence (depth + velocity)")
    ax.text(0.15, 0.95, "Initial random search", transform=ax.transAxes,
            ha="center", va="top")
    ax.text(0.65, 0.95, "Neighbourhood search", transform=ax.transAxes,
            ha="center", va="top")
    ax.legend()




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_007.png
   :alt: Problem 2: NA-I convergence (depth + velocity)
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_007.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcc1193390>



.. GENERATED FROM PYTHON SOURCE LINES 628-631

5.2 Stage II: Appraisal 
~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 631-635

.. code-block:: Python


    # Compute log-PPD for each Stage I sample (same choice as Problem 1)
    log_ppd_dv = -ds_objectives_dv








.. GENERATED FROM PYTHON SOURCE LINES 637-653

.. code-block:: Python


    # NA-II hyperparameters
    inv_options_dv_app = InversionOptions()
    inv_options_dv_app.set_tool("neighpyII")
    inv_options_dv_app.set_params(
        bounds=bounds_dv,
        initial_ensemble=ds_samples_dv,
        log_ppd=log_ppd_dv,
        n_resample=20000,
        n_walkers=10,
    )

    inv_dv_app = Inversion(inv_problem_dv, inv_options_dv_app)
    inv_result_dv_app = inv_dv_app.run()
    inv_result_dv_app.summary()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ============================
    Summary for inversion result
    ============================
    SUCCESS
    ----------------------------
    new_samples: [[ 3.79248677  3.22706369  9.79108    ...  4.18531913 57.20266037
       4.46681503]
     [ 8.3232985   3.22053145 14.96308353 ...  4.18659007 58.11078973
       4.45000291]
     [10.62230822  3.20394371  6.736847   ...  4.2095668  43.51925628
       4.45956822]
     ...
     [11.66904023  3.16648035 15.01326392 ...  4.37806961 48.76035271
       4.49182749]
     [16.96376464  3.29605069 12.21207934 ...  4.43404615 42.62371863
       4.3759516 ]
     [14.44362114  3.12602231  2.08165964 ...  4.48057032 55.05805663
       4.4281363 ]]




.. GENERATED FROM PYTHON SOURCE LINES 655-676

.. code-block:: Python


    appraisal_samples_dv = inv_result_dv_app.new_samples
    appraisal_mean_dv = appraisal_samples_dv.mean(axis=0)

    fig, ax = plt.subplots(figsize=(5, 8))
    for i in range(0, len(appraisal_samples_dv), 10):
        m = get_model_parameters(appraisal_samples_dv[i])
        plot_model_staircase(m, ax, color="k", alpha=0.01, lw=0.5)

    plot_model_staircase(good_model, ax, color="b", lw=2, label="True model")
    plot_model_staircase(get_model_parameters(best_dv), ax, color="g", lw=2,
                         ls="--", label="Best model (NA-I)")
    plot_model_staircase(get_model_parameters(appraisal_mean_dv), ax, color="r",
                         lw=2, ls="--", label="Mean model (NA-II)")
    ax.set_xlim(2.5, 5.0)
    ax.set_ylim(60, 0)
    ax.set_xlabel("Vs (km/s)")
    ax.set_ylabel("Depth (km)")
    ax.set_title("Problem 2: NA-II appraisal ensemble")
    ax.legend(loc="lower left")




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_008.png
   :alt: Problem 2: NA-II appraisal ensemble
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_008.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcb64f8410>



.. GENERATED FROM PYTHON SOURCE LINES 681-684

Corner plot
^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 684-724

.. code-block:: Python


    # 8x8 corner plot for depth+velocity problem
    true_dv = get_inversion_parameters(good_model)
    n_params_dv = 8
    var_names_dv = ["d1", "Vs1", "d2", "Vs2", "d3", "Vs3", "d4", "Vs4"]

    fig, axes = plt.subplots(n_params_dv, n_params_dv, figsize=(12, 12),
                             tight_layout=True,
                             gridspec_kw={"hspace": 0, "wspace": 0})
    for i in range(n_params_dv):
        for j in range(n_params_dv):
            axes[i, j].set_xlim(bounds_dv[j])
            axes[i, j].set_xticks([])
            axes[i, j].set_yticks([])
            if i == j:
                axes[i, j].hist(appraisal_samples_dv[::10, j], bins=100,
                                histtype="step", color="k")
                axes[i, j].axvline(best_dv[j], c="g", ls="--")
                axes[i, j].axvline(true_dv[j], c="b", ls="--")
                axes[i, j].axvline(appraisal_mean_dv[j], c="r", ls="--")
            elif j < i:
                axes[i, j].scatter(appraisal_samples_dv[::10, j],
                                   appraisal_samples_dv[::10, i],
                                   s=1, c="k", alpha=0.01)
                axes[i, j].scatter(best_dv[j], best_dv[i], c="g", s=10, zorder=10)
                axes[i, j].scatter(true_dv[j], true_dv[i], c="b", s=10, zorder=10)
                axes[i, j].scatter(appraisal_mean_dv[j], appraisal_mean_dv[i],
                                   c="r", s=10, zorder=10)
                axes[i, j].set_ylim(bounds_dv[i])
            else:
                axes[i, j].axis("off")

    for ax_, xlabel in zip(axes[-1, :], var_names_dv):
        ax_.set_xlabel(xlabel)
    for ax_, ylabel in zip(axes[:, 0], var_names_dv):
        if ylabel != var_names_dv[0]:
            ax_.set_ylabel(ylabel)

    fig.suptitle("Problem 2: Posterior corner plot (depth + velocity)")




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_009.png
   :alt: Problem 2: Posterior corner plot (depth + velocity)
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_009.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 0.98, 'Problem 2: Posterior corner plot (depth + velocity)')



.. GENERATED FROM PYTHON SOURCE LINES 729-732

Predicted data ensemble
^^^^^^^^^^^^^^^^^^^^^^^


.. GENERATED FROM PYTHON SOURCE LINES 732-755

.. code-block:: Python


    fig, ax = plt.subplots(figsize=(10, 4))

    for i in range(0, len(appraisal_samples_dv), 15):
        m = get_model_parameters(appraisal_samples_dv[i])
        try:
            t_pred, rf_pred = pyrf96.rfcalc(m)
            ax.plot(t_pred, rf_pred, color="k", alpha=0.01, lw=0.5)
        except Exception:
            pass

    ax.plot(t2, observed_data, color="k", lw=1.5, label="Observed", zorder=10)
    t_best, rf_best = pyrf96.rfcalc(get_model_parameters(best_dv))
    ax.plot(t_best, rf_best, color="g", ls="--", lw=1.5, label="Best model (NA-I)",
            zorder=11)
    t_mean, rf_mean = pyrf96.rfcalc(get_model_parameters(appraisal_mean_dv))
    ax.plot(t_mean, rf_mean, color="r", ls="--", lw=1.5, label="Mean model (NA-II)",
            zorder=11)
    ax.set_xlabel("Time (s)")
    ax.set_ylabel("Amplitude")
    ax.set_title("Problem 2: Predicted receiver functions from NA-II ensemble")
    ax.legend()




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_010.png
   :alt: Problem 2: Predicted receiver functions from NA-II ensemble
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_receiver_function_neighpy_010.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7efcca251bd0>



.. GENERATED FROM PYTHON SOURCE LINES 760-762

--------------


.. GENERATED FROM PYTHON SOURCE LINES 765-768

Watermark
---------


.. GENERATED FROM PYTHON SOURCE LINES 768-774

.. code-block:: Python


    watermark_list = ["cofi", "numpy", "scipy", "matplotlib"]
    for pkg in watermark_list:
        pkg_var = __import__(pkg)
        print(pkg, getattr(pkg_var, "__version__"))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    cofi 0.2.11+68.gc319bc9
    numpy 2.2.6
    scipy 1.17.1
    matplotlib 3.10.9




.. GENERATED FROM PYTHON SOURCE LINES 780-780

sphinx_gallery_thumbnail_number = -1


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (2 minutes 21.252 seconds)


.. _sphx_glr_download_examples_generated_scripts_synth_data_receiver_function_neighpy.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: receiver_function_neighpy.ipynb <receiver_function_neighpy.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: receiver_function_neighpy.py <receiver_function_neighpy.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: receiver_function_neighpy.zip <receiver_function_neighpy.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_