Stan Development

Use this skill when working with Stan models for Bayesian inference and probabilistic programming, particularly when integrating with R (cmdstanr) or Python (cmdstanpy).

Modern Stan Syntax

Array Declarations (IMPORTANT)

Use new array syntax introduced in Stan 2.26+:

// Good - Modern syntax (Stan 2.26+) array[10] int x; // Array of 10 integers array[N] real y; // Array of N reals array[N, M] real z; // 2D array array[N] vector[K] v; // Array of vectors

// Avoid - Old syntax (deprecated) int x[10]; // Old style real y[N]; // Old style real z[N, M]; // Old style vector[K] v[N]; // Old style

Why the change?

• More consistent and readable

• Clearer distinction between arrays and matrix/vector types

• Aligns with modern Stan style

Other Key Syntax

Data types:

// Scalar types int n; real x;

// Vector/Matrix types vector[N] v; // Column vector row_vector[N] rv; // Row vector matrix[N, M] A; // Matrix

// Arrays of vectors/matrices array[K] vector[N] vectors; array[K] matrix[N, M] matrices;

// Constrained types real<lower=0> sigma; // Lower bound real<lower=0, upper=1> theta; // Bounded simplex[K] pi; // Simplex constraint ordered[N] sorted_values; // Ordered constraint

Stan File Organization

Standard Structure

project/ ├── inst/stan/ # Stan models (R packages) │ ├── model.stan # Main model file │ ├── functions/ # Reusable Stan functions │ │ ├── likelihood.stan │ │ ├── priors.stan │ │ └── utils.stan │ └── chunks/ # Reusable code chunks (optional) └── stan/ # Alternative location (Python projects)

Modular Stan Code

Main model includes functions:

functions { #include functions/likelihood.stan #include functions/priors.stan #include functions/utils.stan }

data { int<lower=0> N; array[N] real y; }

parameters { real mu; real<lower=0> sigma; }

model { // Use functions from included files y ~ custom_likelihood(mu, sigma); mu ~ custom_prior(); }

Benefits of modular approach:

• Reusable functions across models

• Easier testing of individual components

• Better organization for complex models

• Cleaner main model files

Stan Integration with R (cmdstanr)

Basic Workflow

library(cmdstanr)

Compile Stan model

model <- cmdstan_model("path/to/model.stan")

Prepare data (list with names matching Stan data block)

stan_data <- list( N = 100, y = rnorm(100) )

Fit model

fit <- model$sample( data = stan_data, chains = 4, parallel_chains = 4, iter_warmup = 1000, iter_sampling = 1000 )

Extract samples

draws <- fit$draws() summary <- fit$summary()

Model Compilation and Caching

Get model path from package

model_path <- system.file("stan", "model.stan", package = "mypackage")

Compile with cmdstanr

model <- cmdstan_model(model_path)

Package-specific model functions

Many packages provide wrapper functions:

my_package_model <- get_model() # Returns compiled model

Inference Algorithms

NUTS sampling (default, most robust)

fit <- model$sample(data = stan_data)

Variational inference (faster, approximate)

fit <- model$variational(data = stan_data)

Pathfinder (fast, approximate)

fit <- model$pathfinder(data = stan_data)

Laplace approximation (very fast, approximate)

fit <- model$laplace(data = stan_data)

Optimization (MAP estimate)

fit <- model$optimize(data = stan_data)

Testing Stan Functions

Exposing Stan Functions to R

In R packages with cmdstanr:

Expose Stan functions for testing

Typically in a package function

expose_stan_functions <- function() { model_path <- system.file("stan", "functions.stan", package = "mypackage")

Note: This typically works on Linux only

cmdstanr::cmdstan_model(model_path, compile_standalone = TRUE) }

Then in tests

test_that("Stan function works correctly", {

Exposed functions are now available in R

result <- stan_function_name(args) expect_equal(result, expected) })

Common pattern in test setup:

tests/testthat/setup.R

if (on_linux()) { expose_stan_functions() }

tests/testthat/test-stan-functions.R

test_that("likelihood function works", { skip_if_not(on_linux()) # Stan function exposure often Linux-only

result <- stan_likelihood(data, params) expect_gt(result, -Inf) })

Stan-R Interface Patterns

Data Preparation

Convert R data structures to Stan format

stan_data_list <- list( N = nrow(data), K = ncol(X), y = data$outcome, X = as.matrix(X),

Arrays use R vectors/matrices directly

group = as.integer(data$group) )

For complex models, packages often provide conversion functions

stan_data <- prepare_stan_data(preprocessed_data, model_spec)

Prior Specification as Data

Empirical Bayes approach - priors as data:

data { // Data int<lower=0> N; array[N] real y;

// Priors as data (empirical Bayes) real prior_mu_mean; real<lower=0> prior_mu_sd; real<lower=0> prior_sigma_alpha; real<lower=0> prior_sigma_beta; }

parameters { real mu; real<lower=0> sigma; }

model { // Use data-specified priors mu ~ normal(prior_mu_mean, prior_mu_sd); sigma ~ gamma(prior_sigma_alpha, prior_sigma_beta); y ~ normal(mu, sigma); }

Benefits:

• Priors can be updated without recompiling

• Easy to specify different priors for different models

• Enables automated prior selection

Distribution Lookup Systems

For packages with multiple distribution options:

R side - get distribution ID

dist_id <- get_distribution_id("lognormal")

Stan side - use distribution

stan_data <- list( distribution = dist_id, # Integer ID

... other data

)

// Stan functions for distribution dispatch real compute_lpdf(real y, int distribution, vector params) { if (distribution == 1) { // Normal return normal_lpdf(y | params[1], params[2]); } else if (distribution == 2) { // Lognormal return lognormal_lpdf(y | params[1], params[2]); } // ... other distributions }

Common Stan Patterns

Vectorization

// Good - vectorized (much faster) y ~ normal(mu, sigma);

// Avoid - loop (slower) for (n in 1:N) { y[n] ~ normal(mu, sigma); }

Efficient Matrix Operations

// Use built-in matrix operations vector[N] mu = X * beta; // Matrix-vector multiplication

// Avoid loops when vectorization possible

Handling Missing Data

data { int<lower=0> N; int<lower=0> N_obs; array[N_obs] int<lower=1, upper=N> obs_idx; array[N_obs] real y_obs; }

parameters { array[N] real y_latent; real mu; real<lower=0> sigma; }

model { // Likelihood for observed data y_obs ~ normal(mu, sigma);

// Prior/model for all data y_latent ~ normal(mu, sigma);

// Constrain observed values y_latent[obs_idx] ~ normal(y_obs, 0.001); // Small noise }

Debugging Stan Models

Common Issues

Divergences:

Increase adapt_delta

fit <- model$sample( data = stan_data, adapt_delta = 0.95 # Default 0.8 )

Slow mixing:

Use non-centered parameterization

Reparameterize hierarchical models

Model diagnostics:

Check convergence

fit$diagnostic_summary()

Check R-hat, ESS

fit$summary(c("Rhat", "ess_bulk", "ess_tail"))

Pairs plot for problematic parameters

bayesplot::mcmc_pairs(fit$draws(), pars = c("mu", "sigma"))

When to Use This Skill

Activate this skill when:

• Writing Stan models

• Integrating Stan with R packages (cmdstanr)

• Testing Stan functions

• Debugging Stan models

• Working with Bayesian hierarchical models

• Implementing custom likelihoods or priors in Stan

This skill provides Stan-specific development patterns. Project-specific model architecture and domain knowledge should remain in project CLAUDE.md files.