︠52eea0a0-ebad-4eed-a3fd-6d6c520f1b74s︠ # Write file from within Sage worksheet -- open('/tmp/8schools.stan','w').write(""" data { int J; // number of schools real y[J]; // estimated treatment effects real sigma[J]; // s.e. of effect estimates } parameters { real mu; real 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); }""") ︡caf2994b-9fe2-47d5-b84e-494c5a0a5f37︡ ︠28b7b338-8200-477b-9682-0386fc293342s︠ %r library(rstan) ︡c94dd83d-7d2a-4048-a9cb-4561c995fbd5︡{"stdout":"Loading required package: Rcpp\nLoading required package: inline\n\nAttaching package: ‘inline’\n\nThe following object is masked from ‘package:Rcpp’:\n\n registerPlugin\n\nrstan (Version 2.5.0, packaged: 2014-11-24 23:24:49 UTC, GitRev: 498f4db1f270)"}︡ ︠78fe5624-f03f-466c-9681-364c77f10b2c︠ %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) ︡02985829-80ee-436b-9891-85350e0796d6︡{"stdout":"TRANSLATING MODEL '8schools' FROM Stan CODE TO C++ CODE NOW.\nCOMPILING THE C++ CODE FOR MODEL '8schools' NOW.\nIn file included from /usr/local/sage/sage-6.4/local/include/boost/random/detail/large_arithmetic.hpp:19:0,\n from /usr/local/sage/sage-6.4/local/include/boost/random/detail/const_mod.hpp:23,\n from /usr/local/sage/sage-6.4/local/include/boost/random/linear_congruential.hpp:30,\n from /usr/local/sage/sage-6.4/local/lib/R/library/rstan/include//stansrc/stan/model/model_header.hpp:14,\n from file2147c9a4659.cpp:8:\n/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]\n BOOST_RANDOM_DETAIL_CONSTEXPR int integer_log2(T t)\n ^\n\nSAMPLING FOR MODEL '8schools' NOW (CHAIN 1).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.258777 seconds (Warm-up)\n# 0.163569 seconds (Sampling)\n# 0.422346 seconds (Total)\n\n\nSAMPLING FOR MODEL '8schools' NOW (CHAIN 2).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.338461 seconds (Warm-up)\n# 0.408218 seconds (Sampling)\n# 0.746679 seconds (Total)\n\n\nSAMPLING FOR MODEL '8schools' NOW (CHAIN 3).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.248708 seconds (Warm-up)\n# 0.199703 seconds (Sampling)\n# 0.448411 seconds (Total)\n\n\nSAMPLING FOR MODEL '8schools' NOW (CHAIN 4).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.229368 seconds (Warm-up)\n# 0.272838 seconds (Sampling)\n# 0.502206 seconds (Total)"}︡ ︠be564446-50fe-416f-868c-263e4def489bs︠ %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) ︡3a24cafb-3364-4e45-a057-696e10d80dd4︡{"stdout":"Warning message:\nIn readLines(\"/tmp/8schools.stan\") :\n incomplete final line found on '/tmp/8schools.stan'\n\n\nTRANSLATING MODEL 'schools_code' FROM Stan CODE TO C++ CODE NOW.\nCOMPILING THE C++ CODE FOR MODEL 'schools_code' NOW.\nIn file included from /usr/local/sage/sage-6.4/local/include/boost/random/detail/large_arithmetic.hpp:19:0,\n from /usr/local/sage/sage-6.4/local/include/boost/random/detail/const_mod.hpp:23,\n from /usr/local/sage/sage-6.4/local/include/boost/random/linear_congruential.hpp:30,\n from /usr/local/sage/sage-6.4/local/lib/R/library/rstan/include//stansrc/stan/model/model_header.hpp:14,\n from file214441ca7aa.cpp:8:\n/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]\n BOOST_RANDOM_DETAIL_CONSTEXPR int integer_log2(T t)\n ^\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.240743 seconds (Warm-up)\n# 0.251067 seconds (Sampling)\n# 0.49181 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.263597 seconds (Warm-up)\n# 0.186642 seconds (Sampling)\n# 0.450239 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.331187 seconds (Warm-up)\n# 0.272641 seconds (Sampling)\n# 0.603828 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4).\n\nIteration: 1 / 1000 [ 0%] (Warmup)\nIteration: 100 / 1000 [ 10%] (Warmup)\nIteration: 200 / 1000 [ 20%] (Warmup)\nIteration: 300 / 1000 [ 30%] (Warmup)\nIteration: 400 / 1000 [ 40%] (Warmup)\nIteration: 500 / 1000 [ 50%] (Warmup)\nIteration: 501 / 1000 [ 50%] (Sampling)\nIteration: 600 / 1000 [ 60%] (Sampling)\nIteration: 700 / 1000 [ 70%] (Sampling)\nIteration: 800 / 1000 [ 80%] (Sampling)\nIteration: 900 / 1000 [ 90%] (Sampling)\nIteration: 1000 / 1000 [100%] (Sampling)\n# Elapsed Time: 0.283146 seconds (Warm-up)\n# 0.290336 seconds (Sampling)\n# 0.573482 seconds (Total)"}︡ ︠a87fe052-7320-4973-b0ff-b13543edead9s︠ %r fit2 <- stan(fit = fit1, data = schools_dat, iter = 10000, chains = 4) ︡36727ca6-32ab-40da-b030-92d7be7ae36a︡{"stdout":"SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).\n\nIteration: 1 / 10000 [ 0%] (Warmup)\nIteration: 1000 / 10000 [ 10%] (Warmup)\nIteration: 2000 / 10000 [ 20%] (Warmup)\nIteration: 3000 / 10000 [ 30%] (Warmup)\nIteration: 4000 / 10000 [ 40%] (Warmup)\nIteration: 5000 / 10000 [ 50%] (Warmup)\nIteration: 5001 / 10000 [ 50%] (Sampling)\nIteration: 6000 / 10000 [ 60%] (Sampling)\nIteration: 7000 / 10000 [ 70%] (Sampling)\nIteration: 8000 / 10000 [ 80%] (Sampling)\nIteration: 9000 / 10000 [ 90%] (Sampling)\nIteration: 10000 / 10000 [100%] (Sampling)\n# Elapsed Time: 2.13098 seconds (Warm-up)\n# 2.55002 seconds (Sampling)\n# 4.681 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 2).\n\nIteration: 1 / 10000 [ 0%] (Warmup)\nIteration: 1000 / 10000 [ 10%] (Warmup)\nIteration: 2000 / 10000 [ 20%] (Warmup)\nIteration: 3000 / 10000 [ 30%] (Warmup)\nIteration: 4000 / 10000 [ 40%] (Warmup)\nIteration: 5000 / 10000 [ 50%] (Warmup)\nIteration: 5001 / 10000 [ 50%] (Sampling)\nIteration: 6000 / 10000 [ 60%] (Sampling)\nIteration: 7000 / 10000 [ 70%] (Sampling)\nIteration: 8000 / 10000 [ 80%] (Sampling)\nIteration: 9000 / 10000 [ 90%] (Sampling)\nIteration: 10000 / 10000 [100%] (Sampling)\n# Elapsed Time: 1.95314 seconds (Warm-up)\n# 1.67748 seconds (Sampling)\n# 3.63063 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 3).\n\nIteration: 1 / 10000 [ 0%] (Warmup)\nIteration: 1000 / 10000 [ 10%] (Warmup)\nIteration: 2000 / 10000 [ 20%] (Warmup)\nIteration: 3000 / 10000 [ 30%] (Warmup)\nIteration: 4000 / 10000 [ 40%] (Warmup)\nIteration: 5000 / 10000 [ 50%] (Warmup)\nIteration: 5001 / 10000 [ 50%] (Sampling)\nIteration: 6000 / 10000 [ 60%] (Sampling)\nIteration: 7000 / 10000 [ 70%] (Sampling)\nIteration: 8000 / 10000 [ 80%] (Sampling)\nIteration: 9000 / 10000 [ 90%] (Sampling)\nIteration: 10000 / 10000 [100%] (Sampling)\n# Elapsed Time: 1.93508 seconds (Warm-up)\n# 1.57316 seconds (Sampling)\n# 3.50823 seconds (Total)\n\n\nSAMPLING FOR MODEL 'schools_code' NOW (CHAIN 4).\n\nIteration: 1 / 10000 [ 0%] (Warmup)\nIteration: 1000 / 10000 [ 10%] (Warmup)\nIteration: 2000 / 10000 [ 20%] (Warmup)\nIteration: 3000 / 10000 [ 30%] (Warmup)\nIteration: 4000 / 10000 [ 40%] (Warmup)\nIteration: 5000 / 10000 [ 50%] (Warmup)\nIteration: 5001 / 10000 [ 50%] (Sampling)\nIteration: 6000 / 10000 [ 60%] (Sampling)\nIteration: 7000 / 10000 [ 70%] (Sampling)\nIteration: 8000 / 10000 [ 80%] (Sampling)\nIteration: 9000 / 10000 [ 90%] (Sampling)\nIteration: 10000 / 10000 [100%] (Sampling)\n# Elapsed Time: 1.93833 seconds (Warm-up)\n# 1.65158 seconds (Sampling)\n# 3.5899 seconds (Total)"}︡ ︠c11c0f83-e788-4ec2-9903-54be5668c9c9s︠ %r print(fit2) plot(fit2) ︡8e940f61-7544-4bb8-ae8d-c11df8ed672c︡{"stdout":"\n"}︡{"once":false,"file":{"show":true,"uuid":"d3d92873-06ce-46ed-9900-2c8e0fd5b966","filename":"/tmp/96cee31e-9c2e-4569-bd88-a18bcc20517b.svg"}}︡{"stdout":"\n"}︡{"stdout":"Inference for Stan model: schools_code.\n4 chains, each with iter=10000; warmup=5000; thin=1; \npost-warmup draws per chain=5000, total post-warmup draws=20000.\n\n mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\nmu 7.92 0.08 5.17 -2.38 4.67 7.89 11.19 18.29 4697 1\ntau 6.72 0.10 5.61 0.24 2.58 5.45 9.38 20.68 3296 1\neta[1] 0.41 0.01 0.94 -1.50 -0.21 0.43 1.04 2.20 11925 1\neta[2] -0.01 0.01 0.86 -1.70 -0.58 -0.02 0.55 1.70 12148 1\neta[3] -0.19 0.01 0.93 -2.01 -0.81 -0.20 0.42 1.65 12730 1\neta[4] -0.04 0.01 0.88 -1.78 -0.62 -0.05 0.54 1.71 12745 1\neta[5] -0.35 0.01 0.88 -2.04 -0.93 -0.37 0.21 1.48 11157 1\neta[6] -0.21 0.01 0.88 -1.92 -0.79 -0.22 0.36 1.58 11805 1\neta[7] 0.35 0.01 0.89 -1.45 -0.23 0.36 0.94 2.09 11462 1\neta[8] 0.06 0.01 0.94 -1.80 -0.57 0.06 0.69 1.90 13915 1\ntheta[1] 11.61 0.09 8.49 -2.18 6.04 10.40 15.83 31.93 8104 1\ntheta[2] 7.79 0.05 6.25 -4.57 3.89 7.74 11.68 20.56 14928 1\ntheta[3] 6.09 0.08 7.89 -11.51 1.91 6.60 10.86 20.77 10503 1\ntheta[4] 7.53 0.06 6.66 -6.46 3.63 7.58 11.62 20.73 13667 1\ntheta[5] 5.05 0.06 6.40 -9.21 1.23 5.53 9.36 16.41 11307 1\ntheta[6] 6.11 0.06 6.74 -8.66 2.24 6.52 10.37 18.61 11651 1\ntheta[7] 10.76 0.06 6.94 -1.39 6.09 10.10 14.80 26.25 11604 1\ntheta[8] 8.44 0.07 7.98 -7.14 3.83 8.17 12.79 25.71 11828 1\nlp__ -4.84 0.04 2.62 -10.69 -6.41 -4.59 -2.99 -0.45 4585 1\n\nSamples were drawn using NUTS(diag_e) at Mon Nov 24 23:31:10 2014.\nFor each parameter, n_eff is a crude measure of effective sample size,\nand Rhat is the potential scale reduction factor on split chains (at \nconvergence, Rhat=1)."}︡ ︠f3dd5858-d919-4d0f-9307-ab698948ac05s︠ %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) ︡73e60e7d-9928-4f23-b180-e6488d6dd0c9︡ ︠0cf10590-eea0-4bd3-a8e8-755aa88a1473s︠ %r print(fit, digits = 1) ︡a9765069-ee3f-488b-b485-00ee74b109b2︡{"stdout":"Inference for Stan model: 8schools.\n4 chains, each with iter=1000; warmup=500; thin=1; \npost-warmup draws per chain=500, total post-warmup draws=2000.\n\n mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\nmu 7.6 0.2 5.2 -2.4 4.4 7.5 10.7 17.8 512 1\ntau 6.3 0.2 5.5 0.2 2.4 4.8 8.8 20.8 565 1\neta[1] 0.4 0.0 0.9 -1.4 -0.2 0.4 1.0 2.1 1264 1\neta[2] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.6 1.9 1265 1\neta[3] -0.2 0.0 1.0 -2.1 -0.9 -0.2 0.5 1.7 1372 1\neta[4] 0.0 0.0 0.9 -1.8 -0.6 0.0 0.6 1.8 1002 1\neta[5] -0.3 0.0 0.9 -2.0 -0.9 -0.4 0.2 1.4 1240 1\neta[6] -0.2 0.0 0.9 -2.0 -0.8 -0.2 0.4 1.6 1337 1\neta[7] 0.4 0.0 0.9 -1.5 -0.2 0.4 1.0 2.1 1181 1\neta[8] 0.0 0.0 0.9 -1.9 -0.6 0.1 0.7 1.8 1355 1\ntheta[1] 10.9 0.2 7.6 -1.7 5.9 9.9 14.6 28.1 1101 1\ntheta[2] 7.7 0.1 6.3 -5.5 3.8 7.6 11.5 20.7 1754 1\ntheta[3] 6.0 0.2 7.5 -11.3 1.8 6.4 10.4 19.9 1388 1\ntheta[4] 7.3 0.1 6.4 -6.1 3.6 7.4 11.2 20.5 1890 1\ntheta[5] 4.9 0.1 6.1 -8.6 1.3 5.3 9.1 15.9 1760 1\ntheta[6] 6.1 0.2 6.8 -8.5 2.3 6.3 10.5 18.9 1187 1\ntheta[7] 10.6 0.2 7.1 -1.8 6.0 9.9 14.4 26.4 1154 1\ntheta[8] 8.1 0.2 8.1 -8.4 3.6 7.8 12.4 25.9 1286 1\nlp__ -5.1 0.1 2.7 -11.1 -6.7 -4.8 -3.1 -0.5 613 1\n\nSamples were drawn using NUTS(diag_e) at Mon Nov 24 23:29:47 2014.\nFor each parameter, n_eff is a crude measure of effective sample size,\nand Rhat is the potential scale reduction factor on split chains (at \nconvergence, Rhat=1)."}︡ ︠9f3b0f28-2762-49bb-af49-377d32e00a89︠