Hamiltonian Monte Carlo (HMC)

Implementation of the Hybrid Monte Carlo (also known as Hamiltonian Monte Carlo) of Duane [21]. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality. Code is derived from Neal’s implementation [62].

Model-Based Constructors

HMC(params::ElementOrVector{Symbol}, epsilon::Real, L::Integer; dtype::Symbol=:forward)
HMC(params::ElementOrVector{Symbol}, epsilon::Real, L::Integer, Sigma::Matrix{T<:Real}; dtype::Symbol=:forward)

Construct a Sampler object for HMC sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.

Arguments

  • params : stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic unlist() function.

  • epsilon : step size.

  • L : number of steps to take in the Leapfrog algorithm.

  • Sigma : covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.

  • dtype : type of differentiation for gradient calculations. Options are
    • :central : central differencing.
    • :forward : forward differencing.

Value

Returns a Sampler{HMCTune} type object.

Example

See the Dyes and other Examples.

Stand-Alone Functions

hmc!(v::HMCVariate, epsilon::Real, L::Integer, logfgrad::Function)
hmc!(v::HMCVariate, epsilon::Real, L::Integer, SigmaF::Cholesky{Float64}, logfgrad::Function)

Simulate one draw from a target distribution using the HMC sampler. Parameters are assumed to be continuous and unconstrained.

Arguments

  • v : current state of parameters to be simulated.
  • epsilon : step size.
  • L : number of steps to take in the Leapfrog algorithm.
  • SigmaF : Cholesky factorization of the covariance matrix for the multivariate normal proposal distribution. If omitted, the identity matrix is assumed.
  • logfgrad : function that takes a single DenseVector argument at which to compute the log-transformed density (up to a normalizing constant) and gradient vector, and returns the respective results as a tuple.

Value

Returns v updated with simulated values and associated tuning parameters.

Example

The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.

################################################################################
## Linear Regression
##   y ~ N(b0 + b1 * x, s2)
##   b0, b1 ~ N(0, 1000)
##   s2 ~ invgamma(0.001, 0.001)
################################################################################

using Mamba

## Data
data = Dict(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)

## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
logfgrad = function(x::DenseVector)
  b0 = x[1]
  b1 = x[2]
  logs2 = x[3]
  r = data[:y] - b0 - b1 * data[:x]
  logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
           (0.5 * dot(r, r) + 0.001) / exp(logs2) -
           0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
  grad = [
    sum(r) / exp(logs2) - b0 / 1000,
    sum(data[:x] .* r) / exp(logs2) - b1 / 1000,
    -0.5 * length(data[:y]) - 0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2)
  ]
  logf, grad
end

## MCMC Simulation with Hamiltonian Monte Carlo
## Without (1) and with (2) a user-specified proposal covariance matrix
n = 5000
sim1 = Chains(n, 3, names = ["b0", "b1", "s2"])
sim2 = Chains(n, 3, names = ["b0", "b1", "s2"])
theta1 = HMCVariate([0.0, 0.0, 0.0])
theta2 = HMCVariate([0.0, 0.0, 0.0])
epsilon = 0.1
L = 50
SigmaF = cholfact(eye(3))
for i in 1:n
  hmc!(theta1, epsilon, L, logfgrad)
  hmc!(theta2, epsilon, L, SigmaF, logfgrad)
  sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])]
  sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])]
end
describe(sim1)
describe(sim2)

HMCVariate Type

Declaration

typealias HMCVariate SamplerVariate{HMCTune}

Fields

  • value::Vector{Float64} : simulated values.
  • tune::HMCTune : tuning parameters for the sampling algorithm.

Constructors

HMCVariate(x::AbstractVector{T<:Real})
HMCVariate(x::AbstractVector{T<:Real}, tune::HMCTune)

Construct a HMCVariate object that stores simulated values and tuning parameters for HMC sampling.

Arguments

  • x : simulated values.
  • tune : tuning parameters for the sampling algorithm. If not supplied, parameters are set to their defaults.

Value

Returns a HMCVariate type object with fields set to the values supplied to arguments x and tune.

HMCTune Type

Declaration

type HMCTune <: SamplerTune

Fields

  • epsilon::Float64 : step size.
  • L::Int : number of steps to take in the Leapfrog algorithm.
  • SigmaF::Cholesky{Float64} : Cholesky factorization of the covariance matrix for the multivariate normal proposal distribution.