emcee and Bayesian Software Tools

Contents

20. emcee and Bayesian Software Tools#

20.1. Review of emcee#

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

c_true = (5, 0.12, 11/4)
N = 50
x = np.sort(10*np.random.rand(N))
yerr = 1 + 5 * np.random.rand(N)
y = c_true[0] + c_true[1]*x + c_true[2]*x**2
y += yerr * np.random.randn(N)

plt.errorbar(x, y, yerr=yerr, fmt="k.");

../../_images/00eb52ea8e08b976208d4a87e0a00d72517cb47e95b4d67cb20747bc4dda9b6a.png

def log_likelihood(theta, x, y, yerr):
    theory = theta[0] + theta[1] * x + theta[2] * x**2
    return -0.5 * np.sum((y-theory) **2 / yerr**2 + np.log(yerr))

def log_prior(theta):
    if -10 < theta[0] < 10 and -10 < theta[1] < 10 and -10 < theta[2] < 10:
       return 0.0
    return - np.inf

def log_prob(theta, x, y, err):
    lp = log_prior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + log_likelihood(theta, x, y, yerr)

nwalkers, ndim = 32, 3
pos = [0, 0, 0] + 1e-3 * np.random.randn(nwalkers, ndim)

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=(x, y, yerr))
sampler.run_mcmc(pos, 5000, progress=True);

100%|██████████████████████████████████████████████████████████████████████████████| 5000/5000 [00:01<00:00, 2749.87it/s]

20.3. blobs: unlock sample stats, posterior predictive and miscellanea#

We will see how to do this first with emcee, but define these terms with pymce.

def log_prob_blobs(theta, x, y, err):
    lp = log_prior(theta)
    like_vector = log_likelihood(theta, x, y, yerr)
    like = np.sum(like_vector)
    return lp + like, like_vector

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob_blobs, args=(x, y, yerr))
sampler.run_mcmc(pos, 5000, progress=True);

100%|██████████████████████████████████████████████████████████████████████████████| 5000/5000 [00:02<00:00, 2238.72it/s]

20.5. Further Reading#

Inference Data - https://python.arviz.org/en/latest/getting_started/XarrayforArviZ.html
PyMC basics - https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/pymc_overview.html
Simple linear regression - https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/GLM_linear.html#glm-linear

%load_ext watermark
%watermark -untzvm -iv -w

Last updated: Thu Apr 17 2025 11:36:34CDT

Python implementation: CPython
Python version       : 3.11.12
IPython version      : 9.1.0

Compiler    : Clang 18.1.8 
OS          : Darwin
Release     : 24.4.0
Machine     : arm64
Processor   : arm
CPU cores   : 12
Architecture: 64bit

arviz     : 0.21.0
pandas    : 2.2.3
matplotlib: 3.10.1
rich      : 14.0.0
numpy     : 1.25.2
emcee     : 3.1.6
pymc      : 5.22.0

Watermark: 2.5.0