---
title: "Validation and Simulation Study for the Rtwalk Package"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Validation and Simulation Study}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
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.

```{r setup, message=FALSE, warning=FALSE}
library(Rtwalk)
library(mvtnorm)
```


## Test Battery

```{r, eval=TRUE, cache=FALSE}
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$samples, BURN_FRAC, "theta", "Standard Normal")
visualize_results(result_1$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$samples, BURN_FRAC, c("theta1", "theta2"), "Bivariate Normal")
visualize_results(result_2$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$samples, BURN_FRAC, c("theta1", "theta2"), "Funnel")
visualize_results(result_3$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 = N_ITER_VIGNETTE, x0 = c(0,0), xp0 = c(-1,1))
calculate_diagnostics(result_4$samples, BURN_FRAC, c("theta1", "theta2"), "Rosenbrock")
visualize_results(result_4$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_component_1 <- log(weight1) + mvtnorm::dmvnorm(x, mean1, cov1, log = TRUE)
  log_component_2 <- log(weight2) + mvtnorm::dmvnorm(x, mean2, cov2, log = TRUE)
  max_log <- max(log_component_1, log_component_2)

  max_log + log(
    exp(log_component_1 - max_log) + exp(log_component_2 - max_log)
  )
}

result_5 <- twalk(log_posterior_5, n_iter = N_ITER_VIGNETTE, x0 = mean1, xp0 = mean2)
calculate_diagnostics(result_5$samples, BURN_FRAC, c("theta1", "theta2"), "Gaussian Mixture")
visualize_results(result_5$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$samples, BURN_FRAC, paste0("theta", 1:n_dim_6), "10D Normal")
visualize_results(result_6$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 <- plogis(eta); y <- rbinom(n_obs, 1, prob)
log_posterior_7 <- function(beta, X, y) { 
  eta <- X %*% beta
  log_lik <- sum(y * eta - pmax(eta, 0) - log1p(exp(-abs(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$samples, BURN_FRAC, c("beta0", "beta1", "beta2"), "Logistic Regression")
visualize_results(result_7$samples, true_beta, "Test 7: Logistic Regression")
```

## Understanding the returned trajectories

Each t-walk chain evolves two points. The primary, time-ordered trajectory used
by `summary()` and the examples above is stored in `samples`; the auxiliary
trajectory is stored in `companion_samples`. The `all_samples` field is a
backward-compatible row binding of both trajectories and is not itself one
time-ordered Markov chain.

The returned `acceptance_rate` follows the Metropolis--Hastings definition from
the original algorithm. `move_rate` reports the proportion of iterations that
actually changed the state, while `no_move_rate` reports accepted identity
proposals.

For parallel runs, `summary()` applies burn-in independently to each element of
`individual_chains` before combining posterior summaries.

## 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.