0Pricing
R Academy · Lesson

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 — lesson 2 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 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')

Frequently asked questions

Is the “Writing Stan Models in R” lesson free?

Yes — the full text of “Writing Stan Models in R” 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 “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; this is — lesson 2 of 4, so you can start here or from the beginning and move at your own pace.

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

  1. Introduction to Bayesian Thinking
  2. Writing Stan Models in R
  3. MCMC Sampling and Diagnostics
  4. Posterior Predictive Checks
← Back to R Academy