Writing Stan Models in R
Define data blocks, parameters, and the model block in Stan syntax.
Writing Stan Models in R is a free R Academy lesson on CoddyKit. This is lesson 2 of 4. Any 3 lessons of this learning path are free to read in full — after that, CoddyKit PRO unlocks every lesson, plus hands-on practice with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, and your progress syncs across the web and the CoddyKit app. The R Academy course includes 4 lessons in total.
What Is Stan?
Stan is a probabilistic programming language for Bayesian statistical modelling. RStan is the R interface to Stan. You write the model in Stan's language (C++-like), and Stan compiles it to efficient C++ code that runs Hamiltonian Monte Carlo (HMC) sampling to approximate the posterior.
# Stan installation check
library(rstan)
# Check version
cat('RStan version:', as.character(packageVersion('rstan')), '\n')
# Enable parallel chains
options(mc.cores = parallel::detectCores())
# Reuse compiled models across sessions
rstan_options(auto_write = TRUE)
cat('Stan ready. Cores:', parallel::detectCores(), '\n')Stan Model Structure
A Stan model has up to six named blocks: data, transformed data, parameters, transformed parameters, model, and generated quantities. The three essential blocks are data, parameters, and model.
library(rstan)
# Stan model as a character string in R
stan_code <- '
data {
int<lower=0> N; // number of observations
vector[N] x; // predictor
vector[N] y; // response
}
parameters {
real alpha; // intercept
real beta; // slope
real<lower=0> sigma; // noise (must be positive)
}
model {
// Priors
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ exponential(1);
// Likelihood
y ~ normal(alpha + beta * x, sigma);
}
'
cat('Stan model defined as a string in R\n')
cat('Blocks: data, parameters, model\n')data{}: Declaring Input Data
The data block declares all external data the model receives from R. Types include int, real, vector[N], matrix[M,N], and array. Constraints like <lower=0> are checked at runtime.
library(rstan)
# Comprehensive data block examples
data_block_examples <- '
data {
// Scalars
int<lower=1> N; // at least 1 observation
int<lower=2> K; // at least 2 groups
// Constrained scalars
real<lower=0, upper=1> rate; // probability
// Vectors
vector[N] y; // continuous response
array[N] int<lower=0, upper=1> z; // binary outcomes
// Matrix
matrix[N, K] X; // design matrix
// Integer array
array[N] int group; // group membership 1..K
}
'
cat(data_block_examples)parameters{}: Parameter Types
The parameters block declares the unknowns Stan will sample. Constraints in this block define the parameter space: <lower=0> for positive quantities, <lower=0, upper=1> for probabilities. Stan automatically applies log-probability Jacobian corrections for constrained parameters.
library(rstan)
# Common parameter declarations
param_examples <- '
parameters {
// Unconstrained
real mu; // mean
vector[K] beta; // regression coefficients
// Positive (sigma, lambda, variance)
real<lower=0> sigma;
real<lower=0> lambda;
// Probability
real<lower=0, upper=1> theta;
// Simplex (sums to 1, for mixture weights)
simplex[K] pi;
// Correlation matrix
corr_matrix[K] Omega;
// Cholesky factor of covariance
cholesky_factor_cov[K] L_Sigma;
}
'
cat(param_examples)model{}: Priors and Likelihood
The model block accumulates the log-posterior via ~ (tilde) syntax. y ~ normal(mu, sigma) is shorthand for adding the log of the Normal PDF to target. You can write explicit log-probability increments with target += normal_lpdf(y | mu, sigma).
library(rstan)
model_block_example <- '
model {
// --- Priors ---
mu ~ normal(0, 10); // weakly informative
sigma ~ cauchy(0, 2.5); // half-Cauchy for scale
beta ~ normal(0, 1); // standardised coefficients
// --- Likelihood ---
// Tilde notation (most common)
y ~ normal(mu + X * beta, sigma);
// Equivalent explicit notation:
// target += normal_lpdf(y | mu + X * beta, sigma);
// For loop (less common but valid)
// for (i in 1:N)
// target += normal_lpdf(y[i] | mu, sigma);
}
'
cat(model_block_example)Fitting a Simple Normal Model
Here is a complete minimal Stan workflow: define the model string, prepare the data list, call stan(), and inspect results with print(). Stan compiles the model the first time and caches it for subsequent runs.
library(rstan)
# Stan model: estimate mean and SD of a normal distribution
normal_model <- '
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.1);
y ~ normal(mu, sigma);
}
'
# Simulate data
set.seed(42)
y_data <- rnorm(50, mean = 5, sd = 2)
# Fit the model
fit <- stan(
model_code = normal_model,
data = list(N = length(y_data), y = y_data),
chains = 2,
iter = 1000,
warmup = 500,
refresh = 0 # suppress iteration output
)
print(fit, pars = c('mu', 'sigma'))Linear Regression in Stan
Stan makes Bayesian linear regression straightforward: specify normal priors on coefficients and an exponential or half-Cauchy prior on sigma. The posterior gives the full uncertainty distribution for each coefficient, not just point estimates.
library(rstan)
lin_reg_model <- '
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
alpha ~ normal(0, 10);
beta ~ normal(0, 10);
sigma ~ exponential(1);
y ~ normal(alpha + beta * x, sigma);
}
generated quantities {
vector[N] y_rep; // posterior predictive
for (i in 1:N)
y_rep[i] = normal_rng(alpha + beta * x[i], sigma);
}
'
set.seed(7)
n <- 80
x <- rnorm(n); y <- 2 + 3 * x + rnorm(n, 0, 1)
fit <- stan(model_code = lin_reg_model,
data = list(N = n, x = x, y = y),
chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('alpha', 'beta', 'sigma'))transformed parameters{} Block
The transformed parameters block computes derived quantities from sampled parameters. These are useful for parameterising models in a numerically stable way (e.g. Cholesky decompositions, log-transformations) while keeping the primary parameters interpretable.
library(rstan)
# Model with transformed parameters
trans_param_example <- '
data {
int<lower=0> N;
array[N] int<lower=0> y; // counts
}
parameters {
real log_lambda; // work in log space for stability
}
transformed parameters {
real<lower=0> lambda;
lambda = exp(log_lambda); // transform back to original scale
}
model {
log_lambda ~ normal(1, 2); // prior on log scale
y ~ poisson(lambda);
}
'
set.seed(42)
y_counts <- rpois(30, lambda = 5)
fit_poisson <- stan(
model_code = trans_param_example,
data = list(N = length(y_counts), y = y_counts),
chains = 2, iter = 1000, refresh = 0
)
print(fit_poisson, pars = c('log_lambda', 'lambda'))generated quantities{} Block
The generated quantities block runs after sampling to compute additional quantities: posterior predictive samples (y_rep), log-likelihood for LOO-CV, or transformed parameters for summaries. Values here are sampled from the posterior predictive distribution.
library(rstan)
# Use generated quantities for posterior predictive checks
model_with_gq <- '
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.5);
y ~ normal(mu, sigma);
}
generated quantities {
vector[N] y_rep; // replicated datasets
real mean_y_rep; // mean of replicated data
for (i in 1:N)
y_rep[i] = normal_rng(mu, sigma);
mean_y_rep = mean(y_rep);
}
'
set.seed(1)
y_obs <- rnorm(40, 3, 1.5)
fit <- stan(model_code = model_with_gq,
data = list(N = length(y_obs), y = y_obs),
chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('mu', 'sigma', 'mean_y_rep'))Passing Data from R to Stan
The data argument to stan() is a named R list. Names must exactly match the variable names declared in the Stan data{} block. Vectors become Stan vector[N]; integers become int; R matrices become Stan matrix[M,N].
library(rstan)
# Data preparation: names must match Stan data block exactly
set.seed(42)
n <- 60
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1.5 + 2 * x1 - 0.8 * x2 + rnorm(n, 0, 0.5)
# Build the design matrix
X <- cbind(x1, x2) # 60 x 2 matrix
# Named list passed to stan(data = ...)
stan_data <- list(
N = n, # int
K = ncol(X), # int
X = X, # matrix[N, K]
y = y # vector[N]
)
cat('Stan data list elements:\n')
for (nm in names(stan_data)) {
cat(' ', nm, ': class =', class(stan_data[[nm]]),
'dim =', paste(dim(stan_data[[nm]]), collapse = 'x'),
'\n')
}Viewing the Compiled Model
After fitting, use print(fit) to see posterior summaries (mean, se_mean, sd, quantiles, Rhat, n_eff) for all parameters. Use stan_plot(fit) for a visual interval plot and traceplot(fit) to assess chain mixing.
library(rstan)
# Reuse the simple normal model from scene 6
normal_model <- '
data { int<lower=0> N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model {
mu ~ normal(0, 10);
sigma ~ exponential(0.1);
y ~ normal(mu, sigma);
}
'
set.seed(5)
fit <- stan(model_code = normal_model,
data = list(N = 50, y = rnorm(50, 7, 3)),
chains = 2, iter = 1000, refresh = 0)
# Detailed summary table
print(fit)
# Extract as data frame
posterior_df <- as.data.frame(fit)
cat('\nPosterior samples shape:', nrow(posterior_df),
'rows x', ncol(posterior_df), 'cols\n')Quick Check
In a Stan model, you want to constrain a parameter to be strictly positive (e.g. a standard deviation). Which declaration is correct?
Recap: Writing Stan Models
Key takeaways:
- Stan model has three essential blocks:
data{},parameters{},model{} - Constraints:
<lower=0>,<upper=1>,<lower=0, upper=1>for probabilities - Types:
int,real,vector[N],matrix[M,N],simplex[K] - Model block: use tilde syntax
y ~ normal(mu, sigma)for priors and likelihood transformed parameters{}for derived quantities;generated quantities{}for post-sampling- Pass data as a named R list matching Stan block names exactly
rstan_options(auto_write = TRUE)caches compiled models
library(rstan)
# Stan model skeleton
model_skeleton <- '
data { int N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model { mu ~ normal(0,10); sigma ~ exponential(1); y ~ normal(mu, sigma); }
'
cat(model_skeleton)
cat('\nFit with: stan(model_code = model_skeleton, data = list(N=..., y=...), chains=4)\n')
Learn R with an AI tutor — free
Write and run real code in your browser, get instant help from a 24/7 AI tutor, and pick up where you left off on the web or in the app.
- Courses
- 43
- Lessons
- 159
Frequently Asked Questions
Is the “Writing Stan Models in R” lesson free?
Yes — any 3 lessons of the R Academy learning path, including “Writing Stan Models in R”, are free to read in full here on the web. After that, CoddyKit PRO unlocks every lesson, plus interactive practice with a built-in code editor and a 24/7 AI tutor. The R Academy course includes 4 lessons in total.
What will I learn in “Writing Stan Models in R”?
Define data blocks, parameters, and the model block in Stan syntax. 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, so you can start here or from the beginning and move at your own pace. This is lesson 2 of 4.
How long does the “Writing Stan Models in R” 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
- Introduction to Bayesian Thinking
- Writing Stan Models in R
- MCMC Sampling and Diagnostics
- Posterior Predictive Checks