Exercise 3: Metropolis-Hastings algorithm(s) for the historical application (Beta-Bernoulli model)

Using the historical example, program an independant Metropolis-Hastings algorithm to estimate the posterior distribution of parameter\(\theta\) (i.e. the probability of having a girl for a birth). The prior distribution on \(\theta\) will be used as the instrumental proposal, and we will start by using a uniform prior on \(\theta\). We will consider the \(493,472\) births observed in Parisbetween 1745 and 1770, of which \(241,945\) were girls.

  1. Program a function that computes the numerator of the posterior density, which can be written \(p(\theta|n,S)\propto\theta^S(1-\theta)^{n-S}\) with \(S = 241\,945\) and \(n = 493\,472\) (plan for a boolean argument that will allow to return — or not — the logarithm of the posterior instead).

    post_num_hist <- function(theta, log = FALSE) {
    
        n <- 493472
        S <- 241945
    
        if (log) {
            num <- S * log(theta) + (n - S) * log(1 - theta)
        } else {
            num <- theta^S * (1 - theta)^(n - S)
        }
        return(num)
    }
    
    post_num_hist(0.2, log = TRUE)
    ## [1] -445522.1
    post_num_hist(0.6, log = TRUE)
    ## [1] -354063.6
  2. Program the corresponding Metropolis-Hastings algorithm, returning a vector of size \(n\) sampled according to the posterior distribution. Also returns the vector of acceptance probabilities \(\alpha\). What happens if this acceptance probability is NOT computed on the \(log\) scale ?

    myMH <- function(niter, post_num) {
        x_save <- numeric(length = niter)
        alpha <- numeric(length = niter)
        # initialise x0
        x <- runif(n = 1, min = 0, max = 1)
        # accpetance-rejection loop
        for (t in 1:niter) {
            # sample y from the proposal (uniform prior)
            y <- runif(n = 1, min = 0, max = 1)
            # compute the acceptance-rejection probability
            alpha[t] <- min(1, exp(post_num(y, log = TRUE) - post_num(x, log = TRUE)))
            # accept or reject
            u <- runif(1)
            if (u <= alpha[t]) {
                x_save[t] <- y
            } else {
                x_save[t] <- x
            }
            # update the current value
            x <- x_save[t]
        }
        return(list(theta = x_save, alpha = alpha))
    }
  3. Compare the posterior density obtained with this Metropolis-Hastings algorithm over 2000 itérations to the theoretical one (the theoretical density can be obtained with the R fuction dbeta(x, 241945 + 1, 251527 + 1) and represented with the R function curve(..., from = 0, to = 1, n = 10000)). Mindfully discard the first 500 iterations of your Metropolis-Hastings algorithm in order to reach the Markov chain convergence before constructing your Monte Carlo sample. Comment those results, especially in light of the acceptance probabilities computed throughout the algorithm, as well as the different sampled values for \(\theta\).

    sampleMH <- myMH(2000, post_num = post_num_hist)
    
    par(mfrow = c(2, 2))
    plot(density(sampleMH$theta[-c(1:500)]), col = "red", xlim = c(0, 1), ylab = "Posterior probability density", 
        xlab = expression(theta), main = "")
    curve(dbeta(x, 241945 + 1, 251527 + 1), from = 0, to = 1, n = 10000, add = TRUE)
    legend("topright", c("M-H", "theory"), col = c("red", "black"), lty = 1)
    
    plot(density(sampleMH$theta[-c(1:500)]), col = "red", ylab = "Posterior probability density", 
        xlab = expression(theta), main = "Zoom")
    curve(dbeta(x, 241945 + 1, 251527 + 1), from = 0, to = 1, n = 10000, add = TRUE)
    legend("topright", c("M-H", "theory"), col = c("red", "black"), lty = 1)
    
    plot(sampleMH$alpha, type = "h", xlab = "Iteration", ylab = "Probabilité d'acceptation", 
        ylim = c(0, 1), col = "springgreen")
    plot(sampleMH$theta, type = "l", xlab = "Iteration", ylab = expression(paste("Sampled value for ", 
        theta)), ylim = c(0, 1))

  4. Now imagine we only observe \(100\) births, among which \(49\) girls, and use a \(\text{Beta}(\alpha = 3, \beta=3)\) distribution as prior. Program the corresponding M-H algorithm and study the new results (one can do \(10,000\) iterations of this new M-H algorithm for instance, again mindfully discarding the first 500 iterations).

    post_num_beta <- function(theta, a = 3, b = 3, log = TRUE) {
    
        n <- 100  # number of trials (births)
        S <- 49  # number of success (feminine births)
    
        if (log) {
            num <- (a + S - 1) * (log(theta)) + (b + n - S - 1) * log(1 - theta)
        } else {
            num <- theta^(a + S - 1) * (1 - theta)^(b + n - S - 1)
        }
        return(num)
    }
    
    myMH_betaprior <- function(niter, post_num, a = 3, b = 3) {
        theta_save <- numeric(length = niter)  # vector of theta values ready to be saved
        alpha <- numeric(length = niter)  # vector of alpha values ready to be saved
        # initialisation of theta
        theta <- runif(n = 1, min = 0, max = 1)
        for (i in 1:niter) {
            # sample a value from the proposal (beta prior)
            theta_prop <- rbeta(n = 1, a, b)
            # compute acceptance-rejection probability
            alpha[i] <- min(1, exp(post_num(theta_prop, a = a, b = b, log = TRUE) - 
                post_num(theta, a = a, b = b, log = TRUE)))
            # acceptance-rejection step
            u <- runif(1)
            if (u <= alpha[i]) {
                theta <- theta_prop  # acceptance of theta_prop as new current value
            }
            # saving the current value of theta
            theta_save[i] <- theta
        }
        return(list(theta = theta_save, alpha = alpha))
    }
    sampleMH <- myMH_betaprior(10000, post_num = post_num_beta)
    
    
    par(mfrow = c(2, 2))
    plot(density(sampleMH$theta[-c(1:500)]), col = "red", xlim = c(0, 1), ylab = "Posterior probability density", 
        xlab = expression(theta), main = "")
    curve(dbeta(x, 49 + 1, 51 + 1), from = 0, to = 1, add = TRUE)
    legend("topright", c("M-H", "theory"), col = c("red", "black"), lty = 1)
    plot.new()
    plot(sampleMH$alpha, type = "h", xlab = "Iteration", ylab = "Probabilité d'acceptation", 
        ylim = c(0, 1), col = "springgreen")
    plot(sampleMH$theta, type = "l", xlab = "Iteration", ylab = expression(paste("Sampled value for ", 
        theta)), ylim = c(0, 1))

  5. Using the data from the historical example and with a \(\text{Beta}(\alpha = 3, \beta=3)\) prior, program a random-walk Metropolis-Hastings algorithm (with a Gaussian random step with a standard deviation of \(0.02\) for instance). Once again, study the results obtained this way (one can change the size of the random step).

    post_num_beta_hist <- function(theta, a = 3, b = 3, log = TRUE){
    
      n <- 493472 # number of trials (births)
      S <- 241945 # number of success (feminine births)
    
      if(log){
        num <- (a + S - 1)*(log(theta)) + (b + n - S - 1)*log(1-theta)
      }else{
        num <- theta^(a + S - 1)*(1-theta)^(b + n - S - 1)
      }
      return(num)
    }
    
    myMH_betaprior_marchalea <- function(niter, post_num, a=3, b=3){
    
      theta_save <- numeric(length = niter)
      alpha <- numeric(length = niter)
    
      #initialisation
      theta <- runif(n = 1, min = 0, max = 1)
    
      for(i in 1:niter){
        # sample a value from the proposal  (random walk)
        theta_prop <- theta + runif(1, -0.01, 0.01) #rnorm(1, mean = 0, sd = 0.02)
    
        # compute acceptance-rejection probability
        alpha[i] <- min(1, exp(post_num(theta_prop, a=a, b=b, log=TRUE) - 
                               post_num(theta, a=a, b=b, log=TRUE)))
    
        # acceptance-rejection step
        u <- runif(1)
        if(u <= alpha[i]){
          theta <- theta_prop # accept theta_prop and update current value
        }
    
        # save current value
        theta_save[i] <- theta
      }
    
      return(list("theta" = theta_save, "alpha" = alpha))
    }
    
    sampleMH <- myMH_betaprior_marchalea(20000, post_num = post_num_beta_hist)
    
    
    par(mfrow=c(2,2))
    plot(density(sampleMH$theta[-c(1:1000)]), col = "red", #xlim=c(0,1),
         ylab = "Posterior probability density", 
         xlab=expression(theta), main = "")
    curve(dbeta(x, 241945 + 1, 251527+1), from = 0, to = 1, n = 10000, add = TRUE)
    legend("topright", c("M-H", "theory"), col = c("red", "black"), lty = 1)
    plot(sampleMH$alpha, type = "h", xlab = "Iteration", 
         ylab = "Probabilité d'acceptation", ylim = c(0,1), col = "springgreen")
    plot(sampleMH$theta, type="l", xlab = "Iteration", 
         ylab = expression(paste("Sampled value for ", theta)), ylim = c(0,1))
    plot(sampleMH$theta, type="l", xlab = "Iteration", main = "Zoom",
         ylab = expression(paste("Sampled value for ", theta)), ylim = c(0.45, 0.55))