How to debug a model#
Introduction#
There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on.
Because PyMC uses PyTensor expressions to build the model, and not functions, there is no way to place a print statement into a likelihood function. Instead, you can use the pytensor.printing.Print class to print intermediate values.
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
%matplotlib inline
%config InlineBackend.figure_format = "retina"
RANDOM_SEED = 8927
How to print intermediate values of PyTensor functions#
Since PyTensor functions are compiled to C, you have to use pytensor.printing.Print class to print intermediate values (imported below as Print). Python print function will not work. Below is a simple example of using Print. For more information, see Debugging PyTensor.
import pytensor.tensor as pt
from pytensor import function
from pytensor.printing import Print
x = pt.dvector("x")
y = pt.dvector("y")
func = function([x, y], 1 / (x - y))
func([1, 2, 3], [1, 0, -1])
array([ inf, 0.5 , 0.25])
To see what causes the inf value in the output, we can print intermediate values of \((x-y)\) using Print. Print class simply passes along its caller but prints out its value along a user-define message:
x - y = __str__ = [0. 2. 4.]
array([ inf, 0.5 , 0.25])
Print reveals the root cause: \((x-y)\) takes a zero value when \(x=1, y=1\), causing the inf output.
How to capture Print output for further analysis#
When we expect many rows of output from Print, it can be desirable to redirect the output to a string buffer and access the values later on (thanks to Lindley Lentati for inspiring this example). Here is a toy example using Python print function:
import sys
from io import StringIO
old_stdout = sys.stdout
mystdout = sys.stdout = StringIO()
for i in range(5):
print(f"Test values: {i}")
output = mystdout.getvalue().split("\n")
sys.stdout = old_stdout # setting sys.stdout back
output
['Test values: 0',
'Test values: 1',
'Test values: 2',
'Test values: 3',
'Test values: 4',
'']
Troubleshooting a toy PyMC model#
rng = np.random.default_rng(RANDOM_SEED)
x = rng.normal(size=100)
with pm.Model() as model:
# priors
mu = pm.Normal("mu", mu=0, sigma=1)
sd = pm.Normal("sd", mu=0, sigma=1)
# setting out printing for mu and sd
mu_print = Print("mu")(mu)
sd_print = Print("sd")(sd)
# likelihood
obs = pm.Normal("obs", mu=mu_print, sigma=sd_print, observed=x)
pm.model_to_graphviz(model)
with model:
step = pm.Metropolis()
trace = pm.sample(5, step, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)
Only 5 samples in chain.
sd __str__ = 0.0
mu __str__ = 0.0
---------------------------------------------------------------------------
SamplingError Traceback (most recent call last)
Input In [9], in <cell line: 1>()
1 with model:
2 step = pm.Metropolis()
----> 3 trace = pm.sample(5, step, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)
File c:\users\igork\pycharmprojects\pymc\pymc\sampling.py:558, in sample(draws, step, init, n_init, initvals, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, **kwargs)
556 # One final check that shapes and logps at the starting points are okay.
557 for ip in initial_points:
--> 558 model.check_start_vals(ip)
559 _check_start_shape(model, ip)
561 sample_args = {
562 "draws": draws,
563 "step": step,
(...)
573 "discard_tuned_samples": discard_tuned_samples,
574 }
File c:\users\igork\pycharmprojects\pymc\pymc\model.py:1794, in Model.check_start_vals(self, start)
1791 initial_eval = self.point_logps(point=elem)
1793 if not all(np.isfinite(v) for v in initial_eval.values()):
-> 1794 raise SamplingError(
1795 "Initial evaluation of model at starting point failed!\n"
1796 f"Starting values:\n{elem}\n\n"
1797 f"Initial evaluation results:\n{initial_eval}"
1798 )
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'mu': array(0.), 'sd': array(0.)}
Initial evaluation results:
{'mu': -0.92, 'sd': -0.92, 'obs': -inf}
Exception handling of PyMC v4 has improved, so now SamplingError exception prints out the intermediate values of mu and sd which led to likelihood of -inf. However, this technique of printing intermediate values with aeasara.printing.Print can be valuable in more complicated cases.
Bringing it all together#
rng = np.random.default_rng(RANDOM_SEED)
y = rng.normal(loc=5, size=20)
old_stdout = sys.stdout
mystdout = sys.stdout = StringIO()
with pm.Model() as model:
mu = pm.Normal("mu", mu=0, sigma=10)
a = pm.Normal("a", mu=0, sigma=10, initval=0.1)
b = pm.Normal("b", mu=0, sigma=10, initval=0.1)
sd_print = Print("Delta")(a / b)
obs = pm.Normal("obs", mu=mu, sigma=sd_print, observed=y)
# limiting number of samples and chains to simplify output
trace = pm.sample(draws=10, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)
output = mystdout.getvalue()
sys.stdout = old_stdout # setting sys.stdout back
Only 10 samples in chain.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (1 chains in 1 job)
NUTS: [mu, a, b]
Sampling 1 chain for 0 tune and 10 draw iterations (0 + 10 draws total) took 1 seconds.
The chain contains only diverging samples. The model is probably misspecified.
The acceptance probability does not match the target. It is 0, but should be close to 0.8. Try to increase the number of tuning steps.
C:\Users\igork\AppData\Local\Temp\ipykernel_14804\1992602661.py:15: UserWarning: The number of samples is too small to check convergence reliably.
trace = pm.sample(draws=10, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)
output
'Delta __str__ = -85.74093608165128\nDelta __str__ = -9.182002291671038\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.315734173890055\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.312485782438435\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.314669656412736\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.31581619157038\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.315114719133609\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.31511040479387\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.314077394936474\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.313673830463395\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.31561025339713\nDelta __str__ = 0.10737295473067673\nDelta __str__ = -9.31526569370057\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\nDelta __str__ = 0.10737295473067673\n'
Raw output is a bit messy and requires some cleanup and formatting to convert to numpy.ndarray. In the example below regex is used to clean up the output, and then it is evaluated with eval to give a list of floats. Code below also works with higher-dimensional outputs (in case you want to experiment with different models).
import re
# output cleanup and conversion to numpy array
# this is code accepts more complicated inputs
pattern = re.compile("Delta __str__ = ")
output = re.sub(pattern, " ", output)
pattern = re.compile("\\s+")
output = re.sub(pattern, ",", output)
pattern = re.compile(r"\[,")
output = re.sub(pattern, "[", output)
output += "]"
output = "[" + output[1:]
output = eval(output)
output = np.array(output)
output
array([-85.74093608, -9.18200229, 0.10737295, 0.10737295,
0.10737295, -9.31573417, 0.10737295, -9.31248578,
0.10737295, -9.31466966, 0.10737295, -9.31581619,
0.10737295, -9.31511472, 0.10737295, -9.3151104 ,
0.10737295, -9.31407739, 0.10737295, -9.31367383,
0.10737295, -9.31561025, 0.10737295, -9.31526569,
0.10737295, 0.10737295, 0.10737295, 0.10737295,
0.10737295, 0.10737295, 0.10737295, 0.10737295,
0.10737295, 0.10737295])
Notice that we requested 5 draws, but got 34 sets of \(a/b\) values. The reason is that for each iteration, all proposed values are printed (not just the accepted values). Negative values are clearly problematic.
output.shape
(34,)
Watermark#
%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,xarray
Last updated: Tue Aug 02 2022
Python implementation: CPython
Python version : 3.10.5
IPython version : 8.4.0
pytensor: 2.7.5
xarray: 2022.3.0
matplotlib: 3.5.2
pytensor : 2.7.5
numpy : 1.23.0
arviz : 0.12.1
sys : 3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 06:57:19) [MSC v.1929 64 bit (AMD64)]
re : 2.2.1
pandas : 1.4.3
pymc : 4.1.2
Watermark: 2.3.1
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: