19 Metropolis

Previously, we used the rejection algorithm to build intuition about sampling from an unnormalized posterior. We showed that it can be easier to sample from a distribution than to work with closed-form results. The rejection algorithm is simple and intuitive, but requires (1) finding a suitable envelope constant and (2) rejecting many samples in realistic problems (e.g., logistic regression).

We can improve the rejection algorithm by making proposals that depend on the current value. Many algorithms—including Stan!—do this. The Metropolis algorithm is one example. The Metropolis algorithm belongs to a larger class of Markov chain Monte Carlo (MCMC) algorithms,¹ which includes Metropolis-Hastings, Gibbs sampling, and the Hamiltonian Monte Carlo algorithm used by Stan.

¹ MCMC refers to a family of algorithms for generating samples from a target probability distribution by constructing a Markov chain whose stationary distribution is that target. Over many iterations, the samples produced by the chain approximate (dependent) draws from the posterior distribution, allowing us to estimate quantities such as means, variances, and credible intervals even when the posterior cannot be computed analytically.

19.1 Algorithm: Metropolis

Inputs:

The unnormalized log-posterior \(\log f(\theta\mid y)\).
Desired number of iterations \(S\). Note: the samples are dependent so choose large \(S\) (e.g., \(S = 20{,}000\)), perhaps very large (e.g., \(S = 200{,}000\)).
Tuning parameter \(\tau\), which controls the width of the uniform proposal density.

Algorithm:

Initialize: Choose an initial \(\theta^{(1)}\). Set \(s=1\).
Repeat until \(s=S\):
- Propose: Draw \(z\) from a uniform distribution from \([\theta^{(s)} - \tau]\) to \([\theta^{(s)} + \tau]\).
- Compute log-acceptance: \(\Delta = \log f(z)-\log f \big(\theta^{(s)}\big)\)
- Accept/reject: Draw \(u \sim \text{Uniform}(0,1)\). Accept iff \(\log u \le \Delta\). If accept, set \(\theta^{(s+1)}=z\); otherwise \(\theta^{(s+1)}=\theta^{(s)}\).²
- Iterate: \(s\leftarrow s+1\).

² Importantly, this accepts with probability equal to \(\frac{f(\text{proposal})}{f(\text{current})}\) if the ratio is between zero and one and always if the ratio is greater than one. This step is similar to the rejection algorithm, except we keep the current value in place of the rejected proposal in Metropolis.

Output: A Markov chain \(\{\theta^{(s)}\}_{s=1}^S\) with stationary distribution \(f(\theta\mid y)\). Consecutive draws are not independent; in practice, discard a burn-in and tune \(\tau\) to maintain a reasonable acceptance rate (perhaps 0.2-0.5, problem-dependent).

We can write this algorithm in R.

19.1.1 Minimal version

The version below has only the critical pieces.

metrop <- function(logf, theta_start, S = 10000, tau = 0.1, ...) {
  # initialize matrix of samples with starting values
  k <- length(theta_start)
  samples <- matrix(NA_real_, nrow = S, ncol = k)
  samples[1, ] <- theta_start

  # proceed with algorithm
  for (s in 2:S) {
    # extract current location
    current <- samples[s - 1, ]

    # generate symmetric random-walk proposal
    proposed_move <- runif(k, -tau, tau)
    proposal <- current + proposed_move

    # acceptance step
    delta <- logf(proposal, ...) - logf(current, ...)
    if (delta > 0) {
      accept <- TRUE
    } else {
      accept <- (log(runif(1)) <= delta)
    }

    # update samples
    if (accept) {
      samples[s, ] <- proposal
    } else {
      samples[s, ] <- current
    }
  }

  samples
}

19.1.2 Some features added

The version below adds a few features, like a helpful message and an optional progress bar.

metrop <- function(logf, theta_start, S = 10000, tau = 0.1, progress_bar = FALSE, ...) {
  # record start time
  start_time <- Sys.time()

  # initialize matrix of samples with starting values
  k <- length(theta_start)
  samples <- matrix(NA_real_, nrow = S, ncol = k)
  samples[1, ] <- theta_start
  n_accepted <- 0

  # optionally initialize progress bar
  if (progress_bar) {
    pb <- txtProgressBar(min = 0, max = S, style = 3)
  }

  # proceed with algorithm
  for (s in 2:S) {
    # extract current location
    current <- samples[s - 1, ]

    # generate symmetric random-walk proposal
    proposed_move <- runif(k, -tau, tau)
    proposal <- current + proposed_move 

    # acceptance step
    delta  <- logf(proposal, ...) - logf(current, ...)
    if (delta > 0) {
      accept <- TRUE
    } else {
      accept <- (log(runif(1)) <= delta)
    }

    # update samples
    if (accept) {
      samples[s, ] <- proposal
      n_accepted <- n_accepted + 1
    } else {
      samples[s, ] <- current
    }

    # update progress bar occasionally
    if (progress_bar && (s %% (S / 100) == 0 || s == S)) {
      setTxtProgressBar(pb, s)
    }
  }

  # close progress bar if used
  if (progress_bar) close(pb)

  # print summary report
  message(
    paste0(
      "💪 Successfully generated ", scales::comma(S), " dependent samples! 🎉\n\n",
      "✅ Accepted moves: ", scales::comma(n_accepted), "\n",
      "﹪ Acceptance rate: ", scales::percent(n_accepted / (S - 1), accuracy = 1),
      " (tune tau so that acceptance rate is about 20%-50%).\n",
      "⏰ Total time: ", prettyunits::pretty_dt(Sys.time() - start_time),  "\n",
      "🧠 Reminder: if a proposal is not accepted, the previous sample is reused. This makes the samples depend on starting values and each other. Discard initial samples, use R-hat to assess convergence, and use ESS to assess effective sample size."
    )
  )

  list(
    prop_accepted = n_accepted / (S - 1),
    samples = samples
  )
}

19.2 Toothpaste Cap

This example has a known analytical solution. For \(k=8\) and \(N=150\) with a \(\text{beta}(3, 15)\) prior, the posterior is \(\text{beta}(11,157)\).

# create log-target for Beta posterior; -Inf outside (0,1)
logf <- function(p) {
  ifelse(p <= 0 | p >= 1, 
         -Inf, 
         dbeta(p, 11, 157, log = TRUE))
}

# run algorithm starting in the middle of (0,1)
m <- metrop(
  logf,            # log-target
  theta_start = 0.5
)

We can compare the analytical posterior mean with the simulation.

# closed-form posterior mean
11 / (11 + 157)

[1] 0.0655

# simulations (samples is S x k; here k = 1)
mean(m$samples)

[1] 0.0663

We can plot the histogram and time series of the samples. The chain explores the posterior nicely overall; zooming into the first 100 iterations shows the initial transient; later windows show dependence across draws.

par(mfrow = c(2, 2))  # arrange figures
hist(m$samples, main = "All 10,000 samples")
plot(m$samples, type = "l", main = "All 10,000 samples")
plot(m$samples[1:100], type = "l", main = "First 100 samples")
plot(m$samples[1001:1100], type = "l", main = "Samples 1,001 through 1,100")

There are two things worth noting.

The samples are not independent. Many consecutive draws are identical (rejections) or move only a little; the chain explores the target slowly.
It takes a few iterations for the chain to reach the high-probability region (burn-in). Starting at \(\pi=0.5\) is far from the posterior mass; it takes roughly a few dozen iterations to reach \(\approx 0.1\).

19.2.1 Burn-in

The MCMC samples are heavily dependent on the starting values, so discard early samples that depend strongly on the starting value (the “burn-in”). Discarding the first 10% to 50% of the samples is typical.

Stan and {brms} call this period the “warmup” rather than the burn-in. During warmup, Stan’s sampler also tunes itself (e.g., its step size), so the warmup samples are not draws from a fixed algorithm. Our Metropolis sampler keeps \(\tau\) fixed, so “burn-in” fits here.

Rather than use all the samples (including our chosen starting value) to compute the posterior mean, let’s use only the last 5,000 samples.

# simulations, keep only samples 5,001 through 10,000 (discard first half)
mean(m$samples[5001:10000, 1])

[1] 0.0663

19.2.2 \(\hat{R}\) or R-hat

But how do we know how long to run our algorithms? How do we know if the samples have converged to the target distribution? This is still a difficult problem. No diagnostic can show that the chains have converged, but R-hat can show that they have not. We check that R-hat is less than 1.01 (Vehtari et al. 2021). If R-hat is larger than 1.01, the chains have not converged and we need more iterations. If R-hat is less than 1.01, we have found no sign of a problem.

Vats, Dootika, and Christina Knudson. 2021. “Revisiting the Gelman–Rubin Diagnostic.” Statistical Science 36 (4): 518–29. https://doi.org/10.1214/20-STS812.

Margossian, Charles C., Matthew D. Hoffman, Pavel Sountsov, Lionel Riou-Durand, Aki Vehtari, and Andrew Gelman. 2025. “Nested \(\widehat{R}\): Assessing the Convergence of Markov Chain Monte Carlo When Running Many Short Chains.” Bayesian Analysis 20 (4): 1587–1614. https://doi.org/10.1214/24-BA1453.

To compute R-hat, run several chains from overdispersed starting values. It is important that these be high and low relative to the target distribution. If chains have converged, within- and between-chain variation align and \(\hat{R}\approx 1\). The literature and practices on R-hat are evolving rapidly; see ?posterior::rhat for the latest references. Vats and Knudson (2021) connect the R-hat threshold to the effective sample size, and Margossian et al. (2025) develop a version of R-hat for running many short chains in parallel.

{brms} warns only when some R-hat is larger than 1.05. We use the stricter threshold of 1.01 from Vehtari et al. (2021), so check R-hat yourself even when {brms} prints no warning.

Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved \(\widehat{R}\) for Assessing Convergence of MCMC (with Discussion).” Bayesian Analysis 16 (2): 667–718. https://doi.org/10.1214/20-BA1221.

For example, we can run the algorithm for 25 iterations and see that R-hat is much larger than 1.01.

# run algorithm four times, making four "chains"
m1 <- metrop(logf, theta_start = 0.01, S = 25)
m2 <- metrop(logf, theta_start = 0.25, S = 25)
m3 <- metrop(logf, theta_start = 0.75, S = 25)
m4 <- metrop(logf, theta_start = 0.99, S = 25)

# a matrix of chains for first (and only) dimension
matrix_of_chains <- cbind(
  m1$samples[, 1],
  m2$samples[, 1],
  m3$samples[, 1],
  m4$samples[, 1]
)

# compute r-hat using the rhat() function in {posterior}
posterior::rhat(matrix_of_chains)

[1] 2.69

If we increase S to a larger value and discard the initial samples, then R-hat falls below 1.01.

# run algorithm four times, making four "chains"
m1 <- metrop(logf, theta_start = 0.01, S = 25000)
m2 <- metrop(logf, theta_start = 0.25, S = 25000)
m3 <- metrop(logf, theta_start = 0.75, S = 25000)
m4 <- metrop(logf, theta_start = 0.99, S = 25000)

# a matrix of chains for first (and only) dimension
matrix_of_chains <- cbind(
  m1$samples[10001:25000, 1],
  m2$samples[10001:25000, 1],
  m3$samples[10001:25000, 1],
  m4$samples[10001:25000, 1]
)

# compute r-hat using the rhat() function in {posterior}
posterior::rhat(matrix_of_chains)

[1] 1

19.2.3 Parallel

Multiple chains are naturally parallelizable. In fact, these are called “embarrassingly parallel,” because the chains are run entirely independently—nothing from the first chain is needed for the second, and so on.

# count physical cores (leave some free)
parallel::detectCores(logical = FALSE)

[1] 12

Here is how we might use the {foreach} package to run the chains in parallel. This isn’t important here, because each chain runs very quickly. However, it is important in practice, and modern packages make it easy.

# load packages
library(foreach)
library(doParallel)

# overdispersed starting values
starting_values <- seq(0.05, 0.95, length.out = 10)

# register parallel backend
cl <- makeCluster(10)
registerDoParallel(cl)

# run algorithm 10 times in parallel; combine samples with cbind()
matrix_of_chains <- foreach(s = starting_values, .combine = cbind) %dopar% {
  m <- metrop(
    logf,
    theta_start = s,
    S = 25
  )
  m$samples
}

# stop the parallel backend
stopCluster(cl)

# columns are chains; rows are samples
head(matrix_of_chains)

       [,1]   [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] 0.0500 0.1500 0.250 0.350 0.450 0.550 0.650 0.750 0.850 0.950
[2,] 0.0466 0.0792 0.250 0.267 0.416 0.482 0.650 0.683 0.850 0.934
[3,] 0.0640 0.0792 0.257 0.267 0.416 0.459 0.558 0.683 0.850 0.934
[4,] 0.0660 0.0772 0.157 0.211 0.416 0.401 0.550 0.683 0.850 0.934
[5,] 0.0660 0.0544 0.157 0.211 0.376 0.330 0.535 0.675 0.850 0.892
[6,] 0.0776 0.0632 0.128 0.211 0.376 0.330 0.453 0.586 0.809 0.881

# compute R-hat using {posterior}
posterior::rhat(matrix_of_chains)

[1] 2.06

The figure below shows the 25 samples from each of the 10 chains. Notice the strong dependence on starting values! In practice, we’d need to generate many more samples, discard the first samples, and check that R-hat is less than 1.01.

Code

colnames(matrix_of_chains) <- paste("Chain", 1:ncol(matrix_of_chains))
gg_df <- matrix_of_chains |>
  as_tibble() |>
  mutate(Iteration = 1:n()) |>
  pivot_longer(cols = `Chain 1`:`Chain 10`, names_to = "Chain", values_to = "Sample") |>
  separate(Chain, into = c("tmp", "Chain Number"), remove = FALSE) |>
  mutate(`Chain Number` = as.numeric(`Chain Number`),
         Chain = reorder(Chain, `Chain Number`, ordered = TRUE))

ggplot(gg_df, aes(x = Iteration, y = Sample, color = Chain)) +
  geom_line() +
  theme_ipsum(base_family = "Source Sans 3") +
  labs(title = "Samples generated by Metropolis algorithm",
       subtitle = "10 separate chains")

19.2.4 ESS

Dependent Metropolis samples carry less information than the same number of independent samples. The Effective Sample Size (ESS) summarizes the amount of information in our dependent samples—“our dependent samples are like ______ independent samples.” The literature and practices on ESS are evolving rapidly; see ?posterior::ess_bulk for the latest references.

# the four long chains from above, minus the first 10,000 samples of each
matrix_of_chains <- cbind(
  m1$samples[10001:25000, 1],
  m2$samples[10001:25000, 1],
  m3$samples[10001:25000, 1],
  m4$samples[10001:25000, 1]
)
dim(matrix_of_chains)  # 15,000 x 4 = 60,000 samples

[1] 15000     4

# are like
posterior::ess_bulk(matrix_of_chains)  # ESS for the center of the distribution

[1] 13084

posterior::ess_tail(matrix_of_chains)  # ESS for tail estimands

[1] 13339

19.2.5 Tuning for intuition

The tuning parameter matters a lot. I’ve chosen a smart default for many problems, but the value of tau should be something like the SD of the posterior. Below, we can see how “bad” values of tau work.

# varying the tuning parameter tau
m1 <- metrop(logf, theta_start = 0.99, S = 1000, tau = 0.1)
m2 <- metrop(logf, theta_start = 0.99, S = 1000, tau = 0.01)
m3 <- metrop(logf, theta_start = 0.99, S = 1000, tau = 0.001)
m4 <- metrop(logf, theta_start = 0.99, S = 1000, tau = 3)

par(mfrow = c(2, 2))
plot(m1$samples, type = "l", main = "tau = 0.1\nNice!", ylim = c(0, 1))
plot(m2$samples, type = "l", main = "tau = 0.01\nToo small! (Accept too often.)", ylim = c(0, 1))
plot(m3$samples, type = "l", main = "tau = 0.001\nMuch too small! (Accept much too often.)", ylim = c(0, 1))
plot(m4$samples, type = "l", main = "tau = 3\nToo wide! (Accept too rarely.)", ylim = c(0, 1))

19.3 Sawtooth Prior

We can use the Metropolis algorithm for our “weird” sawtooth posterior. We’ll use the log-posterior in this case.

# sawtooth log-prior (unnormalized): log of ((n_teeth*p) %% 1) on (0,1)
log_prior_saw <- function(p, n_teeth = 5) {
  ifelse(p <= 0 | p >= 1, 
  -Inf, 
  log((n_teeth * p) %% 1))
}

# likelihood for 10 tosses w/ one success
log_lik <- function(p) {
  ifelse(p <= 0 | p >= 1, 
  -Inf, 
  dbinom(1, size = 10, prob = p, log = TRUE))
}

# log-posterior
log_post <- function(p) log_lik(p) + log_prior_saw(p)

# run metropolis algorithm
m <- metrop(log_post, theta_start = 0.5)

# plot samples
par(mfrow = c(1, 2))
plot(m$samples, type = "l")
hist(m$samples, breaks = 100)

19.4 Logistic Regression

19.4.1 Data

Recall the logit model from a previous chapter. This time, let’s rescale the numeric variables to have mean of zero and SD of 0.5. This keeps coefficients on a similar scale, which helps us with the computation. For more on this rescaling strategy and its advantages, see Gelman (2008).

Gelman, Andrew. 2008. “Scaling Regression Inputs by Dividing by Two Standard Deviations.” Statistics in Medicine 27 (15): 2865–73. https://doi.org/10.1002/sim.3107.

# load only the turnout data frame
turnout <- ZeligData::turnout  # see ?ZeligData::turnout for details

# alias for arm::rescale(): mean 0, SD 0.5
rs <- function(x) { arm::rescale(x) }

# create the rescaled variables as columns; rescaling inside the formula
# (e.g., rs(age)) breaks predictions on new data, because rescale() would
# then use the new data's mean and SD
turnout <- turnout |>
  mutate(across(age:income, rs, .names = "rs_{.col}")) |>
  glimpse()

Rows: 2,000
Columns: 8
$ race       <fct> white, white, white, white, white, white, white, white, whi…
$ age        <int> 60, 51, 24, 38, 25, 67, 40, 56, 32, 75, 46, 52, 22, 60, 24,…
$ educate    <dbl> 14, 10, 12, 8, 12, 12, 12, 10, 12, 16, 15, 12, 12, 12, 14, …
$ income     <dbl> 3.346, 1.856, 0.630, 3.418, 2.785, 2.387, 4.286, 9.320, 3.8…
$ vote       <int> 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1,…
$ rs_age     <dbl> 0.4197, 0.1628, -0.6079, -0.2082, -0.5793, 0.6195, -0.1512,…
$ rs_educate <dbl> 0.2868, -0.3067, -0.0099, -0.6034, -0.0099, -0.0099, -0.009…
$ rs_income  <dbl> -0.09566, -0.35915, -0.57595, -0.08284, -0.19482, -0.26532,…

f <- vote ~ rs_age + rs_educate + rs_income + race

# fit logit model
fit <- glm(f, family = binomial, data = turnout)

# print estimates
arm::display(fit, digits = 4)

glm(formula = f, family = binomial, data = turnout)
            coef.est coef.se
(Intercept) 1.0578   0.1396 
rs_age      0.9934   0.1212 
rs_educate  1.1837   0.1370 
rs_income   1.0013   0.1535 
racewhite   0.2508   0.1465 
---
  n = 2000, k = 5
  residual deviance = 2024.0, null deviance = 2266.7 (difference = 242.8)

19.4.2 Log-Posterior

Now define the unnormalized log-posterior. For this example, we imagine an improper, constant prior.

# make X and y
mf <- model.frame(f, data = turnout)
X <- model.matrix(f, data = mf)
y <- model.response(mf)

# ❌ correct, but unstable, log unnormalized posterior 
log_posterior <- function(beta, y, X) {
  linpred <- X %*% beta
  p <- plogis(linpred)
  # for |linpred| large, plogis() rounds to exactly 0 or 1, so dbinom(log = TRUE) is log(0) = -Inf;
  #   if the current and proposed values are both -Inf, delta = -Inf - (-Inf) = NaN
  sum(dbinom(y, size = 1, prob = plogis(linpred), log = TRUE))
}

# ✅ same unnormalized log posterior, but computed on the log scale to avoid the rounding
log_posterior <- function(beta, y, X) {
  linpred <- drop(X %*% beta)
  log1pexp <- ifelse(linpred > 0, linpred + log1p(exp(-linpred)), log1p(exp(linpred)))
  sum(y * linpred - log1pexp)
}

19.4.3 Running the Algorithm

After some experimenting, I found that 100,000 iterations with a 20,000 iteration burn-in works well.

# set seed
set.seed(1234)

# sample with metropolis
S <- 100000
m1 <- metrop(log_posterior, S = S, tau = 0.1, theta_start = rep(-2, ncol(X)), y = y, X = X)
m2 <- metrop(log_posterior, S = S, tau = 0.1, theta_start = rep(-1, ncol(X)), y = y, X = X)
m3 <- metrop(log_posterior, S = S, tau = 0.1, theta_start = rep(1, ncol(X)), y = y, X = X)
m4 <- metrop(log_posterior, S = S, tau = 0.1, theta_start = rep(2, ncol(X)), y = y, X = X)

Now we can verify that the samplers have converged (R-hat less than about 1.01) and that we have enough samples (ESS more than about 2,000).

# first parameter (intercept); discard the first 20,000 samples as burn-in
start <- .2*S + 1
end <-   S
matrix_of_chains <- cbind(
  m1$samples[start:end, 1],
  m2$samples[start:end, 1],
  m3$samples[start:end, 1],
  m4$samples[start:end, 1]
)

# compute r-hat (for intercept)
posterior::rhat(matrix_of_chains)

[1] 1

# compute ess (for intercept)
posterior::ess_bulk(matrix_of_chains)

[1] 3188

posterior::ess_tail(matrix_of_chains)

[1] 7472

Code

colnames(matrix_of_chains) <- paste("Chain", 1:ncol(matrix_of_chains))
gg_df <- matrix_of_chains |>
  as_tibble() |>
  mutate(Iteration = 1:n()) |>
  pivot_longer(cols = `Chain 1`:`Chain 4`, names_to = "Chain", values_to = "Sample") |>
  separate(Chain, into = c("tmp", "Chain Number"), remove = FALSE) |>
  mutate(`Chain Number` = as.numeric(`Chain Number`),
         Chain = reorder(Chain, `Chain Number`, ordered = TRUE))

gg1 <- ggplot(gg_df, aes(x = Iteration, y = Sample, color = Chain)) +
  geom_line() +
  facet_wrap(vars(Chain))

gg2 <- ggplot(gg_df, aes(x = Sample, color = Chain)) +
  geom_density()

library(patchwork)
gg1 + gg2

The code above computes the R-hat and ESS for only the first parameter (intercept)—we should check the rest as well. The table below shows that R-hat is less than 1.01 for every parameter, so we have found no sign that the chains have failed to converge.

Coefficient	R-hat	Converged (R-hat < 1.01)	ESS (bulk)	ESS (tail)
(Intercept)	1	✅	3,188	7,472
rs_age	1	✅	6,803	14,584
rs_educate	1	✅	4,889	10,768
rs_income	1	✅	4,445	9,648
racewhite	1	✅	3,157	7,773

19.4.4 Comparing

Finally, we can compare the posterior mean and 95% credible intervals to the ML estimate and 95% Wald confidence intervals. They are almost identical.

19.4.5 Quantities of Interest

To compute a quantity of interest, we need to use our Bayesian invariance property. We need to transform, then summarize, the simulations.

Here’s an example of the first difference for a “typical” case.

# make X_lo
X_lo <- cbind(
  "constant" = 1, # intercept
  "rs_age" = -0.5, # 1 SD below avg -- see ?arm::rescale
  "rs_educate" = 0,
  "rs_income" = 0,
  "white" = 1 # white indicator = 1 
)

# make X_hi by modifying the relevant value of X_lo
X_hi <- X_lo
X_hi[, "rs_age"] <- 0.5  # 1 SD above avg

# function to compute first difference
fd_fn <- function(beta, hi, lo) {
  plogis(hi%*%beta) - plogis(lo%*%beta)
}

# put the simulations of the coefficients into a matrix
# note 1: each row is one beta-tilde
# note 2: we're discarding the first samples (start and end defined above)
# note 3: we're just stacking the chains on top of each other
beta_tilde <- rbind(
  m1$samples[start:end, ], # chain 1, minus burn-in
  m2$samples[start:end, ], # chain 2, minus burn-in
  m3$samples[start:end, ], # chain 3, minus burn-in
  m4$samples[start:end, ]  # chain 4, minus burn-in
)

# transform simulations of coefficients into simulations of first-difference
# note 1: for clarity, just do this one simulation at a time,
# note 2: i indexes the simulations
fd_tilde <- numeric(nrow(beta_tilde))  # container
for (i in 1:nrow(beta_tilde)) {
  fd_tilde[i] <- fd_fn(beta_tilde[i, ], hi = X_hi, lo = X_lo)
}

# posterior mean
mean(fd_tilde)

[1] 0.167

And we can replicate this with {marginaleffects}

# load packages
library(marginaleffects)

# compute qi with the rescaled columns (see the note in the data chunk)
comparisons(fit, variables = list(rs_age = c(-0.5, 0.5)), 
            newdata = datagrid(grid_type = "mean_or_mode"))


 Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
    0.166       0.02 8.32   <0.001 53.4 0.127  0.205

Term: rs_age
Type: response
Comparison: 0.5 - -0.5

19.5 Summary

There are three important points in this chapter.

MCMC

MCMC offers a powerful and general way to sample from an unnormalized (log-)posterior. The Metropolis algorithm is one such method. Hamiltonian Monte Carlo (HMC) via Stan is even better.

R-hat and burn-in

The first MCMC samples are highly dependent on the starting values. Because of this, you need to:

Discard the samples from a burn-in period.
Run multiple chains and check that R-hat is less than 1.01.

ESS and large samples

The MCMC samples are dependent on the previous samples. Because of this, you need to:

Generate more samples than you would need if they were independent.
Use ESS to understand your effective sample size.