Chapter 9 - Markov Chain Monte Carlo

This chapter 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. Problems are labeled Easy (E), Medium (M), and Hard(H).

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

9E1. Which of the following is a requirement of the simple Metropolis algorithm?

  1. The parameters must be discrete.
  2. The likelihood function must be Gaussian.
  3. The proposal distribution must be symmetric.
# The proposal distribution must be symmetric is a requirement of the simple Metropolis algorithm.

9E2. Gibbs sampling is more efficient than the Metropolis algorithm. How does it achieve this extra efficiency? Are there any limitations to the Gibbs sampling strategy?

# The proposed parameter values distribution can be adjusted to current parameter values. 
# Gibbs sampling:Pairs of priors and likelihoods have solutions for the posterior of an individual parameter.

9E3. Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?

# can't handle discrete parameters, because there is no slope.

9E4. Explain the difference between the effective number of samples, n_eff as calculated by Stan, and the actual number of samples.

# n_eff : estimation of the effective number of independent samples 
# the effective number of samples : estimation of the number of independent samples from the posterior distribution

9E5. Which value should Rhat approach, when a chain is sampling the posterior distribution correctly?

# approach to 1

9E6. Sketch a good trace plot for a Markov chain, one that is effectively sampling from the posterior distribution. What is good about its shape? Then sketch a trace plot for a malfunctioning Markov chain. What about its shape indicates malfunction?

library(rethinking)
data("rugged")
data <- rugged
data$log_gdp <- log(data$rgdppc_2000)
data_1 <- data[complete.cases(data$rgdppc_2000),]
data_1$log_gdp_std <- data_1$log_gdp / mean(data_1$log_gdp)
data_1$rugged_std <- data_1$rugged/max(data_1$rugged)
data_1$cid <- ifelse(data_1$cont_africa==1,1,2)
data_2 <- list(log_gdp_std = data_1$log_gdp_std,
                 rugged_std = data_1$rugged_std,
                 cid = as.integer(data_1$cid))
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 = data_2, chains = 4, cores = 4)
summary(model)
## Inference for Stan model: f3314e777e4c586121dcc9de98266129.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## a[1]    0.89    0.00 0.02   0.86   0.88   0.89   0.90   0.92  3040    1
## a[2]    1.05    0.00 0.01   1.03   1.04   1.05   1.06   1.07  2689    1
## b[1]    0.14    0.00 0.08  -0.01   0.09   0.14   0.18   0.28  2493    1
## b[2]   -0.14    0.00 0.06  -0.25  -0.18  -0.14  -0.11  -0.03  2260    1
## sigma   0.11    0.00 0.01   0.10   0.11   0.11   0.12   0.12  2567    1
## lp__  285.11    0.05 1.60 281.10 284.27 285.44 286.31 287.24  1001    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Sep 21 22:52:06 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
traceplot(model)
## [1] 1000
## [1] 1
## [1] 1000
y <- c(-1, 1)
model_1 <- ulam(alist(y ~ dnorm(mu, sigma),
                   mu <- alpha,
                   alpha ~ dnorm(0, 1000),
                   sigma ~ dexp(0.0001)), 
             data = list(y=y), chains = 3)
## 
## SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 
## Chain 1:  Elapsed Time: 0.024845 seconds (Warm-up)
## Chain 1:                0.008426 seconds (Sampling)
## Chain 1:                0.033271 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 3e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## 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 2: 
## Chain 2:  Elapsed Time: 0.029354 seconds (Warm-up)
## Chain 2:                0.038145 seconds (Sampling)
## Chain 2:                0.067499 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## 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 3: 
## Chain 3:  Elapsed Time: 0.016998 seconds (Warm-up)
## Chain 3:                0.007305 seconds (Sampling)
## Chain 3:                0.024303 seconds (Total)
## Chain 3:
## Warning: There were 93 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.07, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
summary(model_1)
## Inference for Stan model: 726d002e27cec1633082261fcfedb813.
## 3 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=1500.
## 
##         mean se_mean      sd    2.5%    25%   50%    75%   97.5% n_eff Rhat
## alpha  -5.14   19.33  298.70 -660.93 -40.27 -0.08  33.99  588.49   239 1.01
## sigma 477.00   87.53 1191.83    8.00  23.00 88.51 411.94 3329.73   185 1.02
## lp__   -5.18    0.29    1.97   -9.39  -6.62 -5.04  -3.56   -2.16    47 1.07
## 
## Samples were drawn using NUTS(diag_e) at Mon Sep 21 22:52:25 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
traceplot(model_1)

## [1] 1000
## [1] 1
## [1] 1000

9E7. Repeat the problem above, but now for a trace rank plot.

trankplot(model)
trankplot(model_1)

9M1. 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). Use ulam to estimate the posterior. Does the different prior have any detectible influence on the posterior distribution of sigma? Why or why not?

data("rugged")
data <- rugged
data$log_gdp <- log(data$rgdppc_2000)
data_1 <- data[complete.cases(data$rgdppc_2000),]
data_1$log_gdp_std <- data_1$log_gdp / mean(data_1$log_gdp)
data_1$rugged_std <- data_1$rugged/max(data_1$rugged)
data_1$cid <- ifelse(data_1$cont_africa==1,1,2)
data_2 <- list(log_gdp_std = data_1$log_gdp_std,
                 rugged_std = data_1$rugged_std,
                 cid = as.integer(data_1$cid))
model_1 <- 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 = data_2, chains = 1)
## 
## SAMPLING FOR MODEL 'ea0634026dfe9e2f56f6695bc98c30a0' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 
## Chain 1:  Elapsed Time: 0.050081 seconds (Warm-up)
## Chain 1:                0.032895 seconds (Sampling)
## Chain 1:                0.082976 seconds (Total)
## Chain 1:

9M2. Modify the terrain ruggedness model again. This time, change the prior for b[cid] to dexp(0.3). What does this do to the posterior distribution? Can you explain it?

data("rugged")
data <- rugged
data$log_gdp <- log(data$rgdppc_2000)
data_1 <- data[complete.cases(data$rgdppc_2000),]
data_1$log_gdp_std <- data_1$log_gdp / mean(data_1$log_gdp)
data_1$rugged_std <- data_1$rugged/max(data_1$rugged)
data_1$cid <- ifelse(data_1$cont_africa==1,1,2)
data_2 <- list(log_gdp_std = data_1$log_gdp_std,
                 rugged_std = data_1$rugged_std,
                 cid = as.integer(data_1$cid))
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] ~ dexp(0.3),
                 sigma ~ dexp(1)), 
           data = data_2)
## 
## SAMPLING FOR MODEL '2776e20029c9097e438c12597ff62d98' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 
## Chain 1:  Elapsed Time: 0.164153 seconds (Warm-up)
## Chain 1:                0.051027 seconds (Sampling)
## Chain 1:                0.21518 seconds (Total)
## Chain 1:
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
precis(model, depth = 2)
##             mean          sd         5.5%      94.5%    n_eff     Rhat4
## a[1]  0.88713411 0.015670857 0.8630258135 0.91278803 558.0143 0.9983758
## a[2]  1.04833690 0.011347894 1.0303199019 1.06761429 462.1599 1.0014825
## b[1]  0.14665898 0.067603912 0.0416923674 0.26272848 328.4541 0.9991979
## b[2]  0.01860442 0.017081580 0.0005062262 0.04795292 285.0583 0.9995922
## sigma 0.11402797 0.005949751 0.1051475912 0.12412851 386.1218 0.9995428

9M3. Re-estimate one of the Stan models from the chapter, but at different numbers 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?

data("rugged")
data <- rugged
data$log_gdp <- log(data$rgdppc_2000)
data_1 <- data[complete.cases(data$rgdppc_2000),]
data_1$log_gdp_std <- data_1$log_gdp / mean(data_1$log_gdp)
data_1$rugged_std <- data_1$rugged/max(data_1$rugged)
data_1$cid <- ifelse(data_1$cont_africa==1,1,2)
data_2 <- list(log_gdp_std = data_1$log_gdp_std,
                 rugged_std = data_1$rugged_std,
                 cid = as.integer(data_1$cid))
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 = data_2, chains = 1, warmup = 600)
## recompiling to avoid crashing R session
## 
## SAMPLING FOR MODEL 'f3314e777e4c586121dcc9de98266129' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 600 / 1000 [ 60%]  (Warmup)
## Chain 1: Iteration: 601 / 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: 
## Chain 1:  Elapsed Time: 0.058457 seconds (Warm-up)
## Chain 1:                0.026428 seconds (Sampling)
## Chain 1:                0.084885 seconds (Total)
## Chain 1:
precis(model, depth = 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8869126 0.016228411  0.85942522  0.91258992 520.3251 0.9991403
## a[2]   1.0510263 0.009860642  1.03436564  1.06677300 632.6724 0.9975102
## b[1]   0.1350327 0.075431956  0.01772231  0.25379827 722.5457 1.0006389
## b[2]  -0.1427423 0.053241863 -0.22959632 -0.05794828 614.4194 0.9985187
## sigma  0.1116093 0.006301957  0.10171395  0.12132624 523.0683 0.9975002
model_1 <- 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 = data_2, chains = 1, warmup = 400)
## recompiling to avoid crashing R session
## 
## SAMPLING FOR MODEL 'f3314e777e4c586121dcc9de98266129' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000323 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.23 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 401 / 1000 [ 40%]  (Sampling)
## Chain 1: Iteration: 500 / 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: 
## Chain 1:  Elapsed Time: 0.044227 seconds (Warm-up)
## Chain 1:                0.039821 seconds (Sampling)
## Chain 1:                0.084048 seconds (Total)
## Chain 1:
precis(model_1, depth = 2)
##             mean          sd         5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8870140 0.015913811  0.861575243  0.91112497 625.0983 1.0055942
## a[2]   1.0506078 0.010508017  1.032849658  1.06773774 974.1634 0.9986636
## b[1]   0.1290332 0.080765301 -0.003312213  0.26080605 725.1856 1.0075243
## b[2]  -0.1425303 0.054511302 -0.235860812 -0.05589022 793.4539 0.9987639
## sigma  0.1112966 0.006206491  0.102273461  0.12211443 699.5145 1.0022826

9H1. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.

model <- ulam(
 alist(
   a ~ dnorm(0,1),
   b ~ dcauchy(0,1)
 ), data=list(y=1) , chains=1 )
## 
## SAMPLING FOR MODEL '3bd3f4d287e9cccab124308e5415245c' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 9e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## 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: 
## Chain 1:  Elapsed Time: 0.026127 seconds (Warm-up)
## Chain 1:                0.019544 seconds (Sampling)
## Chain 1:                0.045671 seconds (Total)
## Chain 1:
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
traceplot(model)
## [1] 1000
## [1] 1
## [1] 1000

Compare the samples for the parameters a and b. 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?

#Plot a is a normal distribution as the prior is aroung 0 and spread in between 2 and -2. Plot b is Cauchy distribution which #contains some extreme value go up to over 30 and -50.