mcmckit

Documentation

Minimal, plug-and-play MCMC samplers for Python.

Built for Bayesian model updating and inverse problems in engineering, where the expensive part is your forward model and you want the sampler to stay out of the way. Depends only on NumPy and SciPy.

Install

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]"

Quick start — 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.

1. Metropolis-Hastings — the whole idea, start to finish

The simplest sampler, with a fixed proposal covariance and no adaptation. If you read one example, read this one:

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. That is the entire interface. Everything else in this package is a variation on it.

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.

2. RAM — the same loop, tuning itself

Plain MH needs a good cov to work well, and you rarely know one in advance. RAM learns it while sampling. The only change is that its adaptation state S is threaded through alongside everything else, and it needs the step number:

import numpy as np
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, your rules
        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, which is the whole point of using it.

One 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 (result.ess()), not by acceptance alone. See Choosing a sampler for what that looks like in practice.

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 mh_step, ram_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

Sign convention

Every function takes a log-posterior: bigger is a better fit. If your code carries a negative log-posterior, wrap it once:

log_post = lambda theta: -my_nll(theta)

Keeping your forward model's output

If your callable returns (log_post, aux), the extra payload rides along and comes back attached to the accepted sample. Natural frequencies, mode shapes, residuals — whatever your model already computed, without re-running it:

def log_post(theta):
    freqs = surrogate(theta)
    return -0.5 * np.sum(((freqs - measured) / sigma)**2), freqs

x, logp, S, accepted, freqs = ram_step(log_post, x, logp, S, i, aux=freqs)

Return a plain float instead and no aux comes back, so the signature stays its usual width.

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(title="Posterior")

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

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 works on a population of particles per stage rather than one chain position, so it does not have a single-step form. It 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 draws.

Stateful sampler objects

The original class interface is still there for stop-and-inspect workflows, and is unchanged:

sampler = mc.RAM(n_samples=10_000, initial_cov=np.eye(2))
sampler.initialize(problem, x0=[0.0, 0.0])
sampler.step()
result = sampler.get_result()

These take a Problem, which bundles a prior and likelihood:

problem = mc.Problem(prior=log_prior, likelihood=log_likelihood,
                     param_names=["E", "zeta"])

The full-run helpers accept either a Problem or a bare log_post callable.

TMCMC (Transitional MCMC)

TMCMC bridges the prior to the posterior through tempered distributions and provides a log-evidence (log marginal likelihood) estimate:

prior_samples = np.random.uniform(-10, 10, size=(1000, 2))

tmcmc = mc.TMCMC(n_particles=1000, n_mcmc_steps=3)
result = tmcmc.run(problem, prior_samples=prior_samples)

print(f"log-evidence: {result.log_evidence:.3f}")
tmcmc.plot_stages(max_stages=6, title="Prior → Posterior")

Parallel evaluation

Opt-in, off by default. The cost in model updating is the forward model inside your likelihood, so mcmckit can spread those calls over cores:

if __name__ == "__main__":                      # required: workers re-import the module
    result = mc.TMCMC(n_particles=1000, n_workers=4).run(problem, prior_samples=ps)
    multi = mc.run_chains(sampler, problem, x0, n_chains=4, n_workers=4)

n_workers=-1 uses one worker per core. TMCMC particles and independent chains are parallelised; a single chain is sequential by construction and is not.

Process workers pickle your likelihood, so it must be a module-level function rather than a lambda or closure; mcmckit says so clearly instead of failing deep inside the executor. TMCMC gives bit-identical results with and without workers.

On a 14-core machine with a ~4 ms likelihood: 2.60x on 4 workers, 3.36x on 8.

Black-box solvers work — examples/openseespy_parallel.py updates a 6-storey shear building through OpenSeesPy and gets bit-identical results at 1, 4 and 8 workers. Use processes, not threads, for any solver with global state: OpenSeesPy keeps one global model domain, and threading it raises OpenSeesError or segfaults. See the docs.

Visualisation

Result objects have built-in plotting:

result.plot_trace()
result.plot_marginals()
result.plot_corner(style="corner", true_values=[2.0, -1.0])
# styles: "corner" | "scatter" | "full" | "kde"

Driving the loop yourself gives you a plain array, which you can wrap when you want those plots:

from mcmckit import Result

result = Result(samples=chain, param_names=["E", "zeta"])
result.plot_corner()

Examples

Script	Demonstrates
examples/own_loop.py	Start here. Step functions, your own loop, early stopping, aux output
examples/simple_gaussian.py	MH, all corner styles, burn-in, raw samples
examples/mala_vs_mh.py	MH vs MALA side-by-side
examples/ram_example.py	RAM self-tuning from bad initial covariance
examples/adaptive_samplers.py	MH / RAM / DRAM / AdaptiveMALA comparison
examples/tmcmc_and_gibbs.py	TMCMC stages + Gibbs scalar/block
examples/diagnostics_multichain.py	Multi-chain runs, R-hat and ESS diagnostics
examples/noise_estimation.py	Inferring measurement noise alongside parameters
examples/structural_identification.py	Stiffness identification from modal data
examples/hierarchical_updating.py	Hierarchical models across multiple structures
examples/model_comparison.py	Bayes factors from TMCMC log-evidence
examples/model_averaging.py	Posterior model averaging
examples/sequential.py	Sequential / online updating with PosteriorPrior
examples/sequential_10dof.py	10-parameter sequential structural identification
examples/openseespy_parallel.py	Black-box OpenSeesPy forward model, run in parallel

Citing

If you use mcmckit in academic work, please cite it:

Caglio, L. mcmckit: minimal, plug-and-play MCMC samplers for Python. https://doi.org/10.5281/zenodo.22300143

That DOI always resolves to the newest release. To cite the exact version you ran, use its own DOI from the Zenodo record; v0.3.0 is 10.5281/zenodo.22300144. CITATION.cff carries the same metadata, and GitHub's "Cite this repository" button renders it as BibTeX or APA.

Requirements

• Python ≥ 3.9
• numpy, scipy
• matplotlib (optional, for plots)