0

New to Python. I am trying to port the R code for the Metropolis Hastings Algorithm found here over to Python. I have successfully duplicated the results in R, but am struggling with the Python version. I just can't seem to obtain anything close to the R results.

R code:

target <- function(x) {
  return(ifelse(x < 0, 0, exp(-x)))
}

x <- rep(0, 10000)
x[1] <- 3    #initialize; I've set arbitrarily set this to 3
for (i in 2:10000){
  current_x <- x[i - 1]
  proposal <- rnorm(n = 1, mean = 0, sd = 1)
  proposed_x <- current_x + proposal
  A <- target(proposed_x) / target(current_x) 
  if (runif(1) < A){
    x[i] <- proposed_x       # accept move with probabily min(1,A)
    } 
    else 
    {
    x[i] <- current_x        # otherwise "reject" move, and stay where we are
    }

}

hist(x, xlim = c(0, 10), probability = TRUE, main = "Histogram of values of x visited by MH algorithm")

Python code:

import math
import numpy as np
import matplotlib.pyplot as plt
import random as random

def target(x):
  return(np.where(x < 0, 0, math.exp(-x)))

list_x = [3] # start with current_x = 3

for i in range(1,10000):
  current_x = list_x[-1] # pull last item from list_x
  proposal = np.random.normal(0,1)
  proposed_x = current_x + proposal
  alpha = target(proposed_x)/target(current_x)
  if min(1,alpha) < 1:
  #if np.random.uniform(0,1) < alpha:  #min = 0, max = 1, sample size = 1
    list_x.append(proposed_x)   # accept and append proposed_x to list_x
  else: 
    list_x.append(current_x)    # reject and append current_x to list_x
 
# print(list_x)
# plt.xlim([0,10])
plt.hist(list_x, edgecolor='black')
plt.show()
user1884367
  • 425
  • 2
  • 7
  • 15

1 Answers1

0

I believe you have a bug in the sampling

Corrected code

import numpy as np
import matplotlib.pyplot as plt

def target(x):
    return(np.where(x < 0, 0, np.exp(-x)))

list_x = [3] # start with current_x = 3

for i in range(1,10000):
    current_x = list_x[-1] # pull last item from list_x
    proposal = np.random.normal(0,1)
    proposed_x = current_x + proposal
    alpha = target(proposed_x)/target(current_x)
    if np.random.uniform() < alpha:
        list_x.append(proposed_x)   # accept and append proposed_x to list_x
    else:
        list_x.append(current_x)    # reject and append current_x to list_x

plt.hist(list_x, edgecolor='black')
plt.show()

display something like graph below, looks exactly like exponential distribution

enter image description here

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64