Metropolis-Hastings算法从R迁移至Python的结果不符问题求助

我是Python新手，正在将Metropolis-Hastings算法的R代码移植到Python。R代码能成功复现预期结果，但Python版本始终无法得到相近结果。以下是原R代码和我的Python代码：

原R代码 #

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代码 #

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

问题原因与修正代码 #

你的Python代码核心错误在于接受-拒绝逻辑完全错误：你注释掉了正确的np.random.uniform(0,1) < alpha判断，换成了if min(1,alpha) < 1——这个逻辑会在alpha小于1时直接接受，alpha大于等于1时直接拒绝，完全违背了Metropolis-Hastings算法的随机接受规则。

此外还有两个小细节需要调整：

1. target函数中用np.exp代替math.exp，避免标量/数组兼容问题；
2. 直方图设置density=True（对应R的probability=TRUE）和xlim=(0,10)，和R的可视化保持一致。

修正后的Python代码：

import numpy as np
import matplotlib.pyplot as plt

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

list_x = [3]  # 初始值设为3

for _ in range(9999):  # 循环9999次，得到10000个样本
    current_x = list_x[-1]
    proposal = np.random.normal(0, 1)
    proposed_x = current_x + proposal
    alpha = target(proposed_x) / target(current_x)
    # 正确的接受-拒绝判断：生成0-1随机数，与min(1, alpha)比较
    if np.random.uniform(0, 1) < min(1, alpha):
        list_x.append(proposed_x)
    else:
        list_x.append(current_x)

# 绘制与R一致的直方图
plt.hist(list_x, edgecolor='black', density=True, range=(0,10))
plt.title("Histogram of values of x visited by MH algorithm")
plt.show()

解释 #

修正后的代码遵循了Metropolis-Hastings算法的核心规则：

• 计算接受概率alpha = target(proposed_x)/target(current_x)，取min(1, alpha)确保概率不超过1；
• 生成0到1之间的均匀随机数，若该数小于接受概率，则接受新样本，否则保留当前样本；
• 可视化参数与R对齐，方便对比结果。

内容的提问来源于stack exchange，提问作者user1884367