Implementing a model in Stan

The preceding pages developed Hamiltonian Monte Carlo and the No-U-Turn Sampler. Stan is a probabilistic programming language that uses these methods after the analyst declares observed data, unknown parameters, and density contributions. The current reference manual describes how those declarations enter the target. We will read the program by asking, in turn, what is observed, what is unknown, and which terms define the target density.

Consider a constructed study of overt subject pronouns. Code y_i=1 when a sampled clause contains an overt subject pronoun and y_i=0 otherwise. An overt subject pronoun is a pronounced pronominal expression in subject position, as opposed to an unpronounced subject licensed in a language that permits such omission. This binary code does not distinguish pronoun forms or discourse functions. The statistical model is

Y_i\mid\Pi=\pi \sim \operatorname{Bernoulli}(\pi)

with prior

\Pi\sim\operatorname{Beta}(2,2).

The Stan program is

data {
  int<lower=0> n;
  array[n] int<lower=0, upper=1> y;
}

parameters {
  real<lower=0, upper=1> pi;
}

model {
  pi ~ beta(2, 2);
  y ~ bernoulli(pi);
}

Declaring the observed data

data {
  int<lower=0> n;
  array[n] int<lower=0, upper=1> y;
}

The integer n is the number of clauses. The array y contains n observed integers constrained to zero or one. Values in the data block are supplied before sampling and remain fixed during posterior computation.

The constraints check allowable data values. If y contains a two, Stan rejects the input rather than silently treating it as a Bernoulli outcome.

Declaring the unknown parameter

parameters {
  real<lower=0, upper=1> pi;
}

The parameter pi is unknown and must lie between zero and one. Stan internally transforms constrained parameters to an unconstrained scale for sampling and adjusts the density for that transformation.

The support constraint is not the prior, but instead says which values are allowable. The beta statement in the model block says how prior density varies across those values.

Adding density contributions

model {
  pi ~ beta(2, 2);
  y ~ bernoulli(pi);
}

The first statement adds the beta prior log density for pi. The second adds one Bernoulli log likelihood contribution for every element of y.

Stan’s tilde notation is log probability notation in the model block. It does not simulate pi and then simulate y. The program is equivalent, up to constants not involving pi, to

model {
  target += beta_lpdf(pi | 2, 2);
  target += bernoulli_lpmf(y | pi);
}

The variable target accumulates the model’s log probability contributions. Stan samples an unconstrained internal version of pi. Mapping that value back to the interval (0,1) stretches some regions more than others, so Stan adds a Jacobian adjustment, a change-of-variables correction that preserves the intended density. The model contributions and this correction together define the unnormalized target density on the internal scale.

Preparing a matching data object

Suppose the constructed sample has 12 overt pronouns among 20 clauses.

Code
y <- c(rep(1L, 12), rep(0L, 8))
stan_data <- list(n = length(y), y = y)

stopifnot(stan_data$n == 20L)
stopifnot(sum(stan_data$y) == 12L)
stopifnot(all(stan_data$y %in% 0:1))

analytic_posterior <- c(alpha = 2 + sum(y),
                        beta = 2 + length(y) - sum(y))
analytic_posterior

The conjugate update derived earlier gives a \operatorname{Beta}(14,10) benchmark. A successful Stan fit should approximate this posterior within Monte Carlo error. The benchmark checks the implementation for this simple model; most Stan models will not have an analytic comparison.

Matching code to factorization

Stan vectorizes y ~ bernoulli(pi), but the statement still represents

\sum_{i=1}^n \log p(y_i\mid\pi).

If several clauses come from one speaker, concise vectorization does not remove that dependence. The program has asserted conditional independence by adding separate Bernoulli terms with one common probability.

Translate each density statement back into the product or sum it represents, then ask whether that factorization matches the linguistic sampling process.

Separating program validity from model adequacy

Stan can parse a program, calculate gradients, and return draws from the declared posterior. None of those outcomes establishes that independent Bernoulli observations are a suitable model for the clauses.

A completed sampling run is not a validated analysis. The analyst still must check computational behavior, prior implications, posterior predictions, and the linguistic correspondence between data rows and modeled units.

Stan supplies the sampling machinery. The analyst supplies the model and the inferential target. The generated quantities page next separates calculations that define this target from calculations performed after each retained draw.

The program should be stored with the exact data transformation and interface call used to fit it. Identical Stan code can target a different posterior when variables are recoded or rows are filtered.

Check your understanding

  1. Why does y belong in the data block rather than the parameters block?
  2. What is the difference between <lower=0, upper=1> and pi ~ beta(2, 2)?
  3. What does a tilde statement add in the model block?
  4. Why can a program run without representing speaker or document dependence?