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(tidybayes)


install.packages(c('devtools','coda','mvtnorm'))
## Error in contrib.url(repos, "source"): trying to use CRAN without setting a mirror
library(devtools)
install_github("rmcelreath/rethinking")

library(rethinking)
rstan_options(auto_write = TRUE)
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)


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)

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 ...
m8.3 <- 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(m8.3 , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865637 0.015676124  0.86151027  0.91161722
## a[2]   1.0505719 0.009936882  1.03469085  1.06645297
## b[1]   0.1325070 0.074206399  0.01391089  0.25110320
## b[2]  -0.1425815 0.054750886 -0.23008401 -0.05507903
## sigma  0.1094972 0.005935716  0.10001077  0.11898362
pairs(m8.3)

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?

m9.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)
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499c67b11344.o] Error 127
## Error in sink(type = "output"): invalid connection
precis(m9.2 , depth=2)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'precis': object 'm9.2' not found

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?

# use fixed start values for comparability of runs
start <- list(a = c(1, 1), b = c(0, 0), sigma = 1)

m9.3 <- 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,
  start = start,
  chains = 4, 
  cores = 4,
  iter = 100
)
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499cab82caa.o] Error 127
## Error in sink(type = "output"): invalid connection
# define warmup values to run through
warm_list <- c(5, 10, 100, 500, 1000) 

# first make matrix to hold n_eff results
n_eff <- matrix(NA, nrow = length(warm_list), ncol = 5) 

for (i in 1:length(warm_list)) { # loop over warm_list and collect n_eff
  w <- warm_list[i]
  m_temp <- ulam(m9.3, chains = 4, cores = 4, iter = 1000 + w, warmup = w, start = start)
  n_eff[i, ] <- precis(m_temp, 2)$n_eff
}
## Error in ulam(m9.3, chains = 4, cores = 4, iter = 1000 + w, warmup = w, : object 'm9.3' not found
colnames(n_eff) <- rownames(precis(m_temp, 2))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'precis': object 'm_temp' not found
rownames(n_eff) <- warm_list
n_eff # columns show parameters, rows show n_eff
##      [,1] [,2] [,3] [,4] [,5]
## 5      NA   NA   NA   NA   NA
## 10     NA   NA   NA   NA   NA
## 100    NA   NA   NA   NA   NA
## 500    NA   NA   NA   NA   NA
## 1000   NA   NA   NA   NA   NA

9-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 )
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499c664c58f.o] Error 127
## Error in sink(type = "output"): invalid connection
# The prior is about 0 and the spread is between 2 and -2, thus plot 'a' should be a normal distribution. Plot 'b' #depicts a Cauchy distribution with several extreme values ranging from 30 to -40.

post <- extract.samples(mp)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'extract.samples': object 'mp' not found
par(mfrow = c(1, 2))

dens(post$a)
## Error in dens(post$a): object 'post' not found
curve(dnorm(x, 0, 1), from = -4, to = 4, add = T, lty = 2)
## Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet
legend("topright", lty = c(1, 2), legend = c("Sample", "Exact density"), bty = "n")
## Error in strwidth(legend, units = "user", cex = cex, font = text.font): plot.new has not been called yet
mtext("Normal")
## Error in mtext("Normal"): plot.new has not been called yet
dens(post$b, col = "red", xlim = c(-10, 10))
## Error in dens(post$b, col = "red", xlim = c(-10, 10)): object 'post' not found
curve(dcauchy(x, 0, 1),
  from = -10, to = 10, add = T, lty = 2,
  col = "red"
)
## Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet
mtext("Cauchy")
## Error in mtext("Cauchy"): plot.new has not been called yet

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?

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

# load data and copy
data(WaffleDivorce)
d <- WaffleDivorce

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

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

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_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499c3a31748b.o] Error 127
## Error in sink(type = "output"): invalid connection
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_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499c68592fb.o] Error 127
## Error in sink(type = "output"): invalid connection
m5.3 <- ulam(
  alist(
    D ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  data = d_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)
## Error in compileCode(f, code, language = language, verbose = verbose): sh: line 1: C:/rtools40/usr/mingw_64/bin/g++: No such file or directorymake: *** [C:/PROGRA~1/R/R-40~1.5/etc/x64/Makeconf:229: file499c6bf63ebb.o] Error 127
## Error in sink(type = "output"): invalid connection
set.seed(77)
compare( m5.1 , m5.2 , m5.3 , func=WAIC )
## Error in compare(m5.1, m5.2, m5.3, func = WAIC): object 'm5.1' not found