R Academy · Lesson

MCMC Sampling and Diagnostics

Run sampling, inspect chains, and interpret Rhat and ESS diagnostics.

MCMC Sampling and Diagnostics is a free R Academy lesson on CoddyKit — lesson 3 of 4. You can read the complete lesson below for free — then practise it hands-on in the browser with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, one of 4 lessons in the course, and your progress syncs across the web and the CoddyKit app.

What Is MCMC?

Markov Chain Monte Carlo (MCMC) is a family of algorithms for drawing samples from a probability distribution when direct sampling is impossible. In Bayesian statistics, MCMC samples from the posterior distribution P(parameters | data).

RStan implements the No-U-Turn Sampler (NUTS), a state-of-the-art MCMC algorithm.

A Simple Stan Model

A Stan model is a text block defining data types, parameters, and the log-posterior. The simplest model estimates the mean of a normal distribution with known variance.

# library(rstan)
#
# stan_code <- '
# data {
#   int<lower=0> N;
#   vector[N] y;
# }
# parameters {
#   real mu;
#   real<lower=0> sigma;
# }
# model {
#   mu    ~ normal(0, 10);   // prior
#   sigma ~ exponential(1);   // prior
#   y     ~ normal(mu, sigma); // likelihood
# }
# '

Calling stan() to Sample

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) compiles the model (once), then draws samples. With iter=2000 and warmup=1000, each chain produces 1000 post-warmup samples — 4000 total across 4 chains.

# library(rstan)
# options(mc.cores = parallel::detectCores())
#
# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# stan_data <- list(N = length(y), y = y)
#
# fit <- stan(
#   model_code = stan_code,
#   data       = stan_data,
#   chains     = 4,
#   iter       = 2000,
#   warmup     = 1000,
#   seed       = 42
# )

print(fit) — The Summary Table

print(fit) shows a summary table for each parameter with posterior mean, standard deviation, quantiles, Rhat, and n_eff. These two diagnostics are the first things to check.

# print(fit)
#
# Example output:
#       mean se_mean   sd  2.5%   25%   50%   75%  97.5%  n_eff Rhat
# mu    2.05    0.00 0.15  1.76  1.95  2.05  2.15   2.34   3842    1
# sigma 0.22    0.00 0.06  0.13  0.18  0.21  0.25   0.37   3521    1
# lp__  4.38    0.02 1.01  1.60  3.91  4.71  5.18   5.50   2148    1

The Rhat Convergence Criterion

Rhat (potential scale reduction factor) compares variance within chains to variance between chains. Values close to 1.0 indicate that all chains have converged to the same distribution.

Rhat < 1.01 — converged (current standard)
Rhat > 1.01 — chains have not mixed; run more iterations
Rhat > 1.1 — serious convergence problem

# Check Rhat for all parameters:
# s <- summary(fit)$summary
# rhat_vals <- s[, 'Rhat']
# cat('Max Rhat:', max(rhat_vals, na.rm = TRUE), '
')
# if (any(rhat_vals > 1.01, na.rm = TRUE)) {
#   warning('Convergence issue detected!')
# } else {
#   cat('All Rhat < 1.01 — chains converged
')
# }

n_eff — Effective Sample Size

n_eff (effective sample size) accounts for autocorrelation between successive MCMC samples. Correlated samples carry less information than independent ones.

n_eff close to total iterations — near-independent samples, excellent
n_eff / total_samples > 0.1 — generally acceptable
Very low n_eff — high autocorrelation; consider reparameterizing the model

# s <- summary(fit)$summary
# n_eff_vals <- s[, 'n_eff']
# total_samples <- 4 * 1000   # chains * post-warmup iter
# ratio <- n_eff_vals / total_samples
# cat('n_eff ratio (mu)   :', round(ratio['mu'], 2), '
')
# cat('n_eff ratio (sigma):', round(ratio['sigma'], 2), '
')

traceplot() for Visual Convergence

traceplot(fit, pars = 'mu') plots the sampled values of mu across iterations for each chain. Converged chains look like a fuzzy caterpillar — all chains overlapping with no trends or drifts. Divergent chains wander or stay separate.

# library(rstan)
#
# traceplot(fit, pars = c('mu', 'sigma'), inc_warmup = FALSE)
#
# Good traceplot characteristics:
# - All 4 chains overlapping completely (same range)
# - No visible drift or trend
# - Rapid mixing (values jump around quickly)
# - No flat regions (stuck sampler)
cat('A healthy traceplot looks like a fuzzy caterpillar
')

pairs() for Posterior Correlations

pairs(fit, pars = c('mu', 'sigma')) shows a scatter plot matrix of posterior samples. It reveals correlations between parameters and highlights divergent transitions (plotted in red) which indicate regions the sampler struggles with.

# pairs(fit, pars = c('mu', 'sigma'))
#
# What to look for:
# - Elliptical clouds: mild correlation (OK)
# - Banana / funnel shapes: reparameterization needed
# - Red dots (divergences): geometry problem in posterior
#   => increase adapt_delta: stan(..., control=list(adapt_delta=0.95))
cat('Red dots in pairs() indicate divergent transitions — investigate!
')

Extracting Posterior Samples

extract(fit, pars = 'mu')$mu returns a numeric vector of all post-warmup samples for mu. Use these samples to compute any posterior summary: mean, credible intervals, probability of a condition.

# mu_samples <- extract(fit, pars = 'mu')$mu
# cat('Posterior mean :', mean(mu_samples), '
')
# cat('95% CI:', quantile(mu_samples, c(0.025, 0.975)), '
')
# cat('P(mu > 2):', mean(mu_samples > 2), '
')
# hist(mu_samples, main = 'Posterior of mu', xlab = 'mu', col = 'steelblue')

Launching ShinyStan for Interactive Diagnostics

shinystan::launch_shinystan(fit) opens an interactive Shiny app with traceplots, posterior distributions, pairs plots, and NUTS diagnostics all in one place. It is the most comprehensive tool for exploring an RStan fit.

# install.packages('shinystan')
# library(shinystan)
#
# shinystan::launch_shinystan(fit)
#
# ShinyStan tabs:
# - Diagnose: Rhat, n_eff, divergences, energy
# - Explore:  marginal posteriors, scatter plots
# - Model:    Stan code, data
# - NUTS:     step size, tree depth per chain

Common Convergence Fixes

When Rhat > 1.01 or you see divergences:

Increase iter and warmup
Increase adapt_delta towards 1.0 (e.g., 0.95) in control
Reparameterize — use non-centered parameterization for hierarchical models
Tighten priors if they are too diffuse
Check the data — outliers or scale differences cause sampler problems

# Re-run with higher adapt_delta to reduce divergences:
# fit2 <- stan(
#   model_code = stan_code,
#   data       = stan_data,
#   chains     = 4,
#   iter       = 4000,
#   warmup     = 2000,
#   control    = list(adapt_delta = 0.95, max_treedepth = 12),
#   seed       = 42
# )

Quick Check: Rhat Threshold

What is the current standard threshold for Rhat that indicates a Stan model has converged?

MCMC Sampling and Diagnostics Recap

Key RStan MCMC workflow:

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) fits the model
print(fit) shows Rhat and n_eff — the primary convergence diagnostics
Rhat < 1.01 and n_eff / total > 0.1 indicate a well-behaved sample
traceplot() — visual mixing check; pairs() — reveals posterior geometry issues
extract(fit, pars='mu')$mu — access raw posterior samples
shinystan::launch_shinystan(fit) — comprehensive interactive diagnostics

Frequently asked questions

Is the “MCMC Sampling and Diagnostics” lesson free?

Yes — the full text of “MCMC Sampling and Diagnostics” is free to read here on the web, and the R Academy course includes 4 lessons in total. To practise it interactively (a built-in code editor and a 24/7 AI tutor) and unlock the rest of the R Academy course, upgrade to CoddyKit PRO.

What will I learn in “MCMC Sampling and Diagnostics”?

Run sampling, inspect chains, and interpret Rhat and ESS diagnostics. You practise R Academy with hands-on code you run directly in the browser, and a 24/7 AI tutor answers your questions as you work through the lesson.

Do I need any experience to start R Academy?

No prior experience is required. R Academy on CoddyKit is structured for beginners through advanced learners; this is — lesson 3 of 4, so you can start here or from the beginning and move at your own pace.

How long does the “MCMC Sampling and Diagnostics” lesson take?

Most CoddyKit lessons take about 5–10 minutes. Each one is bite-sized and interactive, so you make steady progress and pick up exactly where you left off across the web and the app.

Can I write and run code in this R Academy lesson?

Yes. Every R Academy lesson includes a built-in code editor, so you write and run real code right in your browser and get instant AI feedback — no local setup required.

All lessons in this course

← Back to R Academy