MCMC via Langevin diffusions.

Definition and Syntax

bool mala(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);
bool mala(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, 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.
• 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 mala_n_draws number of posterior draws to keep.
• int mala_n_burnin number of burnin draws.
• double mala_step_size step size $\epsilon$ below.
• arma::mat mala_precond_mat preconditioning matrix $\mathbf{M}$.

Details

Let $\theta^{(i)}$ denote a $d$-dimensional vector of stored values at stage $i$ of the algorithm. The basic MALA algorithm is as follows.

• Proposal Step. Let $$\mu(\theta^{(i)}) := \theta^{(i)} + \frac{\epsilon^2}{2} \mathbf{M} \nabla_\theta \ln K(\theta^{(i)} | X)$$ where $K$ is the posterior kernel, $\mathbf{M}$ is set by mala_precond_mat, and $\epsilon = $ mala_step_size.
  
  Generate a proposal draw: $$\theta^{(*)} = \mu(\theta^{(i)}) + \epsilon W$$ where $W \sim N(0,\mathbf{M})$.
• Accept/Reject Step. Denote the proposal density by $q(\theta^{(*)} | \theta^{(i)}) := \phi(\theta^{(*)}; \mu(\theta^{(i)}), \epsilon^2 \mathbf{M})$ and let $$\alpha = \min \left\{ 1, [ K(\theta^{(*)} | X) q(\theta^{(i)} | \theta^{(*)})] / [ K(\theta^{(i)} | X) q(\theta^{(*)} | \theta^{(i)})] \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;
}

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 = 1000;
    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; // mu
    initial_val(1) = sigma + 1; // sigma

    mcmc::algo_settings_t settings;

    settings.mala_step_size = 0.10;
    settings.mala_n_burnin = 1000;
    settings.mala_n_draws = 1000;

    arma::mat draws_out;
    mcmc::mala(initial_val,draws_out,log_target_dens,&dta,settings);

    std::cout << "mala: mcmc mean = " << arma::mean(draws_out) << std::endl;

    return 0;
}