Validation and Simulation Study for the Rtwalk Package

This vignette runs the t-walk on targets whose answers are known and compares the estimates with the truth. It also shows the basic workflow of the package: running chains, checking convergence, and summarizing the output. Two further vignettes go deeper: vignette("reproducing_christen_fox") reproduces the experiments of the original article, and vignette("sampler_comparison") compares the t-walk with other general-purpose samplers.

library(Rtwalk)
set.seed(2024)

A small helper puts the posterior summaries of summary() next to the true values:

compare <- function(fit, truth, burnin_frac = 0.2) {
  s <- summary(fit, burnin_frac = burnin_frac)
  data.frame(
    parameter = s$Parameter, true_mean = truth$mean, estimate = s$Mean,
    mcse = s$MCSE, true_sd = truth$sd, sd = s$SD, ess_bulk = s$ESS_bulk,
    rhat = s$Rhat
  )
}

1. A first example: the standard normal

fit <- twalk(function(x) dnorm(x, log = TRUE), n_iter = 20000,
             x0 = -2, xp0 = 2, show_progress = FALSE)
fit
#> t-walk MCMC output
#> Iterations per chain: 20000 
#> Chains: 1 
#> Dimension: 1 
#> 
#> SAMPLER PERFORMANCE
#> Acceptance rate: 62.17%
#> Actual move rate: 62.17%
#> Accepted no-move rate: 0.00%
#> Kernel statistics: available in `$kernel_statistics`
compare(fit, list(mean = 0, sd = 1))
#>   parameter true_mean estimate   mcse true_sd     sd ess_bulk  rhat
#> 1    theta1         0  -0.0283 0.0379       1 1.0569      789 1.003

The estimate is within a few Monte Carlo standard errors (mcse) of the true value, the bulk effective sample size is large and R-hat is close to 1.

2. Correlated bivariate normal, four chains

Several chains started from dispersed points make the R-hat diagnostic informative. Starting points are given as matrices with one row per chain.

mu <- c(1, -0.5)
Sigma <- matrix(c(4, 2.1, 2.1, 2.25), 2)
log_bvn <- function(x) mvtnorm::dmvnorm(x, mu, Sigma, log = TRUE)

starts <- function() matrix(rnorm(8, sd = 5), nrow = 4)
fit_bvn <- twalk(log_bvn, n_iter = 10000, x0 = starts(), xp0 = starts(),
                 n_chains = 4, n_cores = 1, show_progress = FALSE)
compare(fit_bvn, list(mean = mu, sd = sqrt(diag(Sigma))))
#>   parameter true_mean estimate   mcse true_sd     sd ess_bulk  rhat
#> 1    theta1       1.0   1.2458 0.0734     2.0 1.8922      667 1.005
#> 2    theta2      -0.5  -0.3539 0.0535     1.5 1.4122      699 1.009

plot() shows the trace of the log posterior, which is the first thing to inspect for burn-in, and of each parameter, with one colour per chain:

plot(fit_bvn)

3. A constrained parameter: Gamma posterior

Return -Inf outside the support. Named starting points give named parameters, both inside log_posterior and in the output.

y <- c(2.1, 0.7, 1.9, 3.4, 0.4, 1.2, 2.8, 1.5)
# Exponential likelihood with rate lambda and a Gamma(2, 1) prior:
# the posterior is Gamma(2 + n, 1 + sum(y)).
log_post_rate <- function(theta) {
  if (theta[["lambda"]] <= 0) return(-Inf)
  sum(dexp(y, theta[["lambda"]], log = TRUE)) + dgamma(theta[["lambda"]], 2, 1, log = TRUE)
}

fit_rate <- twalk(log_post_rate, n_iter = 20000, x0 = c(lambda = 0.5),
                  xp0 = c(lambda = 1), show_progress = FALSE)
shape <- 2 + length(y); rate <- 1 + sum(y)
compare(fit_rate, list(mean = shape / rate, sd = sqrt(shape) / rate))
#>   parameter true_mean estimate   mcse   true_sd     sd ess_bulk rhat
#> 1    lambda 0.6666667   0.6631 0.0079 0.2108185 0.2107      814    1

4. Ten correlated dimensions

d <- 10
Sigma10 <- 0.7^abs(outer(1:d, 1:d, "-"))
mu10 <- 1:d
log_mvn <- function(x) mvtnorm::dmvnorm(x, mu10, Sigma10, log = TRUE)

fit10 <- twalk(log_mvn, n_iter = 1e5, x0 = rep(0, d), xp0 = rep(2, d),
               show_progress = FALSE, thin = 10)
compare(fit10, list(mean = mu10, sd = rep(1, d)))
#>    parameter true_mean estimate   mcse true_sd     sd ess_bulk  rhat
#> 1     theta1         1   1.2170 0.1053       1 0.9906       89 1.005
#> 2     theta2         2   2.1741 0.1277       1 1.0001       64 1.008
#> 3     theta3         3   3.1228 0.1324       1 0.9865       56 1.004
#> 4     theta4         4   4.1371 0.1530       1 0.9788       44 1.008
#> 5     theta5         5   5.1238 0.1579       1 0.9792       42 1.002
#> 6     theta6         6   6.1319 0.1246       1 1.0547       79 1.002
#> 7     theta7         7   7.0902 0.1419       1 1.0357       56 1.003
#> 8     theta8         8   8.0093 0.1329       1 0.9749       57 1.015
#> 9     theta9         9   8.9690 0.0980       1 0.9010       86 1.026
#> 10   theta10        10  10.0047 0.0875       1 0.8701      100 1.017

Thinning (thin = 10) stores one iteration in ten, which reduces memory use in long runs; diagnostics are computed on the stored draws.

5. Bayesian logistic regression

Additional arguments of twalk() are passed to the log posterior.

n_obs <- 1000
X <- cbind(1, rnorm(n_obs), rnorm(n_obs, sd = 1.5))
beta_true <- c(0.5, -1.2, 0.8)
y_bin <- rbinom(n_obs, 1, plogis(X %*% beta_true))

log_post_logit <- function(beta, X, y) {
  eta <- drop(X %*% beta)
  sum(y * eta - pmax(eta, 0) - log1p(exp(-abs(eta)))) +
    sum(dnorm(beta, 0, 5, log = TRUE))
}

fit_logit <- twalk(log_post_logit, n_iter = 30000,
                   x0 = c(b0 = 0, b1 = 0, b2 = 0), xp0 = c(0.2, -0.2, 0.1),
                   X = X, y = y_bin, show_progress = FALSE)

glm_fit <- glm(y_bin ~ X - 1, family = binomial)
compare(fit_logit, list(mean = unname(coef(glm_fit)),
                        sd = unname(sqrt(diag(vcov(glm_fit))))))
#>   parameter  true_mean estimate   mcse    true_sd     sd ess_bulk  rhat
#> 1        b0  0.5778943   0.5828 0.0054 0.08083740 0.0842      245 1.018
#> 2        b1 -1.1976142  -1.2040 0.0064 0.09831746 0.1011      249 1.023
#> 3        b2  0.7625827   0.7648 0.0037 0.06478374 0.0663      326 1.019

With 1000 observations and a vague prior the posterior is close to the sampling distribution of the maximum likelihood estimate, so the two agree.

6. Continuing a run

If the diagnostics suggest that more iterations are needed, twalk_continue() extends the same Markov chain from where it stopped:

fit_more <- twalk_continue(fit_logit, log_post_logit, n_iter = 30000,
                           X = X, y = y_bin, show_progress = FALSE)
fit_more$n_iter
#> [1] 60000
summary(fit_more)$ESS_bulk
#> [1] 531 616 663

7. Interoperability

Results convert to the formats of the coda and posterior packages:

coda::effectiveSize(coda::as.mcmc(fit_logit, burnin_frac = 0.2))
#>       b0       b1       b2 
#> 278.0470 346.2790 300.5871

if (requireNamespace("posterior", quietly = TRUE)) {
  draws <- posterior::as_draws_df(fit_bvn, burnin_frac = 0.2)
  posterior::summarise_draws(draws)
}
#> # A tibble: 3 × 10
#>   variable   mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
#>   <chr>     <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 theta1    1.25   1.21  1.89  1.87  -1.82  4.46  1.01     667.    1024.
#> 2 theta2   -0.354 -0.406 1.41  1.34  -2.63  2.04  1.01     699.     977.
#> 3 lp__     -3.57  -3.24  0.985 0.694 -5.57 -2.65  1.01     737.    1478.

8. How the two points and the rates are reported

Each t-walk chain evolves two points. The primary, time-ordered trajectory used by summary() is stored in samples, the auxiliary one in companion_samples, and the log posterior of the stored points in log_posterior.

The returned acceptance_rate follows the Metropolis–Hastings definition of the original algorithm. move_rate reports the proportion of iterations that actually changed the state, while no_move_rate reports accepted identity proposals; kernel_statistics breaks these rates down by move:

fit10$kernel_statistics
#>     kernel proposed accepted moved acceptance_rate  move_rate
#> 1 traverse    49259     2680  2369      0.05440630 0.04809273
#> 2     walk    49095    15025 14723      0.30603931 0.29988797
#> 3     blow      811       59    56      0.07274969 0.06905055
#> 4      hop      835      317   308      0.37964072 0.36886228