Research Methods in Political Science II
Day 2-2: Introduction to Stan
Yuki Yanai
November 5, 2016
Preparation
Load some packages.
library('coefplot')
library('rethinking')
library('ggmcmc')Fitting a Normal Linear Model in Different Ways
Data generation
First, let’s create a fake data set we are analyzing today.
set.seed(2016-11-05)
n <- 50
true_beta <- c(2, 0.2, 1.5, -2)
true_sigma <- 3
x2 <- runif(n, min = 1, max = 10)
x3 <- rnorm(n, mean = 0, sd = 4)
x4 <- runif(n, min = -10, max = 10)
y <- true_beta[1] + true_beta[2] * x2 + true_beta[3] * x3 +
true_beta[4] * x4 + rnorm(n, mean = 0, sd = true_sigma)
dd <- data.frame(y = y, x2 = x2, x3 = x3, x4 = x4)Statistical Model
We have the outcome \(y\) and three explanatory variables \(x_2\), \(x_3\), and \(x_4\). Suppose we believe for theoretical reasons that the outcome \(y\) is generated by the following mechanism. \[y_i \sim \mbox{Normal}(\mu_i, \sigma^2),\] \[\mu_i = \beta_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4,\] \[\beta_p \sim \mbox{Normal}(0, 100^2),\] \[\sigma \sim \mbox{Uniform}(0, 100).\]
Fit the model by lm()
Since the prior are non-informative in the model, we can use lm() to estimate \(\beta\)s.
fit_1 <- lm(y ~ x2 + x3 + x4, data = dd)
summary(fit_1)##
## Call:
## lm(formula = y ~ x2 + x3 + x4, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.877 -1.934 -0.183 2.434 5.976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.11688 1.05534 2.953 0.00494
## x2 0.07788 0.16540 0.471 0.63995
## x3 1.50223 0.13465 11.156 1.12e-14
## x4 -2.05099 0.08643 -23.729 < 2e-16
##
## Residual standard error: 3.24 on 46 degrees of freedom
## Multiple R-squared: 0.9304, Adjusted R-squared: 0.9259
## F-statistic: 204.9 on 3 and 46 DF, p-value: < 2.2e-16Let’s visualize the result with coefplot().
coefplot(fit_1, intercept = TRUE)
Fit the model by map()
Alternatively, we can use rethinking::map(), which has been our favorite method so far.
map_model <- alist(
y ~ dnorm(mu, sigma),
mu <- b1 + b2*x2 + b3*x3 + b4*x4,
b1 ~ dnorm(0, 100),
b2 ~ dnorm(0, 100),
b3 ~ dnorm(0, 100),
b4 ~ dnorm(0, 100),
sigma ~ dunif(0, 100)
)
fit_2 <- map(map_model, data = dd)
precis(fit_2)## Mean StdDev 5.5% 94.5%
## b1 3.12 1.01 1.50 4.73
## b2 0.08 0.16 -0.18 0.33
## b3 1.50 0.13 1.30 1.71
## b4 -2.05 0.08 -2.18 -1.92
## sigma 3.11 0.31 2.61 3.60#plot(precis(fit_2))Thesre estimates were obtained by quadratic approximation.
Fit the model by Stan
Indirectly call RStan with map2sta()
Now, let’s fit the model with RStan. We can feed the same model we used for map().
Here, we run 4 chains. We sample 2,000 values in each chain, and discard 1,000 values in each as warm-up. As a result, we will get 4,000 smaple values. Since each chain samples values independently from each other, we can simulate chains in parallel if the computer has multiple cores. Below, we specify the number of cores we use in parallel by core argument.
fit_3 <- map2stan(map_model, data = dd, cores = parallel::detectCores(),
chains = 4, iter = 2000, warmup = 1000)## In file included from file21c61c89ec02.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:42:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file21c61c89ec02.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:43:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
## static void set_zero_all_adjoints_nested() {
## ^
## In file included from file21c61c89ec02.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:54:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## In file included from file21c61c89ec02.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:235:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint.hpp:61:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array.hpp:21:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/base.hpp:28:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## 7 warnings generated.
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 5e-05 seconds (Sampling)
## 5.5e-05 seconds (Total)
##
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]precis(fit_3)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## b1 3.11 1.05 1.43 4.78 2100 1
## b2 0.08 0.17 -0.21 0.32 2066 1
## b3 1.50 0.14 1.28 1.71 2842 1
## b4 -2.05 0.09 -2.19 -1.92 2560 1
## sigma 3.33 0.34 2.77 3.84 2483 1Compare this result to the result of fit_2. In fit_2, the fourth and fifth columns show the values for 5.5 and 94.5 percentiles, respectively, In contrast, this result shows the lower and upper bound of the 89% HPDI (highest posterior density interval). In this example, they don’t differ much. But if the posterior is far from normal, quadratic approximation doesn’t work well, and the percentiles can differ a lot from the HPDI.
Use RStan directly
We have learned map2stan() is very convenient. For a simple model like this, we can instantly get the estimates by map2stan(). However, it is better to know how to use RStan itself without rethinking packages because it gives us more felxibility.
When use RStan, we should specify a model in a seprate file with the extention .stan. Every Stan model has to have three blocks: data, parameters, and model. In addition to them, it can have transformed parameters and generated quantities blocks.
Let’s translate our model into Stan and save it as the file named d2-model1.stan in the stan-models/ folder. RStudio helps you write a Stan model by highliting functions.
data {
int<lower=0> N;
real Y[N];
real X2[N];
real X3[N];
real X4[N];
}
parameters {
real beta[4];
real<lower=0> sigma;
}
transformed parameters {
real mu[N];
for (n in 1:N) {
mu[n] = beta[1] + beta[2]*X2[n] + beta[3]*X3[n] + beta[4]*X4[n];
}
}
model {
for (n in 1:N) {
Y[n] ~ normal(mu[n], sigma);
}
for (p in 1:4) {
beta[p] ~ normal(0, 100);
}
sigma ~ uniform(0, 100);
}Could you tell the differece between the Stan model and the map model?
Now, let’s fit the model. We use rstan::stan() function. But before that, if the computer has a multi-core processor, it is faster to run the simulation in parallel. To do so, set
options(mc.cores = parallel::detectCores())And we need to pass a list of data.
dd1 <- list(N = nrow(dd), Y = dd$y, X2 = dd$x2, X3 = dd$x3, X4 = dd$x4)Now, we’re rady.
fit_4 <- stan(file = 'stan-models/d2-model1.stan', data = dd1,
chains = 4, iter = 2000, warmup = 1000)## In file included from file21c6b232743.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:42:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file21c6b232743.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:43:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
## static void set_zero_all_adjoints_nested() {
## ^
## In file included from file21c6b232743.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:54:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## In file included from file21c6b232743.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:235:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint.hpp:61:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array.hpp:21:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/base.hpp:28:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## 7 warnings generated.As you probably noticed, it takes time to compile to a Stan model. We sometimes wanna re-fit the model without compiling the same model again. For that purpose, we can separate the compiling phrse from the sampling phase. First, comile the model by `rsta::stan_model().
stan_m1 <- stan_model('stan-models/d2-model1.stan')Then, sample the posterior by rstan::sampling().
fit_4 <- sampling(stan_m1, data = dd1, seed = 12345,
chains = 4, iter = 2000, warmup = 1000)Let’s look at the result.
print(fit_4)## Inference for Stan model: d2-model1.
## 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
## beta[1] 3.14 0.02 1.04 1.11 2.43 3.12 3.84 5.19 1786 1
## beta[2] 0.07 0.00 0.16 -0.24 -0.03 0.08 0.18 0.39 1816 1
## beta[3] 1.50 0.00 0.13 1.24 1.41 1.50 1.59 1.76 2995 1
## beta[4] -2.05 0.00 0.09 -2.23 -2.11 -2.05 -1.99 -1.88 2875 1
## sigma 3.32 0.01 0.34 2.72 3.08 3.29 3.53 4.10 3070 1
## mu[1] -7.33 0.02 0.81 -8.90 -7.88 -7.34 -6.81 -5.69 2387 1
## mu[2] -0.95 0.01 0.87 -2.62 -1.55 -0.93 -0.36 0.76 4000 1
## mu[3] 2.66 0.02 0.98 0.72 2.01 2.68 3.31 4.61 3484 1
## mu[4] -4.22 0.01 0.62 -5.43 -4.64 -4.21 -3.82 -2.99 3348 1
## mu[5] 6.31 0.02 0.94 4.50 5.68 6.31 6.96 8.15 3565 1
## mu[6] -0.62 0.01 0.67 -1.90 -1.08 -0.62 -0.17 0.73 3628 1
## mu[7] -10.43 0.02 0.96 -12.30 -11.04 -10.44 -9.78 -8.52 3086 1
## mu[8] 2.78 0.01 0.65 1.52 2.34 2.78 3.22 4.04 4000 1
## mu[9] -14.34 0.02 0.99 -16.28 -14.99 -14.34 -13.68 -12.34 3163 1
## mu[10] -5.93 0.01 0.70 -7.27 -6.41 -5.93 -5.47 -4.55 4000 1
## mu[11] 7.98 0.02 0.80 6.40 7.45 7.98 8.50 9.56 2049 1
## mu[12] 4.70 0.02 0.85 3.05 4.14 4.70 5.27 6.38 2237 1
## mu[13] 12.89 0.02 0.99 10.96 12.24 12.89 13.53 14.85 1954 1
## mu[14] 18.84 0.02 1.05 16.76 18.15 18.86 19.51 20.96 2135 1
## mu[15] 5.97 0.02 0.84 4.37 5.41 5.98 6.52 7.62 3010 1
## mu[16] 20.01 0.02 0.86 18.28 19.45 20.01 20.60 21.72 2978 1
## mu[17] -9.47 0.01 0.82 -11.05 -10.02 -9.48 -8.93 -7.86 3450 1
## mu[18] 25.56 0.02 1.08 23.40 24.85 25.56 26.25 27.73 2430 1
## mu[19] -5.30 0.01 0.85 -6.95 -5.89 -5.32 -4.74 -3.62 4000 1
## mu[20] 8.19 0.02 1.10 6.10 7.45 8.20 8.93 10.36 2806 1
## mu[21] 8.17 0.02 0.96 6.28 7.53 8.17 8.82 10.07 3052 1
## mu[22] -2.70 0.02 1.02 -4.67 -3.41 -2.72 -2.00 -0.69 3619 1
## mu[23] -18.25 0.02 1.07 -20.35 -18.96 -18.27 -17.53 -16.12 3196 1
## mu[24] -0.50 0.01 0.73 -1.91 -1.01 -0.50 -0.01 0.97 4000 1
## mu[25] 9.33 0.02 0.84 7.68 8.78 9.33 9.89 10.98 3106 1
## mu[26] 4.34 0.01 0.69 3.00 3.88 4.35 4.80 5.69 4000 1
## mu[27] 9.03 0.01 0.61 7.83 8.63 9.05 9.42 10.24 4000 1
## mu[28] 1.36 0.01 0.81 -0.18 0.81 1.37 1.90 2.93 4000 1
## mu[29] -13.82 0.02 0.87 -15.54 -14.41 -13.83 -13.24 -12.11 3241 1
## mu[30] 5.91 0.02 1.05 3.84 5.21 5.93 6.60 8.02 2727 1
## mu[31] 11.18 0.02 0.89 9.41 10.57 11.16 11.80 12.92 3235 1
## mu[32] -7.12 0.02 1.10 -9.27 -7.84 -7.15 -6.38 -4.91 2880 1
## mu[33] -8.00 0.02 0.96 -9.85 -8.63 -8.00 -7.39 -6.08 2377 1
## mu[34] 12.59 0.02 0.92 10.79 11.97 12.60 13.22 14.38 3391 1
## mu[35] -10.01 0.02 1.04 -11.98 -10.72 -10.03 -9.33 -7.95 3902 1
## mu[36] 21.09 0.02 0.99 19.14 20.44 21.08 21.77 23.06 2901 1
## mu[37] 12.83 0.02 0.86 11.12 12.27 12.84 13.38 14.53 2066 1
## mu[38] 6.69 0.01 0.82 5.10 6.13 6.69 7.23 8.29 3010 1
## mu[39] -2.65 0.01 0.59 -3.80 -3.06 -2.65 -2.26 -1.48 4000 1
## mu[40] -22.63 0.02 1.22 -24.99 -23.43 -22.65 -21.82 -20.19 2655 1
## mu[41] 1.05 0.01 0.60 -0.12 0.64 1.06 1.45 2.24 4000 1
## mu[42] 12.20 0.02 1.01 10.19 11.53 12.20 12.86 14.19 2107 1
## mu[43] -5.88 0.02 1.04 -7.86 -6.56 -5.89 -5.20 -3.82 2246 1
## mu[44] 4.40 0.02 1.21 1.98 3.58 4.41 5.18 6.81 3211 1
## mu[45] 24.03 0.02 1.02 22.09 23.33 24.04 24.73 26.04 3078 1
## mu[46] 11.76 0.01 0.65 10.44 11.33 11.77 12.18 13.06 2626 1
## mu[47] 15.58 0.01 0.89 13.88 14.97 15.57 16.20 17.30 3658 1
## mu[48] -19.49 0.03 1.32 -22.07 -20.36 -19.49 -18.65 -16.81 2378 1
## mu[49] 5.72 0.01 0.57 4.59 5.33 5.72 6.10 6.81 4000 1
## mu[50] 21.19 0.02 1.17 18.92 20.39 21.20 21.98 23.44 3081 1
## lp__ -83.02 0.04 1.54 -86.85 -83.84 -82.70 -81.88 -80.97 1788 1
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 4 16:03:13 2016.
## 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).Since Rhat is 1 for all paramters, it seems that we reached at convergence. However, if there is a parameter with Rhat > 1.01, we should sample the chains longer.
We can create some plots as follows.
plot(fit_4, pars = 'beta')## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
Let’s make a trace plot. We can use ggmcmc::ggmcmc(). The output is saved in the file named (d2-model1-traceplot.pdf)[output/d2-model1-traceplots] in the output folder.
ggmcmc(ggs(fit_4, inc_warmup = TRUE), plot = 'traceplot',
file = 'output/d2-model1-traceplots.pdf')## Plotting traceplots## Time taken to generate the report: 15 seconds.You can extract the sampled values by rstan::extract(). Since dplyr package also contains the function called extract(), it is safer to call it with the package name.
stan_sample <- rstan::extract(fit_4)Let’s check some posterior distributions.
df <- data.frame(b1 = stan_sample$beta[,1],
b2 = stan_sample$beta[,2],
b3 = stan_sample$beta[,3],
b4 = stan_sample$beta[,4])
hist_beta <- ggplot(df, aes(y = ..density..))
hist_b1 <- hist_beta + geom_histogram(aes(x = b1), color = 'black') +
geom_vline(aes(xintercept = 0), color = 'red')
plot(hist_b1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
hist_b2 <- hist_beta + geom_histogram(aes(x = b2), color = 'black') +
geom_vline(aes(xintercept = 0), color = 'red')
plot(hist_b2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
hist_b3 <- hist_beta + geom_histogram(aes(x = b3), color = 'black') +
geom_vline(aes(xintercept = 0), color = 'red')
plot(hist_b3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
hist_b4 <- hist_beta + geom_histogram(aes(x = b4), color = 'black') +
geom_vline(aes(xintercept = 0), color = 'red')
plot(hist_b4)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Now let’s make a figure presenting the predicted values of y. Since we have to sample of estimates, we can calculate the predicted values. However, it is easier to ask Stan to compute the values. Let’s modify the model as follows.
data {
int<lower=0> N;
real Y[N];
real X2[N];
real X3[N];
real X4[N];
}
parameters {
real beta[4];
real<lower=0> sigma;
}
transformed parameters {
real mu[N];
for (n in 1:N) {
mu[n] = beta[1] + beta[2]*X2[n] + beta[3]*X3[n] + beta[4]*X4[n];
}
}
model {
for (n in 1:N) {
Y[n] ~ normal(mu[n], sigma);
}
for (p in 1:4) {
beta[p] ~ normal(0, 100);
}
sigma ~ uniform(0, 100);
}
generated quantities {
real y_pred[N];
for (n in 1:N) {
y_pred[n] = normal_rng(mu[n], sigma);
}
}Save this model as d2-model1b.stan in the stan-models folder.
We need to compile the model and sample the values.
stan_m1b <- stan_model('stan-models/d2-model1b.stan')## In file included from file21c6770f60bd.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:42:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
## In file included from file21c6770f60bd.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:43:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
## static void set_zero_all_adjoints_nested() {
## ^
## In file included from file21c6770f60bd.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:54:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
## ^
## In file included from file21c6770f60bd.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:235:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint.hpp:61:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array.hpp:21:
## In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/base.hpp:28:
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: unused typedef 'index_range' [-Wunused-local-typedef]
## typedef typename Array::index_range index_range;
## ^
## /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: unused typedef 'index' [-Wunused-local-typedef]
## typedef typename Array::index index;
## ^
## 7 warnings generated.fit_5 <- sampling(stan_m1b, data = dd1, seed = 777,
chains = 4, iter = 2000, warmup = 1000)Check the result.
print(fit_5)## Inference for Stan model: d2-model1b.
## 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
## beta[1] 3.14 0.02 1.07 1.12 2.43 3.11 3.84 5.30 1999
## beta[2] 0.08 0.00 0.17 -0.26 -0.03 0.08 0.19 0.40 1914
## beta[3] 1.50 0.00 0.14 1.23 1.41 1.50 1.60 1.77 2938
## beta[4] -2.05 0.00 0.09 -2.22 -2.11 -2.05 -2.00 -1.88 2866
## sigma 3.32 0.01 0.35 2.73 3.08 3.30 3.53 4.11 2326
## mu[1] -7.33 0.02 0.80 -8.82 -7.87 -7.34 -6.81 -5.70 2499
## mu[2] -0.95 0.02 0.89 -2.75 -1.53 -0.94 -0.35 0.83 2933
## mu[3] 2.65 0.02 1.01 0.68 1.97 2.66 3.32 4.64 2898
## mu[4] -4.21 0.01 0.60 -5.36 -4.62 -4.22 -3.81 -2.99 3132
## mu[5] 6.35 0.02 0.96 4.51 5.72 6.35 6.99 8.23 3381
## mu[6] -0.62 0.01 0.68 -1.97 -1.08 -0.61 -0.15 0.73 3302
## mu[7] -10.42 0.02 0.97 -12.36 -11.05 -10.43 -9.76 -8.53 2350
## mu[8] 2.79 0.01 0.66 1.47 2.36 2.79 3.22 4.08 3479
## mu[9] -14.33 0.02 0.98 -16.35 -14.99 -14.35 -13.67 -12.42 2520
## mu[10] -5.93 0.01 0.70 -7.35 -6.39 -5.93 -5.45 -4.53 3027
## mu[11] 7.98 0.02 0.82 6.42 7.43 7.97 8.52 9.63 2275
## mu[12] 4.72 0.02 0.87 3.06 4.14 4.71 5.28 6.52 2393
## mu[13] 12.91 0.02 1.01 10.96 12.24 12.88 13.57 14.90 2126
## mu[14] 18.85 0.02 1.06 16.78 18.14 18.84 19.55 20.95 2325
## mu[15] 5.99 0.02 0.85 4.28 5.43 6.00 6.56 7.67 2646
## mu[16] 20.03 0.01 0.84 18.36 19.47 20.02 20.59 21.67 3226
## mu[17] -9.48 0.02 0.81 -11.12 -10.01 -9.48 -8.93 -7.84 2838
## mu[18] 25.59 0.02 1.06 23.53 24.88 25.57 26.30 27.63 2812
## mu[19] -5.28 0.02 0.86 -7.00 -5.84 -5.28 -4.70 -3.57 2732
## mu[20] 8.23 0.02 1.13 6.04 7.48 8.22 8.96 10.46 2727
## mu[21] 8.21 0.02 0.98 6.28 7.55 8.20 8.87 10.16 2807
## mu[22] -2.67 0.02 1.04 -4.67 -3.35 -2.67 -1.96 -0.63 2820
## mu[23] -18.24 0.02 1.05 -20.32 -18.92 -18.25 -17.55 -16.19 2626
## mu[24] -0.47 0.01 0.74 -1.88 -0.95 -0.48 0.02 0.99 4000
## mu[25] 9.37 0.02 0.85 7.68 8.80 9.36 9.93 11.05 2951
## mu[26] 4.37 0.01 0.70 3.02 3.90 4.35 4.83 5.77 3772
## mu[27] 9.04 0.01 0.61 7.84 8.65 9.04 9.44 10.20 3634
## mu[28] 1.37 0.02 0.83 -0.27 0.82 1.38 1.91 2.99 2920
## mu[29] -13.82 0.02 0.84 -15.51 -14.39 -13.82 -13.27 -12.16 2778
## mu[30] 5.90 0.02 1.08 3.81 5.18 5.90 6.61 8.00 2552
## mu[31] 11.22 0.02 0.91 9.44 10.63 11.19 11.81 12.99 3236
## mu[32] -7.10 0.02 1.11 -9.23 -7.85 -7.10 -6.38 -4.86 2738
## mu[33] -7.99 0.02 0.95 -9.78 -8.63 -8.00 -7.37 -6.04 2461
## mu[34] 12.61 0.02 0.92 10.80 12.00 12.61 13.22 14.39 3032
## mu[35] -9.99 0.02 1.04 -11.95 -10.69 -10.00 -9.30 -7.99 3183
## mu[36] 21.13 0.02 0.97 19.27 20.48 21.13 21.78 23.09 3177
## mu[37] 12.84 0.02 0.87 11.16 12.24 12.83 13.43 14.56 2292
## mu[38] 6.71 0.02 0.83 5.07 6.17 6.71 7.26 8.37 2695
## mu[39] -2.64 0.01 0.58 -3.73 -3.03 -2.64 -2.24 -1.49 3316
## mu[40] -22.64 0.02 1.18 -24.95 -23.40 -22.63 -21.86 -20.32 2632
## mu[41] 1.08 0.01 0.60 -0.09 0.68 1.07 1.48 2.29 4000
## mu[42] 12.22 0.02 1.03 10.22 11.54 12.19 12.88 14.31 2219
## mu[43] -5.87 0.02 1.04 -7.83 -6.58 -5.88 -5.20 -3.73 2353
## mu[44] 4.38 0.02 1.25 1.91 3.54 4.38 5.21 6.83 2732
## mu[45] 24.06 0.02 0.99 22.11 23.40 24.05 24.72 26.01 3321
## mu[46] 11.77 0.01 0.65 10.52 11.34 11.77 12.20 13.03 3055
## mu[47] 15.60 0.02 0.88 13.86 15.03 15.61 16.19 17.31 3447
## mu[48] -19.52 0.03 1.31 -22.10 -20.38 -19.52 -18.67 -16.96 2480
## mu[49] 5.73 0.01 0.57 4.59 5.35 5.74 6.10 6.86 4000
## mu[50] 21.21 0.02 1.16 18.92 20.44 21.22 21.99 23.49 2984
## y_pred[1] -7.36 0.06 3.37 -13.87 -9.57 -7.38 -5.11 -0.56 3721
## y_pred[2] -0.90 0.05 3.45 -7.70 -3.18 -0.88 1.44 5.78 3946
## y_pred[3] 2.60 0.05 3.45 -4.28 0.33 2.64 4.90 9.33 3986
## y_pred[4] -4.30 0.05 3.41 -10.89 -6.59 -4.25 -2.05 2.34 4000
## y_pred[5] 6.39 0.06 3.49 -0.62 4.13 6.42 8.75 13.17 4000
## y_pred[6] -0.57 0.06 3.49 -7.37 -2.84 -0.56 1.71 6.37 3924
## y_pred[7] -10.38 0.05 3.44 -17.15 -12.70 -10.45 -8.10 -3.62 4000
## y_pred[8] 2.82 0.05 3.39 -3.96 0.59 2.84 5.01 9.59 4000
## y_pred[9] -14.33 0.06 3.55 -21.37 -16.70 -14.32 -12.02 -7.54 3995
## y_pred[10] -5.93 0.05 3.41 -12.62 -8.19 -5.92 -3.66 0.63 4000
## y_pred[11] 7.94 0.05 3.44 1.19 5.60 7.95 10.20 14.75 4000
## y_pred[12] 4.80 0.06 3.43 -2.12 2.55 4.77 7.08 11.66 3634
## y_pred[13] 12.88 0.06 3.47 6.10 10.66 12.87 15.16 19.62 3745
## y_pred[14] 18.93 0.05 3.46 12.23 16.59 18.86 21.25 25.59 4000
## y_pred[15] 5.98 0.06 3.46 -0.99 3.70 5.94 8.29 12.77 3632
## y_pred[16] 19.93 0.05 3.42 13.16 17.71 19.97 22.22 26.76 4000
## y_pred[17] -9.42 0.06 3.40 -15.91 -11.70 -9.41 -7.11 -2.79 3525
## y_pred[18] 25.63 0.06 3.51 18.58 23.35 25.65 27.93 32.70 3845
## y_pred[19] -5.27 0.05 3.42 -11.97 -7.55 -5.27 -3.07 1.60 4000
## y_pred[20] 8.29 0.06 3.55 1.29 5.90 8.24 10.63 15.42 3884
## y_pred[21] 8.18 0.06 3.40 1.57 5.94 8.17 10.49 14.97 3715
## y_pred[22] -2.71 0.06 3.55 -9.86 -5.06 -2.70 -0.27 4.32 3834
## y_pred[23] -18.23 0.06 3.56 -25.36 -20.54 -18.27 -15.84 -11.16 4000
## y_pred[24] -0.49 0.05 3.47 -7.26 -2.81 -0.48 1.85 6.31 4000
## y_pred[25] 9.43 0.06 3.44 2.74 7.10 9.40 11.72 16.17 3902
## y_pred[26] 4.35 0.05 3.32 -2.10 2.15 4.27 6.54 10.84 4000
## y_pred[27] 9.03 0.06 3.42 2.39 6.73 9.00 11.29 15.68 3770
## y_pred[28] 1.37 0.05 3.41 -5.54 -0.83 1.40 3.65 7.97 4000
## y_pred[29] -13.83 0.05 3.40 -20.66 -16.09 -13.78 -11.53 -7.11 3828
## y_pred[30] 5.81 0.06 3.50 -1.07 3.49 5.76 8.10 12.72 3935
## y_pred[31] 11.19 0.06 3.48 4.50 8.86 11.26 13.46 17.94 3841
## y_pred[32] -7.01 0.05 3.47 -13.75 -9.29 -7.06 -4.77 -0.01 3998
## y_pred[33] -7.99 0.06 3.50 -14.94 -10.27 -7.99 -5.64 -1.29 3780
## y_pred[34] 12.66 0.06 3.45 5.85 10.27 12.70 15.02 19.48 3840
## y_pred[35] -9.98 0.06 3.52 -16.80 -12.29 -9.94 -7.62 -3.06 3844
## y_pred[36] 21.13 0.06 3.48 14.28 18.85 21.12 23.48 27.86 3753
## y_pred[37] 12.92 0.05 3.35 6.46 10.69 12.83 15.16 19.69 3833
## y_pred[38] 6.69 0.06 3.45 -0.14 4.35 6.66 8.87 13.66 3469
## y_pred[39] -2.61 0.06 3.42 -9.47 -4.82 -2.59 -0.32 3.97 3760
## y_pred[40] -22.54 0.06 3.51 -29.41 -24.88 -22.53 -20.26 -15.51 3654
## y_pred[41] 1.09 0.05 3.41 -5.81 -1.11 1.03 3.28 7.97 3955
## y_pred[42] 12.19 0.06 3.45 5.47 9.90 12.22 14.51 19.05 3927
## y_pred[43] -5.91 0.06 3.48 -12.88 -8.23 -5.90 -3.63 0.80 3858
## y_pred[44] 4.39 0.06 3.53 -2.39 1.97 4.43 6.66 11.34 3510
## y_pred[45] 24.02 0.06 3.54 16.88 21.73 24.01 26.34 31.04 3551
## y_pred[46] 11.75 0.06 3.46 4.95 9.48 11.70 14.05 18.42 3838
## y_pred[47] 15.69 0.05 3.48 8.71 13.40 15.68 17.96 22.57 4000
## y_pred[48] -19.43 0.06 3.63 -26.48 -21.82 -19.51 -17.06 -12.19 3674
## y_pred[49] 5.68 0.05 3.37 -0.88 3.40 5.75 7.92 12.36 4000
## y_pred[50] 21.25 0.06 3.55 14.28 18.92 21.18 23.63 28.25 3730
## lp__ -83.04 0.04 1.58 -86.92 -83.84 -82.71 -81.89 -80.95 1759
## Rhat
## beta[1] 1
## beta[2] 1
## beta[3] 1
## beta[4] 1
## sigma 1
## mu[1] 1
## mu[2] 1
## mu[3] 1
## mu[4] 1
## mu[5] 1
## mu[6] 1
## mu[7] 1
## mu[8] 1
## mu[9] 1
## mu[10] 1
## mu[11] 1
## mu[12] 1
## mu[13] 1
## mu[14] 1
## mu[15] 1
## mu[16] 1
## mu[17] 1
## mu[18] 1
## mu[19] 1
## mu[20] 1
## mu[21] 1
## mu[22] 1
## mu[23] 1
## mu[24] 1
## mu[25] 1
## mu[26] 1
## mu[27] 1
## mu[28] 1
## mu[29] 1
## mu[30] 1
## mu[31] 1
## mu[32] 1
## mu[33] 1
## mu[34] 1
## mu[35] 1
## mu[36] 1
## mu[37] 1
## mu[38] 1
## mu[39] 1
## mu[40] 1
## mu[41] 1
## mu[42] 1
## mu[43] 1
## mu[44] 1
## mu[45] 1
## mu[46] 1
## mu[47] 1
## mu[48] 1
## mu[49] 1
## mu[50] 1
## y_pred[1] 1
## y_pred[2] 1
## y_pred[3] 1
## y_pred[4] 1
## y_pred[5] 1
## y_pred[6] 1
## y_pred[7] 1
## y_pred[8] 1
## y_pred[9] 1
## y_pred[10] 1
## y_pred[11] 1
## y_pred[12] 1
## y_pred[13] 1
## y_pred[14] 1
## y_pred[15] 1
## y_pred[16] 1
## y_pred[17] 1
## y_pred[18] 1
## y_pred[19] 1
## y_pred[20] 1
## y_pred[21] 1
## y_pred[22] 1
## y_pred[23] 1
## y_pred[24] 1
## y_pred[25] 1
## y_pred[26] 1
## y_pred[27] 1
## y_pred[28] 1
## y_pred[29] 1
## y_pred[30] 1
## y_pred[31] 1
## y_pred[32] 1
## y_pred[33] 1
## y_pred[34] 1
## y_pred[35] 1
## y_pred[36] 1
## y_pred[37] 1
## y_pred[38] 1
## y_pred[39] 1
## y_pred[40] 1
## y_pred[41] 1
## y_pred[42] 1
## y_pred[43] 1
## y_pred[44] 1
## y_pred[45] 1
## y_pred[46] 1
## y_pred[47] 1
## y_pred[48] 1
## y_pred[49] 1
## y_pred[50] 1
## lp__ 1
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 4 16:04:24 2016.
## 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).Looking at Rhats, it looks OK.
Extract the sample.
stan_sample <- rstan::extract(fit_5)Now this list contains the predicted values of y. Using these values, let’s make a plot.
d_quantiles <- t(apply(stan_sample$y_pred, 2, quantile, prob = c(.1, .5, .9)))
colnames(d_quantiles) <- c('p10', 'p50', 'p90')
d_quantiles <- data.frame(dd, d_quantiles)
p <- ggplot(d_quantiles, aes(x = y, y = p50, ymin = p10, ymax = p90)) +
coord_fixed(ratio = 1) +
geom_pointrange() +
geom_abline(aes(intercept = 0, slope = 1), color = 'black', linetype = 'dashed') +
labs(x = 'observed outcome', y = 'predicted outcome')
print(p)
One big advantage of Bayesian estimation is that we can calculate a variety of probabilities because we have distributions of parameters. For instance, the probability that \(\beta_2\) is larger than 0 is
sum(stan_sample$beta[,2] > 0) / length(stan_sample$beta[,2])## [1] 0.67975Similarly, the probability that \(\sigma\) is in the interval \([2.5, 3.5]\) is
sum(stan_sample$sigma >= 2.5 & stan_sample$sigma <= 3.5) / length(stan_sample$sigma)## [1] 0.711