Assignment #8
Nischal Bondalapati
5/24/2022
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?
data(rugged)
data1 <- rugged
data1$log_gdp <- log(data1$rgdppc_2000)
data2 <- data1[ complete.cases(data1$rgdppc_2000) , ]
data2$log_gdp_std <- data2$log_gdp / mean(data2$log_gdp)
data2$rugged_std <- data2$rugged / max(data2$rugged)
data2$cid <- ifelse( data2$cont_africa == 1,1,2)
model1 <- 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=data2)
precis(model1,depth=2)## mean sd 5.5% 94.5%
## a[1] 0.8865493 0.015674673 0.86149813 0.91160044
## a[2] 1.0505701 0.009935964 1.03469048 1.06644966
## b[1] 0.1323451 0.074200042 0.01375909 0.25093109
## b[2] -0.1426724 0.054745851 -0.23016682 -0.05517793
## sigma 0.1094870 0.005934329 0.10000277 0.11897118pairs(model1)
model2 <- 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=data2)
precis(model2,depth=2)## mean sd 5.5% 94.5%
## a[1] 0.8865629 0.015680472 0.86150249 0.91162333
## a[2] 1.0505660 0.009939688 1.03468050 1.06645159
## b[1] 0.1325016 0.074226237 0.01387377 0.25112949
## b[2] -0.1425731 0.054765973 -0.23009974 -0.05504654
## sigma 0.1095284 0.005939946 0.10003521 0.11902157pairs(model2)
# Different priors did not have any detectible influence on the posterior distribution of sigma.
# It's probably because the data size is large enough.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?
model3 <- 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=data2)
precis(model3,depth=2)## mean sd 5.5% 94.5%
## a[1] 0.8865658 0.015679157 0.86150749 0.91162414
## a[2] 1.0505694 0.009938835 1.03468522 1.06645357
## b[1] 0.1325014 0.074220217 0.01388318 0.25111966
## b[2] -0.1425724 0.054761396 -0.23009165 -0.05505308
## sigma 0.1095189 0.005938662 0.10002779 0.11901005pairs(model3)
# Again, there is no apparent difference in the posterior distribution.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?
data(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ complete.cases(d$rgdppc_2000) , ]
dd.trim <- dd[ , c("log_gdp","rugged","cont_africa") ]
formula8m3 <- alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged + bA*cont_africa + bAR*rugged*cont_africa ,
a ~ dnorm(0,100),
bR ~ dnorm(0,10),
bA ~ dnorm(0,10),
bAR ~ dnorm(0,10),
sigma ~ dcauchy(0,2)
)
#Try different number of warm-ups
m8m3.w1 <- map2stan(formula8m3, data=dd.trim, iter=1001, warmup=1)##
## SAMPLING FOR MODEL 'a5e2b039ebe62a861814af5bf3ed4a79' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 6.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.65 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: No variance estimation is
## Chain 1: performed for num_warmup < 20
## Chain 1:
## Chain 1: Iteration: 1 / 1001 [ 0%] (Warmup)
## Chain 1: Iteration: 2 / 1001 [ 0%] (Sampling)
## Chain 1: Iteration: 101 / 1001 [ 10%] (Sampling)
## Chain 1: Iteration: 201 / 1001 [ 20%] (Sampling)
## Chain 1: Iteration: 301 / 1001 [ 30%] (Sampling)
## Chain 1: Iteration: 401 / 1001 [ 40%] (Sampling)
## Chain 1: Iteration: 501 / 1001 [ 50%] (Sampling)
## Chain 1: Iteration: 601 / 1001 [ 60%] (Sampling)
## Chain 1: Iteration: 701 / 1001 [ 70%] (Sampling)
## Chain 1: Iteration: 801 / 1001 [ 80%] (Sampling)
## Chain 1: Iteration: 901 / 1001 [ 90%] (Sampling)
## Chain 1: Iteration: 1001 / 1001 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.000109 seconds (Warm-up)
## Chain 1: 0.045239 seconds (Sampling)
## Chain 1: 0.045348 seconds (Total)
## Chain 1:precis(m8m3.w1) # awful results, n_eff=1## mean sd 5.5% 94.5% n_eff Rhat4
## a 1.276760 0 1.276760 1.276760 0.5007511 0.9989995
## bR -1.060461 0 -1.060461 -1.060461 0.5007511 0.9989995
## bA -1.564990 0 -1.564990 -1.564990 0.5007511 0.9989995
## bAR -1.917977 0 -1.917977 -1.917977 0.5007511 0.9989995
## sigma 1.696225 0 1.696225 1.696225 0.5007511 0.9989995m8m3.w10 <- map2stan(formula8m3, data=dd.trim, iter=1010, warmup=10)##
## SAMPLING FOR MODEL '80d712d269753490de9b45c0ce477697' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: No variance estimation is
## Chain 1: performed for num_warmup < 20
## Chain 1:
## Chain 1: Iteration: 1 / 1010 [ 0%] (Warmup)
## Chain 1: Iteration: 11 / 1010 [ 1%] (Sampling)
## Chain 1: Iteration: 111 / 1010 [ 10%] (Sampling)
## Chain 1: Iteration: 212 / 1010 [ 20%] (Sampling)
## Chain 1: Iteration: 313 / 1010 [ 30%] (Sampling)
## Chain 1: Iteration: 414 / 1010 [ 40%] (Sampling)
## Chain 1: Iteration: 515 / 1010 [ 50%] (Sampling)
## Chain 1: Iteration: 616 / 1010 [ 60%] (Sampling)
## Chain 1: Iteration: 717 / 1010 [ 70%] (Sampling)
## Chain 1: Iteration: 818 / 1010 [ 80%] (Sampling)
## Chain 1: Iteration: 919 / 1010 [ 90%] (Sampling)
## Chain 1: Iteration: 1010 / 1010 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.000821 seconds (Warm-up)
## Chain 1: 0.947471 seconds (Sampling)
## Chain 1: 0.948292 seconds (Total)
## Chain 1:precis(m8m3.w10) # not so awful results, n_eff=~100..200, troubles with estimates bA & sigma## mean sd 5.5% 94.5% n_eff Rhat4
## a 9.0090636 1.1958344 8.9125080 9.42053455 35.43273 1.024415
## bR -0.1738851 0.2410068 -0.3251808 -0.05030203 75.18333 1.003800
## bA -1.8533731 0.5297900 -2.3129273 -1.46361917 38.86823 1.007217
## bAR 0.3266993 0.3533128 0.1055254 0.59536911 34.09648 1.017763
## sigma 1.0705281 0.7862609 0.8733901 1.05900783 43.36793 1.023173plot(m8m3.w10)
m8m3.w100 <- map2stan(formula8m3, data=dd.trim, iter=1100, warmup=100)##
## SAMPLING FOR MODEL '0487dde54e11598f0c22d7f82052b248' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## 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:
## Chain 1: Iteration: 1 / 1100 [ 0%] (Warmup)
## Chain 1: Iteration: 101 / 1100 [ 9%] (Sampling)
## Chain 1: Iteration: 210 / 1100 [ 19%] (Sampling)
## Chain 1: Iteration: 320 / 1100 [ 29%] (Sampling)
## Chain 1: Iteration: 430 / 1100 [ 39%] (Sampling)
## Chain 1: Iteration: 540 / 1100 [ 49%] (Sampling)
## Chain 1: Iteration: 650 / 1100 [ 59%] (Sampling)
## Chain 1: Iteration: 760 / 1100 [ 69%] (Sampling)
## Chain 1: Iteration: 870 / 1100 [ 79%] (Sampling)
## Chain 1: Iteration: 980 / 1100 [ 89%] (Sampling)
## Chain 1: Iteration: 1090 / 1100 [ 99%] (Sampling)
## Chain 1: Iteration: 1100 / 1100 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.022009 seconds (Warm-up)
## Chain 1: 0.242995 seconds (Sampling)
## Chain 1: 0.265004 seconds (Total)
## Chain 1:precis(m8m3.w100) # enough of warmup## mean sd 5.5% 94.5% n_eff Rhat4
## a 9.2232467 0.14370464 8.9981501 9.45392125 262.0319 1.0023984
## bR -0.2052828 0.07893688 -0.3293788 -0.08320893 330.7064 1.0049679
## bA -1.9493293 0.23418168 -2.3150927 -1.60850712 239.3908 1.0063767
## bAR 0.3967211 0.13645549 0.1819950 0.61734075 294.1457 1.0035792
## sigma 0.9506511 0.05230445 0.8714026 1.03608490 1272.8603 0.9992234plot(m8m3.w100)
m8m3.w100 <- map2stan(formula8m3, data=dd.trim, iter=1100, warmup=150)##
## SAMPLING FOR MODEL '402abcd70a82b0ca09b90b227a1b1c4f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 1100 [ 0%] (Warmup)
## Chain 1: Iteration: 110 / 1100 [ 10%] (Warmup)
## Chain 1: Iteration: 151 / 1100 [ 13%] (Sampling)
## Chain 1: Iteration: 260 / 1100 [ 23%] (Sampling)
## Chain 1: Iteration: 370 / 1100 [ 33%] (Sampling)
## Chain 1: Iteration: 480 / 1100 [ 43%] (Sampling)
## Chain 1: Iteration: 590 / 1100 [ 53%] (Sampling)
## Chain 1: Iteration: 700 / 1100 [ 63%] (Sampling)
## Chain 1: Iteration: 810 / 1100 [ 73%] (Sampling)
## Chain 1: Iteration: 920 / 1100 [ 83%] (Sampling)
## Chain 1: Iteration: 1030 / 1100 [ 93%] (Sampling)
## Chain 1: Iteration: 1100 / 1100 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.039755 seconds (Warm-up)
## Chain 1: 0.209053 seconds (Sampling)
## Chain 1: 0.248808 seconds (Total)
## Chain 1:precis(m8m3.w100) # enough of warmup## mean sd 5.5% 94.5% n_eff Rhat4
## a 9.2296506 0.14134970 9.0141702 9.47189189 389.2066 1.001830
## bR -0.2062499 0.07780766 -0.3355150 -0.08091629 447.8372 1.002765
## bA -1.9607112 0.22052385 -2.2967485 -1.61660111 313.0264 1.000181
## bAR 0.3985882 0.12801653 0.2062488 0.59950519 434.9286 1.001383
## sigma 0.9498095 0.04839802 0.8768174 1.02614815 895.0560 1.002066plot(m8m3.w100)
m8m3.w100 <- map2stan(formula8m3, data=dd.trim, iter=1100, warmup=200)##
## SAMPLING FOR MODEL '965da3ac8732709907b43f2b5450ba77' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 1100 [ 0%] (Warmup)
## Chain 1: Iteration: 110 / 1100 [ 10%] (Warmup)
## Chain 1: Iteration: 201 / 1100 [ 18%] (Sampling)
## Chain 1: Iteration: 310 / 1100 [ 28%] (Sampling)
## Chain 1: Iteration: 420 / 1100 [ 38%] (Sampling)
## Chain 1: Iteration: 530 / 1100 [ 48%] (Sampling)
## Chain 1: Iteration: 640 / 1100 [ 58%] (Sampling)
## Chain 1: Iteration: 750 / 1100 [ 68%] (Sampling)
## Chain 1: Iteration: 860 / 1100 [ 78%] (Sampling)
## Chain 1: Iteration: 970 / 1100 [ 88%] (Sampling)
## Chain 1: Iteration: 1080 / 1100 [ 98%] (Sampling)
## Chain 1: Iteration: 1100 / 1100 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.046844 seconds (Warm-up)
## Chain 1: 0.199775 seconds (Sampling)
## Chain 1: 0.246619 seconds (Total)
## Chain 1:precis(m8m3.w100) # enough of warmup## mean sd 5.5% 94.5% n_eff Rhat4
## a 9.2264886 0.13864395 8.9941284 9.43960358 434.1535 0.9997483
## bR -0.2071663 0.07542822 -0.3200264 -0.08219634 446.3326 0.9995511
## bA -1.9360720 0.21518321 -2.2677956 -1.58667685 481.7326 0.9991942
## bAR 0.3912298 0.12539861 0.1896259 0.59060616 527.4999 0.9997030
## sigma 0.9516791 0.05072172 0.8737147 1.03609017 591.4396 1.0036831plot(m8m3.w100)
# Looks like 100-150 warmup are enough. 200 warmup performs worse than 150.9-4. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.
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?
mp <- 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 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 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.013041 seconds (Warm-up)
## Chain 1: 0.011437 seconds (Sampling)
## Chain 1: 0.024478 seconds (Total)
## Chain 1:traceplot(mp)
precis(mp)## mean sd 5.5% 94.5% n_eff Rhat4
## a 0.005441882 0.9411759 -1.427882 1.505640 246.5845 0.997998
## b -0.184235375 5.8987412 -5.564838 4.861189 142.9511 1.008421plot(mp, n_cols = 1, col = "blue")
# Plot a is more normally distributed while plot b tends to be a Cauchy distribution.
# Most of the values in a range from -2 to 2 while b has much more extreme values / outliers, which can go from -60 to 60.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.
library(tidybayes)
data(WaffleDivorce)
data = WaffleDivorce
data$Divorce_sd=standardize(data$Divorce)
data$Marriage_sd=standardize(data$Marriage)
data$MedianAgeMarriage_sd=standardize(data$MedianAgeMarriage)
d_trim = list(D = data$Divorce_sd, M = data$Marriage_sd, A = data$MedianAgeMarriage_sd)
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
)
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
)
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
)
set.seed(77)
compare( m5.1 , m5.2 , m5.3 , func=WAIC )## WAIC SE dWAIC dSE pWAIC weight
## m5.1 125.3762 12.499661 0.000000 NA 3.442654 0.732068740
## m5.3 127.3928 12.566473 2.016577 0.8233072 4.602175 0.267090096
## m5.2 138.9139 9.707619 13.537686 9.1471820 2.761965 0.000841164# Model m5.1 with only median age at marriage as a predictor performs best, but is not really distinguishable from model m5.3.
# However, the model with marriage rate only, m5.2 clearly performs worse than both.