Metropolis–Hastings for the logit-normal posterior

Topic: Markov chains, Gibbs sampling, and Metropolis–Hastings · more questions on this topic
‹ A bivariate-normal Gibbs sampler Random-walk Metropolis for a correlated normal ›

Use the posterior kernel in the logit-normal prior exercise.

  1. For n=40 and x=16, use a Metropolis–Hastings algorithm to retain 10{,}000 draws after a burn-in of 1{,}000. Plot a kernel-density estimate and describe it.

  2. Estimate the posterior mean, compare it with the prior mean, and explain how the data have changed the distribution.

  3. Estimate and interpret an equal-tailed 95\% credible interval.

  4. Repeat the numerical summaries for n=400 and x=160. Explain why ten times as much data changes the spread far more than the centre.

Question. For 0<\theta<1, take the posterior target

\pi(\theta\mid x,n)\propto \theta^{x-1}(1-\theta)^{n-x-1} \exp\!\left[-\frac12\{\operatorname{logit}(\theta)-1\}^2\right].

  1. For n=40 and x=16, use a Metropolis–Hastings algorithm to retain 10{,}000 draws after a burn-in of 1{,}000. Plot a kernel-density estimate and describe it.

  2. Estimate the posterior mean, compare it with the prior mean, and explain how the data have changed the distribution.

  3. Estimate and interpret an equal-tailed 95\% credible interval.

  4. Repeat the numerical summaries for n=400 and x=160. Explain why ten times as much data changes the spread far more than the centre.

Solution. For numerical stability define the log of the unnormalised target by

\ell(\theta;x,n) =(x-1)\log\theta+(n-x-1)\log(1-\theta) -\frac12\{\operatorname{logit}(\theta)-1\}^2.

Use an independent \operatorname{Beta}(A,B) proposal with density g. Since this is not a random walk, the proposal terms do not cancel. For current state \theta and proposal y, \begin{align*} r(\theta,y) &=\frac{\pi(y)g(\theta)}{\pi(\theta)g(y)},\\ \log r(\theta,y) &=\ell(y;x,n)-\ell(\theta;x,n) +\log g(\theta)-\log g(y),\\ a(\theta,y)&=\min\{1,r(\theta,y)\}. \end{align*} The following implementation uses \operatorname{Beta}(4,6) for (x,n)=(16,40) and \operatorname{Beta}(40,60) for (160,400).

log_posterior <- function(theta, x, n) {
  ifelse(theta <= 0 | theta >= 1, -Inf,
         (x - 1) * log(theta) +
         (n - x - 1) * log1p(-theta) -
         0.5 * (qlogis(theta) - 1)^2)
}
independent_mh <- function(x, n, A, B, keep = 10000, burn = 1000) {
  total <- keep + burn
  out <- numeric(total)
  out[1] <- x / n
  for (t in 2:total) {
    proposed <- rbeta(1, A, B)
    log_ratio <-
      log_posterior(proposed, x, n) -
      log_posterior(out[t - 1], x, n) +
      dbeta(out[t - 1], A, B, log = TRUE) -
      dbeta(proposed, A, B, log = TRUE)
    if (log(runif(1)) < min(0, log_ratio)) {
      out[t] <- proposed
    } else {
      out[t] <- out[t - 1]
    }
  }
  out[(burn + 1):total]
}
set.seed(3500)
draw40  <- independent_mh(16,  40,  4,  6)
draw400 <- independent_mh(160, 400, 40, 60)
mean(draw40)
quantile(draw40, c(0.025, 0.975))
mean(draw400)
quantile(draw400, c(0.025, 0.975))
plot(density(draw40, from = 0, to = 1),
     xlab = expression(theta), main = "n = 40")

For the displayed seed, the n=40 run gives

\widehat{\mathbb{E}(\theta\mid x)}=0.432, \qquad \text{$95\%$ interval }(0.294,0.578).

Direct numerical integration gives 0.4320 and (0.291,0.579), so the Monte Carlo discrepancy is of the expected size. The prior mean is approximately 0.6967: observing only 16/40=0.4 heads moves the posterior sharply towards 0.4 and makes it appreciably more concentrated. The interval means that, under the model and prior, the posterior probability that \theta lies between the reported endpoints is approximately 0.95.

For n=400 and x=160, the same seed gives

\widehat{\mathbb{E}(\theta\mid x)}=0.404, \qquad \text{$95\%$ interval }(0.355,0.453),

while numerical integration gives 0.4035 and (0.356,0.452). Both data sets have empirical proportion 0.4, so the centre changes little. The ten-fold increase in sample size makes the likelihood much more concentrated and leaves the prior with much less leverage.

The function dbeta in the Hastings ratio is essential: g is a proposal density. Replacing it by the distribution function pbeta changes the invariant distribution and produces incorrect posterior summaries.