Metropolis Hastings in R and Rcpp. Normal model with known variance

metropolis.cpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector f_MHC(NumericVector y, double delta2, double S) {
  double s2 = 1; // known data variance
  double mu = 5; // prior mean
  double t2 = 10; // prior variance

  // Initialize storage
  NumericVector res_theta(S);

  // Starting values
  NumericVector theta(1);
  theta[0] = 0;
  res_theta[0] = theta[0];

  // Loop
  for (int s = 1; s < S; s++) {
    NumericVector theta_star = rnorm(1, theta[0], sqrt(delta2));
    double log_MH_ratio = sum(dnorm(y, theta_star[0], sqrt(s2), true)) +
                          dnorm(theta_star, mu, sqrt(t2), true)[0] -
                          sum(dnorm(y, theta[0], sqrt(s2), true)) -
                          dnorm(theta, mu, sqrt(t2), true)[0];

    if (log(runif(1)[0]) < log_MH_ratio) {
      theta = theta_star;
    }
    res_theta[s] = theta[0];
  }

  return res_theta;
}

/*** R
y <- c(9.37, 10.18, 9.16, 11.60, 10.33)

f_MHR <- function(y, delta2, S) {
  s2 <- 1
  t2 <- 10
  mu <- 5
  theta <- 0
  res_theta <- vector("numeric", length = S)
  res_theta[1] <- theta
  
  for (s in 2:S) {
    theta.star <- rnorm(1, theta, sqrt(delta2))
    log.r <- sum(dnorm(y, theta.star, sqrt(s2), log = TRUE)) +
      dnorm(theta.star, mu, sqrt(t2), log = TRUE) - 
      sum(dnorm(y, theta, sqrt(s2), log = TRUE)) -
      dnorm(theta, mu, sqrt(t2), log = TRUE)
    if (log(runif(1)) < log.r) {
      theta <- theta.star
    }
    res_theta[s] <- theta
  }
  res_theta
}

set.seed(1)
resR <- f_MHR(y = y, delta2 = 2, S = 10000)
plot(resR, type = "l")

set.seed(1)
resC <- f_MHC(y = y, delta2 = 2, S = 10000)

par(mfrow = c(1, 2))
plot(resR, type = "l")
plot(resC, type = "l")

head(resR)
head(resC)
*/