1
$\begingroup$

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:

enter image description here

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

Thanks!

$\endgroup$

1 Answer 1

3
+50
$\begingroup$

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

$\endgroup$
5
  • 1
    $\begingroup$ @ Xi'an: Thank you so much for your answer! I will correct the Posterior and make sure that the likelihood only appears once $\endgroup$ Jul 30, 2023 at 7:08
  • $\begingroup$ Can you please recommend a way to fix this mistake without adding more mistakes to the code? I am worried that if I try to change this, I might add more mistakes to the code :( Thank you so much for all your help! $\endgroup$ Jul 30, 2023 at 7:14
  • $\begingroup$ @ Xi'an: I am also attempting to learn how to simulate (Continuous Time) Markov Chains by hand: stats.stackexchange.com/questions/622408/… $\endgroup$ Aug 1, 2023 at 2:46
  • $\begingroup$ You do not have to modify anything because the target distribution is equal to zero at both $p=0,1$, hence the proposed value is never accepted. $\endgroup$
    – Xi'an
    Aug 1, 2023 at 10:59
  • $\begingroup$ Thank you for your reply! Is there a way I can improve my approach to make it closer to the standard approach? $\endgroup$ Aug 1, 2023 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.