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. 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

9-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. Visualize the posteriors for both models. Does the different prior have any detectible influence on the posterior distribution of sigma? Why or why not?

library(rethinking)
data(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ complete.cases(d$rgdppc_2000),]
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)

m1 <- quap(
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=dd)
precis(m1 , depth=2)
##             mean          sd       5.5%       94.5%
## a[1]   0.8865628 0.015675343  0.8615106  0.91161501
## a[2]   1.0505707 0.009936380  1.0346905  1.06645097
## b[1]   0.1325055 0.074202847  0.0139150  0.25109596
## b[2]  -0.1425762 0.054748190 -0.2300744 -0.05507801
## sigma  0.1094916 0.005934959  0.1000064  0.11897682
m1_unif <- quap(
  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=dd)

precis(m1_unif , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865655 0.015680676  0.86150477  0.91162627
## a[2]   1.0505693 0.009939815  1.03468358  1.06645507
## b[1]   0.1325015 0.074227148  0.01387216  0.25113079
## b[2]  -0.1425732 0.054766665 -0.23010089 -0.05504547
## sigma  0.1095298 0.005940141  0.10003633  0.11902332
pairs(m1)

pairs(m1_unif)

##### There is no visual difference between resulting models, at least difference is indistinguishable from several runs of the same model.

9-2. 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?

m2 <- quap(
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(0.3)
) ,
data=dd )
precis(m2 , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865649 0.015679201  0.86150655  0.91162333
## a[2]   1.0505696 0.009938864  1.03468534  1.06645379
## b[1]   0.1325029 0.074220422  0.01388433  0.25112147
## b[2]  -0.1425741 0.054761552 -0.23009368 -0.05505461
## sigma  0.1095192 0.005938706  0.10002805  0.11901045
pairs(m2)

##### There’s no visual difference between resulting models.

9-3. 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?

m3 <- quap(
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(0.3)
) ,
data=dd)
precis(m3 , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865778 0.015677828  0.86152165  0.91163404
## a[2]   1.0505617 0.009937970  1.03467891  1.06644450
## b[1]   0.1325023 0.074214056  0.01389389  0.25111068
## b[2]  -0.1425729 0.054756712 -0.23008469 -0.05506108
## sigma  0.1095092 0.005937351  0.10002020  0.11899827