Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

RStan Getting Started

maverickg edited this page · 65 revisions

 

Latest Version: 2.6.0   (06 Feb 2015)

All install instructions below are for the latest version of RStan.

Introduction

RStan is the R (http://www.r-project.org/) interface for Stan. This how-to includes:

For more information on Stan and its modeling language, see:

Prerequisites

R

R version 2.15.1 or later is required (if you are on a Mac, you may need to reinstall if you upgraded Xcode); the latest version of R is available from

http://www.r-project.org/

Follow the download link, then choose a mirror (we recommend http://cran.rstudio.com/ because it redirects to the closest reliable mirror), then click on the link for your platform (Windows, Linux, or Mac). For Windows, there is an additional step of choosing the "base" package before the download.

The Linux and Mac versions of the R command-line and GUI should install and work with the default configurations.

Toolchain

How to Install RStan

  • Open R (either the R GUI, in the terminal using command R, or by opening the RStudio application).

  • Set the number of processes to use for the build to the number of cores on your machine you want to devote to the build. For example, to use 4 processes, execute the following in R.

Sys.setenv(MAKEFLAGS = "-j4") 
  • Run the following install script in R, which may ask you to select a CRAN mirror.
source('http://mc-stan.org/rstan/install.R', echo = TRUE, max.deparse.length = 2000)
install_rstan()
  • Restart You may need to restart R after the installation.

How to Use RStan

Load rstan

The package names is rstan, so we need to use library(rstan) to load the package.

library(rstan) 

Example 1: Eight Schools

This is an example in Section 5.5 of Gelman et al (2003), which studied coaching effects from eight schools. For simplicity, we call this example "eight schools."

First, we specify this model in a file called 8schools.stan as follows (it can be found here):

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}

In this model, we let theta be transformed parameters of mu and eta instead of directly declaring theta as parameters. By parameterizing this way, the sampler will run more efficiently. Assuming we have 8schools.rstan file in our working directory, we can prepare the data and fit the model as the following R code shows.

schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

fit <- stan(file = '8schools.stan', data = schools_dat, 
            iter = 1000, chains = 4)

In RStan, we can also specify a Stan model using a character string by using argument model_code of function stan instead. The following is a demonstration.

# 
# This is a demo of using model_code argument since 
# we can use this file directly or put the string in R 
# directly. 
#  
schools_code <- paste(readLines('8schools.stan'), collapse = '\n')
fit1 <- stan(model_code = schools_code, data = schools_dat, 
             iter = 1000, chains = 4)

Once a model is fitted, we can use the fitted result as an input to fit the model with other data or settings. This would save us time of compiling the C++ code for the model. By specifying parameter fit for function stan, we can fit the model again. For example, if we want to sample more iterations, we can proceed as follows.

fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4)

The object fit2, returned from function stan is an S4 object of class stanfit. Methods such as print and plot are associated with the fitted result so we can use the following code to check out the results in fit2. print provides a summary for the parameter of the model as well as the log-posterior with name lp__ (see the following example output). For more methods and details of class stanfit, see the help of class stanfit.

In particular, we can use extract function on stanfit objects to obtain the samples. extract extracts samples from the stanfit object as a list of arrays for parameters of interest, or just an array. In addition, S3 functions as.array and as.matrix are defined for stanfit object (using help("as.array.stanfit") to check out the help document in R).

print(fit2)
plot(fit2)

la <- extract(fit2, permuted = TRUE) # return a list of arrays 
mu <- la$mu 

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit2, permuted = FALSE) 

### use S3 functions as.array (or as.matrix) on stanfit objects
a2 <- as.array(fit2)
m <- as.matrix(fit2)
> print(fit, digits = 1)
Inference for Stan model: schools_code.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

         mean se_mean  sd  2.5%  25%  50%  75% 97.5% n_eff Rhat
mu        7.9     0.2 4.9  -2.1  4.5  7.9 11.0  17.8   422    1
tau       6.3     0.3 5.0   0.2  2.5  5.2  8.9  18.7   214    1
eta[1]    0.4     0.0 0.9  -1.5 -0.2  0.4  1.0   2.1   928    1
eta[2]    0.0     0.0 0.9  -1.8 -0.6  0.0  0.5   1.8  1640    1
eta[3]   -0.2     0.0 1.0  -2.1 -0.8 -0.2  0.4   1.8  1243    1
eta[4]    0.0     0.0 0.9  -1.7 -0.6  0.0  0.6   1.7  1421    1
eta[5]   -0.3     0.0 0.9  -2.0 -1.0 -0.4  0.3   1.5   883    1
eta[6]   -0.2     0.0 0.9  -2.0 -0.8 -0.2  0.4   1.6   926    1
eta[7]    0.4     0.0 0.9  -1.4 -0.2  0.4  0.9   2.1   969    1
eta[8]    0.1     0.0 1.0  -1.8 -0.6  0.1  0.7   2.0  1365    1
theta[1] 11.4     0.3 8.1  -1.4  5.9 10.3 15.2  30.6   574    1
theta[2]  7.7     0.2 6.1  -3.7  3.9  7.8 11.4  19.5   762    1
theta[3]  5.8     0.3 7.9 -12.1  1.8  6.5 10.5  19.9   715    1
theta[4]  8.0     0.2 6.5  -5.4  3.9  8.1 12.3  20.2   977    1
theta[5]  5.0     0.3 6.7 -10.3  1.3  5.7  9.5  16.5   667    1
theta[6]  6.0     0.2 6.6  -8.4  2.0  6.2 10.2  18.6   976    1
theta[7] 10.8     0.3 6.8  -1.1  6.2 10.2 14.9  26.0   596    1
theta[8]  8.6     0.3 7.9  -6.1  4.0  8.1 12.6  27.7   629    1
lp__     -5.0     0.1 2.6 -10.7 -6.6 -4.8 -3.1  -0.5   367    1

Samples were drawn using NUTS2 at Fri Apr 12 22:09:54 2013.
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).

In addition, as in BUGS (or JAGS), CmdStan (the command line interface to Stan) needs all the data to be in an R dump file. In the case we have this file, rstan provides function read_rdump to read all the data into an R list. For example, we have a file named 8schools.rdump that contains the following text in our working directory.

J <- 8
y <- c(28,  8, -3,  7, -1,  1, 18, 12)
sigma_y <- c(15, 10, 16, 11,  9, 11, 10, 18)

Then we can read the data from "8schools.rdump" as follows.

schools_dat <- read_rdump('8schools.rdump')

The R dump file actually can be sourced using function source in R into the global environment. In this case, we can specify the data for function stan using object names. That is,

source('8schools.rdump') 
fit <- stan(file = '8schools.stan', data = c("J", "y", "sigma_y"), 
            iter = 1000, chains = 4) 

Example 2: Rats

The Rats example is also a popular example. For example, we can find the OpenBUGS version from here, which originally is from Gelfand et al (1990). The data are about the growth of 30 rats weekly for five weeks. In the following table, we list the data, in which we use x to denote the dates the data were collected. We can try this example using the linked data rats.txt and model code rats.stan.

Rat x=8 x=15 x=22 x=29 x=36 Rat x=8 x=15 x=22 x=29 x=36
1 151 199 246 283 320 16 160 207 248 288 324
2 145 199 249 293 354 17 142 187 234 280 316
3 147 214 263 312 328 18 156 203 243 283 317
4 155 200 237 272 297 19 157 212 259 307 336
5 135 188 230 280 323 20 152 203 246 286 321
6 159 210 252 298 331 21 154 205 253 298 334
7 141 189 231 275 305 22 139 190 225 267 302
8 159 201 248 297 338 23 146 191 229 272 302
9 177 236 285 350 376 24 157 211 250 285 323
10 134 182 220 260 296 25 132 185 237 286 331
11 160 208 261 313 352 26 160 207 257 303 345
12 143 188 220 273 314 27 169 216 261 295 333
13 154 200 244 289 325 28 157 205 248 289 316
14 171 221 270 326 358 29 137 180 219 258 291
15 163 216 242 281 312 30 153 200 244 286 324
library(rstan)  
y <- read.table('rats.txt', header = TRUE)
x <- c(8, 15, 22, 29, 36)
rats_dat <- list(N = nrow(y), T = ncol(y), 
                 x = x, y = y, xbar = mean(x))
rats_fit <- stan(file = 'rats.stan', data = rats_dat, verbose = FALSE)

Sample multiple chains in parallel

To facilitate running multiple chains in parallel, RStan provides function sflist2stanfit to consolidate a list of stanfit objects for the same model to one stanfit object. But it is the responsibility of users to run multiple chains on a local computer or on a cluster (RStan itself does not implement running multiple chains in parallel: in function stan, multiple chains are run sequentially). Here we give a simple example using package parallel for Unix-like operating systems (such as Mac OS X and Linux) and Windows, respectively. Also see the example of function `sflist2stanfit' in rstan.

library(parallel)
library(rstan)
example(sflist2stanfit, run.dontrun = TRUE)
Unix-like OS (use function mclapply in package parallel)

For example, to run four chains for the model in file foo.stan in parallel on two cores, use:

library(parallel)
library(rstan)

foo_data <- ...;
rng_seed <- ...;
foo <- stan("foo.stan", data = foo_data, chains = 0)
sflist <- 
  mclapply(1:4, mc.cores = 2, 
           function(i) stan(fit = foo, data = foo_data, 
                            seed = rng_seed, 
                            chains = 1, chain_id = i, 
                            refresh = -1))

The 1:4 indicates how many chains to run, and the mc.cores indicates how many cores to use; as in this example, these don't have to match. It's important that you use different values for chain ID and that you use chains=1 in each call.

If you already have a fit object (here foo) for foo.stan, you can specify fit=foo as an argument to stan() to avoid the overhead of compiling the model again for each chain. By specifying chains=0 originally, the model compiles but does not draw any samples.

To collect the results back into a single Stan fit object for inference, use:

fit <- sflist2stanfit(sflist)
Windows (use function parLapply in package parallel)

On Windows, to run four chains for the model in file foo.stan in parallel on two cores, follow the demo code below.

library(parallel)
library(rstan)

foo_data <- ...
rng_seed <- ...
foo <- stan(file = 'foo.stan', data = foo_data, chains = 0)

CL = makeCluster(2, outfile = 'parallel.log')
clusterExport(cl = CL, c("foo_data", "foo", "rng_seed")) 
sflist1 <- parLapply(CL, 1:4, fun = function(cid) {  
  require(rstan)
  stan(fit = foo, data = foo_data, chains = 1, 
       iter = 2000, seed = rng_seed, chain_id = cid)
})

fit <- sflist2stanfit(sflist1)
print(fit)
stopCluster(CL)

In the above code, we pass a file to makeCluster so we can look at the output from the child processes. This output file is helpful for debugging.

More Help

More details about RStan can be found in the documentation including the vignette of package rstan. For example, using help(stan) and help("stanfit-class") to check out the help for function stan and S4 class stanfit.
And see Stan's modeling language manual for details about Stan's samplers, optimizers, and the Stan modeling language.

In addition, the Stan User's Mailing list can be used to discuss the use of Stan, post examples or ask questions about (R)Stan. When help is needed, it is important to provide enough information such as the following:

  • model code in Stan modeling language
  • data
  • necessary R code
  • dump of error message using verbose=TRUE for calling function stan
  • version of the C++ compiler, for example, using g++ -v to obtain this if gcc is used
  • information about R by using function sessionInfo in R

References

  • Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2003). Bayesian Data Analysis, CRC Press, London, 2nd Edition.
  • The Stan Development Team (2014). Stan Modeling Language User's Guide and Reference Manual.
  • Gelfand, A. E., Hills S. E., Racine-Poon, A., and Smith A. F. M. (1990). "Illustration of Bayesian Inference in Normal Data Models Using Gibbs Sampling", Journal of the American Statistical Association, 85, 972-985.
  • Stan
  • R
  • BUGS
  • OpenBUGS
  • JAGS
  • Rcpp
Something went wrong with that request. Please try again.