mcmckit

mcmckit is a minimal, plug-and-play set of MCMC samplers for Python, built for Bayesian model updating and inverse problems in engineering: structural health monitoring, finite element model updating, structural dynamics.

The expensive part of model updating is your forward model. mcmckit is designed to stay out of its way.

---

You own the loop

Every sampler is a plain function that advances the chain by one step. State goes in as arguments and comes back as return values. Nothing is hidden on an object, so the recursion is yours to write, stop, inspect and modify.

Start with the simplest sampler, Metropolis-Hastings, which has a fixed proposal and no adaptation state at all:

import numpy as np
from mcmckit import mh_step

def log_post(theta):                  # your model goes here
    return -0.5 * np.sum(theta**2)    # a standard normal, for illustration

x = np.zeros(2)                       # current position
logp = log_post(x)                    # its log posterior
cov = np.eye(2) * 0.5                 # proposal covariance

chain = []                            # a list, so the run can be any length
n_accepted = 0

for i in range(10_000):
    x, logp, accepted = mh_step(log_post, x, logp, cov)
    chain.append(x)
    n_accepted += accepted

chain = np.array(chain)               # (n_iterations, n_parameters)
posterior = chain[1000:]              # burn-in is just a slice
print(posterior.mean(0), posterior.std(0))
print(f"acceptance rate: {n_accepted / len(chain):.2f}")

mh_step proposes a move, accepts or rejects it, and hands back the new position, its log posterior, and whether the move was taken. Everything else in this package is a variation on that.

Because the chain is a plain list you append to, nothing needs to know the run length in advance. Break out of the loop whenever your own criterion says so, and np.array whatever you collected.

The acceptance rate is your tuning signal. This one prints about 0.67, which is too high: the proposal is too small, so the chain accepts nearly everything and inches along. Widening cov would fix it. Roughly 0.2 to 0.4 is healthy for a random walk.

Plain MH needs a good cov, which you rarely know in advance. RAM learns one while sampling; the only change is that its adaptation state S is threaded through too, and it needs the step number:

from mcmckit import ram_step

x = np.zeros(2)
logp = log_post(x)
S = np.linalg.cholesky(np.eye(2) * 0.1**2)   # a rough guess is fine

chain = []
n_accepted = 0

for i in range(1, 10_001):                   # 1-indexed: drives adaptation
    x, logp, S, accepted = ram_step(log_post, x, logp, S, i)

    chain.append(x)
    n_accepted += accepted
    if i % 1000 == 0:                        # your convergence check
        print(i, x, logp)

chain = np.array(chain)
print(f"acceptance rate: {n_accepted / len(chain):.2f}")

Same shape, one extra threaded value. This prints about 0.25, against RAM's 0.234 target: starting from a proposal 5x too small, it found a sensible one on its own.

A caveat worth knowing early: a high acceptance rate is not a good sign. It usually means the steps are too small and the chain is crawling. Judge a run by effective sample size, not acceptance alone. Choosing a sampler shows what that looks like.

Everything is a plain module-level function, so import what you use and call it bare, or keep the package namespace if you prefer. Both are the same call:

from mcmckit import ram_step, dram_step      # bare
import mcmckit as mc                          # namespaced: mc.ram_step(...)

Function	Threaded state	Returns
mh_step	— (fixed cov)	x, logp, accepted
ram_step	S, step index i	x, logp, S, accepted
dram_step	DRAMState	x, logp, state, accepted
mala_step	grad	x, logp, grad, accepted
adaptive_mala_step	grad, log_step, i	x, logp, grad, log_step, accepted
gibbs_step	— (blocks, proposal_std)	x, logp, accepted_per_block

See Step functions for the full reference, and Your own loop for a worked structural example.

---

Or hand over the loop

When you do not need control of the recursion, the full-run helpers are thin loops over exactly the same step functions:

from mcmckit import ram

result = ram(log_post, x0=[0.0, 0.0], n_samples=10_000)

print(result.mean(), result.std())
result.discard(1000).plot_corner()

metropolis, ram, dram, mala, adaptive_mala and gibbs all follow this shape and return a Result with statistics and plots.

---

Samplers

Step function	Full run	Method	Notes
mh_step	metropolis	Random-walk MH	Fixed proposal covariance
mala_step	mala	Langevin	Needs the gradient
ram_step	ram	Robust Adaptive MH	Self-tunes covariance (Vihola 2012)
dram_step	dram	Delayed Rejection + AM	Best general-purpose adaptive sampler
adaptive_mala_step	adaptive_mala	Adaptive Langevin	Log-space step-size tuning
gibbs_step	gibbs	Metropolis-within-Gibbs	Block updates, per-block rates
—	TMCMC	Transitional MCMC	Prior→posterior bridge, log-evidence

TMCMC advances a whole population of particles per stage rather than a single chain position, so it has no single-step form and stays a class.

Not sure which to use? Choosing a sampler compares them all on the same problem, with traces, effective sample sizes and a hard target that separates them. Plots shows everything the package can draw.

---

Installation

pip install mcmckit            # core (numpy + scipy)
pip install mcmckit[plot]      # + matplotlib for plots

For development:

git clone https://github.com/LuigiCaglio/mcmckit
cd mcmckit
pip install -e ".[dev,plot]"

---

Package layout

mcmckit/
├── steps.py                 ← single-step functions: you own the loop
├── runners.py               ← full-run helpers: thin loops over steps.py
├── core/
│   ├── problem.py           ← Problem (prior + likelihood interface)
│   ├── result.py            ← Result (samples, plots, statistics)
│   └── noise.py             ← GaussianNoiseLikelihood
└── samplers/                ← stateful sampler classes, and TMCMC