Manifold Markov chain Monte Carlo methods in Python
MIT License
Mici is a Python package providing implementations of Markov chain Monte Carlo (MCMC) methods for approximate inference in probabilistic models, with a particular focus on MCMC methods based on simulating Hamiltonian dynamics on a manifold.
Key features include
To install and use Mici the minimal requirements are a Python 3.10+ environment with NumPy and SciPy installed. The latest Mici release on PyPI (and its dependencies) can be installed in the current Python environment by running
pip install mici
To instead install the latest development version from the main
branch on Github run
pip install git+https://github.com/matt-graham/mici
If available in the installed Python environment the following additional packages provide extra functionality and features
arviz.InferenceData
container object usingarviz.convert_to_inference_data
or implicitly converted by passing thedata
argumentautograd.numpy
and autograd.scipy
interfaces). To sample chains inautograd
functions you also need to installmultiprocess.Pool
to be used in preference to the in-builtmutiprocessing.Pool
for parallelisation as multiprocess supportspip install mici[autograd]
.jax
pip install mici[jax]
.symnum.numpy
pip install mici[symnum]
.Mici is named for Augusta 'Mici' Teller, who along with Arianna Rosenbluth developed the code for the MANIAC I computer used in the seminal paper Equations of State Calculations by Fast Computing Machines which introduced the first example of a Markov chain Monte Carlo method.
Other Python packages for performing MCMC inference include PyMC, PyStan (the Python interface to Stan), Pyro / NumPyro, TensorFlow Probability, emcee, Sampyl and BlackJAX.
Unlike PyMC, PyStan, (Num)Pyro and TensorFlow Probability which are complete probabilistic programming frameworks including functionality for defining a probabilistic model / program, but like emcee, Sampyl and BlackJAX, Mici is solely focused on providing implementations of inference algorithms, with the user expected to be able to define at a minimum a function specifying the negative log (unnormalized) density of the distribution of interest.
Further while PyStan, (Num)Pyro and TensorFlow Probability all push the sampling loop into external compiled non-Python code, in Mici the sampling loop is run directly within Python. This has the consequence that for small models in which the negative log density of the target distribution and other model functions are cheap to evaluate, the interpreter overhead in iterating over the chains in Python can dominate the computational cost, making sampling much slower than packages which outsource the sampling loop to a efficient compiled implementation.
API documentation for the package is available
here. The three main user-facing
modules within the mici
package are the systems
, integrators
and
samplers
modules and you will generally need to create an instance of one
class from each module.
mici.systems
-
Hamiltonian systems encapsulating model functions and their derivatives
EuclideanMetricSystem
- systems with a metric on the position space withGaussianEuclideanMetricSystem
- systems in which the target distributionRiemannianMetricSystem
- systems with a metric on the position spaceSoftAbsRiemannianMetricSystem
- system with SoftAbsDenseConstrainedEuclideanMetricSystem
- Euclidean-metric system subjectmici.integrators
-
symplectic integrators for Hamiltonian dynamics
LeapfrogIntegrator
- explicit leapfrog (Strmer-Verlet) integrator forImplicitLeapfrogIntegrator
- implicit (or generalized) leapfrogImplicitMidpointIntegrator
- implicit midpointSymmetricCompositionIntegrator
- family of symplectic integrators for HamiltoniansConstrainedLeapfrogIntegrator
- constrained leapfrog integrator formici.samplers
- MCMC
samplers for peforming inference
StaticMetropolisHMC
- static integration time Hamiltonian Monte CarloRandomMetropolisHMC
- random integration time Hamiltonian Monte CarloDynamicSliceHMC
- dynamic integration time Hamiltonian Monte CarloDynamicMultinomialHMC
- dynamic integration time Hamiltonian Monte CarloThe manifold MCMC methods implemented in Mici have been used in several research projects. Below links are provided to a selection of Jupyter notebooks associated with these projects as demonstrations of how to use Mici and to illustrate some of the settings in which manifold MCMC methods can be computationally advantageous.
A simple complete example of using the package to compute approximate samples from a
distribution on a two-dimensional torus embedded in a three-dimensional space is given
below. The computed samples are visualized in the animation above. Here we use
SymNum to automatically construct functions to
calculate the required derivatives (gradient of negative log density of target
distribution and Jacobian of constraint function), sample four chains in parallel using
multiprocessing
, use ArviZ to calculate
diagnostics and use Matplotlib to plot the samples.
import mici
import numpy as np
import symnum
import symnum.numpy as snp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import arviz
# Define fixed model parameters
R = 1.0 # toroidal radius (0, )
r = 0.5 # poloidal radius (0, R)
= 0.9 # density fluctuation amplitude [0, 1)
# State dimension
dim_q = 3
# Define constraint function such that the set {q : constr(q) == 0} is a torus
@symnum.numpify(dim_q)
def constr(q):
x, y, z = q
return snp.array([((x**2 + y**2) ** 0.5 - R) ** 2 + z**2 - r**2])
# Define negative log density for the target distribution on torus
# (with respect to 2D 'area' measure for torus)
@symnum.numpify(dim_q)
def neg_log_dens(q):
x, y, z = q
= snp.arctan2(y, x)
= snp.arctan2(z, x / snp.cos() - R)
return snp.log1p(r * snp.cos() / R) - snp.log1p(snp.sin(4 * ) * snp.cos() * )
# Specify constrained Hamiltonian system with default identity metric
system = mici.systems.DenseConstrainedEuclideanMetricSystem(
neg_log_dens,
constr,
backend="symnum",
)
# System is constrained therefore use constrained leapfrog integrator
integrator = mici.integrators.ConstrainedLeapfrogIntegrator(system)
# Seed a random number generator
rng = np.random.default_rng(seed=1234)
# Use dynamic integration-time HMC implementation as MCMC sampler
sampler = mici.samplers.DynamicMultinomialHMC(system, integrator, rng)
# Sample initial positions on torus using parameterisation (, ) [0, 2)
# x, y, z = (R + r * cos()) * cos(), (R + r * cos()) * sin(), r * sin()
n_chain = 4
_init, _init = rng.uniform(0, 2 * np.pi, size=(2, n_chain))
q_init = np.stack(
[
(R + r * np.cos(_init)) * np.cos(_init),
(R + r * np.cos(_init)) * np.sin(_init),
r * np.sin(_init),
],
-1,
)
# Define function to extract variables to trace during sampling
def trace_func(state):
x, y, z = state.pos
return {"x": x, "y": y, "z": z}
# Sample 4 chains in parallel with 500 adaptive warm up iterations in which the
# integrator step size is tuned, followed by 2000 non-adaptive iterations
final_states, traces, stats = sampler.sample_chains(
n_warm_up_iter=500,
n_main_iter=2000,
init_states=q_init,
n_process=4,
trace_funcs=[trace_func],
)
# Print average accept probability and number of integrator steps per chain
for c in range(n_chain):
print(f"Chain {c}:")
print(f" Average accept prob. = {stats['accept_stat'][c].mean():.2f}")
print(f" Average number steps = {stats['n_step'][c].mean():.1f}")
# Print summary statistics and diagnostics computed using ArviZ
print(arviz.summary(traces))
# Visualize concatentated chain samples as animated 3D scatter plot
fig, ax = plt.subplots(
figsize=(4, 4),
subplot_kw={"projection": "3d", "proj_type": "ortho"},
)
(points_3d,) = ax.plot(*(np.concatenate(traces[k]) for k in "xyz"), ".", ms=0.5)
ax.axis("off")
for set_lim in [ax.set_xlim, ax.set_ylim, ax.set_zlim]:
set_lim((-1, 1))
def update(i):
angle = 45 * (np.sin(2 * np.pi * i / 60) + 1)
ax.view_init(elev=angle, azim=angle)
return (points_3d,)
anim = animation.FuncAnimation(fig, update, frames=60, interval=100)
plt.show()