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?

data(rugged)
df <- rugged

df$log_gdp <- log(df$rgdppc_2000)
df1 <- df[ complete.cases(df$rgdppc_2000) , ]

df1$log_gdp_std <- df1$log_gdp/ mean(df1$log_gdp)
df1$rugged_std<- df1$rugged/max(df1$rugged)

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

dat_slim <- list(
log_gdp_std = df1$log_gdp_std,
rugged_std = df1$rugged_std,
cid = as.integer(df1$cid )
)
#uniform(0,1)
m1 <- 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 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 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 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: 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 2 finished in 0.1 seconds.
## 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 3 finished in 0.1 seconds.
## 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 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
#exponential(1)
m2 <- 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 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 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 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 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.2 seconds.
## 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 2 finished in 0.1 seconds.
## 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 3 finished in 0.2 seconds.
## 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 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
Plot_df <- data.frame(
  Posteriors = c(
    extract.samples(m1, n = 1000)$sigma,
    extract.samples(m2, n = 1000)$sigma
  ),
  Name = rep(c("Uni","Exp"), each = 1000),
  Model = rep(c("m1", "m2"), each = 1000)
)
ggplot(Plot_df, aes(y = Model, x = Posteriors))+
  stat_halfeye()

#The different priors do not have much influence on posterior distributions of sigma.

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?

m3 <- 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 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 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: 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: 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: 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 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 0.3 seconds.
## Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (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.3 seconds.
## Chain 2 finished in 0.3 seconds.
## 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.
precis(m3, depth = 2)
##             mean          sd        5.5%      94.5%    n_eff     Rhat4
## a[1]  0.88648754 0.015766443 0.861581295 0.91250462 1619.763 1.0002958
## a[2]  1.04798303 0.009953436 1.031757250 1.06336110 1885.610 0.9993342
## b[1]  0.14995758 0.069836406 0.042989104 0.26845495 1326.162 1.0007094
## b[2]  0.01844249 0.017093107 0.001366559 0.05185391 1756.127 0.9997642
## sigma 0.11386486 0.006183228 0.104411140 0.12386783 1746.403 0.9997244

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 warmup is enough?

#Try five warmup values of 5, 10, 100, 1000, 10000, the n_eff does not change much when warmup reaches around 100. So a warmup of 100 is enough. 
warmup=5
warmup=10
warmup=100
warmup=1000
warmup=10000
m4 <- ulam(m1, chains = 4, cores = 4, iter = 1000 + warmup, warmup =warmup)
## Running MCMC with 4 parallel chains, with 1 thread(s) per chain...
## 
## Chain 1 Iteration:     1 / 11000 [  0%]  (Warmup) 
## Chain 1 Iteration:   100 / 11000 [  0%]  (Warmup) 
## Chain 1 Iteration:   200 / 11000 [  1%]  (Warmup) 
## Chain 1 Iteration:   300 / 11000 [  2%]  (Warmup) 
## Chain 1 Iteration:   400 / 11000 [  3%]  (Warmup) 
## Chain 1 Iteration:   500 / 11000 [  4%]  (Warmup) 
## Chain 2 Iteration:     1 / 11000 [  0%]  (Warmup) 
## Chain 2 Iteration:   100 / 11000 [  0%]  (Warmup) 
## Chain 2 Iteration:   200 / 11000 [  1%]  (Warmup) 
## Chain 2 Iteration:   300 / 11000 [  2%]  (Warmup) 
## Chain 2 Iteration:   400 / 11000 [  3%]  (Warmup) 
## Chain 2 Iteration:   500 / 11000 [  4%]  (Warmup) 
## Chain 3 Iteration:     1 / 11000 [  0%]  (Warmup) 
## Chain 3 Iteration:   100 / 11000 [  0%]  (Warmup) 
## Chain 3 Iteration:   200 / 11000 [  1%]  (Warmup) 
## Chain 3 Iteration:   300 / 11000 [  2%]  (Warmup) 
## Chain 4 Iteration:     1 / 11000 [  0%]  (Warmup) 
## Chain 4 Iteration:   100 / 11000 [  0%]  (Warmup) 
## Chain 4 Iteration:   200 / 11000 [  1%]  (Warmup) 
## Chain 4 Iteration:   300 / 11000 [  2%]  (Warmup) 
## Chain 4 Iteration:   400 / 11000 [  3%]  (Warmup) 
## Chain 1 Iteration:   600 / 11000 [  5%]  (Warmup) 
## Chain 1 Iteration:   700 / 11000 [  6%]  (Warmup) 
## Chain 1 Iteration:   800 / 11000 [  7%]  (Warmup) 
## Chain 1 Iteration:   900 / 11000 [  8%]  (Warmup) 
## Chain 1 Iteration:  1000 / 11000 [  9%]  (Warmup) 
## Chain 1 Iteration:  1100 / 11000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  1200 / 11000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  1300 / 11000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  1400 / 11000 [ 12%]  (Warmup) 
## Chain 1 Iteration:  1500 / 11000 [ 13%]  (Warmup) 
## Chain 1 Iteration:  1600 / 11000 [ 14%]  (Warmup) 
## Chain 1 Iteration:  1700 / 11000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  1800 / 11000 [ 16%]  (Warmup) 
## Chain 1 Iteration:  1900 / 11000 [ 17%]  (Warmup) 
## Chain 1 Iteration:  2000 / 11000 [ 18%]  (Warmup) 
## Chain 1 Iteration:  2100 / 11000 [ 19%]  (Warmup) 
## Chain 1 Iteration:  2200 / 11000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  2300 / 11000 [ 20%]  (Warmup) 
## Chain 2 Iteration:   600 / 11000 [  5%]  (Warmup) 
## Chain 2 Iteration:   700 / 11000 [  6%]  (Warmup) 
## Chain 2 Iteration:   800 / 11000 [  7%]  (Warmup) 
## Chain 2 Iteration:   900 / 11000 [  8%]  (Warmup) 
## Chain 2 Iteration:  1000 / 11000 [  9%]  (Warmup) 
## Chain 2 Iteration:  1100 / 11000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  1200 / 11000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  1300 / 11000 [ 11%]  (Warmup) 
## Chain 2 Iteration:  1400 / 11000 [ 12%]  (Warmup) 
## Chain 2 Iteration:  1500 / 11000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  1600 / 11000 [ 14%]  (Warmup) 
## Chain 2 Iteration:  1700 / 11000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  1800 / 11000 [ 16%]  (Warmup) 
## Chain 2 Iteration:  1900 / 11000 [ 17%]  (Warmup) 
## Chain 2 Iteration:  2000 / 11000 [ 18%]  (Warmup) 
## Chain 2 Iteration:  2100 / 11000 [ 19%]  (Warmup) 
## Chain 2 Iteration:  2200 / 11000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  2300 / 11000 [ 20%]  (Warmup) 
## Chain 3 Iteration:   400 / 11000 [  3%]  (Warmup) 
## Chain 3 Iteration:   500 / 11000 [  4%]  (Warmup) 
## Chain 3 Iteration:   600 / 11000 [  5%]  (Warmup) 
## Chain 3 Iteration:   700 / 11000 [  6%]  (Warmup) 
## Chain 3 Iteration:   800 / 11000 [  7%]  (Warmup) 
## Chain 3 Iteration:   900 / 11000 [  8%]  (Warmup) 
## Chain 3 Iteration:  1000 / 11000 [  9%]  (Warmup) 
## Chain 3 Iteration:  1100 / 11000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  1200 / 11000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  1300 / 11000 [ 11%]  (Warmup) 
## Chain 3 Iteration:  1400 / 11000 [ 12%]  (Warmup) 
## Chain 3 Iteration:  1500 / 11000 [ 13%]  (Warmup) 
## Chain 3 Iteration:  1600 / 11000 [ 14%]  (Warmup) 
## Chain 3 Iteration:  1700 / 11000 [ 15%]  (Warmup) 
## Chain 3 Iteration:  1800 / 11000 [ 16%]  (Warmup) 
## Chain 3 Iteration:  1900 / 11000 [ 17%]  (Warmup) 
## Chain 3 Iteration:  2000 / 11000 [ 18%]  (Warmup) 
## Chain 3 Iteration:  2100 / 11000 [ 19%]  (Warmup) 
## Chain 4 Iteration:   500 / 11000 [  4%]  (Warmup) 
## Chain 4 Iteration:   600 / 11000 [  5%]  (Warmup) 
## Chain 4 Iteration:   700 / 11000 [  6%]  (Warmup) 
## Chain 4 Iteration:   800 / 11000 [  7%]  (Warmup) 
## Chain 4 Iteration:   900 / 11000 [  8%]  (Warmup) 
## Chain 4 Iteration:  1000 / 11000 [  9%]  (Warmup) 
## Chain 4 Iteration:  1100 / 11000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  1200 / 11000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  1300 / 11000 [ 11%]  (Warmup) 
## Chain 4 Iteration:  1400 / 11000 [ 12%]  (Warmup) 
## Chain 4 Iteration:  1500 / 11000 [ 13%]  (Warmup) 
## Chain 4 Iteration:  1600 / 11000 [ 14%]  (Warmup) 
## Chain 4 Iteration:  1700 / 11000 [ 15%]  (Warmup) 
## Chain 4 Iteration:  1800 / 11000 [ 16%]  (Warmup) 
## Chain 4 Iteration:  1900 / 11000 [ 17%]  (Warmup) 
## Chain 4 Iteration:  2000 / 11000 [ 18%]  (Warmup) 
## Chain 4 Iteration:  2100 / 11000 [ 19%]  (Warmup) 
## Chain 4 Iteration:  2200 / 11000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  2400 / 11000 [ 21%]  (Warmup) 
## Chain 1 Iteration:  2500 / 11000 [ 22%]  (Warmup) 
## Chain 1 Iteration:  2600 / 11000 [ 23%]  (Warmup) 
## Chain 1 Iteration:  2700 / 11000 [ 24%]  (Warmup) 
## Chain 1 Iteration:  2800 / 11000 [ 25%]  (Warmup) 
## Chain 1 Iteration:  2900 / 11000 [ 26%]  (Warmup) 
## Chain 1 Iteration:  3000 / 11000 [ 27%]  (Warmup) 
## Chain 1 Iteration:  3100 / 11000 [ 28%]  (Warmup) 
## Chain 1 Iteration:  3200 / 11000 [ 29%]  (Warmup) 
## Chain 1 Iteration:  3300 / 11000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  3400 / 11000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  3500 / 11000 [ 31%]  (Warmup) 
## Chain 2 Iteration:  2400 / 11000 [ 21%]  (Warmup) 
## Chain 2 Iteration:  2500 / 11000 [ 22%]  (Warmup) 
## Chain 2 Iteration:  2600 / 11000 [ 23%]  (Warmup) 
## Chain 2 Iteration:  2700 / 11000 [ 24%]  (Warmup) 
## Chain 2 Iteration:  2800 / 11000 [ 25%]  (Warmup) 
## Chain 2 Iteration:  2900 / 11000 [ 26%]  (Warmup) 
## Chain 2 Iteration:  3000 / 11000 [ 27%]  (Warmup) 
## Chain 2 Iteration:  3100 / 11000 [ 28%]  (Warmup) 
## Chain 2 Iteration:  3200 / 11000 [ 29%]  (Warmup) 
## Chain 2 Iteration:  3300 / 11000 [ 30%]  (Warmup) 
## Chain 2 Iteration:  3400 / 11000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  2200 / 11000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  2300 / 11000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  2400 / 11000 [ 21%]  (Warmup) 
## Chain 3 Iteration:  2500 / 11000 [ 22%]  (Warmup) 
## Chain 3 Iteration:  2600 / 11000 [ 23%]  (Warmup) 
## Chain 3 Iteration:  2700 / 11000 [ 24%]  (Warmup) 
## Chain 3 Iteration:  2800 / 11000 [ 25%]  (Warmup) 
## Chain 3 Iteration:  2900 / 11000 [ 26%]  (Warmup) 
## Chain 3 Iteration:  3000 / 11000 [ 27%]  (Warmup) 
## Chain 3 Iteration:  3100 / 11000 [ 28%]  (Warmup) 
## Chain 3 Iteration:  3200 / 11000 [ 29%]  (Warmup) 
## Chain 4 Iteration:  2300 / 11000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  2400 / 11000 [ 21%]  (Warmup) 
## Chain 4 Iteration:  2500 / 11000 [ 22%]  (Warmup) 
## Chain 4 Iteration:  2600 / 11000 [ 23%]  (Warmup) 
## Chain 4 Iteration:  2700 / 11000 [ 24%]  (Warmup) 
## Chain 4 Iteration:  2800 / 11000 [ 25%]  (Warmup) 
## Chain 4 Iteration:  2900 / 11000 [ 26%]  (Warmup) 
## Chain 4 Iteration:  3000 / 11000 [ 27%]  (Warmup) 
## Chain 4 Iteration:  3100 / 11000 [ 28%]  (Warmup) 
## Chain 4 Iteration:  3200 / 11000 [ 29%]  (Warmup) 
## Chain 4 Iteration:  3300 / 11000 [ 30%]  (Warmup) 
## Chain 4 Iteration:  3400 / 11000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  3600 / 11000 [ 32%]  (Warmup) 
## Chain 1 Iteration:  3700 / 11000 [ 33%]  (Warmup) 
## Chain 1 Iteration:  3800 / 11000 [ 34%]  (Warmup) 
## Chain 1 Iteration:  3900 / 11000 [ 35%]  (Warmup) 
## Chain 1 Iteration:  4000 / 11000 [ 36%]  (Warmup) 
## Chain 1 Iteration:  4100 / 11000 [ 37%]  (Warmup) 
## Chain 1 Iteration:  4200 / 11000 [ 38%]  (Warmup) 
## Chain 1 Iteration:  4300 / 11000 [ 39%]  (Warmup) 
## Chain 1 Iteration:  4400 / 11000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  4500 / 11000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  3500 / 11000 [ 31%]  (Warmup) 
## Chain 2 Iteration:  3600 / 11000 [ 32%]  (Warmup) 
## Chain 2 Iteration:  3700 / 11000 [ 33%]  (Warmup) 
## Chain 2 Iteration:  3800 / 11000 [ 34%]  (Warmup) 
## Chain 2 Iteration:  3900 / 11000 [ 35%]  (Warmup) 
## Chain 2 Iteration:  4000 / 11000 [ 36%]  (Warmup) 
## Chain 2 Iteration:  4100 / 11000 [ 37%]  (Warmup) 
## Chain 2 Iteration:  4200 / 11000 [ 38%]  (Warmup) 
## Chain 2 Iteration:  4300 / 11000 [ 39%]  (Warmup) 
## Chain 2 Iteration:  4400 / 11000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  4500 / 11000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  3300 / 11000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  3400 / 11000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  3500 / 11000 [ 31%]  (Warmup) 
## Chain 3 Iteration:  3600 / 11000 [ 32%]  (Warmup) 
## Chain 3 Iteration:  3700 / 11000 [ 33%]  (Warmup) 
## Chain 3 Iteration:  3800 / 11000 [ 34%]  (Warmup) 
## Chain 3 Iteration:  3900 / 11000 [ 35%]  (Warmup) 
## Chain 3 Iteration:  4000 / 11000 [ 36%]  (Warmup) 
## Chain 3 Iteration:  4100 / 11000 [ 37%]  (Warmup) 
## Chain 3 Iteration:  4200 / 11000 [ 38%]  (Warmup) 
## Chain 3 Iteration:  4300 / 11000 [ 39%]  (Warmup) 
## Chain 4 Iteration:  3500 / 11000 [ 31%]  (Warmup) 
## Chain 4 Iteration:  3600 / 11000 [ 32%]  (Warmup) 
## Chain 4 Iteration:  3700 / 11000 [ 33%]  (Warmup) 
## Chain 4 Iteration:  3800 / 11000 [ 34%]  (Warmup) 
## Chain 4 Iteration:  3900 / 11000 [ 35%]  (Warmup) 
## Chain 4 Iteration:  4000 / 11000 [ 36%]  (Warmup) 
## Chain 4 Iteration:  4100 / 11000 [ 37%]  (Warmup) 
## Chain 4 Iteration:  4200 / 11000 [ 38%]  (Warmup) 
## Chain 4 Iteration:  4300 / 11000 [ 39%]  (Warmup) 
## Chain 4 Iteration:  4400 / 11000 [ 40%]  (Warmup) 
## Chain 4 Iteration:  4500 / 11000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  4600 / 11000 [ 41%]  (Warmup) 
## Chain 1 Iteration:  4700 / 11000 [ 42%]  (Warmup) 
## Chain 1 Iteration:  4800 / 11000 [ 43%]  (Warmup) 
## Chain 1 Iteration:  4900 / 11000 [ 44%]  (Warmup) 
## Chain 1 Iteration:  5000 / 11000 [ 45%]  (Warmup) 
## Chain 1 Iteration:  5100 / 11000 [ 46%]  (Warmup) 
## Chain 1 Iteration:  5200 / 11000 [ 47%]  (Warmup) 
## Chain 1 Iteration:  5300 / 11000 [ 48%]  (Warmup) 
## Chain 1 Iteration:  5400 / 11000 [ 49%]  (Warmup) 
## Chain 1 Iteration:  5500 / 11000 [ 50%]  (Warmup) 
## Chain 1 Iteration:  5600 / 11000 [ 50%]  (Warmup) 
## Chain 2 Iteration:  4600 / 11000 [ 41%]  (Warmup) 
## Chain 2 Iteration:  4700 / 11000 [ 42%]  (Warmup) 
## Chain 2 Iteration:  4800 / 11000 [ 43%]  (Warmup) 
## Chain 2 Iteration:  4900 / 11000 [ 44%]  (Warmup) 
## Chain 2 Iteration:  5000 / 11000 [ 45%]  (Warmup) 
## Chain 2 Iteration:  5100 / 11000 [ 46%]  (Warmup) 
## Chain 2 Iteration:  5200 / 11000 [ 47%]  (Warmup) 
## Chain 2 Iteration:  5300 / 11000 [ 48%]  (Warmup) 
## Chain 2 Iteration:  5400 / 11000 [ 49%]  (Warmup) 
## Chain 2 Iteration:  5500 / 11000 [ 50%]  (Warmup) 
## Chain 3 Iteration:  4400 / 11000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  4500 / 11000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  4600 / 11000 [ 41%]  (Warmup) 
## Chain 3 Iteration:  4700 / 11000 [ 42%]  (Warmup) 
## Chain 3 Iteration:  4800 / 11000 [ 43%]  (Warmup) 
## Chain 3 Iteration:  4900 / 11000 [ 44%]  (Warmup) 
## Chain 3 Iteration:  5000 / 11000 [ 45%]  (Warmup) 
## Chain 3 Iteration:  5100 / 11000 [ 46%]  (Warmup) 
## Chain 3 Iteration:  5200 / 11000 [ 47%]  (Warmup) 
## Chain 3 Iteration:  5300 / 11000 [ 48%]  (Warmup) 
## Chain 4 Iteration:  4600 / 11000 [ 41%]  (Warmup) 
## Chain 4 Iteration:  4700 / 11000 [ 42%]  (Warmup) 
## Chain 4 Iteration:  4800 / 11000 [ 43%]  (Warmup) 
## Chain 4 Iteration:  4900 / 11000 [ 44%]  (Warmup) 
## Chain 4 Iteration:  5000 / 11000 [ 45%]  (Warmup) 
## Chain 4 Iteration:  5100 / 11000 [ 46%]  (Warmup) 
## Chain 4 Iteration:  5200 / 11000 [ 47%]  (Warmup) 
## Chain 4 Iteration:  5300 / 11000 [ 48%]  (Warmup) 
## Chain 4 Iteration:  5400 / 11000 [ 49%]  (Warmup) 
## Chain 4 Iteration:  5500 / 11000 [ 50%]  (Warmup) 
## Chain 1 Iteration:  5700 / 11000 [ 51%]  (Warmup) 
## Chain 1 Iteration:  5800 / 11000 [ 52%]  (Warmup) 
## Chain 1 Iteration:  5900 / 11000 [ 53%]  (Warmup) 
## Chain 1 Iteration:  6000 / 11000 [ 54%]  (Warmup) 
## Chain 1 Iteration:  6100 / 11000 [ 55%]  (Warmup) 
## Chain 1 Iteration:  6200 / 11000 [ 56%]  (Warmup) 
## Chain 1 Iteration:  6300 / 11000 [ 57%]  (Warmup) 
## Chain 1 Iteration:  6400 / 11000 [ 58%]  (Warmup) 
## Chain 1 Iteration:  6500 / 11000 [ 59%]  (Warmup) 
## Chain 1 Iteration:  6600 / 11000 [ 60%]  (Warmup) 
## Chain 1 Iteration:  6700 / 11000 [ 60%]  (Warmup) 
## Chain 2 Iteration:  5600 / 11000 [ 50%]  (Warmup) 
## Chain 2 Iteration:  5700 / 11000 [ 51%]  (Warmup) 
## Chain 2 Iteration:  5800 / 11000 [ 52%]  (Warmup) 
## Chain 2 Iteration:  5900 / 11000 [ 53%]  (Warmup) 
## Chain 2 Iteration:  6000 / 11000 [ 54%]  (Warmup) 
## Chain 2 Iteration:  6100 / 11000 [ 55%]  (Warmup) 
## Chain 2 Iteration:  6200 / 11000 [ 56%]  (Warmup) 
## Chain 2 Iteration:  6300 / 11000 [ 57%]  (Warmup) 
## Chain 2 Iteration:  6400 / 11000 [ 58%]  (Warmup) 
## Chain 2 Iteration:  6500 / 11000 [ 59%]  (Warmup) 
## Chain 2 Iteration:  6600 / 11000 [ 60%]  (Warmup) 
## Chain 3 Iteration:  5400 / 11000 [ 49%]  (Warmup) 
## Chain 3 Iteration:  5500 / 11000 [ 50%]  (Warmup) 
## Chain 3 Iteration:  5600 / 11000 [ 50%]  (Warmup) 
## Chain 3 Iteration:  5700 / 11000 [ 51%]  (Warmup) 
## Chain 3 Iteration:  5800 / 11000 [ 52%]  (Warmup) 
## Chain 3 Iteration:  5900 / 11000 [ 53%]  (Warmup) 
## Chain 3 Iteration:  6000 / 11000 [ 54%]  (Warmup) 
## Chain 3 Iteration:  6100 / 11000 [ 55%]  (Warmup) 
## Chain 3 Iteration:  6200 / 11000 [ 56%]  (Warmup) 
## Chain 3 Iteration:  6300 / 11000 [ 57%]  (Warmup) 
## Chain 3 Iteration:  6400 / 11000 [ 58%]  (Warmup) 
## Chain 4 Iteration:  5600 / 11000 [ 50%]  (Warmup) 
## Chain 4 Iteration:  5700 / 11000 [ 51%]  (Warmup) 
## Chain 4 Iteration:  5800 / 11000 [ 52%]  (Warmup) 
## Chain 4 Iteration:  5900 / 11000 [ 53%]  (Warmup) 
## Chain 4 Iteration:  6000 / 11000 [ 54%]  (Warmup) 
## Chain 4 Iteration:  6100 / 11000 [ 55%]  (Warmup) 
## Chain 4 Iteration:  6200 / 11000 [ 56%]  (Warmup) 
## Chain 4 Iteration:  6300 / 11000 [ 57%]  (Warmup) 
## Chain 4 Iteration:  6400 / 11000 [ 58%]  (Warmup) 
## Chain 4 Iteration:  6500 / 11000 [ 59%]  (Warmup) 
## Chain 4 Iteration:  6600 / 11000 [ 60%]  (Warmup) 
## Chain 1 Iteration:  6800 / 11000 [ 61%]  (Warmup) 
## Chain 1 Iteration:  6900 / 11000 [ 62%]  (Warmup) 
## Chain 1 Iteration:  7000 / 11000 [ 63%]  (Warmup) 
## Chain 1 Iteration:  7100 / 11000 [ 64%]  (Warmup) 
## Chain 1 Iteration:  7200 / 11000 [ 65%]  (Warmup) 
## Chain 1 Iteration:  7300 / 11000 [ 66%]  (Warmup) 
## Chain 1 Iteration:  7400 / 11000 [ 67%]  (Warmup) 
## Chain 1 Iteration:  7500 / 11000 [ 68%]  (Warmup) 
## Chain 1 Iteration:  7600 / 11000 [ 69%]  (Warmup) 
## Chain 1 Iteration:  7700 / 11000 [ 70%]  (Warmup) 
## Chain 1 Iteration:  7800 / 11000 [ 70%]  (Warmup) 
## Chain 2 Iteration:  6700 / 11000 [ 60%]  (Warmup) 
## Chain 2 Iteration:  6800 / 11000 [ 61%]  (Warmup) 
## Chain 2 Iteration:  6900 / 11000 [ 62%]  (Warmup) 
## Chain 2 Iteration:  7000 / 11000 [ 63%]  (Warmup) 
## Chain 2 Iteration:  7100 / 11000 [ 64%]  (Warmup) 
## Chain 2 Iteration:  7200 / 11000 [ 65%]  (Warmup) 
## Chain 2 Iteration:  7300 / 11000 [ 66%]  (Warmup) 
## Chain 2 Iteration:  7400 / 11000 [ 67%]  (Warmup) 
## Chain 2 Iteration:  7500 / 11000 [ 68%]  (Warmup) 
## Chain 2 Iteration:  7600 / 11000 [ 69%]  (Warmup) 
## Chain 2 Iteration:  7700 / 11000 [ 70%]  (Warmup) 
## Chain 3 Iteration:  6500 / 11000 [ 59%]  (Warmup) 
## Chain 3 Iteration:  6600 / 11000 [ 60%]  (Warmup) 
## Chain 3 Iteration:  6700 / 11000 [ 60%]  (Warmup) 
## Chain 3 Iteration:  6800 / 11000 [ 61%]  (Warmup) 
## Chain 3 Iteration:  6900 / 11000 [ 62%]  (Warmup) 
## Chain 3 Iteration:  7000 / 11000 [ 63%]  (Warmup) 
## Chain 3 Iteration:  7100 / 11000 [ 64%]  (Warmup) 
## Chain 3 Iteration:  7200 / 11000 [ 65%]  (Warmup) 
## Chain 3 Iteration:  7300 / 11000 [ 66%]  (Warmup) 
## Chain 3 Iteration:  7400 / 11000 [ 67%]  (Warmup) 
## Chain 3 Iteration:  7500 / 11000 [ 68%]  (Warmup) 
## Chain 4 Iteration:  6700 / 11000 [ 60%]  (Warmup) 
## Chain 4 Iteration:  6800 / 11000 [ 61%]  (Warmup) 
## Chain 4 Iteration:  6900 / 11000 [ 62%]  (Warmup) 
## Chain 4 Iteration:  7000 / 11000 [ 63%]  (Warmup) 
## Chain 4 Iteration:  7100 / 11000 [ 64%]  (Warmup) 
## Chain 4 Iteration:  7200 / 11000 [ 65%]  (Warmup) 
## Chain 4 Iteration:  7300 / 11000 [ 66%]  (Warmup) 
## Chain 4 Iteration:  7400 / 11000 [ 67%]  (Warmup) 
## Chain 4 Iteration:  7500 / 11000 [ 68%]  (Warmup) 
## Chain 4 Iteration:  7600 / 11000 [ 69%]  (Warmup) 
## Chain 4 Iteration:  7700 / 11000 [ 70%]  (Warmup) 
## Chain 1 Iteration:  7900 / 11000 [ 71%]  (Warmup) 
## Chain 1 Iteration:  8000 / 11000 [ 72%]  (Warmup) 
## Chain 1 Iteration:  8100 / 11000 [ 73%]  (Warmup) 
## Chain 1 Iteration:  8200 / 11000 [ 74%]  (Warmup) 
## Chain 1 Iteration:  8300 / 11000 [ 75%]  (Warmup) 
## Chain 1 Iteration:  8400 / 11000 [ 76%]  (Warmup) 
## Chain 1 Iteration:  8500 / 11000 [ 77%]  (Warmup) 
## Chain 1 Iteration:  8600 / 11000 [ 78%]  (Warmup) 
## Chain 1 Iteration:  8700 / 11000 [ 79%]  (Warmup) 
## Chain 1 Iteration:  8800 / 11000 [ 80%]  (Warmup) 
## Chain 2 Iteration:  7800 / 11000 [ 70%]  (Warmup) 
## Chain 2 Iteration:  7900 / 11000 [ 71%]  (Warmup) 
## Chain 2 Iteration:  8000 / 11000 [ 72%]  (Warmup) 
## Chain 2 Iteration:  8100 / 11000 [ 73%]  (Warmup) 
## Chain 2 Iteration:  8200 / 11000 [ 74%]  (Warmup) 
## Chain 2 Iteration:  8300 / 11000 [ 75%]  (Warmup) 
## Chain 2 Iteration:  8400 / 11000 [ 76%]  (Warmup) 
## Chain 2 Iteration:  8500 / 11000 [ 77%]  (Warmup) 
## Chain 2 Iteration:  8600 / 11000 [ 78%]  (Warmup) 
## Chain 2 Iteration:  8700 / 11000 [ 79%]  (Warmup) 
## Chain 3 Iteration:  7600 / 11000 [ 69%]  (Warmup) 
## Chain 3 Iteration:  7700 / 11000 [ 70%]  (Warmup) 
## Chain 3 Iteration:  7800 / 11000 [ 70%]  (Warmup) 
## Chain 3 Iteration:  7900 / 11000 [ 71%]  (Warmup) 
## Chain 3 Iteration:  8000 / 11000 [ 72%]  (Warmup) 
## Chain 3 Iteration:  8100 / 11000 [ 73%]  (Warmup) 
## Chain 3 Iteration:  8200 / 11000 [ 74%]  (Warmup) 
## Chain 3 Iteration:  8300 / 11000 [ 75%]  (Warmup) 
## Chain 3 Iteration:  8400 / 11000 [ 76%]  (Warmup) 
## Chain 3 Iteration:  8500 / 11000 [ 77%]  (Warmup) 
## Chain 4 Iteration:  7800 / 11000 [ 70%]  (Warmup) 
## Chain 4 Iteration:  7900 / 11000 [ 71%]  (Warmup) 
## Chain 4 Iteration:  8000 / 11000 [ 72%]  (Warmup) 
## Chain 4 Iteration:  8100 / 11000 [ 73%]  (Warmup) 
## Chain 4 Iteration:  8200 / 11000 [ 74%]  (Warmup) 
## Chain 4 Iteration:  8300 / 11000 [ 75%]  (Warmup) 
## Chain 4 Iteration:  8400 / 11000 [ 76%]  (Warmup) 
## Chain 4 Iteration:  8500 / 11000 [ 77%]  (Warmup) 
## Chain 4 Iteration:  8600 / 11000 [ 78%]  (Warmup) 
## Chain 4 Iteration:  8700 / 11000 [ 79%]  (Warmup) 
## Chain 1 Iteration:  8900 / 11000 [ 80%]  (Warmup) 
## Chain 1 Iteration:  9000 / 11000 [ 81%]  (Warmup) 
## Chain 1 Iteration:  9100 / 11000 [ 82%]  (Warmup) 
## Chain 1 Iteration:  9200 / 11000 [ 83%]  (Warmup) 
## Chain 1 Iteration:  9300 / 11000 [ 84%]  (Warmup) 
## Chain 1 Iteration:  9400 / 11000 [ 85%]  (Warmup) 
## Chain 1 Iteration:  9500 / 11000 [ 86%]  (Warmup) 
## Chain 1 Iteration:  9600 / 11000 [ 87%]  (Warmup) 
## Chain 1 Iteration:  9700 / 11000 [ 88%]  (Warmup) 
## Chain 1 Iteration:  9800 / 11000 [ 89%]  (Warmup) 
## Chain 1 Iteration:  9900 / 11000 [ 90%]  (Warmup) 
## Chain 2 Iteration:  8800 / 11000 [ 80%]  (Warmup) 
## Chain 2 Iteration:  8900 / 11000 [ 80%]  (Warmup) 
## Chain 2 Iteration:  9000 / 11000 [ 81%]  (Warmup) 
## Chain 2 Iteration:  9100 / 11000 [ 82%]  (Warmup) 
## Chain 2 Iteration:  9200 / 11000 [ 83%]  (Warmup) 
## Chain 2 Iteration:  9300 / 11000 [ 84%]  (Warmup) 
## Chain 2 Iteration:  9400 / 11000 [ 85%]  (Warmup) 
## Chain 2 Iteration:  9500 / 11000 [ 86%]  (Warmup) 
## Chain 2 Iteration:  9600 / 11000 [ 87%]  (Warmup) 
## Chain 2 Iteration:  9700 / 11000 [ 88%]  (Warmup) 
## Chain 2 Iteration:  9800 / 11000 [ 89%]  (Warmup) 
## Chain 3 Iteration:  8600 / 11000 [ 78%]  (Warmup) 
## Chain 3 Iteration:  8700 / 11000 [ 79%]  (Warmup) 
## Chain 3 Iteration:  8800 / 11000 [ 80%]  (Warmup) 
## Chain 3 Iteration:  8900 / 11000 [ 80%]  (Warmup) 
## Chain 3 Iteration:  9000 / 11000 [ 81%]  (Warmup) 
## Chain 3 Iteration:  9100 / 11000 [ 82%]  (Warmup) 
## Chain 3 Iteration:  9200 / 11000 [ 83%]  (Warmup) 
## Chain 3 Iteration:  9300 / 11000 [ 84%]  (Warmup) 
## Chain 3 Iteration:  9400 / 11000 [ 85%]  (Warmup) 
## Chain 3 Iteration:  9500 / 11000 [ 86%]  (Warmup) 
## Chain 3 Iteration:  9600 / 11000 [ 87%]  (Warmup) 
## Chain 4 Iteration:  8800 / 11000 [ 80%]  (Warmup) 
## Chain 4 Iteration:  8900 / 11000 [ 80%]  (Warmup) 
## Chain 4 Iteration:  9000 / 11000 [ 81%]  (Warmup) 
## Chain 4 Iteration:  9100 / 11000 [ 82%]  (Warmup) 
## Chain 4 Iteration:  9200 / 11000 [ 83%]  (Warmup) 
## Chain 4 Iteration:  9300 / 11000 [ 84%]  (Warmup) 
## Chain 4 Iteration:  9400 / 11000 [ 85%]  (Warmup) 
## Chain 4 Iteration:  9500 / 11000 [ 86%]  (Warmup) 
## Chain 4 Iteration:  9600 / 11000 [ 87%]  (Warmup) 
## Chain 4 Iteration:  9700 / 11000 [ 88%]  (Warmup) 
## Chain 4 Iteration:  9800 / 11000 [ 89%]  (Warmup) 
## Chain 1 Iteration: 10000 / 11000 [ 90%]  (Warmup) 
## Chain 1 Iteration: 10001 / 11000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 10100 / 11000 [ 91%]  (Sampling) 
## Chain 1 Iteration: 10200 / 11000 [ 92%]  (Sampling) 
## Chain 1 Iteration: 10300 / 11000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 10400 / 11000 [ 94%]  (Sampling) 
## Chain 1 Iteration: 10500 / 11000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 10600 / 11000 [ 96%]  (Sampling) 
## Chain 1 Iteration: 10700 / 11000 [ 97%]  (Sampling) 
## Chain 1 Iteration: 10800 / 11000 [ 98%]  (Sampling) 
## Chain 1 Iteration: 10900 / 11000 [ 99%]  (Sampling) 
## Chain 2 Iteration:  9900 / 11000 [ 90%]  (Warmup) 
## Chain 2 Iteration: 10000 / 11000 [ 90%]  (Warmup) 
## Chain 2 Iteration: 10001 / 11000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 10100 / 11000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 10200 / 11000 [ 92%]  (Sampling) 
## Chain 2 Iteration: 10300 / 11000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 10400 / 11000 [ 94%]  (Sampling) 
## Chain 2 Iteration: 10500 / 11000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 10600 / 11000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 10700 / 11000 [ 97%]  (Sampling) 
## Chain 2 Iteration: 10800 / 11000 [ 98%]  (Sampling) 
## Chain 3 Iteration:  9700 / 11000 [ 88%]  (Warmup) 
## Chain 3 Iteration:  9800 / 11000 [ 89%]  (Warmup) 
## Chain 3 Iteration:  9900 / 11000 [ 90%]  (Warmup) 
## Chain 3 Iteration: 10000 / 11000 [ 90%]  (Warmup) 
## Chain 3 Iteration: 10001 / 11000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 10100 / 11000 [ 91%]  (Sampling) 
## Chain 3 Iteration: 10200 / 11000 [ 92%]  (Sampling) 
## Chain 3 Iteration: 10300 / 11000 [ 93%]  (Sampling) 
## Chain 3 Iteration: 10400 / 11000 [ 94%]  (Sampling) 
## Chain 3 Iteration: 10500 / 11000 [ 95%]  (Sampling) 
## Chain 3 Iteration: 10600 / 11000 [ 96%]  (Sampling) 
## Chain 3 Iteration: 10700 / 11000 [ 97%]  (Sampling) 
## Chain 4 Iteration:  9900 / 11000 [ 90%]  (Warmup) 
## Chain 4 Iteration: 10000 / 11000 [ 90%]  (Warmup) 
## Chain 4 Iteration: 10001 / 11000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 10100 / 11000 [ 91%]  (Sampling) 
## Chain 4 Iteration: 10200 / 11000 [ 92%]  (Sampling) 
## Chain 4 Iteration: 10300 / 11000 [ 93%]  (Sampling) 
## Chain 4 Iteration: 10400 / 11000 [ 94%]  (Sampling) 
## Chain 4 Iteration: 10500 / 11000 [ 95%]  (Sampling) 
## Chain 4 Iteration: 10600 / 11000 [ 96%]  (Sampling) 
## Chain 4 Iteration: 10700 / 11000 [ 97%]  (Sampling) 
## Chain 4 Iteration: 10800 / 11000 [ 98%]  (Sampling) 
## Chain 4 Iteration: 10900 / 11000 [ 99%]  (Sampling) 
## Chain 1 Iteration: 11000 / 11000 [100%]  (Sampling) 
## Chain 1 finished in 1.4 seconds.
## Chain 2 Iteration: 10900 / 11000 [ 99%]  (Sampling) 
## Chain 2 Iteration: 11000 / 11000 [100%]  (Sampling) 
## Chain 2 finished in 1.4 seconds.
## Chain 4 Iteration: 11000 / 11000 [100%]  (Sampling) 
## Chain 4 finished in 1.3 seconds.
## Chain 3 Iteration: 10800 / 11000 [ 98%]  (Sampling) 
## Chain 3 Iteration: 10900 / 11000 [ 99%]  (Sampling) 
## Chain 3 Iteration: 11000 / 11000 [100%]  (Sampling) 
## Chain 3 finished in 1.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.4 seconds.
## Total execution time: 1.6 seconds.
precis(m4, 2)$n_eff
## [1] 5778.502 6268.835 6330.144 5023.706 5460.597

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.0 seconds.
precis(mp , depth=2)
##          mean       sd      5.5%    94.5%    n_eff    Rhat4
## a 0.003269664 1.014250 -1.564439 1.661123 181.7171 1.005592
## b 0.059529999 5.530565 -5.288713 3.877110 202.3918 1.004288
traceplot(mp)
Posteriors <- extract.samples(mp, n = 1000)
dens(Posteriors$a)

dens(Posteriors$b)
#The cauchy distribution has been approximated well, the normal distribution looks strange.

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?

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)
df <- WaffleDivorce
dat_slim <- list(Divorce = df$Divorce,
       Marriage = df$Marriage,
       MedianAgeMarriage = df$MedianAgeMarriage)

m5.1 <- ulam(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bA * MedianAgeMarriage,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  data = dat_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 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 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 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 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: 800 / 1000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## 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 3 finished in 0.1 seconds.
## 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 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
m5.2 <- ulam(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bM * Marriage,
    a ~ dnorm(0, 0.2),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  data = dat_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 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 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 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 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: 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 2 finished in 0.1 seconds.
## 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 3 finished in 0.1 seconds.
## 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 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
m5.3<- ulam(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bA * MedianAgeMarriage + bM * Marriage,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  data = dat_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 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
## Chain 2 Iteration: 200 / 1000 [ 20%]  (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 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: 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 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 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 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.4 seconds.
## Total execution time: 0.5 seconds.
compare(m5.1, m5.2, m5.3, func=WAIC)
##          WAIC        SE    dWAIC      dSE    pWAIC       weight
## m5.3 204.9603 10.890266  0.00000       NA 3.470172 0.9961153448
## m5.2 216.3817 14.657685 11.42147 8.134441 2.939908 0.0032973845
## m5.1 219.8325  9.226623 14.87226 9.713900 1.629097 0.0005872708
compare(m5.1, m5.2, m5.3, func=PSIS)
##          PSIS        SE    dPSIS      dSE    pPSIS       weight
## m5.3 205.0456 11.027829  0.00000       NA 3.512864 0.9960358220
## m5.2 216.4324 14.832891 11.38681 8.133208 2.965270 0.0033547572
## m5.1 219.8437  9.322831 14.79806 9.737065 1.634685 0.0006094208
#From the results, model 5.3 with predictor median age at marriage and marriage has the best performance, while model 5.1 and model 5.2 are having similar performance,