MCMC Sampling with Stan

Overview

This skill provides procedural guidance for implementing Markov Chain Monte Carlo (MCMC) sampling using Stan and RStan. It covers the complete workflow from environment setup through model specification to convergence validation, with emphasis on avoiding common pitfalls in hierarchical Bayesian modeling.

When to Use This Skill

• Installing and configuring RStan in R environments
• Writing Stan models for Bayesian inference
• Implementing hierarchical/multilevel models
• Configuring MCMC sampling parameters
• Diagnosing and validating MCMC convergence

Pre-Implementation Checklist

Before writing any code, verify the following in order:

1. System Dependencies First

Critical: Always check system dependencies before attempting R package installation.

# Check for C++ compiler and build tools
which g++ || which clang++
dpkg -l | grep -E "(build-essential|r-base-dev)" # Debian/Ubuntu

Required system packages for RStan:

• build-essential or equivalent C++ toolchain
• r-base-dev for R development headers
• Sufficient memory (Stan compilation is memory-intensive)

2. RStan Installation Strategy

Install RStan with explicit version pinning when specified:

# If specific version required (e.g., 2.32.7)
install.packages("rstan", repos = "https://cloud.r-project.org/")

# Verify installation and version
library(rstan)
packageVersion("rstan")

Common installation failures:

• Missing C++ compiler → Install build-essential first
• Compilation errors → Check g++ version compatibility
• Memory exhaustion → Close other applications, increase swap

3. Environment Verification

Before proceeding with model development:

# Verify RStan loads without errors
library(rstan)

# Test Stan compilation with minimal model
stan_model(model_code = "parameters { real x; } model { x ~ normal(0,1); }")

Stan Model Specification

Hierarchical Model Structure

For hierarchical Bayesian models, follow this block structure:

data {
  // Declare all input data with types and constraints
  int<lower=0> N;           // Number of observations
  int<lower=0> y[N];        // Observed counts
  int<lower=0> n[N];        // Trial sizes
}

parameters {
  // Hyperparameters (population-level)
  real<lower=0> alpha;
  real<lower=0> beta;

  // Group-level parameters
  real<lower=0, upper=1> theta[N];
}

model {
  // Hyperprior (if using custom/improper priors)
  // Example: p(α,β) ∝ (α+β)^(-5/2)
  target += -2.5 * log(alpha + beta);

  // Prior for group parameters
  theta ~ beta(alpha, beta);

  // Likelihood
  y ~ binomial(n, theta);
}

Custom Prior Implementation

When implementing custom or improper priors:

Document the mathematical transformation: Always comment the relationship between the mathematical prior and the target += statement.

MCMC Sampling Configuration

Recommended Control Parameters

For hierarchical models with potential sampling difficulties:

fit <- sampling(
  model,
  data = stan_data,
  chains = 4,                    # Multiple chains for convergence assessment
  iter = 100000,                 # Total iterations (including warmup)
  seed = 1,                      # Reproducibility
  control = list(
    adapt_delta = 0.95,          # Increase for divergent transitions
    max_treedepth = 15           # Increase for complex posteriors
  )
)

Parameter Selection Rationale

Warmup Considerations

• Default warmup is 50% of iter
• For 100,000 iterations: 50,000 warmup + 50,000 sampling per chain
• Explicitly specify if different warmup proportion needed:

sampling(..., iter = 100000, warmup = 25000)  # 25% warmup

Convergence Diagnostics

Required Validation Steps

Always verify these after sampling:

# 1. Check Rhat (potential scale reduction factor)
summary(fit)$summary[, "Rhat"]
# Target: All Rhat < 1.01 (ideally < 1.005)

# 2. Check effective sample size
summary(fit)$summary[, "n_eff"]
# Target: n_eff > 400 for reliable inference

# 3. Check for divergent transitions
get_num_divergent(fit)
# Target: 0 divergent transitions

# 4. Check treedepth saturation
get_num_max_treedepth(fit)
# Target: 0 or minimal saturation

Diagnostic Interpretation

Common Pitfalls and Solutions

1. Installation Order Errors

Wrong: Install RStan → Installation fails → Check dependencies Right: Check dependencies → Install prerequisites → Install RStan

2. Missing Convergence Checks

Wrong: Run sampling → Extract means → Report results Right: Run sampling → Check Rhat/n_eff/divergences → Validate → Extract results

3. Improper Prior Issues

When using improper priors like p(α,β) ∝ (α+β)^(-5/2):

• Ensure parameters are constrained (<lower=0>)
• Verify posterior is proper (converges)
• Check for boundary behavior near zero

4. Command Syntax in Shell

Avoid piping R output with shell redirection that causes parsing errors:

# Problematic
Rscript -e "code" 2>&1 | grep ...

# Safer
Rscript -e "code" > output.txt 2>&1
grep ... output.txt

Verification Strategy

Stepwise Verification Approach

1. Environment: Verify RStan installation before model development
2. Compilation: Test Stan model compiles without errors
3. Sampling: Run with reduced iterations first to catch issues
4. Diagnostics: Check all convergence metrics before trusting results
5. Results: Extract posteriors only after validation passes

Sanity Checks for Results

• Posterior means should be plausible given prior and data
• 95% credible intervals should have reasonable width
• Multiple chains should show similar posteriors
• Results should be consistent across different seeds (approximately)

References

For detailed Stan documentation and examples, consult:

• Stan User's Guide: Model specification and best practices
• RStan Getting Started: Installation and basic usage
• Stan Functions Reference: Available distributions and functions