This document details a few ways to use Rstan in knitr documents.
Before getting started, let’s load the rstan and knitr libraries
library("rstan")
library("knitr")To run a model with Rstan we could write the Stan model in a character vector and pass it to stan or stan_model using the option model_code. For example,
scode <- "
parameters {
real y[2];
}
model {
y[1] ~ normal(0, 1);
y[2] ~ double_exponential(0, 2);
}
"
fit1 <- stan(model_code = scode, iter = 10, verbose = FALSE) ##
## TRANSLATING MODEL 'scode' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'scode' NOW.
##
## SAMPLING FOR MODEL 'scode' NOW (CHAIN 1).
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
## # Elapsed Time: 0.000155 seconds (Warm-up)
## # 0.000116 seconds (Sampling)
## # 0.000271 seconds (Total)
##
##
## SAMPLING FOR MODEL 'scode' NOW (CHAIN 2).
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
## # Elapsed Time: 8.4e-05 seconds (Warm-up)
## # 7.5e-05 seconds (Sampling)
## # 0.000159 seconds (Total)
##
##
## SAMPLING FOR MODEL 'scode' NOW (CHAIN 3).
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
## # Elapsed Time: 0.000123 seconds (Warm-up)
## # 0.000128 seconds (Sampling)
## # 0.000251 seconds (Total)
##
##
## SAMPLING FOR MODEL 'scode' NOW (CHAIN 4).
##
## Iteration: 1 / 10 [ 10%] (Warmup)
## Iteration: 2 / 10 [ 20%] (Warmup)
## Iteration: 3 / 10 [ 30%] (Warmup)
## Iteration: 4 / 10 [ 40%] (Warmup)
## Iteration: 5 / 10 [ 50%] (Warmup)
## Iteration: 6 / 10 [ 60%] (Sampling)
## Iteration: 7 / 10 [ 70%] (Sampling)
## Iteration: 8 / 10 [ 80%] (Sampling)
## Iteration: 9 / 10 [ 90%] (Sampling)
## Iteration: 10 / 10 [100%] (Sampling)
## # Elapsed Time: 9.3e-05 seconds (Warm-up)
## # 9.1e-05 seconds (Sampling)
## # 0.000184 seconds (Total)print(fit1)## Inference for Stan model: scode.
## 4 chains, each with iter=10; warmup=5; thin=1;
## post-warmup draws per chain=5, total post-warmup draws=20.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## y[1] -0.07 0.33 1.21 -1.74 -1.07 -0.23 0.62 2.37 13 1.92
## y[2] -0.06 1.17 2.42 -3.37 -2.35 -0.02 2.00 3.49 4 2.51
## lp__ -1.70 0.26 1.17 -4.41 -2.06 -1.53 -0.81 -0.38 20 1.61
##
## Samples were drawn using NUTS(diag_e) at Thu Nov 27 18:55:37 2014.
## 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).This document provides two more elegant altenatives to that approach.
cat engine to write the code within a block and save it to a file to use later.stan engine which will compile Stan code in a block to a stanmodel object which can be used later.Knitr has an engine cat which saves the contents of the chunk to a file. So write the engine='stan', and engine.opts=list(file="nameoffile.stan") to indicate the name of the file which to write the stancode. In this example, we’ll write the Stan model in the following chunk to the file “ex1.stan”.
The Stan code is cleanly printed, and we can use the file “ex1.stan”:
fit2 <- stan("ex1.stan") ##
## TRANSLATING MODEL 'ex1' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'ex1' NOW.
##
## SAMPLING FOR MODEL 'ex1' NOW (CHAIN 1).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.008975 seconds (Warm-up)
## # 0.010661 seconds (Sampling)
## # 0.019636 seconds (Total)
##
##
## SAMPLING FOR MODEL 'ex1' NOW (CHAIN 2).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.0091 seconds (Warm-up)
## # 0.009007 seconds (Sampling)
## # 0.018107 seconds (Total)
##
##
## SAMPLING FOR MODEL 'ex1' NOW (CHAIN 3).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009391 seconds (Warm-up)
## # 0.008061 seconds (Sampling)
## # 0.017452 seconds (Total)
##
##
## SAMPLING FOR MODEL 'ex1' NOW (CHAIN 4).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009009 seconds (Warm-up)
## # 0.006802 seconds (Sampling)
## # 0.015811 seconds (Total)print(fit2)## Inference for Stan model: ex1.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## y[1] -0.01 0.03 1.00 -1.99 -0.69 -0.01 0.66 1.88 1528 1
## y[2] 0.08 0.08 2.75 -5.73 -1.21 0.01 1.31 6.45 1176 1
## lp__ -1.46 0.04 1.20 -4.57 -2.01 -1.16 -0.56 -0.10 953 1
##
## Samples were drawn using NUTS(diag_e) at Thu Nov 27 18:56:12 2014.
## 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).As of now, knitr doesn’t have an engine for Stan but we can define our own. The following engine works similarly to the Rcpp engine; it will compile the code in the chunk.
eng_stan = function(options) {
code = paste(options$code, collapse = '\n')
opts = options$engine.opts
if (!is.environment(opts$env)) env = knit_global()
else env = opts$env
opts$env = NULL
#' name of the modelfit object returned by stan_model
if (is.null(opts$x))
if (!is.null(opts$model_name)) x = model_name
else x = formals(getFromNamespace('stan_model', 'rstan'))$model_name
else x = as.character(opts$x)
opts$x = NULL
if (options$eval) {
message(sprintf('Creating rstan::stanmodel object \'%s\' from Stan code', x))
assign(x,
do.call(getFromNamespace('stan_model', 'rstan'),
c(list(model_code = code), opts)),
envir = env)
}
engine_output(options, code, '')
}
knit_engines$set(stan = eng_stan)To use the engine, set the following options in the chunk:
engine='stan'engine.opts A list with options to pass to the engine. You need to specify x as the name of the stanmodel variable. Any other options are passed to stan_model.The following chunk uses the stan engine and assigns the result of stan_model to ex1.
parameters {
real y[2];
}
model {
y[1] ~ normal(0, 1);
y[2] ~ double_exponential(0, 2);
} Now there is an object of class stanmodel named ex1,
print(ex1)## S4 class stanmodel 'anon_model' coded as follows:
## parameters {
## real y[2];
## }
## model {
## y[1] ~ normal(0, 1);
## y[2] ~ double_exponential(0, 2);
## }We can sample from ex1 using the sampling method:
fit3 <- sampling(ex1) ##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009483 seconds (Warm-up)
## # 0.008584 seconds (Sampling)
## # 0.018067 seconds (Total)
##
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009375 seconds (Warm-up)
## # 0.008633 seconds (Sampling)
## # 0.018008 seconds (Total)
##
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009049 seconds (Warm-up)
## # 0.008701 seconds (Sampling)
## # 0.01775 seconds (Total)
##
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.009122 seconds (Warm-up)
## # 0.008229 seconds (Sampling)
## # 0.017351 seconds (Total)print(fit3)## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## y[1] -0.03 0.03 1.02 -2.01 -0.75 -0.06 0.67 1.93 1517 1.00
## y[2] 0.01 0.08 2.80 -6.13 -1.33 -0.02 1.29 5.85 1269 1.00
## lp__ -1.50 0.04 1.25 -4.70 -2.07 -1.17 -0.60 -0.10 1077 1.01
##
## Samples were drawn using NUTS(diag_e) at Thu Nov 27 18:56:45 2014.
## 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).