Adaptive Metropolis within Gibbs (AMWG)

Implementation of a Metropolis-within-Gibbs sampler [51][64][74] for iteratively simulating autocorrelated draws from a distribution that can be specified up to a constant of proportionality.

Stand-Alone Function

amwg!(v::AMWGVariate, sigma::Vector{Float64}, logf::Function; adapt::Bool=true, batchsize::Integer=50, target::Real=0.44)

Simulate one draw from a target distribution using an adaptive Metropolis-within-Gibbs sampler. Parameters are assumed to be continuous and unconstrained.

Arguments

  • v : current state of parameters to be simulated. When running the sampler in adaptive mode, the v argument in a successive call to the function should contain the tune field returned by the previous call.
  • sigma : initial standard deviations for the univariate normal proposal distributions.
  • logf : function to compute the log-transformed density (up to a normalizing constant) at v.value.
  • adapt : whether to adaptively update the proposal distribution.
  • batchsize : number of samples that must be newly accumulated before applying an adaptive update to the proposal distributions.
  • target : a target acceptance rate for the adaptive algorithm.

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 = [
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
]

## Log-transformed Posterior(b0, b1, log(s2)) + Constant
logf = function(x)
   b0 = x[1]
   b1 = x[2]
   logs2 = x[3]
   r = data[:y] - b0 - b1 * data[:x]
   (-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
end

## MCMC Simulation with Adaptive Metopolis-within-Gibbs Sampling
n = 5000
burnin = 1000
sim = Chains(n, 3, names = ["b0", "b1", "s2"])
theta = AMWGVariate([0.0, 0.0, 0.0])
sigma = ones(3)
for i in 1:n
  amwg!(theta, sigma, logf, adapt = (i <= burnin))
  sim[i,:,1] = [theta[1:2], exp(theta[3])]
end
describe(sim)

AMWGVariate Type

Declaration

AMWGVariate <: VectorVariate

Fields

  • value::Vector{VariateType} : vector of sampled values.
  • tune::AMWGTune : tuning parameters for the sampling algorithm.

Constructors

AMWGVariate(x::Vector{VariateType}, tune::AMWGTune)
AMWGVariate(x::Vector{VariateType}, tune=nothing)

Construct a AMWGVariate object that stores sampled values and tuning parameters for adaptive Metropolis-within-Gibbs sampling.

Arguments

  • x : vector of sampled values.
  • tune : tuning parameters for the sampling algorithm. If nothing is supplied, parameters are set to their defaults.

Value

Returns a AMWGVariate type object with fields pointing to the values supplied to arguments x and tune.

AMWGTune Type

Declaration

type AMWGTune

Fields

  • adapt::Bool : whether the proposal distribution has been adaptively tuned.
  • accept::Vector{Integer} : number of accepted candidate draws generated for each element of the parameter vector during adaptive updating.
  • batchsize::Integer : number of samples that must be accumulated before applying an adaptive update to the proposal distributions.
  • m::Integer : number of adaptive update iterations that have been performed.
  • sigma::Vector{Float64} : updated values of the proposal standard deviations if adapt = true, and the user-defined values otherwise.
  • target::Real : target acceptance rate for the adaptive algorithm.

Sampler Constructor

AMWG(params::Vector{Symbol}, sigma::Vector{T<:Real}; adapt::Symbol=:all, batchsize::Integer=50, target::Real=0.44)

Construct a Sampler object for adaptive Metropolis-within-Gibbs sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.

Arguments

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

  • sigma : initial standard deviations for the univariate normal proposal distributions. The standard deviations are relative to the unconstrained parameter space, where candidate draws are generated.

  • adapt : type of adaptation phase. Options are
    • :all : adapt proposals during all iterations.
    • :burnin : adapt proposals during burn-in iterations.
    • :none : no adaptation (Metropolis-within-Gibbs sampling with fixed proposals).

  • batchsize : number of samples that must be accumulated before applying an adaptive update to the proposal distributions.

  • target : a target acceptance rate for the algorithm.

Value

Returns a Sampler type object.

Example

See the Examples section.