Random Walk Metropolis-Hastings¶

Table of contents

Algorithm Description
Function Declarations
Control Parameters
Examples
- Gaussian Mean

Algorithm Description ¶

The Random Walk Metropolis-Hastings algorithm is the fundamental Markov Chain Monte Carlo method, generating draws from a target posterior distribution via a random walk proposal.

Let \(\theta^{(i)}\) denote a \(d\)-dimensional vector of stored values at stage \(i\) of the algorithm. The RWMH algorithm proceeds in two steps.

(Proposal Step) Generate a proposal draw, \(\theta^{(*)}\), using:

\[\theta^{(*)} = \theta^{(i)} + c \times \Sigma^{1/2} W\]

where \(c\) is a scalar value set via rwmh_settings.par_scale; \(\Sigma\) is a matrix set via rwmh_settings.cov_mat; and \(W \sim N(0,I_d)\).

(Accept/Reject Step) Define

\[\alpha = \min \left\{ 1, K(\theta^{(*)} | X) / K(\theta^{(i)} | X) \right\}\]

where \(K\) denotes the posterior kernel. Then

\[\theta^{(i+1)} = \begin{cases} \theta^{(*)} & \text{ with probability } \alpha \\ \theta^{(i)} & \text{ else } \end{cases}\]

The algorithm stops when the number of draws reaches n_burnin_draws + n_keep_draws, and returns the final n_keep_draws number of draws.

Function Declarations ¶

bool rwmh(const ColVec_t &initial_vals, std::function<fp_t(const ColVec_t &vals_inp, void *target_data)> target_log_kernel, Mat_t &draws_out, void *target_data)¶

The Random Walk Metropolis-Hastings MCMC Algorithm.

Parameters

initial_vals – a column vector of initial values.
target_log_kernel – the log posterior kernel function of the target distribution, taking two arguments:
- vals_inp a vector of inputs; and
- target_data additional data passed to the user-provided function.
draws_out – a matrix of posterior draws, where each row represents one draw.
target_data – additional data passed to the user-provided function.

Returns

a boolean value indicating successful completion of the algorithm.

bool rwmh(const ColVec_t &initial_vals, std::function<fp_t(const ColVec_t &vals_inp, void *target_data)> target_log_kernel, Mat_t &draws_out, void *target_data, algo_settings_t &settings)¶

The Random Walk Metropolis-Hastings MCMC Algorithm.

Parameters

initial_vals – a column vector of initial values.
target_log_kernel – the log posterior kernel function of the target distribution, taking two arguments:
- vals_inp a vector of inputs; and
- target_data additional data passed to the user-provided function.
draws_out – a matrix of posterior draws, where each row represents one draw.
target_data – additional data passed to the user-provided function.
settings – parameters controlling the MCMC routine.

Returns

a boolean value indicating successful completion of the algorithm.

Control Parameters ¶

The basic control parameters are:

size_t rwmh_settings.n_burnin_draws: number of burn-in draws.
size_t rwmh_settings.n_keep_draws: number of draws to keep (post sample burn-in period).
bool vals_bound: whether the search space of the algorithm is bounded. If true, then
- ColVec_t lower_bounds: defines the lower bounds of the search space.
- ColVec_t upper_bounds: defines the upper bounds of the search space.

Additional settings:

int rwmh_settings.omp_n_threads: the number of OpenMP threads to use.
- Default value: -1 (use all available threads divided by 2).
fp_t rwmh_settings.par_scale: scaling parameter for the proposal step.
- Default value: 1.0.
Mat_t rwmh_settings.cov_mat: covariance matrix for the proposal step.
- Default value: diagonal matrix.

Examples ¶

Gaussian Mean ¶

Code to run this example is given below.

Armadillo (Click to show/hide)

#define MCMC_ENABLE_ARMA_WRAPPERS
#include "mcmc.hpp"

struct norm_data_t {
    double sigma;
    arma::vec x;

    double mu_0;
    double sigma_0;
};

double ll_dens(const arma::vec& vals_inp, void* ll_data)
{
    const double pi = arma::datum::pi;

    //

    const double mu = vals_inp(0);

    norm_data_t* dta = reinterpret_cast<norm_data_t*>(ll_data);
    const double sigma = dta->sigma;
    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) );

    //

    return ret;
}

double log_pr_dens(const arma::vec& vals_inp, void* ll_data)
{
    const double pi = arma::datum::pi;

    //

    norm_data_t* dta = reinterpret_cast< norm_data_t* >(ll_data);

    const double mu_0 = dta->mu_0;
    const double sigma_0 = dta->sigma_0;

    const double x = vals_inp(0);

    const double ret = - 0.5*std::log(2*pi) - std::log(sigma_0) - std::pow(x - mu_0,2) / (2*sigma_0*sigma_0);

    return ret;
}

double log_target_dens(const arma::vec& vals_inp, void* ll_data)
{
    return ll_dens(vals_inp,ll_data) + log_pr_dens(vals_inp,ll_data);
}

int main()
{
    const int n_data = 100;
    const double mu = 2.0;

    norm_data_t dta;
    dta.sigma = 1.0;
    dta.mu_0 = 1.0;
    dta.sigma_0 = 2.0;

    arma::vec x_dta = mu + arma::randn(n_data,1);
    dta.x = x_dta;

    arma::vec initial_val(1);
    initial_val(0) = 1.0;

    //

    mcmc::algo_settings_t settings;

    settings.rwmh_settings.par_scale = 0.4;
    settings.rwmh_settings.n_burnin_draws = 2000;
    settings.rwmh_settings.n_keep_draws = 2000;

    //

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

    //

    std::cout << "rwmh mean:\n" << arma::mean(draws_out) << std::endl;
    std::cout << "acceptance rate: " << static_cast<double>(settings.rwmh_settings.n_accept_draws) / settings.rwmh_settings.n_keep_draws << std::endl;

    //

    return 0;
}

Eigen (Click to show/hide)

#define MCMC_ENABLE_EIGEN_WRAPPERS
#include "mcmc.hpp"

inline
Eigen::VectorXd
eigen_randn_colvec(size_t nr)
{
    static std::mt19937 gen{ std::random_device{}() };
    static std::normal_distribution<> dist;

    return Eigen::VectorXd{ nr }.unaryExpr([&](double x) { (void)(x); return dist(gen); });
}

struct norm_data_t {
    double sigma;
    Eigen::VectorXd x;

    double mu_0;
    double sigma_0;
};

double ll_dens(const Eigen::VectorXd& vals_inp, void* ll_data)
{
    const double pi = 3.14159265358979;

    //

    const double mu = vals_inp(0);

    norm_data_t* dta = reinterpret_cast<norm_data_t*>(ll_data);
    const double sigma = dta->sigma;
    const Eigen::VectorXd x = dta->x;

    const int n_vals = x.size();

    //

    const double ret = - n_vals * (0.5 * std::log(2*pi) + std::log(sigma)) - (x.array() - mu).pow(2).sum() / (2*sigma*sigma);

    //

    return ret;
}

double log_pr_dens(const Eigen::VectorXd& vals_inp, void* ll_data)
{
    const double pi = 3.14159265358979;

    //

    norm_data_t* dta = reinterpret_cast< norm_data_t* >(ll_data);

    const double mu_0 = dta->mu_0;
    const double sigma_0 = dta->sigma_0;

    const double x = vals_inp(0);

    const double ret = - 0.5*std::log(2*pi) - std::log(sigma_0) - std::pow(x - mu_0,2) / (2*sigma_0*sigma_0);

    return ret;
}

double log_target_dens(const Eigen::VectorXd& vals_inp, void* ll_data)
{
    return ll_dens(vals_inp,ll_data) + log_pr_dens(vals_inp,ll_data);
}

int main()
{
    const int n_data = 100;
    const double mu = 2.0;

    norm_data_t dta;
    dta.sigma = 1.0;
    dta.mu_0 = 1.0;
    dta.sigma_0 = 2.0;

    Eigen::VectorXd x_dta = mu + eigen_randn_colvec(n_data).array();
    dta.x = x_dta;

    Eigen::VectorXd initial_val(1);
    initial_val(0) = 1.0;

    //

    mcmc::algo_settings_t settings;

    settings.rwmh_settings.par_scale = 0.4;
    settings.rwmh_settings.n_burnin_draws = 2000;
    settings.rwmh_settings.n_keep_draws = 2000;

    //

    Eigen::MatrixXd draws_out;
    mcmc::rwmh(initial_val, log_target_dens, draws_out, &dta, settings);

    //

    std::cout << "hmc mean:\n" << draws_out.colwise().mean() << std::endl;
    std::cout << "acceptance rate: " << static_cast<double>(settings.rwmh_settings.n_accept_draws) / settings.rwmh_settings.n_keep_draws << std::endl;

    //

    return 0;
}

Random Walk Metropolis-Hastings¶

Algorithm Description¶

Function Declarations¶

Control Parameters¶

Examples¶

Gaussian Mean¶

Algorithm Description ¶

Function Declarations ¶

Control Parameters ¶

Examples ¶

Gaussian Mean ¶