CoCalc Shared Filessupport / 2014-11-24-rstan.sagews
Authors: Harald Schilly, ℏal Snyder, William A. Stein
Description: Examples for support purposes.
# Write file from within Sage worksheet --

open('/tmp/8schools.stan','w').write("""
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);
}""")

%r

library(rstan)

Loading required package: Rcpp Loading required package: inline Attaching package: ‘inline’ The following object is masked from ‘package:Rcpp’: registerPlugin rstan (Version 2.5.0, packaged: 2014-11-24 23:24:49 UTC, GitRev: 498f4db1f270)
%r

# WARNING: This takes a while to evaluate.  It's compiling C++ code,
# running it, etc., behind the scenes.  Be patient.

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 = '/tmp/8schools.stan', data = schools_dat,
iter = 1000, chains = 4)


TRANSLATING MODEL '8schools' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL '8schools' NOW. In file included from /usr/local/sage/sage-6.4/local/include/boost/random/detail/large_arithmetic.hpp:19:0, from /usr/local/sage/sage-6.4/local/include/boost/random/detail/const_mod.hpp:23, from /usr/local/sage/sage-6.4/local/include/boost/random/linear_congruential.hpp:30, from /usr/local/sage/sage-6.4/local/lib/R/library/rstan/include//stansrc/stan/model/model_header.hpp:14, from file2147c9a4659.cpp:8: /usr/local/sage/sage-6.4/local/include/boost/random/detail/integer_log2.hpp:71:35: warning: always_inline function might not be inlinable [-Wattributes] BOOST_RANDOM_DETAIL_CONSTEXPR int integer_log2(T t) ^ SAMPLING FOR MODEL '8schools' NOW (CHAIN 1). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.258777 seconds (Warm-up) # 0.163569 seconds (Sampling) # 0.422346 seconds (Total) SAMPLING FOR MODEL '8schools' NOW (CHAIN 2). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.338461 seconds (Warm-up) # 0.408218 seconds (Sampling) # 0.746679 seconds (Total) SAMPLING FOR MODEL '8schools' NOW (CHAIN 3). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.248708 seconds (Warm-up) # 0.199703 seconds (Sampling) # 0.448411 seconds (Total) SAMPLING FOR MODEL '8schools' NOW (CHAIN 4). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.229368 seconds (Warm-up) # 0.272838 seconds (Sampling) # 0.502206 seconds (Total)
%r

#
# 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('/tmp/8schools.stan'), collapse = '\n')
fit1 <- stan(model_code = schools_code, data = schools_dat,
iter = 1000, chains = 4)


Warning message: In readLines("/tmp/8schools.stan") : incomplete final line found on '/tmp/8schools.stan' TRANSLATING MODEL 'schools_code' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'schools_code' NOW. In file included from /usr/local/sage/sage-6.4/local/include/boost/random/detail/large_arithmetic.hpp:19:0, from /usr/local/sage/sage-6.4/local/include/boost/random/detail/const_mod.hpp:23, from /usr/local/sage/sage-6.4/local/include/boost/random/linear_congruential.hpp:30, from /usr/local/sage/sage-6.4/local/lib/R/library/rstan/include//stansrc/stan/model/model_header.hpp:14, from file214441ca7aa.cpp:8: /usr/local/sage/sage-6.4/local/include/boost/random/detail/integer_log2.hpp:71:35: warning: always_inline function might not be inlinable [-Wattributes] BOOST_RANDOM_DETAIL_CONSTEXPR int integer_log2(T t) ^ SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.240743 seconds (Warm-up) # 0.251067 seconds (Sampling) # 0.49181 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.263597 seconds (Warm-up) # 0.186642 seconds (Sampling) # 0.450239 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.331187 seconds (Warm-up) # 0.272641 seconds (Sampling) # 0.603828 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) Iteration: 300 / 1000 [ 30%] (Warmup) Iteration: 400 / 1000 [ 40%] (Warmup) Iteration: 500 / 1000 [ 50%] (Warmup) Iteration: 501 / 1000 [ 50%] (Sampling) Iteration: 600 / 1000 [ 60%] (Sampling) Iteration: 700 / 1000 [ 70%] (Sampling) Iteration: 800 / 1000 [ 80%] (Sampling) Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.283146 seconds (Warm-up) # 0.290336 seconds (Sampling) # 0.573482 seconds (Total)
%r

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


SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1). Iteration: 1 / 10000 [ 0%] (Warmup) Iteration: 1000 / 10000 [ 10%] (Warmup) Iteration: 2000 / 10000 [ 20%] (Warmup) Iteration: 3000 / 10000 [ 30%] (Warmup) Iteration: 4000 / 10000 [ 40%] (Warmup) Iteration: 5000 / 10000 [ 50%] (Warmup) Iteration: 5001 / 10000 [ 50%] (Sampling) Iteration: 6000 / 10000 [ 60%] (Sampling) Iteration: 7000 / 10000 [ 70%] (Sampling) Iteration: 8000 / 10000 [ 80%] (Sampling) Iteration: 9000 / 10000 [ 90%] (Sampling) Iteration: 10000 / 10000 [100%] (Sampling) # Elapsed Time: 2.13098 seconds (Warm-up) # 2.55002 seconds (Sampling) # 4.681 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2). Iteration: 1 / 10000 [ 0%] (Warmup) Iteration: 1000 / 10000 [ 10%] (Warmup) Iteration: 2000 / 10000 [ 20%] (Warmup) Iteration: 3000 / 10000 [ 30%] (Warmup) Iteration: 4000 / 10000 [ 40%] (Warmup) Iteration: 5000 / 10000 [ 50%] (Warmup) Iteration: 5001 / 10000 [ 50%] (Sampling) Iteration: 6000 / 10000 [ 60%] (Sampling) Iteration: 7000 / 10000 [ 70%] (Sampling) Iteration: 8000 / 10000 [ 80%] (Sampling) Iteration: 9000 / 10000 [ 90%] (Sampling) Iteration: 10000 / 10000 [100%] (Sampling) # Elapsed Time: 1.95314 seconds (Warm-up) # 1.67748 seconds (Sampling) # 3.63063 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3). Iteration: 1 / 10000 [ 0%] (Warmup) Iteration: 1000 / 10000 [ 10%] (Warmup) Iteration: 2000 / 10000 [ 20%] (Warmup) Iteration: 3000 / 10000 [ 30%] (Warmup) Iteration: 4000 / 10000 [ 40%] (Warmup) Iteration: 5000 / 10000 [ 50%] (Warmup) Iteration: 5001 / 10000 [ 50%] (Sampling) Iteration: 6000 / 10000 [ 60%] (Sampling) Iteration: 7000 / 10000 [ 70%] (Sampling) Iteration: 8000 / 10000 [ 80%] (Sampling) Iteration: 9000 / 10000 [ 90%] (Sampling) Iteration: 10000 / 10000 [100%] (Sampling) # Elapsed Time: 1.93508 seconds (Warm-up) # 1.57316 seconds (Sampling) # 3.50823 seconds (Total) SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4). Iteration: 1 / 10000 [ 0%] (Warmup) Iteration: 1000 / 10000 [ 10%] (Warmup) Iteration: 2000 / 10000 [ 20%] (Warmup) Iteration: 3000 / 10000 [ 30%] (Warmup) Iteration: 4000 / 10000 [ 40%] (Warmup) Iteration: 5000 / 10000 [ 50%] (Warmup) Iteration: 5001 / 10000 [ 50%] (Sampling) Iteration: 6000 / 10000 [ 60%] (Sampling) Iteration: 7000 / 10000 [ 70%] (Sampling) Iteration: 8000 / 10000 [ 80%] (Sampling) Iteration: 9000 / 10000 [ 90%] (Sampling) Iteration: 10000 / 10000 [100%] (Sampling) # Elapsed Time: 1.93833 seconds (Warm-up) # 1.65158 seconds (Sampling) # 3.5899 seconds (Total)
%r

print(fit2)
plot(fit2)

Inference for Stan model: schools_code. 4 chains, each with iter=10000; warmup=5000; thin=1; post-warmup draws per chain=5000, total post-warmup draws=20000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mu 7.92 0.08 5.17 -2.38 4.67 7.89 11.19 18.29 4697 1 tau 6.72 0.10 5.61 0.24 2.58 5.45 9.38 20.68 3296 1 eta[1] 0.41 0.01 0.94 -1.50 -0.21 0.43 1.04 2.20 11925 1 eta[2] -0.01 0.01 0.86 -1.70 -0.58 -0.02 0.55 1.70 12148 1 eta[3] -0.19 0.01 0.93 -2.01 -0.81 -0.20 0.42 1.65 12730 1 eta[4] -0.04 0.01 0.88 -1.78 -0.62 -0.05 0.54 1.71 12745 1 eta[5] -0.35 0.01 0.88 -2.04 -0.93 -0.37 0.21 1.48 11157 1 eta[6] -0.21 0.01 0.88 -1.92 -0.79 -0.22 0.36 1.58 11805 1 eta[7] 0.35 0.01 0.89 -1.45 -0.23 0.36 0.94 2.09 11462 1 eta[8] 0.06 0.01 0.94 -1.80 -0.57 0.06 0.69 1.90 13915 1 theta[1] 11.61 0.09 8.49 -2.18 6.04 10.40 15.83 31.93 8104 1 theta[2] 7.79 0.05 6.25 -4.57 3.89 7.74 11.68 20.56 14928 1 theta[3] 6.09 0.08 7.89 -11.51 1.91 6.60 10.86 20.77 10503 1 theta[4] 7.53 0.06 6.66 -6.46 3.63 7.58 11.62 20.73 13667 1 theta[5] 5.05 0.06 6.40 -9.21 1.23 5.53 9.36 16.41 11307 1 theta[6] 6.11 0.06 6.74 -8.66 2.24 6.52 10.37 18.61 11651 1 theta[7] 10.76 0.06 6.94 -1.39 6.09 10.10 14.80 26.25 11604 1 theta[8] 8.44 0.07 7.98 -7.14 3.83 8.17 12.79 25.71 11828 1 lp__ -4.84 0.04 2.62 -10.69 -6.41 -4.59 -2.99 -0.45 4585 1 Samples were drawn using NUTS(diag_e) at Mon Nov 24 23:31:10 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).
%r

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)


%r

print(fit, digits = 1)

Inference for Stan model: 8schools. 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.6 0.2 5.2 -2.4 4.4 7.5 10.7 17.8 512 1 tau 6.3 0.2 5.5 0.2 2.4 4.8 8.8 20.8 565 1 eta[1] 0.4 0.0 0.9 -1.4 -0.2 0.4 1.0 2.1 1264 1 eta[2] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.6 1.9 1265 1 eta[3] -0.2 0.0 1.0 -2.1 -0.9 -0.2 0.5 1.7 1372 1 eta[4] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.6 1.8 1002 1 eta[5] -0.3 0.0 0.9 -2.0 -0.9 -0.4 0.2 1.4 1240 1 eta[6] -0.2 0.0 0.9 -2.0 -0.8 -0.2 0.4 1.6 1337 1 eta[7] 0.4 0.0 0.9 -1.5 -0.2 0.4 1.0 2.1 1181 1 eta[8] 0.0 0.0 0.9 -1.9 -0.6 0.1 0.7 1.8 1355 1 theta[1] 10.9 0.2 7.6 -1.7 5.9 9.9 14.6 28.1 1101 1 theta[2] 7.7 0.1 6.3 -5.5 3.8 7.6 11.5 20.7 1754 1 theta[3] 6.0 0.2 7.5 -11.3 1.8 6.4 10.4 19.9 1388 1 theta[4] 7.3 0.1 6.4 -6.1 3.6 7.4 11.2 20.5 1890 1 theta[5] 4.9 0.1 6.1 -8.6 1.3 5.3 9.1 15.9 1760 1 theta[6] 6.1 0.2 6.8 -8.5 2.3 6.3 10.5 18.9 1187 1 theta[7] 10.6 0.2 7.1 -1.8 6.0 9.9 14.4 26.4 1154 1 theta[8] 8.1 0.2 8.1 -8.4 3.6 7.8 12.4 25.9 1286 1 lp__ -5.1 0.1 2.7 -11.1 -6.7 -4.8 -3.1 -0.5 613 1 Samples were drawn using NUTS(diag_e) at Mon Nov 24 23:29:47 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).