An introduction to Stan with R
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
Stan is a probabilistic programming language for specifying statistical models. Stan provides full Bayesian inference for continuousvariable models through Markov Chain Monte Carlo methods such as the NoUTurn sampler, an adaptive form of Hamiltonian Monte Carlo sampling. Penalized maximum likelihood estimates are calculated using optimization methods such as the limited memory BroydenFletcherGoldfarbShanno algorithm. Stan can be called through R using the rstan
package, and through Python using the pystan
package. Both interfaces support sampling and optimizationbased inference with diagnostics and posterior analysis. In this talk it is shown a brief glance about the main properties of Stan. It is also shown a couple of examples: the first one related with a simple Bernoulli model and the second one, about a LotkaVolterra model based on ordinary differential equations.
What is Stan?
 Stan was created by Andrew Gelman and Bob Carpenter, with a development team consisting of 34 members.
 Stan is named in honor of Stanislaw Ulam (19091984), coinventor of the Monte Carlo method.
 Stan is an imperative probabilistic programming language.
 A Stan program defines a probability model.
 It declares data and (constrained) parameter variables.
 It defines log posterior (or penalized likelihood).
 Stan inference: fits model to data and makes predictions.
 It can use Markov Chain Monte Carlo (MCMC) for full Bayesian inference.
 Or Variational Bayesian (VB) for approximate Bayesian inference.
 Or Maximum Likelihood Estimation (MLE) for penalized maximum likelihood estimation.
What does Stan compute?
 Draws from posterior distributions.
 MCMC sampling.
 Sequence of draws \(\theta_{(1)} ,\theta_{(2)}, \ldots, \theta_{(M)}\) where each draw \(\theta_{(i)}\) is marginally distributed according to the posterior \(p(\thetay)\).
 Plot with histograms, kernel density estimates, etc.
You can see a quick look about rstan
in its original wiki page.
Installing rstan
To run Stan in R, it is necessary to install rstan
and a C++ compiler. On Windows, this means that Rtools is required, and you have to check whether the path (in Windows) is correctly fixed for all its binaries. If it is not the case, write in R:
install.packages(devtools) library(devtools) Sys.setenv(PATH = paste("C:\\Rtools\\bin", Sys.getenv("PATH"), sep=";")) Sys.setenv(PATH = paste("C:\\Rtools\\mingw_64\\bin", Sys.getenv("PATH"), sep=";"))
Finally, install rstan
:
install.packages(rstan)
For more information about the frameworks which work with Stan (e.g.: Python) check this link.
Basic syntax in Stan
Defining the model
A Stan model is defined by six program blocks:
 Data (required).
 Transformed data.
 Parameters (required).
 Transformed parameters.
 Model (required).
 Generated quantities.
The data block reads external information.
data { int N; int x[N]; int offset; }
The transformed data block allows for preprocessing of the data.
transformed data { int y[N]; for (n in 1:N) y[n] = x[n]  offset; }
The parameters block defines the sampling space.
parameters { real<lower=0> lambda1; real<lower=0> lambda2; }
The transformed parameters block allows for parameter processing before the posterior is computed.
transformed parameters { real<lower=0> lambda; lambda = lambda1 + lambda2; }
In the model block we define our posterior distributions.
model { y ~ poisson(lambda); lambda1 ~ cauchy(0, 2.5); lambda2 ~ cauchy(0, 2.5); }
Lastly, the generated quantities block allows for postprocessing.
generated quantities { int x_predict; x_predict = poisson_rng(lambda) + offset; }
Types
Stan has two primitive types and both can be bounded.
 int is an integer type.
 real is a floating point precision type.
int<lower=1> N; real<upper=5> alpha; real<lower=1,upper=1> beta; real gamma; real<upper=gamma> zeta;
Reals extend to linear algebra types.
vector[10] a; // Column vector matrix[10, 1] b; row_vector[10] c; // Row vector matrix[1, 10] d;
Arrays of integers, reals, vectors, and matrices are available.
real a[10]; vector[10] b; matrix[10, 10] c;
Stan also implements a variety of constrained types.
simplex[5] theta; // sum(theta) = 1 ordered[5] o; // o[1] < ... < o[5] positive_ordered[5] p; corr_matrix[5] C; // Symmetric and cov_matrix[5] Sigma; // positivedefinite
More about Stan
All the typical control and loop statements are available, too.
if/then/else for (i in 1:I) while (i < I)
There are two ways to modify the posterior.
y ~ normal(0, 1); target += normal_lpdf(y  0, 1); # Deprecated in new versions of Stan: increment_log_posterior(log_normal(y, 0, 1))
And many sampling statements are vectorized.
parameters { real mu[N]; real<lower=0> sigma[N]; } model { // for (n in 1:N) // y[n] ~ normal(mu[n], sigma[n]); y ~ normal(mu, sigma); // Vectorized version }
The Bayesian approach
Probability is epistemic. For instance, John Stuart Mill (Logic 1882, Part III, Ch. 2) said:
[…] the probability of an event is not a quality of the event itself, but a mere name for the degree of ground which we, or someone else, have for expecting it. Every event is in itself certain, not probable; if we knew all, we should either know positively that it will happen, or positively that it will not.
[…] its probability to us means the degree of expectation of its occurrence, which we are warranted in entertaining by our present evidence.
Probabilities quantify uncertainty and we can consider that statistical reasoning is counterfactual.
Bayesian example with Stan: repeated binary trial model
As a first real approach to Stan and its syntax, we will start solving a small example in which the objective is, given a random sample drawn from a Bernoulli population, to estimate the posterior distribution of the missing parameter \(\theta \in \lbrack 0,1]\) (chance of success).
Step 1: Structure of the problem
In this example we will consider the following structure:

Data:
 \(N\in \{0,1,\ldots \}\), number of trials.
 \(y_{n}\in \{0,1\}\), result of trial n (known, modeled data).

Parameter:
\[\theta \in \lbrack 0,1]\]
 Prior distribution
\[p(\theta) = \mathrm{Uniform}(\theta0,1) = 1\]
 Likelihood
\[p(y\theta )=\prod_{n=1}^{N}\mathrm{Bernoulli}(y_{n}\theta) = \prod_{n=1}^{N}\theta ^{y_{n}}(1\theta )^{1y_{n}}\]
 Posterior distribution
\[p(\theta y)\propto p(\theta )p(y\theta )\]
Step 2: Stan Program
We create the Stan program which we will call from R.
bern.stan = " data { int<lower=0> N; // number of trials int<lower=0, upper=1> y[N]; // success on trial n } parameters { real<lower=0, upper=1> theta; // chance of success } model { theta ~ uniform(0, 1); // prior y ~ bernoulli(theta); // likelihood } "
Step 3: The data
In this case, instead of using a given data set, we will simulate a random sample to use in our example.
# Generate data theta = 0.30 N = 20 y = rbinom(N, 1, 0.3) y
## [1] 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
Calculate MLE as the sample mean from data:
sum(y) / N
## [1] 0.25
Step 4: Bayesian Posterior using rstan
The final step is to obtain our estimation using Stan from R.
library(rstan)
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling ## options(mc.cores = parallel::detectCores()). ## To avoid recompilation of unchanged Stan programs, we recommend calling ## rstan_options(auto_write = TRUE)
fit = stan(model_code=bern.stan, data=list(y=y, N=N), iter=5000)
## ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 1). ## Chain 1: ## Chain 1: Gradient evaluation took 7e06 seconds ## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. ## Chain 1: Adjust your expectations accordingly! ## Chain 1: ## Chain 1: ## Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 1: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 1: ## Chain 1: Elapsed Time: 0.012914 seconds (Warmup) ## Chain 1: 0.013376 seconds (Sampling) ## Chain 1: 0.02629 seconds (Total) ## Chain 1: ## ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 2). ## Chain 2: ## Chain 2: Gradient evaluation took 2e06 seconds ## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds. ## Chain 2: Adjust your expectations accordingly! ## Chain 2: ## Chain 2: ## Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 2: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 2: ## Chain 2: Elapsed Time: 0.01289 seconds (Warmup) ## Chain 2: 0.012867 seconds (Sampling) ## Chain 2: 0.025757 seconds (Total) ## Chain 2: ## ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 3). ## Chain 3: ## Chain 3: Gradient evaluation took 3e06 seconds ## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. ## Chain 3: Adjust your expectations accordingly! ## Chain 3: ## Chain 3: ## Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 3: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 3: ## Chain 3: Elapsed Time: 0.012875 seconds (Warmup) ## Chain 3: 0.01241 seconds (Sampling) ## Chain 3: 0.025285 seconds (Total) ## Chain 3: ## ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 4). ## Chain 4: ## Chain 4: Gradient evaluation took 3e06 seconds ## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. ## Chain 4: Adjust your expectations accordingly! ## Chain 4: ## Chain 4: ## Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 4: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 4: ## Chain 4: Elapsed Time: 0.012823 seconds (Warmup) ## Chain 4: 0.014169 seconds (Sampling) ## Chain 4: 0.026992 seconds (Total) ## Chain 4:
print(fit, probs=c(0.1, 0.9))
## Inference for Stan model: 6dcfbccbf2f063595ccc9b93f383e221. ## 4 chains, each with iter=5000; warmup=2500; thin=1; ## postwarmup draws per chain=2500, total postwarmup draws=10000. ## ## mean se_mean sd 10% 90% n_eff Rhat ## theta 0.27 0.00 0.09 0.16 0.39 3821 1 ## lp__ 13.40 0.01 0.73 14.25 12.90 3998 1 ## ## Samples were drawn using NUTS(diag_e) at Tue Jan 22 11:25:43 2019. ## For each parameter, n_eff is a crude measure of effective sample size, ## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
# Extracting the posterior draws theta_draws = extract(fit)$theta # Calculating posterior mean (estimator) mean(theta_draws)
## [1] 0.2715866
# Calculating posterior intervals quantile(theta_draws, probs=c(0.10, 0.90))
## 10% 90% ## 0.1569165 0.3934832
theta_draws_df = data.frame(list(theta=theta_draws)) plotpostre = ggplot(theta_draws_df, aes(x=theta)) + geom_histogram(bins=20, color="gray") plotpostre
RStan: MAP, penalized MLE
Stan’s optimization for estimation; two views:
 max posterior mode, also known as max a posteriori (MAP).
 max penalized likelihood (MLE).
library(rstan) N = 5 y = c(0, 1, 1, 0, 0) model = stan_model(model_code=bern.stan) mle = optimizing(model, data=c("N", "y")) print(mle, digits=2)
## $par ## theta ## 0.4 ## ## $value ## [1] 3.4 ## ## $return_code ## [1] 0
LotkaVolterra’s Model
 Lotka (1925) and Volterra (1926) formulated parametric differential equations that characterize the oscillating populations of predators and preys.
 A statistical model to account for measurement error and unexplained variation uses the deterministic solutions to the LotkaVolterra equations as expected population sizes.
 Full Bayesian inference may be used to estimate future (or past) populations.
 Stan is used to encode the statistical model and perform full Bayesian inference to solve the inverse problem of inferring parameters from noisy data.
See full reference: Stan and LotkaVolterra models.
In this example, we want to fit the model to Canadian lynx predator and snowshoe hare prey with respective populations between 1900 and 1920, based on the number of pelts collected annually by the Hudson’s Bay Company.
Notation and mathematical model
We denote \(u(t)\) and \(v(t)\) as the prey and predator population respectively. The differential equations that relate them are:
\[\frac{d}{dt}u = (\alpha \beta v)u\]
\[\frac{d}{dt}v=(\gamma +\delta u)v\]
Here:
 \(\alpha\): prey growth, intrinsic.
 \(\beta\): prey shrinkage due to predation.
 \(\gamma\): predator shrinkage, intrinsic.
 \(\delta\): predator growth from predation.
LotkaVolterra in Stan
real[] dz_dt(data real t, // time real[] z, // system state real[] theta, // parameters data real[] x_r, // real data data int[] x_i) // integer data { real u = z[1]; // extract state real v = z[2]; real alpha = theta[1]; real beta = theta[2]; real gamma = theta[3]; real delta = theta[4]; real du_dt = (alpha  beta * v) * u; real dv_dt = (gamma + delta * u) * v; return { du_dt, dv_dt }; }
Known variables are observed:
 \(y_{n,k}\): denotes species \(k\) at times \(t_{n}\) for \(n \in 0:N\)
Unknown variables must be inferred (inverse problem):
 Initial state: \(z_{k}^{init}\): initial population for \(k\).
 Subsequent states \(z_{n,k}\): population \(k\) at time \(t_{n}\).
 Parameters \(\alpha\), \(\beta\), \(\gamma\), \(\delta > 0\).
Likelihood assumes errors are proportional (not additive):
\[y_{n,k}\sim \mathrm{LogNormal}(\hat{z}_{n,k}, \sigma_{k})\]
Equivalently:
\[\log y_{n,k} = \log \widehat{z}_{n,k} + \epsilon_{n,k}\]
with
\[\epsilon_{n,k} \sim \mathrm{Normal}(0, \sigma_{k})\]
Creating the model
Variables for known constants and observed data.
data { int<lower = 0> N; // num measurements real ts[N]; // measurement times > 0 real y0[2]; // initial pelts real<lower=0> y[N,2]; // subsequent pelts }
Variables for unknown parameters.
parameters { real<lower=0> theta[4]; // alpha, beta, gamma, delta real<lower=0> z0[2]; // initial population real<lower=0> sigma[2]; // scale of prediction error }
Prior distribution and likelihood.
model { // priors sigma ~ lognormal(0, 0.5); theta[{1, 3}] ~ normal(1, 0.5); theta[{2, 4}] ~ normal(0.05, 0.05); z0[1] ~ lognormal(log(30), 5); z0[2] ~ lognormal(log(5), 5); // likelihood (lognormal) for (k in 1:2) { y0[k] ~ lognormal(log(z0[k]), sigma[k]); y[ , k] ~ lognormal(log(z[, k]), sigma[k]); } }
Solution to the differential equations. We have to define variables for populations predicted by ode
, given:
 System function (
dz_dt
), initial populations (z0
).  Initial time (
0.0
), solution times (ts
).  Parameters (
theta
), data arrays.  Tolerances (
1e6, 1e4
), max iterations (1e3
).
transformed parameters { real z[N, 2] = integrate_ode_rk45( dz_dt, z0, 0.0, ts, theta, rep_array(0.0, 0), rep_array(0, 0), 1e6, 1e4, 1e3); }
LotkaVolterra Parameter estimation
fit = stan(lotkavolterra.stan, data=lynx_hare_data) print(fit, c("theta", "sigma"), probs=c(0.1, 0.5, 0.9))
It is obtained, for example:
mean se_mean sd 10% 50% 90% n_eff Rhat ## theta[1] 0.55 0 0.07 0.46 0.54 0.64 1168 1 ## theta[2] 0.03 0 0.00 0.02 0.03 0.03 1305 1 ## theta[3] 0.80 0 0.10 0.68 0.80 0.94 1117 1 ## theta[4] 0.02 0 0.00 0.02 0.02 0.03 1230 1 ## sigma[1] 0.29 0 0.05 0.23 0.28 0.36 2673 1 ## sigma[2] 0.29 0 0.06 0.23 0.29 0.37 2821 1
Analysis of the obtained results:
 Rhat near 1 signals convergence; n_eff is effective sample size.
 10%, … posterior quantiles; e.g. \(P[\alpha \in (0.46,0.64)y]=0.8\).
 Posterior mean is a Bayesian point estimate: \(\alpha = 0.55\).
 Standard error in posterior mean estimate is 0 (with rounding).
 Posterior standard deviation of \(\alpha\) is 0.07.
Full program to run
file = "https://raw.githubusercontent.com/standev/examplemodels/master/knitr/lotkavolterra/hudsonbaylynxhare.csv" # https://github.com/bblais/SystemsModelingSpring2015Notebooks/blob/master/data/Lynx%20and%20Hare%20Data/lynxhare.csv lynx_hare_df = read.csv(file, header=TRUE, comment.char="#") # lynx_hare_df N = length(lynx_hare_df$Year)  1 ts = 1:N y_init = c(lynx_hare_df$Hare[1], lynx_hare_df$Lynx[1]) y = as.matrix(lynx_hare_df[2:(N + 1), 2:3]) y = cbind(y[ , 2], y[ , 1]); # hare, lynx order lynx_hare_data = list(N, ts, y_init, y) library(rstan) thething = " functions { real[] dz_dt(real t, real[] z, real[] theta, real[] x_r, int[] x_i) { real u = z[1]; real v = z[2]; real alpha = theta[1]; real beta = theta[2]; real gamma = theta[3]; real delta = theta[4]; real du_dt = (alpha  beta * v) * u; real dv_dt = (gamma + delta * u) * v; return { du_dt, dv_dt }; } } data { int<lower = 0> N; // num measurements real ts[N]; // measurement times > 0 real y_init[2]; // initial measured population real<lower = 0> y[N, 2]; // measured population at measurement times } parameters { real<lower = 0> theta[4]; // theta = {alpha, beta, gamma, delta} real<lower = 0> z_init[2]; // initial population real<lower = 0> sigma[2]; // error scale } transformed parameters { real z[N, 2] = integrate_ode_rk45(dz_dt, z_init, 0, ts, theta, rep_array(0.0, 0), rep_array(0, 0), 1e6, 1e5, 1e3); } model { theta[{1, 3}] ~ normal(1, 0.5); theta[{2, 4}] ~ normal(0.05, 0.05); sigma ~ lognormal(1, 1); z_init ~ lognormal(log(10), 1); for (k in 1:2) { y_init[k] ~ lognormal(log(z_init[k]), sigma[k]); y[ , k] ~ lognormal(log(z[, k]), sigma[k]); } } " fit = stan(model_code=thething, data=lynx_hare_data, iter=5000, chains=3, cores=3) print(fit, pars=c("theta", "sigma", "z_init"), probs=c(0.1, 0.5, 0.9), digits = 3)
Other references and examples of Stan
 Andrew Gelman has a blog about
rstan
that shows examples to run in RStudio Cloud. But I had problems to run these codes in the Cloud… 🙁  Example Session in RStudio Cloud.
 All examples of his blog can be downloaded from GitHub.
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.