Statistical Rethinking - Chapter 09

Markov Chain Monte Carlo

Introduction

These are answers and solutions to the exercises at the end of chapter 9 in Satistical Rethinking 2 by Richard McElreath. I have created these notes as a part of my ongoing involvement in the AU Bayes Study Group. Much of my inspiration for these solutions, where necessary, has been obtained from Taras Svirskyi, William Wolf, and Corrie Bartelheimer as well as the solutions provided to instructors by Richard McElreath himself.

R Environment

For today’s exercise, I load the following packages:

library(rethinking)
library(rstan)
library(ggplot2)
library(tidybayes)

Easy Exercises

Practice E1

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

Answer:

  1. Not a requirement. Metropolis can accommodate continuous and discrete parameters.
  2. Not a requirement. Distribution could be any symmetric distribution. Not just Gaussian.
  3. This is a requirement.

Practice E2

Question: 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?

Answer: Gibbs uses adaptive proposals when considering which location in the posterior to sample next. This makes it more efficient because less proposed steps are rejected.

Practice E3

Question: Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?

Answer: Discrete parameters. HMC depends on gradients which to explore using a physics simulation. Discrete parameters would not allow for the construction of any gradients.

Practice E4

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

Answer: Effective sample number (n_eff) identifies the number of ‘ideal’ (i.e. uncorrelated) samples. Since MCMC algorithms explore the posterior as a chain of samples, each sample is usually correlated with the previous one to some extent. Conclusively, n_eff identifies the number of samples used for estimating the posterior mean/distribution whereas actual number of samples is simply the number of data points we have.

n_eff is usually smaller than the actual number of samples (unless we have anti-correlated MCMC samples).

Practice E5

Question: Which value should Rhat approach, when a chain is sampling the posterior distribution correctly?

Answer: \(\hat{R}\) or Rhat, in R, reflects variance within a chain versus variance between chains. If these are the same, \(\hat{R}\) will be \(1.0\) - i.e.: it does not matter from which chain we would infere parameters and predictions. Values higher than 1.0 can indicate problems in the model. Values much higher than 1 indicate serious issues.

Practice E6

Question: 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?

Answer:

Good trace plot

y <- rnorm(1e4, mean = 1, sd = 2)
m.E6Good <- ulam(
  alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha,
    alpha ~ dnorm(0, 10),
    sigma ~ dcauchy(0, 1)
  ),
  data = list(y = y),
  cores = 2,
  chains = 2,
  start = list(
    alpha = 0,
    sigma = 1
  )
)
traceplot(m.E6Good)
## [1] 1000
## [1] 1
## [1] 1000

These trace plots show that the chains quickly find the region with highest posterior probability and stay there.

Bad trace plot

y <- rnorm(1e4, mean = 1, sd = 2)
m.E6Bad <- ulam(
  alist(
    y ~ dnorm(mu, sigma),
    mu <- a1 + a2,
    a1 ~ dnorm(0, 10),
    a2 ~ dnorm(0, 10),
    sigma ~ dcauchy(0, 1)
  ),
  data = list(y = y),
  chains = 2,
  cores = 2,
  start = list(
    a1 = 0,
    a2 = 0,
    sigma = 1
  ),
)
traceplot(m.E6Bad)
## [1] 1000
## [1] 1
## [1] 1000

This is a problem of unidentifiable parameters as a1 and a2 can cancel each other out to arrive at the correct mu and so we see non-stationary behaviour in the trace plots of a1 and a2 while the trace plot for sigma is doing alright.

Medium Exercises

Practice M1

Question: Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior and an exponential prior for the standard deviation, sigma. The uniform prior should be dunif(0,10) and the exponential should be dexp(1). Do the different priors have any detectable influence on the posterior distribution?

Answer: The ruggedness model in question is m8.3 in the book (or m9.1 in ulam() specification). First, I prepare the data like I did previously.

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

Let’s fit that model with the different priors:

## 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,
)

Now on to inspect the model. Let’s start with the parameter estimates in comparison

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     170

These are strikingly the same. What about the individual model outputs in more detail?

precis(m.M1Exp, depth = 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8865564 0.015510598  0.86200325  0.91138982 3031.800 1.0015043
## a[2]   1.0505815 0.009906150  1.03476280  1.06587747 3009.023 0.9986602
## b[1]   0.1316994 0.073266600  0.01442329  0.24912656 2567.528 0.9994985
## b[2]  -0.1427205 0.054396078 -0.23250879 -0.05435822 2577.099 0.9994622
## sigma  0.1115593 0.006429337  0.10149863  0.12232310 3008.556 0.9983962
precis(m.M1Uni, depth = 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8863745 0.016118776  0.86099192  0.91194583 2611.365 1.0010248
## a[2]   1.0503312 0.009879485  1.03510955  1.06537664 2731.848 0.9985196
## b[1]   0.1317989 0.076187921  0.01015468  0.25606572 2151.178 1.0005439
## b[2]  -0.1428677 0.057582799 -0.23752882 -0.05317706 2172.248 1.0006485
## sigma  0.1115922 0.005954098  0.10270642  0.12160327 3098.717 1.0002062

Again, these 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. Let me visualise this:

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

That really does look the same to me.

Practice M2

Question: The Cauchy and exponential priors from the terrain ruggedness model are very weak. They can be made more informative by reducing their scale. Compare the dcauchy and dexp priors for progressively smaller values of the scaling parameter. As these priors become stronger, how does each influence the posterior distribution?

Answer: I write a for loop here to minimise code needs:

RepTimes <- 4 # how many steps I want to try
ScalingFactor <- 10 # by what factor to make priors stronger
# empty lists to store models in
Explist <- as.list(rep(NA, RepTimes))
Caulist <- as.list(rep(NA, RepTimes))
# Loop over all models
for (Mod_Iter in 0:(RepTimes - 1)) {
  dd.trim$ScalingFactor <- ScalingFactor
  dd.trim$Mod_Iter <- Mod_Iter
  ## Exponential prior for sigma
  m.M2Exp <- 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 * ScalingFactor^Mod_Iter)
    ),
    data = dd.trim,
    chains = 4,
    cores = 4,
  )
  Explist[[Mod_Iter + 1]] <- m.M2Exp
  ## Cauchy prior for sigma
  m.M2Cau <- 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 ~ dcauchy(0, 1 / ScalingFactor^Mod_Iter)
    ),
    data = dd.trim,
    chains = 4,
    cores = 4,
  )
  Caulist[[Mod_Iter + 1]] <- m.M2Cau
}
coeftab(Explist[[1]], Explist[[2]], Explist[[3]], Explist[[4]])
##       Explist[[1]] Explist[[2]] Explist[[3]] Explist[[4]]
## a[1]     0.89         0.89         0.89         0.89     
## a[2]     1.05         1.05         1.05         1.05     
## b[1]     0.13         0.13         0.13         0.13     
## b[2]    -0.14        -0.15        -0.14        -0.14     
## sigma    0.11         0.11         0.11         0.09     
## nobs      170          170          170          170
coeftab(Caulist[[1]], Caulist[[2]], Caulist[[3]], Caulist[[4]])
##       Caulist[[1]] Caulist[[2]] Caulist[[3]] Caulist[[4]]
## a[1]     0.89         0.89         0.89         0.89     
## a[2]     1.05         1.05         1.05         1.05     
## b[1]     0.13         0.13         0.13         0.13     
## b[2]    -0.14        -0.14        -0.14        -0.14     
## sigma    0.11         0.11         0.11         0.11     
## nobs      170          170          170          170

The more restrictive exponential priors decrease the estimate for sigma. On the other hand, the more restrictive cauchy priors have no effect, it seems.

Let’s explore why this is by looking at the priors themselves:

par(mfrow = c(1, 2))
curve(dexp(x, 1),
  from = 0, to = 5, ylab = "Density", xlab = "sigma",
  col = "royalblue4"
)
curve(dexp(x, 10), from = 0, to = 5, add = T)
curve(dexp(x, 100), from = 0, to = 5, add = T, col = col.desat("red"))
curve(dexp(x, 1000), from = 0, to = 5, add = T, col = col.desat("green"))
mtext("Exponential Prior")
legend("topright",
  col = c("royalblue4", "black", col.desat("red"), col.desat("green")),
  lty = c(1, 1, 1), legend = c("Exp(1)", "Exp(10)", "Exp(100)", "Exp(1000)"), bty = "n"
)

curve(2 * dcauchy(x, 0, 1),
  from = 0, to = 5, ylab = "Density", xlab = "sigma",
  col = "royalblue4"
)
curve(2 * dcauchy(x, 0, 0.1), from = 0, to = 5, add = T, col = "black")
curve(2 * dcauchy(x, 0, 0.01), from = 0, to = 5, add = T, col = col.desat("red"))
curve(2 * dcauchy(x, 0, 0.001), from = 0, to = 5, add = T, col = col.desat("green"))
mtext("Cauchy Prior")
legend("topright",
  col = c("royalblue4", "black", col.desat("red"), col.desat("green")),
  lty = c(1, 1, 1), legend = c("Cauchy(0, 1)", "Cauchy(0, 0.1)", "Cauchy(0, 0.01)", "Cauchy(0, 0.001)"), bty = "n"
)

The cauchy distributions show thicker tails while the exponential distributions quickly concentrate. Hence why a concentrated Cauchy prior allow more flexibility that a concentrated exponential prior.

Practice M3

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

Answer: The ruggedness model was fine so far so I continue with that one. Here, I build this model with a fixed run length and fixed starting values for each run with changing warmup values:

start <- list(a = c(1, 1), b = c(0, 0), sigma = 1) # use fixed start values for comparability of runs
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) # first make matrix to hold n_eff results
for (i in 1:length(warm_list)) { # loop over warm_list and collect n_eff
  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 # columns show parameters, rows show n_eff
##           a[1]       a[2]       b[1]       b[2]      sigma
## 5     183.2503   43.06876   26.56609    9.49538   30.68489
## 10   1675.4694 1585.83557  987.57928 1437.60536 1349.10366
## 100  2677.0793 2128.86831  795.72972 1198.74889 1123.86516
## 500  3270.1775 2784.87502 2475.93012 2193.86536 2377.33134
## 1000 2697.5280 2604.32903 2596.03499 2307.15721 2235.93952

As we can see, 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.

Hard Exercises

Practice H1

Question: 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, using what you know about the Cauchy distribution?

Answer: First of all, let’s inspect the posterior:

precis(mp)
##          mean        sd      5.5%    94.5%     n_eff    Rhat4
## a -0.02047831  1.001006 -1.641422 1.563926 13542.438 0.999995
## b -0.14880387 12.056711 -5.619659 5.696048  2565.544 1.001984

Oof. Those uncertainties don’t look good at all! So what does the model even do? It simply just samples a from a normal distribution with mean 0 and standard deviation 1. b is sampled from a cauchy distribution. Let’s look at the traceplot for this:

plot(mp, n_cols = 1, col = "royalblue4")

As we can see, there are quite some outliers in the sampling of the cauchy distribution (b). Why is that? 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. Note that this also decreases n_eff. lp in the above is the log-posterior.

Now let’s see how the samples we drew measure up against the underlying functions of a and b, respectively:

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")

As we can see, the normal distribution has been reconstructed well. The cauchy distributions hasn’t.

Practice H2

Question: 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. Explain the results.

Answer: First, I need to load the data and prepare it for ulam():

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)

Now I fit the models with ulam():

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
)

Now we compare the models:

compare(m5.1_stan, m5.2_stan, m5.3_stan, func = PSIS)
##               PSIS        SE     dPSIS       dSE    pPSIS      weight
## m5.1_stan 126.2612 13.117772  0.000000        NA 3.964718 0.748742859
## m5.3_stan 128.4547 13.082707  2.193473 0.8111425 5.159298 0.250049605
## m5.2_stan 139.1208  9.856804 12.859628 9.6374970 2.834047 0.001207536
compare(m5.1_stan, m5.2_stan, m5.3_stan, func = WAIC)
##               WAIC        SE     dWAIC       dSE    pWAIC      weight
## m5.1_stan 126.1634 12.929713  0.000000        NA 3.915828 0.717680850
## m5.3_stan 128.0375 12.679415  1.874051 0.8754981 4.950698 0.281181179
## m5.2_stan 139.0570  9.728896 12.893556 9.5812746 2.802120 0.001137971

WAIC 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. Let’s call these equal in performance and investigate why:

precis(m5.3_stan)
##               mean         sd       5.5%      94.5%    n_eff    Rhat4
## a     -0.002532436 0.10071982 -0.1608884  0.1556612 1798.064 1.003274
## bA    -0.604340251 0.16532280 -0.8622220 -0.3403890 1135.186 1.001120
## bM    -0.061063639 0.16503227 -0.3267091  0.2006084 1079.065 1.000238
## sigma  0.830525323 0.08932759  0.7003602  0.9876803 1260.195 1.000651

While m5.3_stan contains the marriage predictor, it is very unsure of it’s influence. In practical terms, this means that m5.1_stan and m5.3_stan make basically the same predictions

Practice H3

Question: Sometimes changing a prior for one parameter has unanticipated effects on other parameters. This is because when a parameter is highly correlated with another parameter in the posterior, the prior influences both parameters. Here’s an example to work and think through.
Go back to the leg length example in Chapter 5. Here is the code again, which simulates height and leg lengths for 100 imagined individuals:

N <- 100 # number of individuals
height <- rnorm(N, 10, 2) # sim total height of each
leg_prop <- runif(N, 0.4, 0.5) # leg as proportion of height
leg_left <- leg_prop * height + rnorm(N, 0, 0.02) # sim left leg as proportion + error
leg_right <- leg_prop * height + rnorm(N, 0, 0.02) # sim right leg as proportion + error
d <- data.frame(height, leg_left, leg_right) # combine into data frame

And below is the model you fit before, resulting in a highly correlated posterior for the two beta parameters. This time, fit the model using ulam():

m5.8s <- ulam(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + bl * leg_left + br * leg_right,
    a ~ dnorm(10, 100),
    bl ~ dnorm(2, 10),
    br ~ dnorm(2, 10),
    sigma ~ dexp(1)
  ),
  data = d,
  chains = 4,
  cores = 4,
  start = list(
    a = 10,
    bl = 0,
    br = 0.1,
    sigma = 1
  )
)

Compare the posterior distribution produced by the code above to the posterior distribution produced when you change the prior for br so that it is strictly positive:

m5.8s2 <- ulam(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + bl * leg_left + br * leg_right,
    a ~ dnorm(10, 100),
    bl ~ dnorm(2, 10),
    br ~ dnorm(2, 10),
    sigma ~ dexp(1)
  ),
  data = d,
  chains = 4,
  cores = 4,
  constraints = list(br = "lower=0"),
  start = list(
    a = 10,
    bl = 0,
    br = 0.1,
    sigma = 1
  )
)

Note the constraints list. What this does is constrain the prior distribution of br so that it has positive probability only above zero. In other words, that prior ensures that the posterior distribution for br will have no probability mass below zero.
Compare the two posterior distributions for m5.8s and m5.8s2. What has changed in the posterior distribution of both beta parameters? Can you explain the change induced by the change in prior?

Answer: It’s probably easiest to just look at the posterior distributions of the beta prameters through the pairs() function:

pairs(m5.8s, main = "Model 1")

pairs(m5.8s2, main = "Model 2")

As we can see, the beta distributions have shifted drastically between the different models. Interestingly, bl and br were perfectly symmetric in m5.8s, but are skewed in m5.8s2. Given how the height of a person is approximated in both models (a + bl*leg_left + br*leg_right), the distributions of leg lengths are necessarily negatively correlated (you can be of the same height with a short right leg and long left leg, long left leg and short right leg, or two medium-length legs). Thus, by setting br to be strictly positive in m5.8s2 and made it skewed, we have forced bl to be equally skewed in a mirror image of br.

Practice H4

Question: For the two models fit in the previous problem, use WAIC or PSIS to compare the effective numbers of parameters for each model. You will need to use log_lik=TRUE to instruct ulam() to compute the terms that both WAIC and PSIS need. Which model has more effective parameters? Why?

Answer: Let’s run the models:

m.H4_1 <- ulam(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + bl * leg_left + br * leg_right,
    a ~ dnorm(10, 100),
    bl ~ dnorm(2, 10),
    br ~ dnorm(2, 10),
    sigma ~ dexp(1)
  ),
  data = d,
  chains = 4,
  cores = 4,
  start = list(
    a = 10,
    bl = 0,
    br = 0.1,
    sigma = 1
  ),
  log_lik = TRUE
)
m.H4_2 <- ulam(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + bl * leg_left + br * leg_right,
    a ~ dnorm(10, 100),
    bl ~ dnorm(2, 10),
    br ~ dnorm(2, 10),
    sigma ~ dexp(1)
  ),
  data = d,
  chains = 4,
  cores = 4,
  constraints = list(br = "lower=0"),
  start = list(
    a = 10,
    bl = 0,
    br = 0.1,
    sigma = 1
  ),
  log_lik = TRUE
)

Now we compare them with WAIC:

compare(m.H4_1, m.H4_2)
##            WAIC       SE     dWAIC       dSE    pWAIC  weight
## m.H4_2 184.1422 11.79268 0.0000000        NA 3.150165 0.57075
## m.H4_1 184.7120 11.69593 0.5698236 0.3362031 3.410374 0.42925

The models are pretty much tied. The model with truncated priors (m.H4_2) is less flexible as indicated by pWAIC. This is because the prior is more informative and the variance in the posterior distribution is smaller as a result.

Practice H5

Question: Modify the Metropolis algorithm code from the chapter to handle the case that the island populations have a different distribution than the island labels. This means the island’s number will not be the same as its population.

Answer: First of all, we need our 10 islands with population sizes of 1-10, but in random order:

pop_size <- sample(1:10)

Now we can use the code from the chapter almost unaltered safe for one exception - we need to use indexing to translate island location into population size:

num_weeks <- 1e5
positions <- rep(NA, num_weeks)
current <- 10
for (i in 1:num_weeks) {
  positions[i] <- current # record current position
  proposal <- current + sample(c(-1, 1), size = 1) # flip coin to generate proposal
  # now make sure he loops around the archipelago
  if (proposal < 1) proposal <- 10
  if (proposal > 10) proposal <- 1
  prob_move <- pop_size[proposal] / pop_size[current] # move?
  current <- ifelse(runif(1) < prob_move, proposal, current)
}

To see if this works, we can plot population size against frequency of visit by the king:

f <- table(positions) # compute frequencies
plot(as.vector(f), pop_size,
  type = "n", # plot frequencies against relative population sizes
  xlab = "frequency", ylab = "population size"
) # empty plot
text(x = f, y = pop_size, labels = names(f)) # add names of islands / their positions

Practice H6

Question: Modify the Metropolis algorithm code from the chapter to write your own simple MCMC estimator for globe tossing data and model from Chapter 2.

Answer: We want to fit the following model:

\[w∼Binom(θ,n)\] \[θ∼Unif(0,1)\] Our Metropolis algorithm looks like this:

set.seed(42)
# the globe tossing data
w <- 6
n <- 9
# prior on p
p_prior <- function(p) dunif(p, min = 0, max = 1)
# initializing MCMC
iter <- 1e4
p_sample <- rep(0, iter)
p_current <- 0.5 # start value
for (i in 1:iter) {
  p_sample[i] <- p_current # # record current p
  p_proposal <- runif(1, min = 0, max = 1) # generate proposal
  # compute likelihood for current and proposal
  lkhd_current <- dbinom(w, n, p_current)
  lkhd_proposal <- dbinom(w, n, p_proposal)
  prob_proposal <- lkhd_proposal * p_prior(p_proposal)
  prob_current <- lkhd_current * p_prior(p_current)
  prob_accept <- prob_proposal / prob_current
  p_current <- ifelse(runif(1) < prob_accept, p_proposal, p_current)
}

Let’s visualise what happened here:

plot(p_sample, type = "l", col = "royalblue4")

Finally, let’s plot the posterior distribution:

dens(p_sample, col = "royalblue4", adj = 1)
curve(dbeta(x, w + 1, n - w + 1), from = 0, to = 1, add = T, lty = 2)
abline(v = median(p_sample))

Session Info

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.1252    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] tidybayes_2.3.1      rethinking_2.13      rstan_2.21.2         ggplot2_3.3.3        StanHeaders_2.21.0-7
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6           mvtnorm_1.1-1        lattice_0.20-41      tidyr_1.1.3          prettyunits_1.1.1    ps_1.6.0             assertthat_0.2.1     digest_0.6.27        utf8_1.2.1          
## [10] V8_3.4.0             plyr_1.8.6           R6_2.5.0             backports_1.2.1      stats4_4.0.2         evaluate_0.14        coda_0.19-4          highr_0.9            blogdown_1.3        
## [19] pillar_1.6.0         rlang_0.4.10         curl_4.3             callr_3.7.0          jquerylib_0.1.3      R.utils_2.10.1       R.oo_1.24.0          rmarkdown_2.7        styler_1.4.1        
## [28] labeling_0.4.2       stringr_1.4.0        loo_2.4.1            munsell_0.5.0        compiler_4.0.2       xfun_0.22            pkgconfig_2.0.3      pkgbuild_1.2.0       shape_1.4.5         
## [37] htmltools_0.5.1.1    tidyselect_1.1.0     tibble_3.1.1         gridExtra_2.3        bookdown_0.21        arrayhelpers_1.1-0   codetools_0.2-16     matrixStats_0.58.0   fansi_0.4.2         
## [46] crayon_1.4.1         dplyr_1.0.5          withr_2.4.2          MASS_7.3-51.6        R.methodsS3_1.8.1    distributional_0.2.2 ggdist_2.4.0         grid_4.0.2           jsonlite_1.7.2      
## [55] gtable_0.3.0         lifecycle_1.0.0      DBI_1.1.1            magrittr_2.0.1       scales_1.1.1         KernSmooth_2.23-17   RcppParallel_5.1.2   cli_2.4.0            stringi_1.5.3       
## [64] farver_2.1.0         bslib_0.2.4          ellipsis_0.3.1       generics_0.1.0       vctrs_0.3.7          rematch2_2.1.2       forcats_0.5.1        tools_4.0.2          svUnit_1.0.6        
## [73] R.cache_0.14.0       glue_1.4.2           purrr_0.3.4          processx_3.5.1       yaml_2.2.1           inline_0.3.17        colorspace_2.0-0     knitr_1.32           sass_0.3.1
Erik Kusch
Erik Kusch
PhD Student

In my research, I focus on statistical approaches to understanding complex processes and patterns in biology using a variety of data banks.

Related