MCMC via Hamiltonian dynamics.

Definition and Syntax

bool rmhmc(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, std::function<arma::mat (const arma::vec& vals_inp, arma::cube* tensor_deriv_out, void* tensor_data)> tensor_fn, void* tensor_data);
bool rmhmc(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, std::function<arma::mat (const arma::vec& vals_inp, arma::cube* tensor_deriv_out, void* tensor_data)> tensor_fn, void* tensor_data, algo_settings_t& settings);

Function arguments:

• initial_vals a column vector of initial values.
• draws_out a two-dimensional array containing the posterior draws.
• target_log_kernel the target log-posterior kernel function, taking three arguments:
  • vals_inp a vector of input values;
  • grad_out a vector to store the gradient values; and
  • target_data additional parameters passed to the function.
• target_data additional parameters passed to the posterior kernel.
• tensor_fn metric tensor function, taking three arguments:
  • vals_inp a vector of input values;
  • tensor_deriv_out a three-dimensional array to store derivative values of the metric tensor function; and
  • tensor_data additional parameters passed to the function.
• tensor_data additional parameters passed to the tensor function.
• settings parameters controlling the MCMC routine; see below.

MCMC control parameters:

• bool vals_bound whether the search space is bounded. If true, then
• arma::vec lower_bounds this defines the lower bounds.
• arma::vec upper_bounds this defines the upper bounds.
• int hmc_n_draws number of posterior draws to keep.
• int hmc_n_burnin number of burnin draws.
• double hmc_step_size step size $\epsilon$ below.
• int hmc_leap_steps number of leapfrog steps.
• int rmhmc_fp_steps number of fixed-point iterations.

Details

Let $\theta^{(i)}$ denote a $d$-dimensional vector of stored values at stage $i$ of the algorithm. The basic HMC algorithm consists of three steps.

• Initializaton. Denote the Hamiltonian by $$ H \left\{ \theta, p \right\} := - \ln K(\theta | X) + \frac{1}{2} \log \left\{ (2 \pi)^d | \mathbf{G}(\theta) | \right\} + \frac{1}{2} p^\top \mathbf{G}^{-1}(\theta) p $$ where $\mathbf{G}$ is the metric tensor function. Sample $p^{(i)} \sim N(0,\mathbf{G}(\theta^{(i)}))$, and set $\theta^{(*)} = \theta^{(i)}$ and $p^{(*)} = p^{(i)}$.
• Leapfrog Steps.
  For $k = 1, \ldots,$ hmc_leap_steps do
  • Initializaton. Set $\theta_o^{(*)} = \theta^{(*)}$ and $p_o^{(*)} = p^{(*)}$.
  • Momentum Update Step.
    For $j = 1, \ldots,$ rmhmc_fp_steps do $$p_h^{(*)} = p_o^{(*)} - \dfrac{\epsilon}{2} \times \nabla_\theta H \left\{ \theta_o^{(*)},p_h^{(*)} \right\}$$ where the subscript $h$ denotes a half-step. Notice that $p_h^{(*)}$ appears on both sides of this expression.
  • Position (Proposal) Update Step.
    For $j = 1, \ldots,$ rmhmc_fp_steps do $$\theta^{(*)} = \theta_o^{(*)} + \dfrac{\epsilon}{2} \times \left[ \nabla_p H \left\{ \theta_o^{(*)},p_h^{(*)} \right\} + \nabla_p H \left\{ \theta^{(*)},p_h^{(*)} \right\} \right] $$ Notice that $\theta^{(*)}$ appears on both sides of this expression.
  • Momentum Update Step. $$p^{(*)} = p_h^{(*)} + \epsilon \times \nabla_\theta H \left\{ \theta^{(*)},p_h^{(*)} \right\}$$
• Accept/Reject Step. Define $$ \alpha = \min \left\{ 1, \exp( H(\theta^{(i)}, p^{(i)}) - H(\theta^{(*)}, p^{(*)}) ) \right\} $$ Then $$ \theta^{(i+1)} = \begin{cases} \theta^{(*)} & \text{ if } Z < \alpha \\ \theta^{(i)} & \text{ else } \end{cases} $$ where $Z \sim \text{Unif}(0,1)$.

Examples

Normal likelihood.

  #include "mcmc.hpp"
 
 struct norm_data {
    arma::vec x;
 };
 
 double ll_dens(const arma::vec& vals_inp, arma::vec* grad_out, void* ll_data)
 {
    const double mu    = vals_inp(0);
    const double sigma = vals_inp(1);
 
    const double pi = arma::datum::pi;
 
    norm_data* dta = reinterpret_cast<norm_data*>(ll_data);
    const arma::vec x = dta->x;
    const int n_vals = x.n_rows;
 
    //
 
    const double ret = - ((double) n_vals) * (0.5*std::log(2*pi) + std::log(sigma)) - arma::accu( arma::pow(x - mu,2) / (2*sigma*sigma) );
 
    //
 
    if (grad_out)
    {
        grad_out->set_size(2,1);
 
        //
 
        double m_1 = arma::accu(x - mu);
        double m_2 = arma::accu( arma::pow(x - mu,2) );
 
        (*grad_out)(0,0) = m_1 / (sigma*sigma);
        (*grad_out)(1,0) = (m_2 / (sigma*sigma*sigma)) - ((double) n_vals) / sigma;
    }
 
    //
 
    return ret;
 }
 
arma::mat tensor_fn(const arma::vec& vals_inp, arma::cube* tensor_deriv_out, void* tensor_data)
{
    // const double mu    = vals_inp(0);
    const double sigma = vals_inp(1);
 
    norm_data* dta = reinterpret_cast<norm_data*>(ll_data);
    
    const int n_vals = dta->x.n_rows;
 
    //

    const double sigma_sq = sigma*sigma;
 
    arma::mat tensor_out = arma::zeros(2,2);

    tensor_out(0,0) = ((double) n_vals) / sigma_sq;
    tensor_out(1,1) = 2.0 * ((double) n_vals) / sigma_sq;
    
    //
 
    if (tensor_deriv_out) {
       tensor_deriv_out->set_size(2,2,2);
 
       //
 
       tensor_deriv_out->slice(0).zeros();
 
       tensor_deriv_out->slice(1) = -2.0*tensor_out / sigma;
    }
 
    //
 
    return tensor_out;
}

 double log_target_dens(const arma::vec& vals_inp, arma::vec* grad_out, void* ll_data)
 {
    return ll_dens(vals_inp,grad_out,ll_data);
 }
 
 int main()
 {
    const int n_data = 200;
    const double mu = 2.0;
    const double sigma = 2.0;
 
    norm_data dta;
 
    arma::vec x_dta = mu + sigma*arma::randn(n_data,1);
    dta.x = x_dta;
 
    arma::vec initial_val(2);
    initial_val(0) = mu + 1.0; // mu
    initial_val(1) = sigma + 1.0; // sigma
 
    mcmc::algo_settings_t settings;
 
    settings.hmc_step_size = 0.1;
    settings.hmc_n_burnin = 1000;
    settings.hmc_n_draws = 5000;
 
    printf("begin run\n");
    arma::mat draws_out;
    mcmc::rmhmc(initial_val,draws_out,log_target_dens,&dta,tensor_fn,&dta,settings);
 
    arma::cout << "draws:\n" << draws_out.rows(0,9) << arma::endl;
 
    std::cout << "acceptance rate = " << settings.hmc_accept_rate << arma::endl;
 
    std::cout << "mcmc mean = " << arma::mean(draws_out) << std::endl;
    std::cout << "mcmc var:\n" << arma::cov(draws_out) << std::endl;
 
    return 0;
 }