# Writing MCMC Sampling Code by Hand

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!

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

• @ Xi'an: Thank you so much for your answer! I will correct the Posterior and make sure that the likelihood only appears once Commented Jul 30, 2023 at 7:08
• 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! Commented Jul 30, 2023 at 7:14
• @ Xi'an: I am also attempting to learn how to simulate (Continuous Time) Markov Chains by hand: stats.stackexchange.com/questions/622408/… Commented Aug 1, 2023 at 2:46
• 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. Commented Aug 1, 2023 at 10:59
• Thank you for your reply! Is there a way I can improve my approach to make it closer to the standard approach? Commented Aug 1, 2023 at 15:11