BarkerMCMC.jl

Implements an adaptive Monte Carlo Markov Chain sampler that makes use of gradient information. It was the proposed by Livingstone et al. (2021).

The adaptative preconditioning is based on Andrieu and Thoms (2008), Algorithm 4 in Section 5. We followed the Algorithm 7.2 of the supporting information of Livingstone et al. (2021) with slight modifications.

You can find the repository with the source code here.

Usage

Here we sample from the 'banana-shaped' Rosenbruck function:

using BarkerMCMC

# --- Define target distribution and it's gradient
#     (or use automatic differentation)

function log_p_rosebruck_2d(x; k=1/200)
    -k*(100*(x[2] - x[1]^2)^2 + (1 - x[1])^2)
end

function ∇log_p_rosebruck_2d(x; k=1/200)
    [-2*k*(200*x[1]^3 - 200*x[1]*x[2] + x[1] -1),   # d/dx[1]
     -200*k*(x[2] - x[1]^2)]                        # d/dx[2]
end

# --- Generate samples

res = barker_mcmc(log_p_rosebruck_2d,
                  ∇log_p_rosebruck_2d,
                  [5.0, 5.0];
                  n_iter = 1_000,
                  target_acceptance_rate=0.4)

res.samples
res.log_p

# --- Visualize results

# acceptance rate
length(unique(res.samples[:,1])) / size(res.samples, 1)

# You may want to use `MCMCChains.jl` for plots and diagonstics
# (must be installed separately)

using MCMCChains
using StatsPlots

chain = Chains(res.samples, [:x1, :x2])
chains[200:10:end]                 # remove burn-in and apply thinning

plot(chains)
meanplot(chain)
histogram(chain)
autocorplot(chain)
corner(chain)

API

BarkerMCMC.barker_mcmcFunction

Adaptive MCMC sampler that makes use of gradient information. Based on Livingstone et al. (2020) The adaptation is based on Andrieu and Thoms (2008), Algorithm 4 in Section 5.

barker_mcmc(log_p::Function, ∇log_p::Function,
            inits::AbstractVector;
            n_iter = 100::Int,
            σ = 2.4/(length(inits)^(1/6)),
            target_acceptance_rate = 0.4,
            κ::Float64 = 0.6,
            n_iter_adaptation = Inf,
            preconditioning::Function = BarkerMCMC.precond_eigen)

Arguments

  • log_p::Function: function returning the log of the (non-normalized) density
  • ∇log_p::Function: function returning the gradient of log of the (non-normalized) density
  • inits::Vector: initial starting values
  • n_iter = 100: number of iterations
  • σ = 2.4/(length(inits)^(1/6)): global scale of proposal distribution
  • target_acceptance_rate = 0.4: desired acceptance rate
  • κ = 0.6: controls adaptation speed, κ ∈ (0.5, 1). Larger values lead to slower adaptation, see Section 6.1 in Livingstone et al. (2020).
  • n_iter_adaptation = Inf: number of iterations with adaptation
  • preconditioning::Function = BarkerMCMC.precond_eigen: Either BarkerMCMC.precond_eigen or BarkerMCMC.precond_cholesky. Calculating the preconditioning matrix with a cholesky decomposition is slighly cheaper, however, the eigen value decomposition allows for a proper rotation of the proposal distribution.

Return Value

A named tuple with fields:

  • samples: array containing the samples
  • log_p: vector containing the value of log_p for each sample.

References

Andrieu, C., Thoms, J., 2008. A tutorial on adaptive MCMC. Statistics and computing 18, 343–373.

Livingstone, S., Zanella, G., 2021. The Barker proposal: Combining robustness and efficiency in gradient-based MCMC. Journal of the Royal Statistical Society: Series B (Statistical Methodology). https://doi.org/10.1111/rssb.12482

BarkerMCMC.precond_choleskyFunction

Given a covariance matrix Σ, computes the preconditioning matrix M based on cholesky decomposition.

For M holds that cov(M * z) == Σ, where z a uncorrelated vector of random variables with zero mean.

BarkerMCMC.precond_eigenFunction

Given a covariance matrix Σ, computes the preconditioning matrix M based on eigen value decomposition.

For M holds that cov(M * z) == Σ, where z a uncorrelated vector of random variables with zero mean.

Hamiltonian Monte Carlo (gradient based)

Adaptive MCMC (without gradient)

Literature

Andrieu, C., Thoms, J., 2008. A tutorial on adaptive MCMC. Statistics and computing 18, 343–373.

Livingstone, S., Zanella, G., 2021. The Barker proposal: Combining robustness and efficiency in gradient-based MCMC. Journal of the Royal Statistical Society: Series B (Statistical Methodology). https://doi.org/10.1111/rssb.12482