Here we will check that RStan is installed on your machine. First we’ll use a basic model, check that it compiles with stan_model(), then check that it correctly performs MCMC with sampling(). Then we’ll do the same with a slightly more complicated baseball model.

If you get any errors, you should read them carefully!

Quick check

To run a quick check at the command line, navigate to CLASS_MATERIALS/Tutorials/ in this repository, and run

Rscript check_rstan_install.R

If everything is working correctly you should see some sampling for both models. Otherwise you will see some error messages.

Interactive check

Here’s how you can do the same checks, by hand, so you can cut-and-paste into Rstudio, for instance.

First, load the rstan package:

library(rstan)

8 schools

First, run this code to define the variables with the data and the Stan model block:

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))

schools_block <- "
data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}
"

Ok, now run

mod <- stan_model(model_code = schools_block)

and then, if that succeeds,

fit <- sampling(mod, data = schools_dat, chains = 1, iter = 2000)

This should print out some stuff like this:


SAMPLING FOR MODEL '32af4e878d1455a0a170e2073248307e' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.6e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1:
Chain 1:  Elapsed Time: 0.054897 seconds (Warm-up)
Chain 1:                0.047577 seconds (Sampling)
Chain 1:                0.102474 seconds (Total)
Chain 1:
> print(fit)
Inference for Stan model: 32af4e878d1455a0a170e2073248307e.
1 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=1000.

If all is good, you should now be able to do

print(fit)

and see a nice summary of the output, something like this:

Inference for Stan model: 32af4e878d1455a0a170e2073248307e.
1 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=1000.

           mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu         8.07    0.21 4.79  -1.32   5.06   7.91  10.91  17.93   524 1.00
tau        6.30    0.26 5.33   0.18   2.40   5.05   8.64  19.98   418 1.00
eta[1]     0.38    0.03 0.87  -1.37  -0.16   0.41   0.97   2.00   629 1.01
eta[2]     0.00    0.03 0.93  -1.90  -0.57   0.01   0.59   1.81   954 1.00
eta[3]    -0.22    0.03 0.92  -1.91  -0.82  -0.22   0.38   1.53   921 1.00
eta[4]     0.00    0.03 0.97  -1.85  -0.63  -0.03   0.65   1.90   945 1.00
eta[5]    -0.35    0.03 0.84  -1.99  -0.95  -0.36   0.20   1.32   868 1.00
eta[6]    -0.21    0.03 0.91  -1.97  -0.82  -0.24   0.38   1.65   895 1.00
eta[7]     0.38    0.03 0.88  -1.41  -0.22   0.42   0.99   2.03   865 1.00
eta[8]     0.09    0.04 0.92  -1.90  -0.44   0.12   0.68   1.86   603 1.00
theta[1]  11.12    0.30 7.69  -1.06   6.16  10.14  15.28  29.09   643 1.00
theta[2]   8.15    0.21 6.57  -4.24   3.89   7.84  12.19  21.55   967 1.00
theta[3]   6.26    0.28 7.49 -10.09   2.22   6.56  10.82  20.39   723 1.00
theta[4]   7.66    0.19 6.66  -6.75   3.77   7.96  11.58  20.51  1284 1.00
theta[5]   5.35    0.20 6.23  -9.26   1.52   5.96   9.51  16.14   976 1.00
theta[6]   6.36    0.24 7.02  -9.04   2.45   6.84  10.66  19.09   881 1.00
theta[7]  10.77    0.20 6.33   0.35   6.41   9.96  14.64  25.40  1054 1.00
theta[8]   9.12    0.29 7.81  -5.40   4.63   8.67  13.35  26.44   739 1.00
lp__     -39.66    0.14 2.61 -45.40 -41.19 -39.46 -37.85 -35.10   371 1.00

Samples were drawn using NUTS(diag_e) at Fri Nov 20 10:09:24 2020.
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).

Baseball

This is a somewhat more complex model. First, the data and Stan block (note: you may have to change the location of BattingAverage.csv below):

batting <- read.csv("../Lectures/data/BattingAverage.csv", header=TRUE, stringsAsFactors=TRUE)
batting_dat <- list(N=nrow(batting),
                    hits=batting$Hits,
                    at_bats=batting$AtBats,
                    npos=nlevels(batting$PriPos),
                    position=as.numeric(batting$PriPos))
batting_block <- "
data {
    int N;   // number of players
    int hits[N];
    int at_bats[N];
    int npos; // number of positions
    int position[N];
}
parameters {
    real<lower=0, upper=1> theta[N];
    real<lower=0, upper=1> mu[npos];
    real<lower=0> kappa[npos];
}
model {
    real alpha;
    real beta;
    hits ~ binomial(at_bats, theta);
    for (i in 1:N) {
        alpha = mu[position[i]] * kappa[position[i]];
        beta = (1 - mu[position[i]]) * kappa[position[i]];
        theta[i] ~ beta(alpha, beta);
    }
    mu ~ beta(1,1);
    kappa ~ gamma(0.1,0.1);
}
"

Now, let’s see if compilation and sampling work:

mod <- stan_model(model_code=batting_block)
fit <- sampling(mod, data = batting_dat, chains = 1, iter = 2000)

Again, if all is good, you should be able to do

print(fit, pars=c("mu", "kappa"))

and see something like this:

Inference for Stan model: 1e1b401f92e65870a3dbf7e947638a0e.
1 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

           mean se_mean    sd  2.5%    25%    50%    75%  97.5% n_eff Rhat
mu[1]      0.25    0.00  0.01  0.24   0.25   0.25   0.26   0.26   698 1.00
mu[2]      0.25    0.00  0.01  0.24   0.24   0.25   0.25   0.26   622 1.00
mu[3]      0.26    0.00  0.01  0.24   0.25   0.26   0.26   0.27   806 1.00
mu[4]      0.24    0.00  0.01  0.23   0.24   0.24   0.24   0.25   512 1.00
mu[5]      0.26    0.00  0.01  0.24   0.25   0.26   0.26   0.27   660 1.00
mu[6]      0.25    0.00  0.01  0.24   0.25   0.25   0.25   0.26   447 1.00
mu[7]      0.13    0.00  0.01  0.12   0.12   0.13   0.13   0.14   186 1.00
mu[8]      0.26    0.00  0.01  0.25   0.25   0.26   0.26   0.27   790 1.00
mu[9]      0.25    0.00  0.01  0.23   0.24   0.25   0.25   0.26   654 1.00
kappa[1] 111.53    0.90 23.12 73.68  95.61 108.39 124.96 166.97   666 1.00
kappa[2] 103.32    1.14 23.82 59.75  86.91 101.28 118.25 154.18   437 1.00
kappa[3]  97.62    0.84 21.14 63.39  81.69  95.57 111.66 143.43   628 1.00
kappa[4]  96.56    0.92 20.34 61.82  82.45  93.79 107.81 141.29   489 1.01
kappa[5] 101.10    1.15 22.38 63.73  84.45  99.67 115.62 151.91   377 1.00
kappa[6] 106.94    1.18 21.90 70.68  91.04 105.07 119.21 153.82   342 1.00
kappa[7]  54.44    1.40 13.14 34.40  45.29  52.31  61.42  85.32    88 1.00
kappa[8] 118.92    1.28 26.48 74.16 100.05 117.58 136.22 173.49   426 1.00
kappa[9]  95.74    0.86 22.40 61.19  79.52  92.44 109.71 145.52   680 1.01

Samples were drawn using NUTS(diag_e) at Fri Nov 20 10:12:37 2020.
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).