Chapter 9 - Markov Chain Monte Carlo

This week has been an informal introduction to Markov chain Monte Carlo (MCMC) estimation. The goal has been to introduce the purpose and approach MCMC algorithms. The major algorithms introduced were the Metropolis, Gibbs sampling, and Hamiltonian Monte Carlo algorithms. Each has its advantages and disadvantages. The ulam function in the rethinking package was introduced. It uses the Stan (mc-stan.org) Hamiltonian Monte Carlo engine to fit models as they are defined in this book. General advice about diagnosing poor MCMC fits was introduced by the use of a couple of pathological examples.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.

Questions

8-1. Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior for the standard deviation, sigma. The uniform prior should be dunif(0,1). Visualize the priors. Use ulam to estimate the posterior distribution of sigma. Visualize the posterior distribution of sigma for both models. Do not use a ‘pairs’ plot. Does the different prior have any detectable influence on the posterior distribution of sigma? Why or why not?

library(rethinking)
library(ggplot2)

# From page 426
data(rugged)
d <- rugged

# Convert outcome to log
d$log_gdp <- log(d$rgdppc_2000)

# Extract countries with GDP data
dd <- d[ complete.cases(d$rgdppc_2000) , ]

# Re-scale variables
dd$log_gdp_std <- dd$log_gdp / mean(dd$log_gdp)
dd$rugged_std <- dd$rugged / max(dd$rugged)

dd$cid <- ifelse( dd$cont_africa==1 , 1 , 2 )

# Create a sl;m list of the variables we will use
dat_slim <- list(
    log_gdp_std = dd$log_gdp_std,
    rugged_std = dd$rugged_std,
    cid = as.integer(dd$cid)
)
str(dat_slim)
## List of 3
##  $ log_gdp_std: num [1:170] 0.88 0.965 1.166 1.104 0.915 ...
##  $ rugged_std : num [1:170] 0.138 0.553 0.124 0.125 0.433 ...
##  $ cid        : int [1:170] 1 2 2 2 2 2 2 2 2 1 ...
# Original model - exponential sigma
original_model <- ulam(
    alist(
        log_gdp_std ~ dnorm(mu , sigma) ,
        mu <- a[cid] + b[cid]*(rugged_std - 0.215) ,
        a[cid] ~ dnorm(1 , 0.1) ,
        b[cid] ~ dnorm(0 , 0.3) ,
        sigma ~ dexp(1)
    ) , data=dat_slim , chains=4, cores=4)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
# Show result
precis(original_model, 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8870902 0.016150827  0.86159893  0.91275620 3037.640 0.9988386
## a[2]   1.0506090 0.010294473  1.03392835  1.06622000 2913.346 0.9989322
## b[1]   0.1335428 0.078115160  0.01127229  0.25501198 2508.168 1.0010718
## b[2]  -0.1424169 0.056966574 -0.23073058 -0.05159199 2361.610 1.0002380
## sigma  0.1117761 0.006027175  0.10270351  0.12181394 2252.754 1.0007087
# New model - uniform sigma
new_model <- ulam(
    alist(
        log_gdp_std ~ dnorm(mu , sigma) ,
        mu <- a[cid] + b[cid]*(rugged_std - 0.215) ,
        a[cid] ~ dnorm(1 , 0.1) ,
        b[cid] ~ dnorm(0 , 0.3) ,
        sigma ~ dunif(0, 1)
    ) , data=dat_slim , chains=4, cores=4)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
# Show result
precis(new_model, 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8864086 0.015875246  0.86107267  0.91202077 4070.257 0.9993466
## a[2]   1.0506908 0.009874301  1.03518890  1.06688055 3106.583 0.9983578
## b[1]   0.1352215 0.073435468  0.01454241  0.24995550 2725.804 0.9993366
## b[2]  -0.1414425 0.057439764 -0.23454007 -0.04883132 2487.495 1.0002813
## sigma  0.1116231 0.006083988  0.10221314  0.12152027 2423.674 1.0001604
# Visualize priors of sigma
exponential <- rexp(1e3, 1)
uniform <- runif(1e3, 1, 8)

# Create dataframes
exponential_df <- with(density(exponential), data.frame(x,y))
uniform_df <- with(density(uniform), data.frame(x,y))
colors <- c("Exponential" = "red",
            "Uniform" = "blue")

sigma_prior_plot <- ggplot() +
  geom_line(data=exponential_df, aes(x=x, y=y, color="Exponential"), size=1.5) +
  geom_line(data=uniform_df, aes(x=x, y=y, color="Uniform"), size=1.5) +
  ggtitle("Sigma Prior Distributions for The Two Models") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Prior", y="", color="Sigma Prior Distribution") +
  scale_color_manual(values = colors)
  
sigma_prior_plot

# Visualize posterior of sigma
original_model_sigma <- extract.samples(original_model)$sigma
new_model_sigma <- extract.samples(new_model)$sigma

# Create dataframes
exponential_post_df <- with(density(original_model_sigma), data.frame(x,y))
uniform_post_df <- with(density(new_model_sigma), data.frame(x,y))
colors <- c("Exponential" = "red",
            "Uniform" = "blue")

sigma_posterior_plot <- ggplot() +
  geom_line(data=exponential_post_df, aes(x=x, y=y, color="Exponential"), size=1.5) +
  geom_line(data=uniform_post_df, aes(x=x, y=y, color="Uniform"), size=1.5) +
  ggtitle("Sigma Posterior Distributions for The Two Models") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Posterior", y="", color="Sigma Prior Distribution") +
  scale_color_manual(values = colors)
  
sigma_posterior_plot

# The different prior does not have any detectable influence of the posterior distribution of sigma. This is evident from the Sigma Posterior Distributions plot where we can see that the posterior for both the exponential and uniform prior distributions look identical.

8-2. Modify the terrain ruggedness model again. This time, change the prior for b[cid] to dexp(0.3). Plot the joint posterior. Do not use a ‘pairs’ plot. What does this do to the posterior distribution? Can you explain it?

library(data.table)

# New model - uniform sigma
model_8.2 <- ulam(
    alist(
        log_gdp_std ~ dnorm(mu , sigma) ,
        mu <- a[cid] + b[cid]*(rugged_std - 0.215) ,
        a[cid] ~ dnorm(1 , 0.1) ,
        b[cid] ~ dexp(0.3) ,
        sigma ~ dexp(1)
    ) , data=dat_slim , chains=4, cores=4)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
# Show result
precis(model_8.2, 2)
##             mean          sd        5.5%      94.5%     n_eff     Rhat4
## a[1]  0.88693280 0.016798110 0.860422425 0.91349771 1614.6363 0.9990618
## a[2]  1.04777338 0.010373793 1.031426150 1.06433055 1297.8974 0.9999609
## b[1]  0.14633715 0.074810347 0.031152361 0.26964049  839.8453 1.0005730
## b[2]  0.01912431 0.016798358 0.001688269 0.05084325 2001.3254 1.0018782
## sigma 0.11401241 0.006323739 0.104457160 0.12478681 1306.4543 0.9984950
# Visualize joint posterior
b <- extract.samples(model_8.2)$b[,1]
sigma <- extract.samples(model_8.2)$sigma

b_df <- data.table(x=b)
sigma_df <- data.table(x=sigma)

b_posterior_plot <- ggplot(b_df, aes(x=x)) +
  geom_density(fill="gray", alpha=0.5) +
  ggtitle("Posterior for b") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Posterior", y="")
  
b_posterior_plot

# When we change the posterior for b to an exponential distribution, we see that the shape of the posterior doesn't change much, however, it does cut off at 0, since exponential distribution is positive definite.

8-3. Re-estimate one of the Stan models from the chapter, but at different numbers (at least 5) of warmup iterations. Be sure to use the same number of sampling iterations in each case. Compare the n_eff values. How much warm up is enough?

warmup_test <- function(N){
  original_model <- ulam(
    alist(
        log_gdp_std ~ dnorm(mu , sigma) ,
        mu <- a[cid] + b[cid]*(rugged_std - 0.215) ,
        a[cid] ~ dnorm(1 , 0.1) ,
        b[cid] ~ dnorm(0 , 0.3) ,
        sigma ~ dexp(1)
    ) , data=dat_slim , chains=4, cores=4, warmup=N, iter=1000+N)
  
  # Result
  result <- precis(original_model, 2)$n_eff
  return(result)
}

# Simulate WAIC for different number of data points
sim_list <- c(1, 5, 10, 50, 100, 500, 1000, 2000)
warmup_sim <- matrix(NA, nrow=length(sim_list), ncol=5)
for (i in 1:length(sim_list)) (warmup_sim[i,] <- warmup_test(sim_list[i]))
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 WARNING: No variance estimation is 
## Chain 1          performed for num_warmup < 20 
## Chain 1 Iteration:    1 / 1001 [  0%]  (Warmup) 
## Chain 1 Iteration:    2 / 1001 [  0%]  (Sampling) 
## Chain 2 WARNING: No variance estimation is 
## Chain 2          performed for num_warmup < 20 
## Chain 2 Iteration:    1 / 1001 [  0%]  (Warmup) 
## Chain 2 Iteration:    2 / 1001 [  0%]  (Sampling) 
## Chain 2 Iteration:  101 / 1001 [ 10%]  (Sampling) 
## Chain 2 Iteration:  201 / 1001 [ 20%]  (Sampling) 
## Chain 2 Iteration:  301 / 1001 [ 30%]  (Sampling) 
## Chain 2 Iteration:  401 / 1001 [ 40%]  (Sampling) 
## Chain 2 Iteration:  501 / 1001 [ 50%]  (Sampling) 
## Chain 2 Iteration:  601 / 1001 [ 60%]  (Sampling) 
## Chain 2 Iteration:  701 / 1001 [ 70%]  (Sampling) 
## Chain 2 Iteration:  801 / 1001 [ 80%]  (Sampling) 
## Chain 2 Iteration:  901 / 1001 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1001 / 1001 [100%]  (Sampling) 
## Chain 3 WARNING: No variance estimation is 
## Chain 3          performed for num_warmup < 20 
## Chain 3 Iteration:    1 / 1001 [  0%]  (Warmup) 
## Chain 3 Iteration:    2 / 1001 [  0%]  (Sampling) 
## Chain 4 WARNING: No variance estimation is 
## Chain 4          performed for num_warmup < 20 
## Chain 4 Iteration:    1 / 1001 [  0%]  (Warmup) 
## Chain 4 Iteration:    2 / 1001 [  0%]  (Sampling) 
## Chain 2 finished in 0.0 seconds.
## Chain 1 Iteration:  101 / 1001 [ 10%]  (Sampling) 
## Chain 3 Iteration:  101 / 1001 [ 10%]  (Sampling) 
## Chain 4 Iteration:  101 / 1001 [ 10%]  (Sampling) 
## Chain 1 Iteration:  201 / 1001 [ 20%]  (Sampling) 
## Chain 3 Iteration:  201 / 1001 [ 20%]  (Sampling) 
## Chain 4 Iteration:  201 / 1001 [ 20%]  (Sampling) 
## Chain 1 Iteration:  301 / 1001 [ 30%]  (Sampling) 
## Chain 3 Iteration:  301 / 1001 [ 30%]  (Sampling) 
## Chain 4 Iteration:  301 / 1001 [ 30%]  (Sampling) 
## Chain 1 Iteration:  401 / 1001 [ 40%]  (Sampling) 
## Chain 3 Iteration:  401 / 1001 [ 40%]  (Sampling) 
## Chain 4 Iteration:  401 / 1001 [ 40%]  (Sampling) 
## Chain 1 Iteration:  501 / 1001 [ 50%]  (Sampling) 
## Chain 3 Iteration:  501 / 1001 [ 50%]  (Sampling) 
## Chain 4 Iteration:  501 / 1001 [ 50%]  (Sampling) 
## Chain 1 Iteration:  601 / 1001 [ 60%]  (Sampling) 
## Chain 3 Iteration:  601 / 1001 [ 60%]  (Sampling) 
## Chain 4 Iteration:  601 / 1001 [ 60%]  (Sampling) 
## Chain 1 Iteration:  701 / 1001 [ 70%]  (Sampling) 
## Chain 3 Iteration:  701 / 1001 [ 70%]  (Sampling) 
## Chain 4 Iteration:  701 / 1001 [ 70%]  (Sampling) 
## Chain 1 Iteration:  801 / 1001 [ 80%]  (Sampling) 
## Chain 3 Iteration:  801 / 1001 [ 80%]  (Sampling) 
## Chain 4 Iteration:  801 / 1001 [ 80%]  (Sampling) 
## Chain 1 Iteration:  901 / 1001 [ 90%]  (Sampling) 
## Chain 3 Iteration:  901 / 1001 [ 90%]  (Sampling) 
## Chain 4 Iteration:  901 / 1001 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1001 / 1001 [100%]  (Sampling) 
## Chain 1 finished in 8.4 seconds.
## Chain 3 Iteration: 1001 / 1001 [100%]  (Sampling) 
## Chain 3 finished in 8.4 seconds.
## Chain 4 Iteration: 1001 / 1001 [100%]  (Sampling) 
## Chain 4 finished in 8.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 6.3 seconds.
## Total execution time: 8.5 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 WARNING: No variance estimation is 
## Chain 1          performed for num_warmup < 20 
## Chain 1 Iteration:    1 / 1005 [  0%]  (Warmup) 
## Chain 1 Iteration:    6 / 1005 [  0%]  (Sampling) 
## Chain 1 Iteration:  105 / 1005 [ 10%]  (Sampling) 
## Chain 1 Iteration:  205 / 1005 [ 20%]  (Sampling) 
## Chain 1 Iteration:  305 / 1005 [ 30%]  (Sampling) 
## Chain 1 Iteration:  405 / 1005 [ 40%]  (Sampling) 
## Chain 1 Iteration:  505 / 1005 [ 50%]  (Sampling) 
## Chain 1 Iteration:  605 / 1005 [ 60%]  (Sampling) 
## Chain 1 Iteration:  705 / 1005 [ 70%]  (Sampling) 
## Chain 1 Iteration:  805 / 1005 [ 80%]  (Sampling) 
## Chain 1 Iteration:  905 / 1005 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1005 / 1005 [100%]  (Sampling) 
## Chain 2 WARNING: No variance estimation is 
## Chain 2          performed for num_warmup < 20 
## Chain 2 Iteration:    1 / 1005 [  0%]  (Warmup) 
## Chain 2 Iteration:    6 / 1005 [  0%]  (Sampling) 
## Chain 2 Iteration:  105 / 1005 [ 10%]  (Sampling) 
## Chain 2 Iteration:  205 / 1005 [ 20%]  (Sampling) 
## Chain 2 Iteration:  305 / 1005 [ 30%]  (Sampling) 
## Chain 2 Iteration:  405 / 1005 [ 40%]  (Sampling) 
## Chain 2 Iteration:  505 / 1005 [ 50%]  (Sampling) 
## Chain 2 Iteration:  605 / 1005 [ 60%]  (Sampling) 
## Chain 2 Iteration:  705 / 1005 [ 70%]  (Sampling) 
## Chain 2 Iteration:  805 / 1005 [ 80%]  (Sampling) 
## Chain 2 Iteration:  905 / 1005 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1005 / 1005 [100%]  (Sampling) 
## Chain 3 WARNING: No variance estimation is 
## Chain 3          performed for num_warmup < 20 
## Chain 3 Iteration:    1 / 1005 [  0%]  (Warmup) 
## Chain 3 Iteration:    6 / 1005 [  0%]  (Sampling) 
## Chain 3 Iteration:  105 / 1005 [ 10%]  (Sampling) 
## Chain 3 Iteration:  205 / 1005 [ 20%]  (Sampling) 
## Chain 3 Iteration:  305 / 1005 [ 30%]  (Sampling) 
## Chain 3 Iteration:  405 / 1005 [ 40%]  (Sampling) 
## Chain 3 Iteration:  505 / 1005 [ 50%]  (Sampling) 
## Chain 3 Iteration:  605 / 1005 [ 60%]  (Sampling) 
## Chain 3 Iteration:  705 / 1005 [ 70%]  (Sampling) 
## Chain 3 Iteration:  805 / 1005 [ 80%]  (Sampling) 
## Chain 3 Iteration:  905 / 1005 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1005 / 1005 [100%]  (Sampling) 
## Chain 4 WARNING: No variance estimation is 
## Chain 4          performed for num_warmup < 20 
## Chain 4 Iteration:    1 / 1005 [  0%]  (Warmup) 
## Chain 4 Iteration:    6 / 1005 [  0%]  (Sampling) 
## Chain 4 Iteration:  105 / 1005 [ 10%]  (Sampling) 
## Chain 4 Iteration:  205 / 1005 [ 20%]  (Sampling) 
## Chain 4 Iteration:  305 / 1005 [ 30%]  (Sampling) 
## Chain 4 Iteration:  405 / 1005 [ 40%]  (Sampling) 
## Chain 4 Iteration:  505 / 1005 [ 50%]  (Sampling) 
## Chain 4 Iteration:  605 / 1005 [ 60%]  (Sampling) 
## Chain 4 Iteration:  705 / 1005 [ 70%]  (Sampling) 
## Chain 4 Iteration:  805 / 1005 [ 80%]  (Sampling) 
## Chain 4 Iteration:  905 / 1005 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1005 / 1005 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 WARNING: No variance estimation is 
## Chain 1          performed for num_warmup < 20 
## Chain 1 Iteration:    1 / 1010 [  0%]  (Warmup) 
## Chain 1 Iteration:   11 / 1010 [  1%]  (Sampling) 
## Chain 1 Iteration:  110 / 1010 [ 10%]  (Sampling) 
## Chain 1 Iteration:  210 / 1010 [ 20%]  (Sampling) 
## Chain 1 Iteration:  310 / 1010 [ 30%]  (Sampling) 
## Chain 1 Iteration:  410 / 1010 [ 40%]  (Sampling) 
## Chain 1 Iteration:  510 / 1010 [ 50%]  (Sampling) 
## Chain 1 Iteration:  610 / 1010 [ 60%]  (Sampling) 
## Chain 1 Iteration:  710 / 1010 [ 70%]  (Sampling) 
## Chain 1 Iteration:  810 / 1010 [ 80%]  (Sampling) 
## Chain 1 Iteration:  910 / 1010 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1010 / 1010 [100%]  (Sampling) 
## Chain 2 WARNING: No variance estimation is 
## Chain 2          performed for num_warmup < 20 
## Chain 2 Iteration:    1 / 1010 [  0%]  (Warmup) 
## Chain 2 Iteration:   11 / 1010 [  1%]  (Sampling) 
## Chain 2 Iteration:  110 / 1010 [ 10%]  (Sampling) 
## Chain 2 Iteration:  210 / 1010 [ 20%]  (Sampling) 
## Chain 2 Iteration:  310 / 1010 [ 30%]  (Sampling) 
## Chain 2 Iteration:  410 / 1010 [ 40%]  (Sampling) 
## Chain 2 Iteration:  510 / 1010 [ 50%]  (Sampling) 
## Chain 2 Iteration:  610 / 1010 [ 60%]  (Sampling) 
## Chain 3 WARNING: No variance estimation is 
## Chain 3          performed for num_warmup < 20 
## Chain 3 Iteration:    1 / 1010 [  0%]  (Warmup) 
## Chain 3 Iteration:   11 / 1010 [  1%]  (Sampling) 
## Chain 3 Iteration:  110 / 1010 [ 10%]  (Sampling) 
## Chain 3 Iteration:  210 / 1010 [ 20%]  (Sampling) 
## Chain 3 Iteration:  310 / 1010 [ 30%]  (Sampling) 
## Chain 3 Iteration:  410 / 1010 [ 40%]  (Sampling) 
## Chain 3 Iteration:  510 / 1010 [ 50%]  (Sampling) 
## Chain 4 WARNING: No variance estimation is 
## Chain 4          performed for num_warmup < 20 
## Chain 4 Iteration:    1 / 1010 [  0%]  (Warmup) 
## Chain 4 Iteration:   11 / 1010 [  1%]  (Sampling) 
## Chain 4 Iteration:  110 / 1010 [ 10%]  (Sampling) 
## Chain 4 Iteration:  210 / 1010 [ 20%]  (Sampling) 
## Chain 4 Iteration:  310 / 1010 [ 30%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration:  710 / 1010 [ 70%]  (Sampling) 
## Chain 2 Iteration:  810 / 1010 [ 80%]  (Sampling) 
## Chain 2 Iteration:  910 / 1010 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1010 / 1010 [100%]  (Sampling) 
## Chain 3 Iteration:  610 / 1010 [ 60%]  (Sampling) 
## Chain 3 Iteration:  710 / 1010 [ 70%]  (Sampling) 
## Chain 3 Iteration:  810 / 1010 [ 80%]  (Sampling) 
## Chain 3 Iteration:  910 / 1010 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1010 / 1010 [100%]  (Sampling) 
## Chain 4 Iteration:  410 / 1010 [ 40%]  (Sampling) 
## Chain 4 Iteration:  510 / 1010 [ 50%]  (Sampling) 
## Chain 4 Iteration:  610 / 1010 [ 60%]  (Sampling) 
## Chain 4 Iteration:  710 / 1010 [ 70%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration:  810 / 1010 [ 80%]  (Sampling) 
## Chain 4 Iteration:  910 / 1010 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1010 / 1010 [100%]  (Sampling) 
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 7 
## Chain 1            adapt_window = 38 
## Chain 1            term_buffer = 5 
## Chain 1 Iteration:    1 / 1050 [  0%]  (Warmup) 
## Chain 1 Iteration:   51 / 1050 [  4%]  (Sampling) 
## Chain 1 Iteration:  150 / 1050 [ 14%]  (Sampling) 
## Chain 1 Iteration:  250 / 1050 [ 23%]  (Sampling) 
## Chain 1 Iteration:  350 / 1050 [ 33%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 7 
## Chain 2            adapt_window = 38 
## Chain 2            term_buffer = 5 
## Chain 2 Iteration:    1 / 1050 [  0%]  (Warmup) 
## Chain 2 Iteration:   51 / 1050 [  4%]  (Sampling) 
## Chain 2 Iteration:  150 / 1050 [ 14%]  (Sampling) 
## Chain 2 Iteration:  250 / 1050 [ 23%]  (Sampling) 
## Chain 2 Iteration:  350 / 1050 [ 33%]  (Sampling) 
## Chain 2 Iteration:  450 / 1050 [ 42%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 7 
## Chain 3            adapt_window = 38 
## Chain 3            term_buffer = 5 
## Chain 3 Iteration:    1 / 1050 [  0%]  (Warmup) 
## Chain 3 Iteration:   51 / 1050 [  4%]  (Sampling) 
## Chain 3 Iteration:  150 / 1050 [ 14%]  (Sampling) 
## Chain 3 Iteration:  250 / 1050 [ 23%]  (Sampling) 
## Chain 3 Iteration:  350 / 1050 [ 33%]  (Sampling) 
## Chain 3 Iteration:  450 / 1050 [ 42%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 7 
## Chain 4            adapt_window = 38 
## Chain 4            term_buffer = 5 
## Chain 4 Iteration:    1 / 1050 [  0%]  (Warmup) 
## Chain 4 Iteration:   51 / 1050 [  4%]  (Sampling) 
## Chain 4 Iteration:  150 / 1050 [ 14%]  (Sampling) 
## Chain 4 Iteration:  250 / 1050 [ 23%]  (Sampling) 
## Chain 4 Iteration:  350 / 1050 [ 33%]  (Sampling) 
## Chain 4 Iteration:  450 / 1050 [ 42%]  (Sampling) 
## Chain 4 Iteration:  550 / 1050 [ 52%]  (Sampling) 
## Chain 1 Iteration:  450 / 1050 [ 42%]  (Sampling) 
## Chain 1 Iteration:  550 / 1050 [ 52%]  (Sampling) 
## Chain 1 Iteration:  650 / 1050 [ 61%]  (Sampling) 
## Chain 1 Iteration:  750 / 1050 [ 71%]  (Sampling) 
## Chain 1 Iteration:  850 / 1050 [ 80%]  (Sampling) 
## Chain 1 Iteration:  950 / 1050 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1050 / 1050 [100%]  (Sampling) 
## Chain 2 Iteration:  550 / 1050 [ 52%]  (Sampling) 
## Chain 2 Iteration:  650 / 1050 [ 61%]  (Sampling) 
## Chain 2 Iteration:  750 / 1050 [ 71%]  (Sampling) 
## Chain 2 Iteration:  850 / 1050 [ 80%]  (Sampling) 
## Chain 2 Iteration:  950 / 1050 [ 90%]  (Sampling) 
## Chain 3 Iteration:  550 / 1050 [ 52%]  (Sampling) 
## Chain 3 Iteration:  650 / 1050 [ 61%]  (Sampling) 
## Chain 3 Iteration:  750 / 1050 [ 71%]  (Sampling) 
## Chain 3 Iteration:  850 / 1050 [ 80%]  (Sampling) 
## Chain 3 Iteration:  950 / 1050 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1050 / 1050 [100%]  (Sampling) 
## Chain 4 Iteration:  650 / 1050 [ 61%]  (Sampling) 
## Chain 4 Iteration:  750 / 1050 [ 71%]  (Sampling) 
## Chain 4 Iteration:  850 / 1050 [ 80%]  (Sampling) 
## Chain 4 Iteration:  950 / 1050 [ 90%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 2 Iteration: 1050 / 1050 [100%]  (Sampling) 
## Chain 4 Iteration: 1050 / 1050 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.5 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:    1 / 1100 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 1100 [  9%]  (Warmup) 
## Chain 1 Iteration:  101 / 1100 [  9%]  (Sampling) 
## Chain 1 Iteration:  200 / 1100 [ 18%]  (Sampling) 
## Chain 1 Iteration:  300 / 1100 [ 27%]  (Sampling) 
## Chain 1 Iteration:  400 / 1100 [ 36%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:    1 / 1100 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 1100 [  9%]  (Warmup) 
## Chain 2 Iteration:  101 / 1100 [  9%]  (Sampling) 
## Chain 2 Iteration:  200 / 1100 [ 18%]  (Sampling) 
## Chain 2 Iteration:  300 / 1100 [ 27%]  (Sampling) 
## Chain 2 Iteration:  400 / 1100 [ 36%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:    1 / 1100 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 1100 [  9%]  (Warmup) 
## Chain 3 Iteration:  101 / 1100 [  9%]  (Sampling) 
## Chain 3 Iteration:  200 / 1100 [ 18%]  (Sampling) 
## Chain 3 Iteration:  300 / 1100 [ 27%]  (Sampling) 
## Chain 3 Iteration:  400 / 1100 [ 36%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:    1 / 1100 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 1100 [  9%]  (Warmup) 
## Chain 4 Iteration:  101 / 1100 [  9%]  (Sampling) 
## Chain 4 Iteration:  200 / 1100 [ 18%]  (Sampling) 
## Chain 4 Iteration:  300 / 1100 [ 27%]  (Sampling) 
## Chain 1 Iteration:  500 / 1100 [ 45%]  (Sampling) 
## Chain 1 Iteration:  600 / 1100 [ 54%]  (Sampling) 
## Chain 1 Iteration:  700 / 1100 [ 63%]  (Sampling) 
## Chain 1 Iteration:  800 / 1100 [ 72%]  (Sampling) 
## Chain 1 Iteration:  900 / 1100 [ 81%]  (Sampling) 
## Chain 2 Iteration:  500 / 1100 [ 45%]  (Sampling) 
## Chain 2 Iteration:  600 / 1100 [ 54%]  (Sampling) 
## Chain 2 Iteration:  700 / 1100 [ 63%]  (Sampling) 
## Chain 2 Iteration:  800 / 1100 [ 72%]  (Sampling) 
## Chain 2 Iteration:  900 / 1100 [ 81%]  (Sampling) 
## Chain 3 Iteration:  500 / 1100 [ 45%]  (Sampling) 
## Chain 3 Iteration:  600 / 1100 [ 54%]  (Sampling) 
## Chain 3 Iteration:  700 / 1100 [ 63%]  (Sampling) 
## Chain 3 Iteration:  800 / 1100 [ 72%]  (Sampling) 
## Chain 4 Iteration:  400 / 1100 [ 36%]  (Sampling) 
## Chain 4 Iteration:  500 / 1100 [ 45%]  (Sampling) 
## Chain 4 Iteration:  600 / 1100 [ 54%]  (Sampling) 
## Chain 4 Iteration:  700 / 1100 [ 63%]  (Sampling) 
## Chain 4 Iteration:  800 / 1100 [ 72%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1100 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1100 / 1100 [100%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1100 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1100 / 1100 [100%]  (Sampling) 
## Chain 3 Iteration:  900 / 1100 [ 81%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1100 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1100 / 1100 [100%]  (Sampling) 
## Chain 4 Iteration:  900 / 1100 [ 81%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1100 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1100 / 1100 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 1500 [  6%]  (Warmup) 
## Chain 1 Iteration:  200 / 1500 [ 13%]  (Warmup) 
## Chain 1 Iteration:  300 / 1500 [ 20%]  (Warmup) 
## Chain 1 Iteration:  400 / 1500 [ 26%]  (Warmup) 
## Chain 1 Iteration:  500 / 1500 [ 33%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration:  600 / 1500 [ 40%]  (Sampling) 
## Chain 1 Iteration:  700 / 1500 [ 46%]  (Sampling) 
## Chain 1 Iteration:  800 / 1500 [ 53%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 1500 [  6%]  (Warmup) 
## Chain 2 Iteration:  200 / 1500 [ 13%]  (Warmup) 
## Chain 2 Iteration:  300 / 1500 [ 20%]  (Warmup) 
## Chain 2 Iteration:  400 / 1500 [ 26%]  (Warmup) 
## Chain 2 Iteration:  500 / 1500 [ 33%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:  600 / 1500 [ 40%]  (Sampling) 
## Chain 2 Iteration:  700 / 1500 [ 46%]  (Sampling) 
## Chain 2 Iteration:  800 / 1500 [ 53%]  (Sampling) 
## Chain 2 Iteration:  900 / 1500 [ 60%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 1500 [  6%]  (Warmup) 
## Chain 3 Iteration:  200 / 1500 [ 13%]  (Warmup) 
## Chain 3 Iteration:  300 / 1500 [ 20%]  (Warmup) 
## Chain 3 Iteration:  400 / 1500 [ 26%]  (Warmup) 
## Chain 3 Iteration:  500 / 1500 [ 33%]  (Warmup) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  600 / 1500 [ 40%]  (Sampling) 
## Chain 3 Iteration:  700 / 1500 [ 46%]  (Sampling) 
## Chain 3 Iteration:  800 / 1500 [ 53%]  (Sampling) 
## Chain 3 Iteration:  900 / 1500 [ 60%]  (Sampling) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 1500 [  6%]  (Warmup) 
## Chain 4 Iteration:  200 / 1500 [ 13%]  (Warmup) 
## Chain 4 Iteration:  300 / 1500 [ 20%]  (Warmup) 
## Chain 4 Iteration:  400 / 1500 [ 26%]  (Warmup) 
## Chain 4 Iteration:  500 / 1500 [ 33%]  (Warmup) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  600 / 1500 [ 40%]  (Sampling) 
## Chain 4 Iteration:  700 / 1500 [ 46%]  (Sampling) 
## Chain 4 Iteration:  800 / 1500 [ 53%]  (Sampling) 
## Chain 4 Iteration:  900 / 1500 [ 60%]  (Sampling) 
## Chain 1 Iteration:  900 / 1500 [ 60%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1500 [ 66%]  (Sampling) 
## Chain 1 Iteration: 1100 / 1500 [ 73%]  (Sampling) 
## Chain 1 Iteration: 1200 / 1500 [ 80%]  (Sampling) 
## Chain 1 Iteration: 1300 / 1500 [ 86%]  (Sampling) 
## Chain 1 Iteration: 1400 / 1500 [ 93%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1500 [ 66%]  (Sampling) 
## Chain 2 Iteration: 1100 / 1500 [ 73%]  (Sampling) 
## Chain 2 Iteration: 1200 / 1500 [ 80%]  (Sampling) 
## Chain 2 Iteration: 1300 / 1500 [ 86%]  (Sampling) 
## Chain 2 Iteration: 1400 / 1500 [ 93%]  (Sampling) 
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1500 [ 66%]  (Sampling) 
## Chain 3 Iteration: 1100 / 1500 [ 73%]  (Sampling) 
## Chain 3 Iteration: 1200 / 1500 [ 80%]  (Sampling) 
## Chain 3 Iteration: 1300 / 1500 [ 86%]  (Sampling) 
## Chain 3 Iteration: 1400 / 1500 [ 93%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1500 [ 66%]  (Sampling) 
## Chain 4 Iteration: 1100 / 1500 [ 73%]  (Sampling) 
## Chain 4 Iteration: 1200 / 1500 [ 80%]  (Sampling) 
## Chain 4 Iteration: 1300 / 1500 [ 86%]  (Sampling) 
## Chain 4 Iteration: 1400 / 1500 [ 93%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 0.3 seconds.
## Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.
## 
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 3000 [  3%]  (Warmup) 
## Chain 1 Iteration:  200 / 3000 [  6%]  (Warmup) 
## Chain 1 Iteration:  300 / 3000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  400 / 3000 [ 13%]  (Warmup) 
## Chain 1 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 1 Iteration:  600 / 3000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  700 / 3000 [ 23%]  (Warmup) 
## Chain 1 Iteration:  800 / 3000 [ 26%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 3000 [  3%]  (Warmup) 
## Chain 2 Iteration:  200 / 3000 [  6%]  (Warmup) 
## Chain 2 Iteration:  300 / 3000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  400 / 3000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 2 Iteration:  600 / 3000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  700 / 3000 [ 23%]  (Warmup) 
## Chain 2 Iteration:  800 / 3000 [ 26%]  (Warmup) 
## Chain 2 Iteration:  900 / 3000 [ 30%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 3000 [  3%]  (Warmup) 
## Chain 3 Iteration:  200 / 3000 [  6%]  (Warmup) 
## Chain 3 Iteration:  300 / 3000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  400 / 3000 [ 13%]  (Warmup) 
## Chain 3 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 3 Iteration:  600 / 3000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  700 / 3000 [ 23%]  (Warmup) 
## Chain 3 Iteration:  800 / 3000 [ 26%]  (Warmup) 
## Chain 3 Iteration:  900 / 3000 [ 30%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 3000 [  3%]  (Warmup) 
## Chain 4 Iteration:  200 / 3000 [  6%]  (Warmup) 
## Chain 4 Iteration:  300 / 3000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  400 / 3000 [ 13%]  (Warmup) 
## Chain 4 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 4 Iteration:  600 / 3000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  700 / 3000 [ 23%]  (Warmup) 
## Chain 4 Iteration:  800 / 3000 [ 26%]  (Warmup) 
## Chain 4 Iteration:  900 / 3000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  900 / 3000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1100 / 3000 [ 36%]  (Warmup) 
## Chain 1 Iteration: 1200 / 3000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 1300 / 3000 [ 43%]  (Warmup) 
## Chain 1 Iteration: 1400 / 3000 [ 46%]  (Warmup) 
## Chain 1 Iteration: 1500 / 3000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1600 / 3000 [ 53%]  (Warmup) 
## Chain 1 Iteration: 1700 / 3000 [ 56%]  (Warmup) 
## Chain 1 Iteration: 1800 / 3000 [ 60%]  (Warmup) 
## Chain 1 Iteration: 1900 / 3000 [ 63%]  (Warmup) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Warmup) 
## Chain 1 Iteration: 2001 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 2100 / 3000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1100 / 3000 [ 36%]  (Warmup) 
## Chain 2 Iteration: 1200 / 3000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 1300 / 3000 [ 43%]  (Warmup) 
## Chain 2 Iteration: 1400 / 3000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 1500 / 3000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1600 / 3000 [ 53%]  (Warmup) 
## Chain 2 Iteration: 1700 / 3000 [ 56%]  (Warmup) 
## Chain 2 Iteration: 1800 / 3000 [ 60%]  (Warmup) 
## Chain 2 Iteration: 1900 / 3000 [ 63%]  (Warmup) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Warmup) 
## Chain 2 Iteration: 2001 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1100 / 3000 [ 36%]  (Warmup) 
## Chain 3 Iteration: 1200 / 3000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 1300 / 3000 [ 43%]  (Warmup) 
## Chain 3 Iteration: 1400 / 3000 [ 46%]  (Warmup) 
## Chain 3 Iteration: 1500 / 3000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1600 / 3000 [ 53%]  (Warmup) 
## Chain 3 Iteration: 1700 / 3000 [ 56%]  (Warmup) 
## Chain 3 Iteration: 1800 / 3000 [ 60%]  (Warmup) 
## Chain 3 Iteration: 1900 / 3000 [ 63%]  (Warmup) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1100 / 3000 [ 36%]  (Warmup) 
## Chain 4 Iteration: 1200 / 3000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 1300 / 3000 [ 43%]  (Warmup) 
## Chain 4 Iteration: 1400 / 3000 [ 46%]  (Warmup) 
## Chain 4 Iteration: 1500 / 3000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1600 / 3000 [ 53%]  (Warmup) 
## Chain 4 Iteration: 1700 / 3000 [ 56%]  (Warmup) 
## Chain 4 Iteration: 1800 / 3000 [ 60%]  (Warmup) 
## Chain 4 Iteration: 1900 / 3000 [ 63%]  (Warmup) 
## Chain 1 Iteration: 2200 / 3000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 2300 / 3000 [ 76%]  (Sampling) 
## Chain 1 Iteration: 2400 / 3000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 2600 / 3000 [ 86%]  (Sampling) 
## Chain 1 Iteration: 2700 / 3000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 2800 / 3000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 2900 / 3000 [ 96%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 Iteration: 2100 / 3000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 2200 / 3000 [ 73%]  (Sampling) 
## Chain 2 Iteration: 2300 / 3000 [ 76%]  (Sampling) 
## Chain 2 Iteration: 2400 / 3000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 2600 / 3000 [ 86%]  (Sampling) 
## Chain 2 Iteration: 2700 / 3000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 2800 / 3000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 2900 / 3000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Warmup) 
## Chain 3 Iteration: 2001 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2100 / 3000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 2200 / 3000 [ 73%]  (Sampling) 
## Chain 3 Iteration: 2300 / 3000 [ 76%]  (Sampling) 
## Chain 3 Iteration: 2400 / 3000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 3 Iteration: 2600 / 3000 [ 86%]  (Sampling) 
## Chain 3 Iteration: 2700 / 3000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 2800 / 3000 [ 93%]  (Sampling) 
## Chain 3 Iteration: 2900 / 3000 [ 96%]  (Sampling) 
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Warmup) 
## Chain 4 Iteration: 2001 / 3000 [ 66%]  (Sampling) 
## Chain 4 Iteration: 2100 / 3000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 2200 / 3000 [ 73%]  (Sampling) 
## Chain 4 Iteration: 2300 / 3000 [ 76%]  (Sampling) 
## Chain 4 Iteration: 2400 / 3000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 4 Iteration: 2600 / 3000 [ 86%]  (Sampling) 
## Chain 4 Iteration: 2700 / 3000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 2800 / 3000 [ 93%]  (Sampling) 
## Chain 4 Iteration: 2900 / 3000 [ 96%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.5 seconds.
colnames(warmup_sim) <- rownames(precis(original_model, 2))
rownames(warmup_sim) <- sim_list
warmup_sim
##             a[1]        a[2]        b[1]        b[2]       sigma
## 1       2.001001         NaN    2.001001    2.001001    2.001001
## 5       2.024715    2.120734    2.194452    2.015201    2.007198
## 10   1239.713473 1112.664311  246.061062  374.951055  469.920667
## 50   4045.097730 4992.377166 1651.954197 2587.207250 2524.558025
## 100  4975.571070 4629.199789 1887.421621 2668.892385 2636.906637
## 500  5827.654909 5770.849542 4588.530600 4887.846566 5490.831282
## 1000 5075.463679 5475.069574 5359.111469 5711.854133 4455.852923
## 2000 4983.762642 5271.294270 5234.484522 4553.822670 5245.791339
# As can be seen on table warmup_sim, as warmup increases, the number of n_effective stabilizes. Based on the simulation, Once we reach warmup = 500, the fluctiations in n_eff reduces. Therefore, a wamup of 500 is sufficient.

8-4. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.

mp <- ulam(
 alist(
   a ~ dnorm(0,1),
   b ~ dcauchy(0,1)
 ), data=list(y=1) , chains=1 )
## Running MCMC with 1 chain, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.

Compare the samples for the parameters a and b. Plot the trace plots. Can you explain the different trace plots? If you are unfamiliar with the Cauchy distribution, you should look it up. The key feature to attend to is that it has no expected value. Can you connect this fact to the trace plot?

# Check result
precis(mp, 2)
##          mean         sd      5.5%    94.5%     n_eff    Rhat4
## a -0.01931477  0.9484099  -1.54985 1.471257  70.81822 1.000411
## b -0.85504406 11.8414154 -11.14322 6.581033 249.18706 1.001899
# Compare parameters a and b
normal <- rnorm(1e5, 0, 1)
cauchy <- rcauchy(1e5, 0, 1)

# Create dataframes
normal_df <- with(density(normal), data.frame(x,y))
cauchy_df <- with(density(cauchy), data.frame(x,y))

# Plot parameters
colors <- c("a" = "red")
param_plot_a <- ggplot() +
  geom_line(data=normal_df, aes(x=x, y=y, color="a"), size=1.5) +
  ggtitle("Parameters a and b Comparison") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Estimate", y="", color="Parameter") +
  scale_color_manual(values = colors)

colors <- c("b" = "blue")
param_plot_b <- ggplot() +
  geom_line(data=cauchy_df, aes(x=x, y=y, color="b"), size=1.5) +
  ggtitle("Parameters a and b Comparison") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Estimate", y="", color="Parameter") +
  scale_color_manual(values = colors)
  
param_plot_a

param_plot_b

# Comparing parameters a and b, we can see that parameter a is normal, while parameter b is Cauchy.This means that while both parameters is center around 0, the normal posterior has a more gradual increase/decrease in the bell curve in comparison to the Cauchy posterior

# Trace plots
traceplot(mp)

# The traceplot for a looks like it's evenly distributed across the values, so this looks good. However, the traceplot for b showed that the chain may diverge to some extreme values at certain time. This somewhat mirrors the Cauchy distribution, which does not have an expecte value, so should be acceptable.

8-5. Recall the divorce rate example from Chapter 5. Repeat that analysis, using ulam this time, fitting models m5.1, m5.2, and m5.3. Use compare to compare the models on the basis of WAIC or PSIS. To use WAIC or PSIS with ulam, you need add the argument log_log=TRUE. Explain the model comparison results.

data(WaffleDivorce)
d <- WaffleDivorce

# Standardize variables
d$A <- scale(d$MedianAgeMarriage)
d$D <- scale(d$Divorce)
d$M <- scale(d$Marriage)

d_slim <- list(M=d$M, D=d$D, A=d$A)

# Model m5.1
m5.1 <- ulam(
    alist(
        D ~ dnorm(mu , sigma) ,
        mu <- a + bA * A ,
        a ~ dnorm(0 , 0.2) ,
        bA ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
# Model m5.2
m5.2 <- ulam(
    alist(
        D ~ dnorm(mu, sigma) ,
        mu <- a + bM * M ,
        a ~ dnorm(0 , 0.2) ,
        bM ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
# Model m5.3
m5.3 <- ulam(
    alist(
        D ~ dnorm(mu, sigma) ,
        mu <- a + bM*M + bA*A ,
        a ~ dnorm(0 , 0.2) ,
        bM ~ dnorm(0 , 0.5) ,
        bA ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
## Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
## Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
## Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
# WAIC
set.seed(42)
compare(m5.1 , m5.2 , m5.3 , func=WAIC)
##          WAIC        SE     dWAIC       dSE    pWAIC      weight
## m5.1 125.8304 12.688194  0.000000        NA 3.668605 0.701545586
## m5.3 127.5465 12.920491  1.716018 0.8437367 4.821620 0.297459130
## m5.2 138.9465  9.675527 13.116026 9.3106015 2.751297 0.000995284
# Based on the comparison, we see that WAIC for model 5.1 is lower, therefore preferred. Model 5.1 only uses the median age at marriage as a predictor. WAIC for model 5.3 is almost the same as the WAIC for model 5.1, and this model includes both median age at marriage and marriage rate as predictors. The least performant model is model 5.2.