Writing MCMC Sampling Code by Hand

Question

I am trying to understand the steps involved in the Metropolis-Hastings Algorithm and trying to learn how to implement it myself.

As an example, suppose I am interested in estimating the "probability of success" (p) in a Binomial Distribution. I decide to use Bayesian Estimation - that is, I can define the Posterior Probability Distribution of "p" as a function (proportional) of the Binomial Likelihood and some Prior Probability Distribution (e.g. the Beta Distribution):

$$\text{Posterior} = \text{Pr}(p | x) = \frac{ p^{\sum_{i=1}^n x_i + \alpha - 1} \cdot (1-p)^{n - \sum_{i=1}^n x_i + \beta - 1}}{\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \cdot \frac{\Gamma(\sum_{i=1}^n x_i + \alpha)\Gamma(n - \sum_{i=1}^n x_i + \beta)}{\Gamma(n + \alpha + \beta)}} = \frac{ p^{\sum_{i=1}^n x_i + \alpha - 1} \cdot (1-p)^{n - \sum_{i=1}^n x_i + \beta - 1}}{\frac{\Gamma(\sum_{i=1}^n x_i + \alpha)\Gamma(n - \sum_{i=1}^n x_i + \beta)}{\Gamma(n + \alpha + \beta)}}$$

I am aware that the Binomial Distribution and the Beta Distribution have the "Conjugacy Property" (https://en.wikipedia.org/wiki/Conjugate_prior) - result in a "closed-form solution" for the Posterior Distribution and not requiring to use the Metropolis-Hastings Algorithm. However, I am interested in learning how to write the steps of the Metropolis-Hastings Algorithm to take random samples from the Posterior Distribution as an educational exercise (e.g. I could then take the average of these samples to obtain the "Bayes Estimator" for "p").

As I understand, here are the steps necessary in the Metropolis-Hastings Algorithm:

• (Step 1) Define P-Current : A single random number generated from a Uniform(0,1) Distribution
• (Step 2): Define P-Proposed: A single random number generated from a Normal Distribution with Mean = P-Current and some Standard Deviation. This Normal Distribution is considered as the "Candidate Distribution" (note: I heard there is some flexibility in defining the Candidate Distribution)
• (Step 3): Define Rule: IF P-Proposed > 1 THEN P-Proposed = 1. IF P-Proposed < 0 THEN P-Proposed = 0. ELSE P-Proposed = P-Proposed
• (Step 4) Define the Acceptance Ratio: (Likelihood(P-Proposed)* Prior(P-Proposed)) / (Likelihood(P-Current) * Prior(P-Current))
• (Step 5) : Generate a single random number from a Binomial Distribution with p = min(1, Acceptance Ratio). (note: The Binomial Distribution used in Step 5 is has nothing to do with the fact that we are interested in estimating the parameter of a Binomial Distribution - this step would be present regardless of the problem)
• (Step 6) Criteria: IF Step 5 = 1 THEN P-Current = P-Proposed ELSE P-Current = P-Current

As a real example, suppose I observe 20 successes out of 50 trial. I assume a Prior Beta Distribution with alpha = beta = 2. I use a Normal Distribution as a Candidate Distribution and define the Standard Distribution as 0.1 (arbitrary choice).

Using the R programming language, I tried to implement the Metropolis-Hastings algorithm myself (5000 random samples):

    # Define the binomial likelihood function
    binom_likelihood <- function(p, y, n) {
      likelihood <- dbinom(y, size = n, prob = p)
      return(likelihood)
    }
    
    # Define the beta prior distribution
    beta_prior <- function(p, alpha, beta) {
      prior <- dbeta(p, shape1 = alpha, shape2 = beta)
      return(prior)
    }
    
    # Define the acceptance probability function
    acceptance_prob <- function(p_current, p_proposed, y, n, alpha, beta) {
      acceptance_prob <- binom_likelihood(p_proposed, y, n) * beta_prior(p_proposed, alpha, beta) / 
        (binom_likelihood(p_current, y, n) * beta_prior(p_current, alpha, beta))
      return(acceptance_prob)
    }
    
    # Define the Metropolis-Hastings algorithm
    metropolis_hastings <- function(y, n, alpha, beta, num_samples, proposal_sd) {
      p_current <- runif(1, 0, 1)
      p_samples <- numeric(num_samples)
      
      for (i in 1:num_samples) {
        p_proposed <- rnorm(1, mean = p_current, sd = proposal_sd)
        p_proposed <- ifelse(p_proposed < 0, 0, p_proposed)
        p_proposed <- ifelse(p_proposed > 1, 1, p_proposed)
        acceptance_ratio <- acceptance_prob(p_current, p_proposed, y, n, alpha, beta)
        acceptance <- rbinom(1, size = 1, prob = min(1, acceptance_ratio))
        if (acceptance == 1) {
          p_current <- p_proposed
        }
        p_samples[i] <- p_current
      }
      
      return(p_samples)
    }
    
    # Example usage
    set.seed(123)
    y <- 20
    n <- 50
    alpha <- 2
    beta <- 2
    num_samples <- 5000
    proposal_sd <- 0.1
    p_samples <- metropolis_hastings(y, n, alpha, beta, num_samples, proposal_sd)
    
    # Plot the posterior distribution of p
    hist(p_samples, main = "Posterior Distribution of p", xlab = "p", col = "gray")

I got the following output:

I am not sure if I have done this correctly - can someone please confirm?

Thanks!

Answer 1

The code appears to be correct. Note however that the step

p_proposed <- ifelse(p_proposed < 0, 0, p_proposed)
p_proposed <- ifelse(p_proposed > 1, 1, p_proposed)

is not correct, at least in spirit. If p_proposed is outside $(0,1)$, this is a zero posterior density value and hence it should be rejected, not set to the boundary value (even though the final outcome is the same since the density is zero at the boundaries).