The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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.
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
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
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.
