Hamiltonian Monte Carlo (HMC)¶
Implementation of the Hybrid Monte Carlo (also known as Hamiltonian Monte Carlo) of Duane [23]. 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 [65].
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 Stochasticunlist()
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 aSampler{HMCTune}
type object.Example
Stand-Alone Function¶
-
sample!
(v::HMCVariate)¶ Draw one sample 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.
Value
Returnsv
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"]) epsilon = 0.1 L = 50 Sigma = Matrix{Float64}(I, 3, 3) theta1 = HMCVariate([0.0, 0.0, 0.0], epsilon, L, logfgrad) theta2 = HMCVariate([0.0, 0.0, 0.0], epsilon, L, Sigma, logfgrad) for i in 1:n sample!(theta1) sample!(theta2) 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¶
const HMCVariate = SamplerVariate{HMCTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::HMCTune
: tuning parameters for the sampling algorithm.
Constructors¶
-
HMCVariate
(x::AbstractVector{T<:Real}, epsilon::Real, L::Integer, logfgrad::Function)¶ -
HMCVariate
(x::AbstractVector{T<:Real}, epsilon::Real, L::Integer, Sigma::Matrix{U<:Real}, logfgrad::Function) Construct an
HMCVariate
object that stores simulated values and tuning parameters for HMC sampling.Arguments
x
: initial values.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.logfgrad
: function that takes a singleDenseVector
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 anHMCVariate
type object with fields set to the suppliedx
and tuning parameter values.
HMCTune Type¶
Declaration¶
type HMCTune <: SamplerTune
Fields¶
logfgrad::Nullable{Function}
: function supplied to the constructor to compute the log-transformed density and gradient vector, or null if not supplied.epsilon::Float64
: step size.L::Int
: number of steps to take in the Leapfrog algorithm.SigmaL::Union{UniformScaling{Int}, LowerTriangular{Float64}}
: Cholesky factorization of the covariance matrix for the multivariate normal proposal distribution.