Sampler Statistics

import arviz as az
import matplotlib.pyplot as plt
import pymc as pm

print(f"Running on PyMC v{pm.__version__}")

Running on PyMC v5.28.0+58.gf58491a3

az.style.use("arviz-variat")
plt.rcParams["figure.constrained_layout.use"] = False

When checking for convergence or when debugging a badly behaving sampler, it is often helpful to take a closer look at what the sampler is doing. For this purpose some samplers export statistics for each generated sample.

As a minimal example we sample from a standard normal distribution:

model = pm.Model()
with model:
    mu1 = pm.Normal("mu1", mu=0, sigma=1, shape=10)

with model:
    step = pm.NUTS()
    idata = pm.sample(2000, tune=1000, init=None, step=step, chains=4)

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu1]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 6 seconds.

• Note: NUTS provides the following statistics (these are internal statistics that the sampler uses, you don’t need to do anything with them when using PyMC, to learn more about them, pymc.NUTS.

idata.sample_stats

<xarray.DataTree 'sample_stats'>
Group: /sample_stats
    Dimensions:                (chain: 4, draw: 2000)
    Coordinates:
      * chain                  (chain) int64 32B 0 1 2 3
      * draw                   (draw) int64 16kB 0 1 2 3 4 ... 1996 1997 1998 1999
    Data variables: (12/18)
        acceptance_rate        (chain, draw) float64 64kB 0.5117 0.8243 ... 0.6885
        n_steps                (chain, draw) float64 64kB 3.0 7.0 7.0 ... 7.0 7.0
        max_energy_error       (chain, draw) float64 64kB 0.9907 0.5716 ... 0.7114
        divergences            (chain, draw) int64 64kB 0 0 0 0 0 0 ... 0 0 0 0 0 0
        index_in_trajectory    (chain, draw) int64 64kB -3 -4 -5 -3 3 ... 2 3 -4 6 4
        energy                 (chain, draw) float64 64kB 22.85 22.25 ... 17.12
        ...                     ...
        tree_depth             (chain, draw) int64 64kB 2 3 3 3 3 3 ... 3 2 3 3 3 3
        reached_max_treedepth  (chain, draw) bool 8kB False False ... False False
        step_size_bar          (chain, draw) float64 64kB 0.8498 0.8498 ... 0.8669
        lp                     (chain, draw) float64 64kB -16.29 -18.61 ... -11.36
        perf_counter_diff      (chain, draw) float64 64kB 0.0003056 ... 0.0003297
        perf_counter_start     (chain, draw) float64 64kB 37.75 37.75 ... 39.12
    Attributes:
        created_at:                 2026-04-25T12:14:28.293981+00:00
        creation_library:           ArviZ
        creation_library_version:   1.1.1dev0
        creation_library_language:  Python
        inference_library:          pymc
        inference_library_version:  5.28.0+58.gf58491a3
        sample_dims:                ['chain', 'draw']
        sampling_time:              5.952891111373901
        tuning_steps:               1000

xarray.DataTree

'sample_stats'

• Dimensions:
  • chain: 4
  • draw: 2000
• Coordinates: (2)
  • chain
    (chain)
    int64
    0 1 2 3

    array([0, 1, 2, 3])

  • draw
    (draw)
    int64
    0 1 2 3 4 ... 1996 1997 1998 1999

    array([   0,    1,    2, ..., 1997, 1998, 1999], shape=(2000,))

• Data variables: (18)
  • acceptance_rate
    (chain, draw)
    float64
    0.5117 0.8243 0.9935 ... 1.0 0.6885

    array([[0.51174226, 0.82433939, 0.99348419, ..., 0.49443367, 0.59231151,
            0.88179185],
           [1.        , 0.90708712, 0.93269626, ..., 0.82379309, 0.93972096,
            0.91513955],
           [0.63420517, 0.50988746, 0.81480311, ..., 0.71814159, 0.82955929,
            1.        ],
           [0.78426092, 1.        , 1.        , ..., 0.96975301, 1.        ,
            0.68849767]], shape=(4, 2000))

  • n_steps
    (chain, draw)
    float64
    3.0 7.0 7.0 7.0 ... 7.0 7.0 7.0 7.0

    array([[3., 7., 7., ..., 3., 3., 7.],
           [3., 3., 7., ..., 3., 7., 3.],
           [7., 7., 3., ..., 3., 7., 3.],
           [3., 7., 7., ..., 7., 7., 7.]], shape=(4, 2000))

  • max_energy_error
    (chain, draw)
    float64
    0.9907 0.5716 ... -0.4437 0.7114

    array([[ 0.990742  ,  0.57160393, -0.58722209, ...,  1.33365222,
             0.83333381, -0.89979905],
           [-0.98864914,  0.15413002, -0.62117296, ...,  0.75209228,
            -0.35072003, -0.48277515],
           [ 1.00156507,  1.15092182,  0.76886967, ...,  0.60626821,
             0.47703135, -0.50223453],
           [ 0.50828547, -0.51555179, -0.28212408, ..., -0.29894456,
            -0.44372241,  0.71136095]], shape=(4, 2000))

  • divergences
    (chain, draw)
    int64
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], shape=(4, 2000))

  • index_in_trajectory
    (chain, draw)
    int64
    -3 -4 -5 -3 3 3 2 ... 1 2 3 -4 6 4

    array([[-3, -4, -5, ...,  0, -2, -2],
           [-1,  2,  2, ...,  2,  2,  3],
           [-4,  1,  2, ...,  3, -2, -1],
           [ 1, -5,  4, ..., -4,  6,  4]], shape=(4, 2000))

  • energy
    (chain, draw)
    float64
    22.85 22.25 23.79 ... 15.92 17.12

    array([[22.85499102, 22.253824  , 23.79159306, ..., 18.36514144,
            17.83965557, 17.54161575],
           [19.03095328, 19.11227085, 18.07097589, ..., 22.62602352,
            19.74696462, 17.9359895 ],
           [19.26840482, 20.35574472, 23.11597786, ..., 16.70460636,
            20.66267519, 16.71046354],
           [23.99288576, 20.18043294, 20.04696409, ..., 21.05442521,
            15.92210829, 17.11867694]], shape=(4, 2000))

  • smallest_eigval
    (chain, draw)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

  • diverging
    (chain, draw)
    bool
    False False False ... False False

    array([[False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • energy_error
    (chain, draw)
    float64
    0.3979 0.4932 ... -0.4437 -0.05516

    array([[ 0.39792148,  0.49318996, -0.58722209, ...,  0.        ,
             0.83333381, -0.7861699 ],
           [-0.56053047,  0.15413002, -0.62117296, ..., -0.22828296,
            -0.23838603,  0.28225004],
           [-0.26734824,  0.80551207,  0.01931269, ...,  0.47273583,
             0.19432755, -0.46040154],
           [ 0.01873907, -0.19016432, -0.10244615, ..., -0.29894456,
            -0.44372241, -0.05515555]], shape=(4, 2000))

  • process_time_diff
    (chain, draw)
    float64
    0.0003153 0.0006517 ... 0.0003378

    array([[0.0003153 , 0.00065172, 0.00064331, ..., 0.00034069, 0.00034008,
            0.00068047],
           [0.00032311, 0.00032391, 0.00063884, ..., 0.00019314, 0.00034312,
            0.00017632],
           [0.00065287, 0.00060657, 0.0003152 , ..., 0.00021818, 0.00035878,
            0.00021013],
           [0.00032644, 0.00067573, 0.00067412, ..., 0.00033095, 0.00032562,
            0.00033776]], shape=(4, 2000))

  • step_size
    (chain, draw)
    float64
    0.8626 0.8626 ... 0.8859 0.8859

    array([[0.86259684, 0.86259684, 0.86259684, ..., 0.86259684, 0.86259684,
            0.86259684],
           [0.81031828, 0.81031828, 0.81031828, ..., 0.81031828, 0.81031828,
            0.81031828],
           [0.93527807, 0.93527807, 0.93527807, ..., 0.93527807, 0.93527807,
            0.93527807],
           [0.88589666, 0.88589666, 0.88589666, ..., 0.88589666, 0.88589666,
            0.88589666]], shape=(4, 2000))

  • largest_eigval
    (chain, draw)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], shape=(4, 2000))

  • tree_depth
    (chain, draw)
    int64
    2 3 3 3 3 3 2 3 ... 3 3 3 2 3 3 3 3

    array([[2, 3, 3, ..., 2, 2, 3],
           [2, 2, 3, ..., 2, 3, 2],
           [3, 3, 2, ..., 2, 3, 2],
           [2, 3, 3, ..., 3, 3, 3]], shape=(4, 2000))

  • reached_max_treedepth
    (chain, draw)
    bool
    False False False ... False False

    array([[False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False]], shape=(4, 2000))

  • step_size_bar
    (chain, draw)
    float64
    0.8498 0.8498 ... 0.8669 0.8669

    array([[0.84980907, 0.84980907, 0.84980907, ..., 0.84980907, 0.84980907,
            0.84980907],
           [0.85807006, 0.85807006, 0.85807006, ..., 0.85807006, 0.85807006,
            0.85807006],
           [0.86529605, 0.86529605, 0.86529605, ..., 0.86529605, 0.86529605,
            0.86529605],
           [0.86693771, 0.86693771, 0.86693771, ..., 0.86693771, 0.86693771,
            0.86693771]], shape=(4, 2000))

  • lp
    (chain, draw)
    float64
    -16.29 -18.61 ... -11.68 -11.36

    array([[-16.29004717, -18.60854869, -15.15277225, ..., -11.00227046,
            -15.72383988, -11.20869783],
           [-14.47616559, -15.12581827, -12.19859569, ..., -15.36282214,
            -14.21024487, -15.8718359 ],
           [-11.86282826, -16.32912968, -16.24914376, ..., -14.60302444,
            -15.2715556 , -12.77974001],
           [-17.2816567 , -16.12980302, -15.49714683, ..., -14.54270332,
            -11.68068426, -11.36290827]], shape=(4, 2000))

  • perf_counter_diff
    (chain, draw)
    float64
    0.0003056 0.0006314 ... 0.0003297

    array([[0.00030564, 0.00063144, 0.00062311, ..., 0.00033143, 0.00033087,
            0.00066195],
           [0.00031316, 0.00031381, 0.00061899, ..., 0.00018785, 0.00033379,
            0.00017152],
           [0.00063268, 0.00058776, 0.00030542, ..., 0.0002126 , 0.00034996,
            0.00020476],
           [0.00031625, 0.0006547 , 0.00065333, ..., 0.00032304, 0.00031783,
            0.00032966]], shape=(4, 2000))

  • perf_counter_start
    (chain, draw)
    float64
    37.75 37.75 37.75 ... 39.12 39.12

    array([[37.75302424, 37.7534525 , 37.75420829, ..., 38.76399762,
            38.76446004, 38.76492932],
           [38.03858025, 38.03902416, 38.03946503, ..., 39.07884391,
            39.07910813, 39.07951781],
           [37.96663349, 37.96739167, 37.96811278, ..., 39.2227741 ,
            39.22308952, 39.22351981],
           [37.71619748, 37.71668577, 37.71747385, ..., 39.11997091,
            39.1203649 , 39.12075404]], shape=(4, 2000))

• Attributes: (9)

  created_at :

  2026-04-25T12:14:28.293981+00:00

  creation_library :

  ArviZ

  creation_library_version :

  1.1.1dev0

  creation_library_language :

  Python

  inference_library :

  pymc

  inference_library_version :

  5.28.0+58.gf58491a3

  sample_dims :

  ['chain', 'draw']

  sampling_time :

  5.952891111373901

  tuning_steps :

  1000

The sample statistics variables are defined as follows:

• process_time_diff: The time it took to draw the sample, as defined by the python standard library time.process_time. This counts all the CPU time, including worker processes in BLAS and OpenMP.

• step_size: The current integration step size.

• diverging: (boolean) Indicates the presence of leapfrog transitions with large energy deviation from starting and subsequent termination of the trajectory. “large” is defined as max_energy_error going over a threshold.

• lp: The joint log posterior density for the model (up to an additive constant).

• energy: The value of the Hamiltonian energy for the accepted proposal (up to an additive constant).

• energy_error: The difference in the Hamiltonian energy between the initial point and the accepted proposal.

• perf_counter_diff: The time it took to draw the sample, as defined by the python standard library time.perf_counter (wall time).

• perf_counter_start: The value of time.perf_counter at the beginning of the computation of the draw.

• n_steps: The number of leapfrog steps computed. It is related to tree_depth with n_steps <= 2^tree_dept.

• max_energy_error: The maximum absolute difference in Hamiltonian energy between the initial point and all possible samples in the proposed tree.

• acceptance_rate: The average acceptance probabilities of all possible samples in the proposed tree.

• step_size_bar: The current best known step-size. After the tuning samples, the step size is set to this value. This should converge during tuning.

• tree_depth: The number of tree doublings in the balanced binary tree.

Some points to Note:

• Some of the sample statistics used by NUTS are renamed when converting to InferenceData to follow ArviZ’s naming convention, while some are specific to PyMC3 and keep their internal PyMC3 name in the resulting InferenceData object.

• InferenceData also stores additional info like the date, versions used, sampling time and tuning steps as attributes.

idata.sample_stats["tree_depth"].plot(col="chain", ls="none", marker=".", alpha=0.3);

az.plot_dist(
    idata,
    group="sample_stats",
    var_names="acceptance_rate",
    kind="hist",
    visuals={"credible_interval": False},
);

We check if there are any divergences, if yes, how many?

idata.sample_stats["diverging"].sum()

<xarray.DataArray 'diverging' ()> Size: 8B
array(0)

xarray.DataArray

'diverging'

• 0

  array(0)

In this case no divergences are found. If there are any, check this notebook for information on handling divergences.

It is often useful to compare the overall distribution of the energy levels with the change of energy between successive samples. Ideally, they should be very similar:

az.plot_energy(idata);

If the overall distribution of energy levels has longer tails, the efficiency of the sampler will deteriorate quickly.

Multiple samplers

If multiple samplers are used for the same model (e.g. for continuous and discrete variables), the exported values are merged or stacked along a new axis.

coords = {"step": ["BinaryMetropolis", "Metropolis"], "obs": ["mu1"]}
dims = {"accept": ["step"]}

with pm.Model(coords=coords) as model:
    mu1 = pm.Bernoulli("mu1", p=0.8)
    mu2 = pm.Normal("mu2", mu=0, sigma=1, dims="obs")

with model:
    step1 = pm.BinaryMetropolis([mu1])
    step2 = pm.Metropolis([mu2])
    idata = pm.sample(
        10000,
        init=None,
        step=[step1, step2],
        chains=4,
        tune=1000,
        idata_kwargs={"dims": dims, "coords": coords},
    )

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>BinaryMetropolis: [mu1]
>Metropolis: [mu2]

Sampling 4 chains for 1_000 tune and 10_000 draw iterations (4_000 + 40_000 draws total) took 17 seconds.

list(idata.sample_stats.data_vars)

['accepted', 'accept', 'scaling', 'p_jump']

Both samplers export accept, so we get one acceptance probability for each sampler:

az.plot_dist(
    idata,
    group="sample_stats",
    var_names="accept",
    kind="ecdf",
);

We notice that accept sometimes takes really high values (jumps from regions of low probability to regions of much higher probability).

# Range of accept values
idata.sample_stats["accept"].max("draw") - idata.sample_stats["accept"].min("draw")

<xarray.DataArray 'accept' (chain: 4, accept_dim_0: 2)> Size: 64B
array([[  3.75      , 153.17870672],
       [  3.75      , 734.08661568],
       [  3.75      , 127.786444  ],
       [  3.75      , 252.56104158]])
Coordinates:
  * chain         (chain) int64 32B 0 1 2 3
  * accept_dim_0  (accept_dim_0) int64 16B 0 1

xarray.DataArray

'accept'

• chain: 4
• accept_dim_0: 2

• 3.75 153.2 3.75 734.1 3.75 127.8 3.75 252.6

  array([[  3.75      , 153.17870672],
         [  3.75      , 734.08661568],
         [  3.75      , 127.786444  ],
         [  3.75      , 252.56104158]])

• Coordinates: (2)
  • chain
    (chain)
    int64
    0 1 2 3

    array([0, 1, 2, 3])

  • accept_dim_0
    (accept_dim_0)
    int64
    0 1

    array([0, 1])

Authors

• Updated by Meenal Jhajharia in April 2021 (pymc-examples#95)

• Updated to v4 by Christian Luhmann in May 2022 (pymc-examples#338)

• Updated by Osvaldo Martin in Dec 2025

• Updated by Osvaldo Martin in Apr 2026

Watermark

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Sat, 25 Apr 2026

Python implementation: CPython
Python version       : 3.14.4
IPython version      : 9.12.0

arviz     : 1.1.0
matplotlib: 3.10.8
pymc      : 5.28.0+58.gf58491a3

Watermark: 2.6.0

License notice

All the notebooks in this example gallery are provided under the MIT License which allows modification, and redistribution for any use provided the copyright and license notices are preserved.

Citing PyMC examples

To cite this notebook, use the DOI provided by Zenodo for the pymc-examples repository.

Important

Many notebooks are adapted from other sources: blogs, books… In such cases you should cite the original source as well.

Also remember to cite the relevant libraries used by your code.

Here is an citation template in bibtex:

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

which once rendered could look like: