Posterior simulation via a random walk.

Definition and Syntax

bool rwmh(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data);
bool rwmh(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, 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 two arguments:
  • vals_inp a vector of input 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 rwmh_n_draws number of posterior draws to keep.
• int rwmh_n_burnin number of burnin draws.
• double rwmh_par_scale scaling parameter for the proposal covariance matrix rwmh_cov_mat.
• arma::mat rwmh_cov_mat covariance matrix of the random walk proposals.

Details

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

• Proposal Step. Generate a proposal draw: $$\theta^{(*)} = \theta^{(i)} + c \Sigma^{1/2} W$$ where $c$ is determined by rwmh_par_scale, $\Sigma$ is determined by rwmh_cov_mat, and $W \sim N(0,I_d)$.
• Accept/Reject Step. Let $$\alpha = \min \left\{ 1, K(\theta^{(*)} | X) / K(\theta^{(i)} | X) \right\}$$ where $K$ is the posterior kernel. 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 with normal prior.

#include "mcmc.hpp"

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

    double mu_0;
    double sigma_0;
};

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

    norm_data* dta = reinterpret_cast<norm_data*>(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)
{
    norm_data* dta = reinterpret_cast< norm_data* >(ll_data);

    const double mu_0 = dta->mu_0;
    const double sigma_0 = dta->sigma_0;
    const double pi = arma::datum::pi;

    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 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;

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

    return 0;
}