Implementation of some simple adaptive MCMC algorithms
Author mvihola
2 Stars
Updated Last
11 Months Ago
Started In
August 2019


This package provides implementations of some general-purpose random-walk based adaptive MCMC algorithms, including the following:

The aim of the package is to provide a simple and modular general-purpose implementation, which may be easily used to sample from a log-target density, but also used in a variety of custom settings.

See also AdaptiveParticleMCMC.jl which uses this package with SequentialMonteCarlo.jl for adaptive particle MCMC.

Getting the package

using Pkg

Quick start

# Load the package
using AdaptiveMCMC

# Define a function which returns log-density values:
log_p(x) = -.5*sum(x.^2)

# Run 10k iterations of the Adaptive Metropolis:
out = adaptive_rwm(zeros(2), log_p, 10_000; algorithm=:am)

# Calculate '95% credible intervals':
using Statistics
mapslices(x->"$(mean(x)) ± $(1.96std(x))", out.X, dims=2)

The function adaptive_rwm ('Adaptive Random walk Metropolis') is a simple implenentation which does sampling for a given log-target density with the chosen method.

Adaptive parallel tempering

The adaptive_rwm also implements tempering, which is used if an optional argument L≥2 (number of temperature levels) is supplied. Here is a simple multimodal distribution sampled with APT:

# Multimodal target of dimension d.
function multimodalTarget(d::Int, sigma2=0.1^2, sigman=sigma2)
    # The means of mixtures
    m = [2.18 5.76; 3.25 3.47; 5.41 2.65; 4.93 1.50; 8.67 9.59;
         1.70 0.50; 2.70 7.88; 1.83 0.09; 4.24 8.48; 4.59 5.60;
         4.98 3.70; 2.26 0.31; 8.41 1.68; 6.91 5.81; 1.14 2.39;
         5.54 6.86; 3.93 8.82; 6.87 5.40; 8.33 9.50; 1.69 8.11]'
    n_m = size(m,2)
    @assert d>=2 "Dimension should be >= 2"
    let m=m, n_m=size(m,2), d=d
        function log_p(x::Vector{Float64})
            l_dens = -0.5*(mapslices(sum, (m.-x[1:2]).^2, dims=1)/sigma2)
            if d>2
                l_dens .-= 0.5*mapslices(sum, x[3:d].^2, dims=1)/sigman
            l_max = maximum(l_dens) # Prevent underflow by log-sum trick
            l_max + log(sum(exp.(l_dens.-l_max)))

using AdaptiveMCMC
n = 100_000; L = 2
rwm = adaptive_rwm(zeros(2), multimodalTarget(2), n; thin=10)
apt = adaptive_rwm(zeros(2), multimodalTarget(2), div(n,L); L = L, thin=10)

# Assuming you have 'Plots' installed:
using Plots
plot(scatter(rwm.X[1,:], rwm.X[2,:], title="w/o tempering", legend=:none),
scatter(apt.X[1,:], apt.X[2,:], title="w/ tempering", legend=:none), layout=(1,2))

Using with Distributions and LabelledArrays

MCMC is often useful with hierarchical models. These may be conveniently built using Distributions and LabelledArrays packages. The following example assumes these packages to be installed.

using Distributions, LabelledArrays, AdaptiveMCMC
# Define convenience log-transform for continuous univariate distributions
struct LogTransformedDistribution{Dist <: ContinuousUnivariateDistribution}
import Distributions.logpdf
logpdf(d::LogTransformedDistribution, x) = logpdf(d.d, exp(x)) + x
import Base.log
log(d::ContinuousUnivariateDistribution) = LogTransformedDistribution(d)

# This example is modified from Turing Getting Started:
function buildModel(x=0.0, y=1.0)
    let x=x, y=y
            p = 0.0
            p += logpdf(log(InverseGamma(2,3)), v.log_s)
            ss = exp(.5*v.log_s)
            p += logpdf(Normal(0, ss), v.m)
            p += logpdf(Normal(v.m, ss), x)
            p += logpdf(Normal(v.m, ss), y)

# Initial state vector (labelled with keys `s` and `m`)
x0 = LVector(log_s=1.0, m=0.0); log_p = buildModel(3.3, 4.14)
# Hint: If you do not have a good guess of the mode of log_p (which is
# a good initial value for MCMC), you may use optimisation:
#using Optim; o = optimize(x -> -log_p(x), x0); x0 = o.minimizer
out = adaptive_rwm(x0, log_p, 1_000_000; thin=20)
using StatsPlots # Assuming installed
corrplot(out.X', labels=[keys(x0)...])

Custom sampler

The package provides also simple building blocks which you can use within a 'custom' MCMC sampler. Here is an example:

using AdaptiveMCMC

# Sampler in R^d
function mySampler(log_p, n, x0)

    # Initialise random walk sampler state: r.x current state, r.y proposal
    r = RWMState(x0)

    # Initialise Adaptive Metropolis state (with default parameters)
    s = AdaptiveMetropolis(x0)
    # Other adaptations are: AdaptiveScalingMetropolis,
    # AdaptiveScalingWithinAdaptiveMetropolis, and RobustAdaptiveMetropolis

    X = zeros(eltype(x0), length(x0), n) # Allocate output storage
    p_x = log_p(r.x)                     # = log_p(x0); the initial log target
    for k = 1:n

        # Draw new proposal r.x -> r.y:
        draw!(r, s)

        p_y = log_p(r.y)                      # Calculate log target at proposal
        alpha = min(one(p_x), exp(p_y - p_x)) # The Metropolis acceptance probability

        if rand() <= alpha
            p_x = p_y

            # This 'accepts', or interchanges r.x <-> r.y:
            # (NB: do not do r.x = r.y; these are (pointers to) vectors!)

        # Do the adaptation update:
        adapt!(s, r, alpha, k)

        X[:,k] = r.x   # Save the current sample

# Standard normal target for testing
normal_log_p(x) = -mapreduce(e->e*e, +, x)/2

# Run 1M iterations of the sampler targetting 30d standard Normal:
X = mySampler(normal_log_p, 1_000_000, zeros(30))

More information & citations

The algorithms implemented in the package are discussed in the following reference:

  • M. Vihola. Ergonomic and reliable Bayesian inference with adaptive Markov chain Monte Carlo. In W. W. Piegrorsch, R. Levine, H. H. Zhang, and T. C. M. Lee, editors, Handbook of Computational Statistics and Data Science, Wiley, to appear.

If you use the package, please cite this publication.

Details of the implementation, such as about arguments and options, there are help fields written to the functions, which you may access, for instance, by typing ? adaptive_mcmc in the Julia REPL.

Required Packages

Used By Packages

No packages found.