Note

Go to the end to download the full example code

2.2 - Including GemPy¶

Complex probabilistic model¶

import os
import numpy as np
import matplotlib.pyplot as plt
import torch
import pyro
import pyro.distributions as dist
from pyro.infer import MCMC, NUTS, Predictive
from pyro.infer.autoguide import init_to_mean

import gempy as gp
import gempy_engine
import gempy_viewer as gpv
from gempy_engine.core.backend_tensor import BackendTensor
import arviz as az
from gempy_probability.plot_posterior import default_red, default_blue
# sphinx_gallery_thumbnail_number = -1

Set the data path

data_path = os.path.abspath('../')

Define a function for plotting geological settings with wells

def plot_geo_setting_well(geo_model):
    """
    This function plots the geological settings along with the well locations.
    It uses gempy_viewer to create 2D plots of the model.
    """
    # Define well and device locations
    device_loc = np.array([[6e3, 0, 3700]])
    well_1 = 3.41e3
    well_2 = 3.6e3

    # Create a 2D plot
    p2d = gpv.plot_2d(geo_model, show_topography=False, legend=False, show=False)

    # Add well and device markers to the plot
    p2d.axes[0].scatter([3e3], [well_1], marker='^', s=400, c='#71a4b3', zorder=10)
    p2d.axes[0].scatter([9e3], [well_2], marker='^', s=400, c='#71a4b3', zorder=10)
    p2d.axes[0].scatter(device_loc[:, 0], device_loc[:, 2], marker='x', s=400, c='#DA8886', zorder=10)

    # Add vertical lines to represent wells
    p2d.axes[0].vlines(3e3, .5e3, well_1, linewidth=4, color='gray')
    p2d.axes[0].vlines(9e3, .5e3, well_2, linewidth=4, color='gray')

    # Show the plot
    p2d.fig.show()

Creating the Geological Model¶

Here we create a geological model using GemPy. The model defines the spatial extent, resolution, and geological information derived from orientations and surface points data.

geo_model = gp.create_geomodel(
    project_name='Wells',
    extent=[0, 12000, -500, 500, 0, 4000],
    resolution=[100, 2, 100],
    refinement=1,
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=data_path + "/data/2-layers/2-layers_orientations.csv",
        path_to_surface_points=data_path + "/data/2-layers/2-layers_surface_points.csv"
    )
)

Configuring the Model¶

We configure the interpolation options for the geological model. These options determine how the model interpolates between data points.

geo_model.interpolation_options.uni_degree = 0
geo_model.interpolation_options.mesh_extraction = False
geo_model.interpolation_options.sigmoid_slope = 1100.

Setting up a Custom Grid¶

A custom grid is set for the model, defining specific points in space where geological formations will be evaluated.

x_loc = 6000
y_loc = 0
z_loc = np.linspace(0, 4000, 100)
xyz_coord = np.array([[x_loc, y_loc, z] for z in z_loc])
gp.set_custom_grid(geo_model.grid, xyz_coord=xyz_coord)

Active grids: ['regular' 'custom']

<gempy.core.data.grid_modules.grid_types.CustomGrid object at 0x7fd4e3c53be0>

Plotting the Initial Geological Setting¶

Before running any probabilistic analysis, we first visualize the initial geological setting. This step ensures that our model is correctly set up with the initial data.

# Plot initial geological settings
plot_geo_setting_well(geo_model=geo_model)

Interpolating the Initial Guess¶

The model interpolates an initial guess for the geological formations. This step is crucial to provide a starting point for further probabilistic analysis.

gp.compute_model(
    gempy_model=geo_model,
    engine_config=gp.data.GemPyEngineConfig(backend=gp.data.AvailableBackends.numpy)
)
plot_geo_setting_well(geo_model=geo_model)

Setting Backend To: AvailableBackends.numpy
A size: (5, 5)
CG iterations: 5

Probabilistic Geomodeling with Pyro¶

In this section, we introduce a probabilistic approach to geological modeling. By using Pyro, a probabilistic programming language, we define a model that integrates geological data with uncertainty quantification.

sp_coords_copy = geo_model.interpolation_input.surface_points.sp_coords.copy()
# Change the backend to PyTorch for probabilistic modeling
BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH)

Setting Backend To: AvailableBackends.PYTORCH

Defining the Probabilistic Model¶

The Pyro model represents the probabilistic aspects of the geological model. It defines a prior distribution for the top layer’s location and computes the thickness of the geological layer as an observed variable.

def model(y_obs_list):
    """
    This Pyro model represents the probabilistic aspects of the geological model.
    It defines a prior distribution for the top layer's location and
    computes the thickness of the geological layer as an observed variable.
    """
    # Define prior for the top layer's location
    prior_mean = sp_coords_copy[0, 2]
    mu_top = pyro.sample(r'$\mu_{top}$', dist.Normal(prior_mean, torch.tensor(0.02, dtype=torch.float64)))

    # Update the model with the new top layer's location
    interpolation_input = geo_model.interpolation_input
    interpolation_input.surface_points.sp_coords = torch.index_put(
        interpolation_input.surface_points.sp_coords,
        (torch.tensor([0]), torch.tensor([2])),
        mu_top
    )

    # Compute the geological model
    geo_model.solutions = gempy_engine.compute_model(
        interpolation_input=interpolation_input,
        options=geo_model.interpolation_options,
        data_descriptor=geo_model.input_data_descriptor,
        geophysics_input=geo_model.geophysics_input,
    )

    # Compute and observe the thickness of the geological layer
    simulated_well = geo_model.solutions.octrees_output[0].last_output_center.custom_grid_values
    thickness = simulated_well.sum()
    pyro.deterministic(r'$\mu_{thickness}$', thickness.detach())
    y_thickness = pyro.sample(r'$y_{thickness}$', dist.Normal(thickness, 50), obs=y_obs_list)

Running Prior Sampling and Visualization¶

Prior sampling is an essential step in probabilistic modeling. It helps in understanding the distribution of our prior assumptions before observing any data.

Prepare observation data

y_obs_list = torch.tensor([200, 210, 190])

Run prior sampling and visualization

prior = Predictive(model, num_samples=50)(y_obs_list)
data = az.from_pyro(prior=prior)
az.plot_trace(data.prior)
plt.show()

$\mu_{top}$, $\mu_{top}$, $\mu_{thickness}$, $\mu_{thickness}$, $y_{thickness}$, $y_{thickness}$

Final Iteration 5, Residual Norm: 1.4089725943911736e-05
Final Iteration 5, Residual Norm: 6.317415998961143e-06
Final Iteration 5, Residual Norm: 1.468542842347646e-07
Final Iteration 5, Residual Norm: 9.099129314637124e-06
Final Iteration 5, Residual Norm: 2.667314595389484e-06
Final Iteration 5, Residual Norm: 1.4855313948945798e-07
Final Iteration 5, Residual Norm: 5.579229931795512e-05
Final Iteration 5, Residual Norm: 0.0002460598490319718
Final Iteration 5, Residual Norm: 4.7170726870053344e-06
Final Iteration 5, Residual Norm: 2.1329508670549636e-05
Final Iteration 5, Residual Norm: 0.00036439299934579183
Final Iteration 5, Residual Norm: 0.00048370251094574413
Final Iteration 5, Residual Norm: 1.0623705901568265e-06
Final Iteration 5, Residual Norm: 6.863468575162967e-05
Final Iteration 5, Residual Norm: 1.0331479933604487e-05
Final Iteration 5, Residual Norm: 1.049956003777276e-05
Final Iteration 5, Residual Norm: 0.0003624536197216722
Final Iteration 5, Residual Norm: 0.0006498274068455044
Final Iteration 5, Residual Norm: 0.00018471277567577628
Final Iteration 5, Residual Norm: 0.00040254824513041377
Final Iteration 5, Residual Norm: 0.0013596079008816946
Final Iteration 5, Residual Norm: 1.648709638209521e-05
Final Iteration 5, Residual Norm: 7.735022324916065e-06
Final Iteration 5, Residual Norm: 0.0001824839154379149
Final Iteration 5, Residual Norm: 0.0029643142652577897
Final Iteration 5, Residual Norm: 8.423194596482126e-06
Final Iteration 5, Residual Norm: 0.0005033355317047751
Final Iteration 5, Residual Norm: 3.466862162591784e-05
Final Iteration 5, Residual Norm: 4.100952832656466e-06
Final Iteration 5, Residual Norm: 5.186807534557172e-06
Final Iteration 5, Residual Norm: 7.225228101225706e-05
Final Iteration 5, Residual Norm: 0.00012718195369019127
Final Iteration 5, Residual Norm: 1.6491035955590957e-05
Final Iteration 5, Residual Norm: 7.977804546947701e-06
Final Iteration 5, Residual Norm: 6.816675195364853e-05
Final Iteration 5, Residual Norm: 2.4168551204555964e-07
Final Iteration 5, Residual Norm: 3.548664905916464e-08
Final Iteration 5, Residual Norm: 0.0022503253931455245
Final Iteration 5, Residual Norm: 0.0036943902745817376
Final Iteration 5, Residual Norm: 8.690035021310801e-05
Final Iteration 5, Residual Norm: 7.142135244993335e-08
Final Iteration 5, Residual Norm: 3.5356345481439774e-06
Final Iteration 5, Residual Norm: 7.538530566253355e-09
Final Iteration 5, Residual Norm: 3.4309134085127624e-07
Final Iteration 5, Residual Norm: 1.948481493782957e-08
Final Iteration 5, Residual Norm: 0.00014888343414485883
Final Iteration 5, Residual Norm: 4.704095188284034e-05
Final Iteration 5, Residual Norm: 0.00011895489197375922
Final Iteration 5, Residual Norm: 2.789264939712866e-06
Final Iteration 5, Residual Norm: 1.4900200944495394e-05
Final Iteration 5, Residual Norm: 3.218010200331566e-06
Final Iteration 5, Residual Norm: 1.0887753762451048e-06
/home/leguark/.virtualenvs/gempy-geotop-pilot/lib/python3.10/site-packages/arviz/utils.py:184: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  numba_fn = numba.jit(**self.kwargs)(self.function)

Sampling from the Posterior using MCMC¶

We use Markov Chain Monte Carlo (MCMC) with the NUTS (No-U-Turn Sampler) algorithm to sample from the posterior distribution. This gives us an understanding of the distribution of our model parameters after considering the observed data.

Run MCMC using NUTS to sample from the posterior

# Magic sauce
from gempy_engine.core.backend_tensor import BackendTensor
# import gempy_engine.config
# config.DEFAULT_PYKEOPS = False
BackendTensor._change_backend(engine_backend=gp.data.AvailableBackends.PYTORCH, dtype="float64", pykeops_enabled=False)

pyro.primitives.enable_validation(is_validate=True)
nuts_kernel = NUTS(model, step_size=0.0085, adapt_step_size=True, target_accept_prob=0.9, max_tree_depth=10, init_strategy=init_to_mean)
mcmc = MCMC(nuts_kernel, num_samples=200, warmup_steps=50)
mcmc.run(y_obs_list)

Setting Backend To: AvailableBackends.PYTORCH

Warmup:   0%|          | 0/250 [00:00, ?it/s]
Warmup:   0%|          | 1/250 [00:00,  3.12it/s, step size=3.22e-01, acc. prob=0.869]
Warmup:   2%|1         | 4/250 [00:00, 10.90it/s, step size=2.55e-03, acc. prob=0.472]
Warmup:   2%|2         | 6/250 [00:01,  4.24it/s, step size=2.14e-03, acc. prob=0.624]
Warmup:   3%|3         | 8/250 [00:02,  2.77it/s, step size=3.46e-03, acc. prob=0.717]
Warmup:   4%|3         | 9/250 [00:02,  3.16it/s, step size=2.47e-03, acc. prob=0.727]
Warmup:   4%|4         | 10/250 [00:03,  2.73it/s, step size=3.36e-03, acc. prob=0.754]
Warmup:   4%|4         | 11/250 [00:03,  2.81it/s, step size=4.63e-03, acc. prob=0.776]
Warmup:   5%|4         | 12/250 [00:03,  2.87it/s, step size=6.41e-03, acc. prob=0.795]
Warmup:   6%|5         | 14/250 [00:04,  3.31it/s, step size=4.38e-03, acc. prob=0.800]
Warmup:   6%|6         | 15/250 [00:04,  3.08it/s, step size=6.13e-03, acc. prob=0.814]
Warmup:   6%|6         | 16/250 [00:04,  3.41it/s, step size=8.56e-03, acc. prob=0.825]
Warmup:   7%|6         | 17/250 [00:04,  4.06it/s, step size=5.01e-03, acc. prob=0.819]
Warmup:   7%|7         | 18/250 [00:05,  4.26it/s, step size=6.53e-03, acc. prob=0.828]
Warmup:   8%|8         | 20/250 [00:05,  5.91it/s, step size=4.68e-03, acc. prob=0.828]
Warmup:   8%|8         | 21/250 [00:05,  5.69it/s, step size=6.44e-03, acc. prob=0.836]
Warmup:   9%|8         | 22/250 [00:05,  5.65it/s, step size=5.43e-03, acc. prob=0.836]
Warmup:   9%|9         | 23/250 [00:05,  5.48it/s, step size=7.42e-03, acc. prob=0.843]
Warmup:  10%|9         | 24/250 [00:06,  5.35it/s, step size=6.77e-03, acc. prob=0.843]
Warmup:  10%|#         | 25/250 [00:06,  5.31it/s, step size=9.25e-03, acc. prob=0.850]
Warmup:  11%|#         | 27/250 [00:06,  5.89it/s, step size=6.68e-03, acc. prob=0.848]
Warmup:  12%|#1        | 29/250 [00:06,  7.31it/s, step size=1.16e-02, acc. prob=0.858]
Warmup:  12%|#2        | 31/250 [00:06,  9.21it/s, step size=8.17e-03, acc. prob=0.856]
Warmup:  13%|#3        | 33/250 [00:07,  9.11it/s, step size=1.46e-02, acc. prob=0.864]
Warmup:  14%|#4        | 35/250 [00:07,  6.52it/s, step size=1.59e-03, acc. prob=0.842]
Warmup:  14%|#4        | 36/250 [00:08,  3.73it/s, step size=2.15e-03, acc. prob=0.846]
Warmup:  15%|#4        | 37/250 [00:09,  2.67it/s, step size=2.90e-03, acc. prob=0.850]
Warmup:  16%|#5        | 39/250 [00:09,  3.41it/s, step size=1.61e-03, acc. prob=0.846]
Warmup:  16%|#6        | 41/250 [00:09,  4.51it/s, step size=4.38e-04, acc. prob=0.835]
Warmup:  17%|#6        | 42/250 [00:09,  4.47it/s, step size=4.32e-04, acc. prob=0.836]
Warmup:  17%|#7        | 43/250 [00:11,  2.21it/s, step size=5.49e-04, acc. prob=0.840]
Warmup:  18%|#7        | 44/250 [00:13,  1.01it/s, step size=5.81e+00, acc. prob=0.843]
Warmup:  19%|#8        | 47/250 [00:13,  1.97it/s, step size=9.76e-02, acc. prob=0.796]
Warmup:  20%|#9        | 49/250 [00:14,  2.58it/s, step size=4.33e-02, acc. prob=0.781]
Warmup:  20%|##        | 50/250 [00:14,  2.29it/s, step size=4.33e-02, acc. prob=0.786]
Warmup:  20%|##        | 51/250 [00:16,  1.26it/s, step size=4.33e-02, acc. prob=1.000]
Sample:  21%|##        | 52/250 [00:20,  1.45s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  21%|##1       | 53/250 [00:21,  1.33s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  22%|##1       | 54/250 [00:23,  1.45s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  22%|##2       | 55/250 [00:24,  1.39s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  22%|##2       | 56/250 [00:25,  1.35s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  23%|##2       | 57/250 [00:27,  1.35s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  23%|##3       | 58/250 [00:27,  1.08s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  24%|##3       | 59/250 [00:29,  1.28s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  24%|##4       | 60/250 [00:29,  1.03s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  24%|##4       | 61/250 [00:30,  1.07s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  25%|##4       | 62/250 [00:31,  1.04s/it, step size=4.33e-02, acc. prob=1.000]
Sample:  25%|##5       | 63/250 [00:32,  1.27it/s, step size=4.33e-02, acc. prob=1.000]
Sample:  26%|##5       | 64/250 [00:32,  1.23it/s, step size=4.33e-02, acc. prob=1.000]
Sample:  26%|##6       | 65/250 [00:33,  1.38it/s, step size=4.33e-02, acc. prob=1.000]
Sample:  26%|##6       | 66/250 [00:33,  1.50it/s, step size=4.33e-02, acc. prob=1.000]
Sample:  27%|##6       | 67/250 [00:34,  1.40it/s, step size=4.33e-02, acc. prob=0.998]
Sample:  27%|##7       | 68/250 [00:36,  1.04s/it, step size=4.33e-02, acc. prob=0.998]
Sample:  28%|##7       | 69/250 [00:37,  1.17it/s, step size=4.33e-02, acc. prob=0.998]
Sample:  28%|##8       | 70/250 [00:37,  1.41it/s, step size=4.33e-02, acc. prob=0.998]
Sample:  28%|##8       | 71/250 [00:39,  1.02s/it, step size=4.33e-02, acc. prob=0.998]
Sample:  29%|##8       | 72/250 [00:40,  1.18s/it, step size=4.33e-02, acc. prob=0.997]
Sample:  29%|##9       | 73/250 [00:42,  1.34s/it, step size=4.33e-02, acc. prob=0.996]
Sample:  30%|##9       | 74/250 [00:42,  1.01it/s, step size=4.33e-02, acc. prob=0.996]
Sample:  30%|###       | 75/250 [00:43,  1.06s/it, step size=4.33e-02, acc. prob=0.996]
Sample:  30%|###       | 76/250 [00:44,  1.08it/s, step size=4.33e-02, acc. prob=0.996]
Sample:  31%|###       | 77/250 [00:45,  1.04s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  31%|###1      | 78/250 [00:47,  1.20s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  32%|###1      | 79/250 [00:49,  1.51s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  32%|###2      | 80/250 [00:49,  1.17s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  32%|###2      | 81/250 [00:51,  1.41s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  33%|###2      | 82/250 [00:52,  1.29s/it, step size=4.33e-02, acc. prob=0.994]
Sample:  33%|###3      | 83/250 [00:54,  1.36s/it, step size=4.33e-02, acc. prob=0.994]
Sample:  34%|###3      | 84/250 [00:54,  1.11s/it, step size=4.33e-02, acc. prob=0.994]
Sample:  34%|###4      | 85/250 [00:55,  1.15it/s, step size=4.33e-02, acc. prob=0.994]
Sample:  34%|###4      | 86/250 [00:55,  1.27it/s, step size=4.33e-02, acc. prob=0.994]
Sample:  35%|###5      | 88/250 [00:56,  1.61it/s, step size=4.33e-02, acc. prob=0.994]
Sample:  36%|###5      | 89/250 [00:57,  1.60it/s, step size=4.33e-02, acc. prob=0.995]
Sample:  36%|###6      | 90/250 [00:58,  1.38it/s, step size=4.33e-02, acc. prob=0.995]
Sample:  36%|###6      | 91/250 [00:58,  1.73it/s, step size=4.33e-02, acc. prob=0.995]
Sample:  37%|###6      | 92/250 [00:59,  1.27it/s, step size=4.33e-02, acc. prob=0.994]
Sample:  37%|###7      | 93/250 [01:01,  1.15s/it, step size=4.33e-02, acc. prob=0.995]
Sample:  38%|###7      | 94/250 [01:03,  1.27s/it, step size=4.33e-02, acc. prob=0.994]
Sample:  38%|###8      | 95/250 [01:03,  1.01s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  38%|###8      | 96/250 [01:05,  1.10s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  39%|###8      | 97/250 [01:06,  1.05s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  39%|###9      | 98/250 [01:06,  1.11it/s, step size=4.33e-02, acc. prob=0.993]
Sample:  40%|###9      | 99/250 [01:08,  1.12s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  40%|####      | 100/250 [01:09,  1.17s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  40%|####      | 101/250 [01:10,  1.24s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  41%|####      | 102/250 [01:11,  1.13s/it, step size=4.33e-02, acc. prob=0.993]
Sample:  41%|####1     | 103/250 [01:12,  1.10it/s, step size=4.33e-02, acc. prob=0.992]
Sample:  42%|####1     | 104/250 [01:13,  1.14it/s, step size=4.33e-02, acc. prob=0.991]
Sample:  42%|####2     | 105/250 [01:14,  1.02s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  42%|####2     | 106/250 [01:16,  1.24s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  43%|####2     | 107/250 [01:17,  1.37s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  43%|####3     | 108/250 [01:19,  1.37s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  44%|####3     | 109/250 [01:20,  1.33s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  44%|####4     | 110/250 [01:20,  1.01it/s, step size=4.33e-02, acc. prob=0.991]
Sample:  44%|####4     | 111/250 [01:21,  1.01it/s, step size=4.33e-02, acc. prob=0.991]
Sample:  45%|####4     | 112/250 [01:22,  1.08it/s, step size=4.33e-02, acc. prob=0.991]
Sample:  45%|####5     | 113/250 [01:23,  1.05s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  46%|####5     | 114/250 [01:26,  1.48s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  46%|####6     | 115/250 [01:29,  1.90s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  46%|####6     | 116/250 [01:31,  2.02s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  47%|####6     | 117/250 [01:32,  1.62s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  47%|####7     | 118/250 [01:32,  1.30s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  48%|####7     | 119/250 [01:33,  1.07s/it, step size=4.33e-02, acc. prob=0.991]
Sample:  48%|####8     | 120/250 [01:33,  1.18it/s, step size=4.33e-02, acc. prob=0.990]
Sample:  48%|####8     | 121/250 [01:33,  1.40it/s, step size=4.33e-02, acc. prob=0.990]
Sample:  49%|####8     | 122/250 [01:34,  1.27it/s, step size=4.33e-02, acc. prob=0.989]
Sample:  49%|####9     | 123/250 [01:35,  1.38it/s, step size=4.33e-02, acc. prob=0.989]
Sample:  50%|####9     | 124/250 [01:36,  1.29it/s, step size=4.33e-02, acc. prob=0.989]
Sample:  50%|#####     | 125/250 [01:37,  1.25it/s, step size=4.33e-02, acc. prob=0.989]
Sample:  50%|#####     | 126/250 [01:37,  1.54it/s, step size=4.33e-02, acc. prob=0.988]
Sample:  51%|#####     | 127/250 [01:37,  1.96it/s, step size=4.33e-02, acc. prob=0.988]
Sample:  51%|#####1    | 128/250 [01:39,  1.14it/s, step size=4.33e-02, acc. prob=0.988]
Sample:  52%|#####1    | 129/250 [01:42,  1.40s/it, step size=4.33e-02, acc. prob=0.989]
Sample:  52%|#####2    | 130/250 [01:42,  1.03s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  52%|#####2    | 131/250 [01:42,  1.27it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  53%|#####2    | 132/250 [01:42,  1.69it/s, step size=4.33e-02, acc. prob=0.983]
Sample:  53%|#####3    | 133/250 [01:44,  1.05it/s, step size=4.33e-02, acc. prob=0.983]
Sample:  54%|#####3    | 134/250 [01:46,  1.21s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  54%|#####4    | 135/250 [01:46,  1.06s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  54%|#####4    | 136/250 [01:49,  1.51s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  55%|#####4    | 137/250 [01:51,  1.66s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  55%|#####5    | 138/250 [01:52,  1.60s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  56%|#####5    | 139/250 [01:55,  1.85s/it, step size=4.33e-02, acc. prob=0.983]
Sample:  56%|#####6    | 140/250 [01:56,  1.59s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  56%|#####6    | 141/250 [01:58,  1.83s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  57%|#####6    | 142/250 [02:00,  1.82s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  57%|#####7    | 143/250 [02:00,  1.39s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  58%|#####7    | 144/250 [02:01,  1.07s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  58%|#####8    | 145/250 [02:02,  1.06s/it, step size=4.33e-02, acc. prob=0.984]
Sample:  58%|#####8    | 146/250 [02:02,  1.22it/s, step size=4.33e-02, acc. prob=0.984]
Sample:  59%|#####8    | 147/250 [02:02,  1.44it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  59%|#####9    | 148/250 [02:03,  1.68it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  60%|#####9    | 149/250 [02:03,  1.64it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  60%|######    | 150/250 [02:05,  1.16it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  60%|######    | 151/250 [02:07,  1.12s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  61%|######    | 152/250 [02:08,  1.08s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  61%|######1   | 153/250 [02:09,  1.16s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  62%|######1   | 154/250 [02:10,  1.20s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  62%|######2   | 155/250 [02:12,  1.36s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  62%|######2   | 156/250 [02:13,  1.17s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  63%|######2   | 157/250 [02:14,  1.26s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  63%|######3   | 158/250 [02:16,  1.56s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  64%|######3   | 159/250 [02:18,  1.45s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  64%|######4   | 160/250 [02:19,  1.55s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  64%|######4   | 161/250 [02:20,  1.34s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  65%|######4   | 162/250 [02:21,  1.14s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  65%|######5   | 163/250 [02:22,  1.03s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  66%|######5   | 164/250 [02:23,  1.10s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  66%|######6   | 165/250 [02:23,  1.07it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  66%|######6   | 166/250 [02:24,  1.10it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  67%|######6   | 167/250 [02:26,  1.02s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  67%|######7   | 168/250 [02:27,  1.26s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  68%|######7   | 169/250 [02:28,  1.06s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  68%|######8   | 170/250 [02:29,  1.09it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  68%|######8   | 171/250 [02:30,  1.19s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  69%|######8   | 172/250 [02:31,  1.06s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  69%|######9   | 173/250 [02:33,  1.32s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  70%|######9   | 174/250 [02:34,  1.32s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  70%|#######   | 175/250 [02:37,  1.72s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  70%|#######   | 176/250 [02:39,  1.91s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  71%|#######   | 177/250 [02:40,  1.53s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  71%|#######1  | 178/250 [02:41,  1.23s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  72%|#######1  | 179/250 [02:41,  1.09it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  72%|#######2  | 180/250 [02:42,  1.05s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  72%|#######2  | 181/250 [02:44,  1.20s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  73%|#######2  | 182/250 [02:44,  1.05s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  73%|#######3  | 183/250 [02:46,  1.30s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  74%|#######3  | 184/250 [02:48,  1.27s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  74%|#######4  | 185/250 [02:48,  1.15s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  74%|#######4  | 186/250 [02:49,  1.06s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  75%|#######4  | 187/250 [02:50,  1.11it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  75%|#######5  | 188/250 [02:51,  1.08it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  76%|#######5  | 189/250 [02:51,  1.29it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  76%|#######6  | 190/250 [02:53,  1.06it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  76%|#######6  | 191/250 [02:54,  1.04it/s, step size=4.33e-02, acc. prob=0.985]
Sample:  77%|#######6  | 192/250 [02:55,  1.17s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  77%|#######7  | 193/250 [02:56,  1.18s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  78%|#######7  | 194/250 [02:57,  1.11s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  78%|#######8  | 195/250 [02:58,  1.10s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  78%|#######8  | 196/250 [03:00,  1.11s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  79%|#######8  | 197/250 [03:01,  1.23s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  79%|#######9  | 198/250 [03:03,  1.31s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  80%|#######9  | 199/250 [03:05,  1.70s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  80%|########  | 200/250 [03:06,  1.39s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  80%|########  | 201/250 [03:06,  1.17s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  81%|########  | 202/250 [03:09,  1.70s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  81%|########1 | 203/250 [03:11,  1.60s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  82%|########1 | 204/250 [03:13,  1.69s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  82%|########2 | 205/250 [03:15,  1.76s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  82%|########2 | 206/250 [03:16,  1.79s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  83%|########2 | 207/250 [03:17,  1.51s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  83%|########3 | 208/250 [03:18,  1.40s/it, step size=4.33e-02, acc. prob=0.986]
Sample:  84%|########3 | 209/250 [03:19,  1.01s/it, step size=4.33e-02, acc. prob=0.985]
Sample:  84%|########4 | 211/250 [03:19,  1.75it/s, step size=4.33e-02, acc. prob=0.977]
Sample:  85%|########4 | 212/250 [03:19,  2.15it/s, step size=4.33e-02, acc. prob=0.975]
Sample:  85%|########5 | 213/250 [03:21,  1.12it/s, step size=4.33e-02, acc. prob=0.976]
Sample:  86%|########5 | 214/250 [03:24,  1.36s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  86%|########6 | 215/250 [03:26,  1.63s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  86%|########6 | 216/250 [03:28,  1.68s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  87%|########6 | 217/250 [03:28,  1.39s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  87%|########7 | 218/250 [03:29,  1.23s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  88%|########7 | 219/250 [03:30,  1.03s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  88%|########8 | 220/250 [03:32,  1.47s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  88%|########8 | 221/250 [03:34,  1.52s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  89%|########8 | 222/250 [03:35,  1.33s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  89%|########9 | 223/250 [03:36,  1.33s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  90%|########9 | 224/250 [03:37,  1.29s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  90%|######### | 225/250 [03:38,  1.24s/it, step size=4.33e-02, acc. prob=0.976]
Sample:  91%|######### | 227/250 [03:39,  1.44it/s, step size=4.33e-02, acc. prob=0.976]
Sample:  91%|#########1| 228/250 [03:39,  1.33it/s, step size=4.33e-02, acc. prob=0.977]
Sample:  92%|#########1| 229/250 [03:41,  1.16it/s, step size=4.33e-02, acc. prob=0.977]
Sample:  92%|#########2| 230/250 [03:42,  1.08it/s, step size=4.33e-02, acc. prob=0.977]
Sample:  92%|#########2| 231/250 [03:43,  1.06it/s, step size=4.33e-02, acc. prob=0.977]
Sample:  93%|#########2| 232/250 [03:44,  1.09s/it, step size=4.33e-02, acc. prob=0.977]
Sample:  93%|#########3| 233/250 [03:46,  1.30s/it, step size=4.33e-02, acc. prob=0.977]
Sample:  94%|#########3| 234/250 [03:47,  1.35s/it, step size=4.33e-02, acc. prob=0.977]
Sample:  94%|#########3| 235/250 [03:50,  1.75s/it, step size=4.33e-02, acc. prob=0.977]
Sample:  94%|#########4| 236/250 [03:52,  1.77s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  95%|#########4| 237/250 [03:53,  1.45s/it, step size=4.33e-02, acc. prob=0.977]
Sample:  95%|#########5| 238/250 [03:54,  1.42s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  96%|#########5| 239/250 [03:56,  1.45s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  96%|#########6| 240/250 [03:57,  1.43s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  96%|#########6| 241/250 [03:57,  1.09s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  97%|#########6| 242/250 [03:59,  1.30s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  97%|#########7| 243/250 [04:00,  1.25s/it, step size=4.33e-02, acc. prob=0.978]
Sample:  98%|#########7| 244/250 [04:00,  1.03it/s, step size=4.33e-02, acc. prob=0.978]
Sample:  98%|#########8| 245/250 [04:01,  1.06it/s, step size=4.33e-02, acc. prob=0.978]
Sample:  98%|#########8| 246/250 [04:02,  1.21it/s, step size=4.33e-02, acc. prob=0.978]
Sample:  99%|#########8| 247/250 [04:02,  1.44it/s, step size=4.33e-02, acc. prob=0.979]
Sample:  99%|#########9| 248/250 [04:04,  1.09it/s, step size=4.33e-02, acc. prob=0.979]
Sample: 100%|#########9| 249/250 [04:04,  1.19it/s, step size=4.33e-02, acc. prob=0.979]
Sample: 100%|##########| 250/250 [04:06,  1.00it/s, step size=4.33e-02, acc. prob=0.979]
Sample: 100%|##########| 250/250 [04:06,  1.02it/s, step size=4.33e-02, acc. prob=0.979]

Posterior Predictive Checks¶

After obtaining the posterior samples, we perform posterior predictive checks. This step is crucial to evaluate the performance and validity of our probabilistic model.

Sample from posterior predictive and visualize

posterior_samples = mcmc.get_samples()
posterior_predictive = Predictive(model, posterior_samples)(y_obs_list)
data = az.from_pyro(posterior=mcmc, prior=prior, posterior_predictive=posterior_predictive)
az.plot_trace(data)
plt.show()

$\mu_{top}$, $\mu_{top}$

/home/leguark/.virtualenvs/gempy-geotop-pilot/lib/python3.10/site-packages/arviz/data/io_pyro.py:157: UserWarning: Could not get vectorized trace, log_likelihood group will be omitted. Check your model vectorization or set log_likelihood=False
  warnings.warn(

Density Plot of Posterior Predictive¶

A density plot provides a visual representation of the distribution of the posterior predictive checks. It helps in comparing the prior and posterior distributions and in assessing the impact of our observed data on the model.

Plot density of posterior predictive and prior predictive

az.plot_density(
    data=[data.posterior_predictive, data.prior_predictive],
    shade=.9,
    var_names=[r'$\mu_{thickness}$'],
    data_labels=["Posterior Predictive", "Prior Predictive"],
    colors=[default_red, default_blue],
)
plt.show()

# sphinx_gallery_thumbnail_number = 2

$\mu_{thickness}$

Total running time of the script: ( 4 minutes 11.687 seconds)

Gallery generated by Sphinx-Gallery