Validation and Simulation Study for the Rtwalk Package"
In Rtwalk: An MCMC Sampler Using the t-Walk Algorithm

knitr::opts_chunk$set(
  echo = TRUE, 
  comment = NA,   
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  warning = FALSE,
  message = FALSE,
  dev = 'png',      
  fig.align = 'center', 
  dpi = 96,     
  out.width = "95%",
  results = "hide"
)

Introduction

This document presents the Monte Carlo study used to validate the implementation of the Rtwalk package. The objective is to demonstrate that the sampler behaves as expected across a variety of challenging scenarios, replicating classic tests from the MCMC literature.

library(Rtwalk)
library(mvtnorm)

Test Battery

N_ITER_VIGNETTE <- 5000
BURN_FRAC <- 0.2

# --- TEST 1: Standard Univariate Normal ---
log_posterior_1 <- function(x) dnorm(x, log = TRUE)
result_1 <- twalk(log_posterior_1, n_iter = N_ITER_VIGNETTE, x0 = -2, xp0 = 2)
calculate_diagnostics(result_1$all_samples, BURN_FRAC, "theta", "Standard Normal")
visualize_results(result_1$all_samples, 0, "Test 1: Standard Normal")

# --- TEST 2: Correlated Bivariate Normal ---
true_mean_2 <- c(1, -0.5)
true_cov_2 <- matrix(c(4, 0.7*2*1.5, 0.7*2*1.5, 2.25), 2, 2)
log_posterior_2 <- function(x) mvtnorm::dmvnorm(x, mean = true_mean_2, sigma = true_cov_2, log = TRUE)
result_2 <- twalk(log_posterior_2, n_iter = N_ITER_VIGNETTE, x0 = c(0,0), xp0 = c(2,-1))
calculate_diagnostics(result_2$all_samples, BURN_FRAC, c("theta1", "theta2"), "Bivariate Normal")
visualize_results(result_2$all_samples, true_mean_2, "Test 2: Bivariate Normal", true_covariance = true_cov_2)

# --- TEST 3: Funnel Distribution ---
log_posterior_3 <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  log_prior_x1 <- dnorm(x1, mean = 0, sd = 3, log = TRUE)
  log_lik_x2 <- dnorm(x2, mean = 0, sd = exp(x1 / 2), log = TRUE)
  return(log_prior_x1 + log_lik_x2)
}
result_3 <- twalk(log_posterior_3, n_iter = N_ITER_VIGNETTE, x0 = c(-1.5, -0.2), xp0 = c(1.5, -0.2))
calculate_diagnostics(result_3$all_samples, BURN_FRAC, c("theta1", "theta2"), "Funnel")
visualize_results(result_3$all_samples, NULL, "Test 3: Funnel")

# --- TEST 4: Rosenbrock Distribution ---
log_posterior_4 <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  k <- 1 / 20 
  return(-k * (100 * (x2 - x1^2)^2 + (1 - x1)^2))
}
result_4 <- twalk(log_posterior_4, n_iter = 100000, x0 = c(0,0), xp0 = c(-1,1))
calculate_diagnostics(result_4$all_samples, BURN_FRAC, c("theta1", "theta2"), "Rosenbrock")
visualize_results(result_4$all_samples, c(1,1), "Test 4: Rosenbrock")

# --- TEST 5: Gaussian Mixture ---
weight1 <- 0.7; mean1 <- c(6, 0); sigma1_1 <- 4; sigma1_2 <- 5; rho1 <- 0.8
cov1 <- matrix(c(sigma1_1^2, rho1*sigma1_1*sigma1_2, rho1*sigma1_1*sigma1_2, sigma1_2^2), nrow=2)

weight2 <- 0.3; mean2 <- c(-3, 10); sigma2_1 <- 1; sigma2_2 <- 1; rho2 <- 0.1
cov2 <- matrix(c(sigma2_1^2, rho2*sigma2_1*sigma2_2, rho2*sigma2_1*sigma2_2, sigma2_2^2), nrow=2)

log_posterior_5 <- function(x) {
  log(weight1 * mvtnorm::dmvnorm(x, mean1, cov1) + weight2 * mvtnorm::dmvnorm(x, mean2, cov2))
}

result_5 <- twalk(log_posterior_5, n_iter = 100000, x0 = mean1, xp0 = mean2)
calculate_diagnostics(result_5$all_samples, BURN_FRAC, c("theta1", "theta2"), "Gaussian Mixture")
visualize_results(result_5$all_samples, NULL, "Test 5: Gaussian Mixture")

# --- TEST 6: High Dimensionality (10D) ---
n_dim_6 <- 10; true_mean_6 <- 1:n_dim_6; rho <- 0.7
true_cov_6 <- matrix(rho^abs(outer(1:n_dim_6, 1:n_dim_6, "-")), n_dim_6, n_dim_6)
log_posterior_6 <- function(x) mvtnorm::dmvnorm(x, mean = true_mean_6, sigma = true_cov_6, log = TRUE)
result_6 <- twalk(log_posterior_6, n_iter = N_ITER_VIGNETTE, x0 = rep(0, n_dim_6), xp0 = rep(2, n_dim_6))
calculate_diagnostics(result_6$all_samples, BURN_FRAC, paste0("theta", 1:n_dim_6), "10D Normal")
visualize_results(result_6$all_samples, true_mean_6, "Test 6: 10D Normal")

# --- TEST 7: Bayesian Logistic Regression ---
set.seed(123); n_obs <- 2000
true_beta <- c(0.5, -1.2, 0.8)
X <- cbind(1, rnorm(n_obs, 0, 1), rnorm(n_obs, 0, 1.5))
eta <- X %*% true_beta; prob <- 1 / (1 + exp(-eta)); y <- rbinom(n_obs, 1, prob)
log_posterior_7 <- function(beta, X, y) { 
  eta <- X %*% beta
  log_lik <- sum(y * eta - log(1 + exp(eta)))
  log_prior <- sum(dnorm(beta, 0, 5, log = TRUE))
  return(log_lik + log_prior) 
}
result_7 <- twalk(log_posterior_7, n_iter = N_ITER_VIGNETTE, x0 = c(0,0,0), xp0 = c(0.2,-0.2,0.1), X = X, y = y)
calculate_diagnostics(result_7$all_samples, BURN_FRAC, c("beta0", "beta1", "beta2"), "Logistic Regression")
visualize_results(result_7$all_samples, true_beta, "Test 7: Logistic Regression")

Study Conclusion

The results from the test battery demonstrate that this implementation of the t-walk is robust and behaves as expected. The sampler was able to converge and efficiently explore distributions with high correlation and multimodality without the need for manual tuning, validating its utility as a general-purpose tool for Bayesian inference.