Assignment #8
Yuman Liang
2021-10-28
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)
d <- d[complete.cases(d$rgdppc_2000), ]
d$log_gdp_std <- d$log_gdp / mean(d$log_gdp)
d$rugged_std <- d$rugged / max(d$rugged)
d$cid <- ifelse(d$cont_africa == 1, 1, 2)
dd.trim <- list(
log_gdp_std = d$log_gdp_std,
rugged_std = d$rugged_std,
cid = as.integer(d$cid)
)
## Exponential prior for sigma
m.M1Exp <- 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 = dd.trim,
chains = 4,
cores = 4,
)
## Uniform prior for sigma
m.M1Uni <- 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 ~ dnorm(0, 10)
),
data = dd.trim,
chains = 4,
cores = 4,
)# comparison on parameter estimations and model outputs
coeftab(m.M1Exp, m.M1Uni)## m.M1Exp m.M1Uni
## a[1] 0.89 0.89
## a[2] 1.05 1.05
## b[1] 0.13 0.13
## b[2] -0.14 -0.14
## sigma 0.11 0.11
## nobs 170 170precis(m.M1Exp, depth = 2)## mean sd 5.5% 94.5% n_eff Rhat4
## a[1] 0.8867687 0.015978819 0.860967142 0.91168559 2293.336 0.9988503
## a[2] 1.0505465 0.010439552 1.033764278 1.06678102 3073.613 0.9991305
## b[1] 0.1313059 0.078164276 0.001651231 0.25747849 2305.559 0.9996640
## b[2] -0.1431840 0.057228491 -0.235794057 -0.05123892 2598.666 0.9986544
## sigma 0.1116388 0.006243656 0.102098574 0.12207067 2423.680 0.9989816precis(m.M1Uni, depth = 2)## mean sd 5.5% 94.5% n_eff Rhat4
## a[1] 0.8862332 0.015574697 0.861608362 0.91138947 2687.423 0.9992697
## a[2] 1.0504867 0.010310671 1.033714205 1.06650752 2594.069 0.9993662
## b[1] 0.1286910 0.076557969 0.007135008 0.25217660 1971.117 1.0027358
## b[2] -0.1418789 0.057097776 -0.233857818 -0.05340888 3053.542 0.9984129
## sigma 0.1115594 0.006166425 0.102451192 0.12119363 2216.874 0.9985992According to comparison results, they are very similar aside from the effective number of samples (n_eff) which is much higher for all parameter estimates in the model with the exponential prior on sigma (m.M1Exp) except for sigma itself, which boasts a higher n_eff in the uniform-prior model (m.M1Uni). As such, we conclude that while the different priors have an impact on n_eff, they do not change the posterior distributions.
library(tidybayes)
Plot_df <- data.frame(
Posteriors = c(
extract.samples(m.M1Exp, n = 1e4)$sigma,
extract.samples(m.M1Uni, n = 1e4)$sigma
),
Name = rep(c("Exp", "Uni"), each = 1e4),
Model = rep(c("m.M1Exp", "m.M1Uni"), each = 1e4)
)
ggplot(Plot_df, aes(y = Model, x = Posteriors)) +
stat_halfeye() +
labs(x = "Parameter Estimate", y = "Model") +
theme_bw()
#### 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?
m.dexp <- 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 = dd.trim,
chains = 4,
cores = 4,
)pairs(m.dexp)
##### 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?
start <- list(a = c(1, 1), b = c(0, 0), sigma = 1)
m.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] ~ dnorm(0, 0.3),
sigma ~ dexp(1)
),
data = dd.trim,
start = start,
chains = 2, cores = 2,
iter = 100
)
warm_list <- c(5, 10, 100, 500, 1000) # define warmup values to run through
n_eff <- matrix(NA, nrow = length(warm_list), ncol = 5)
for (i in 1:length(warm_list)) {
w <- warm_list[i]
m_temp <- ulam(m.M3, chains = 2, cores = 2, iter = 1000 + w, warmup = w, start = start)
n_eff[i, ] <- precis(m_temp, 2)$n_eff
}colnames(n_eff) <- rownames(precis(m_temp, 2))
rownames(n_eff) <- warm_list
n_eff ## a[1] a[2] b[1] b[2] sigma
## 5 49.52332 2.051028 1.116597 1.133869 2.997175
## 10 1772.33944 2097.039429 733.235732 1072.699095 1128.508213
## 100 2277.76378 2325.402344 907.028035 1093.852277 1242.376340
## 500 2193.97472 2785.835168 2372.957128 2896.383914 2626.673708
## 1000 2694.02300 2958.470890 2466.190900 2641.928713 2124.069310According to above, past just 10 warmup samples, n_eff does not change much (in terms of how useful our samples are). In this case, we could be quite happy with a warmup of 10.
9-4. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.
mp <- map2stan(
alist(
a ~ dnorm(0, 1),
b ~ dcauchy(0, 1)
),
data = list(y = 1),
start = list(a = 0, b = 0),
iter = 1e4,
chains = 2, cores = 2,
warmup = 100,
WAIC = FALSE
)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?
precis(mp)## mean sd 5.5% 94.5% n_eff Rhat4
## a -0.004296876 1.003318 -1.604828 1.596984 13280.724 1.000050
## b 0.208847407 14.043479 -4.688634 4.574779 3680.756 1.000632plot(mp, n_cols = 1, col = "royalblue4")
#### As we can see, there are quite some outliers in the sampling of the cauchy distribution (b). Because the cauchy distribution has very heavy tails thus making it more likely to jump to a value that is far out there in terms of posterior probability.
post <- extract.samples(mp)
par(mfrow = c(1, 2))
dens(post$a)
curve(dnorm(x, 0, 1), from = -4, to = 4, add = T, lty = 2)
legend("topright", lty = c(1, 2), legend = c("Sample", "Exact density"), bty = "n")
mtext("Normal")
dens(post$b, col = "royalblue4", xlim = c(-10, 10))
curve(dcauchy(x, 0, 1),
from = -10, to = 10, add = T, lty = 2,
col = "royalblue4"
)
mtext("Cauchy")
#### The normal distribution has been reconstructed well. The cauchy distributions hasn’t.
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(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce
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_stan <- 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
)
m5.2_stan <- 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
)
m5.3_stan <- 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
)compare(m5.1_stan, m5.2_stan, m5.3_stan, func = PSIS)## PSIS SE dPSIS dSE pPSIS weight
## m5.1_stan 125.7335 12.545387 0.000000 NA 3.583250 0.6937262851
## m5.3_stan 127.3740 12.758069 1.640476 0.7486677 4.631002 0.3054662353
## m5.2_stan 139.2454 9.824084 13.511830 8.9499484 2.883569 0.0008074795compare(m5.1_stan, m5.2_stan, m5.3_stan, func = WAIC)## WAIC SE dWAIC dSE pWAIC weight
## m5.1_stan 125.7686 12.460310 0.000000 NA 3.600797 0.6756072561
## m5.3_stan 127.2410 12.563653 1.472378 0.7310699 4.564500 0.3235723620
## m5.2_stan 139.1958 9.713362 13.427194 9.0017032 2.858798 0.0008203819WAIC tells a similar story as PSIS, but the model only containing age (m5.1_stan) wins. The model with both predictors (m5.3_stan) does almost as well. However, their respective PSIS and WAIC values are nearly identical. Furthermore, both models get assigned all of the WAIC weight.
precis(m5.3_stan)## mean sd 5.5% 94.5% n_eff Rhat4
## a 0.0002427025 0.1024926 -0.1606875 0.1688413 1855.109 0.9982351
## bA -0.6086085365 0.1520980 -0.8538119 -0.3690931 1251.173 1.0010587
## bM -0.0620155686 0.1573815 -0.3094343 0.1919550 1353.417 1.0009371
## sigma 0.8262231406 0.0846412 0.7016010 0.9695795 1505.727 0.9998936