我在R到python的翻译中做错了什么？

Question

我在R中有一个代码，它实现了大都会黑斯廷斯算法：

trials <- 100000
sim <- numeric(trials)
sim[1] <- 2
for (i in 2: trials) {
   old <- sim[i-1]
   prop <- runif(1,0,5)
   acc <- (exp(-(prop-1)^2/2) + exp(-(prop-4)^2/2)) /
     ( (exp(-(old-1)^2/2) + exp(-(old-4)^2/2)) )
   if (runif(1) < acc) 
     sim[i] <- prop
   else
     sim[i] <- old }
mean(sim)
var(sim)

结果是正确的。

但是当我用Python翻译它时，结果是不同的。

trials = 100000
sim = np.repeat(0,trials+1)
sim[0] = 2
for i in range (2, trials):
    old = sim[i-1]
    prop = np.random.uniform(0,5,1)
    acc = (np.exp(-(prop-1)**2/2) + np.exp(-(prop-4)**2/2)) /( (np.exp(-(old-1)**2/2) + np.exp(-(old-4)**2/2)) )
    if np.random.uniform(1) < acc:
        sim[i] = prop
    else:
        sim[i] = old 

为什么？我在这里做错了什么？

Answer 1

首先，您希望在i=1在range(1, trials)中启动python循环，因为0是从0开始的。

其次，在if np.random.uniform(1)中，您只是在生成1.0，所以这需要更改，例如，如果您只想要一个0到1之间的统一数字，则需要更改np.random.uniform(min, max, size=1)或np.random.uniform(size=1)。如果不清楚，请查看文档 for np.random.uniform。

更新

第三，你是在不知不觉中转换成整数，所以这也需要处理。当习惯于R时，这很容易被忘记(我只是做了自己的事情)。我在下面重构了您的代码，这个解决方案应该为您提供一个与您在R中看到的类似的结果。

这里，我将sim转换成一个float向量，并确保使用exp-functions中的浮点数减去和除以。希望这对你有用。

trials = 100000
sim = np.repeat(0,trials).astype(np.float64)
sim[0] = 2.0
for i in range (1, trials):
    old = sim[i-1]
    prop = np.random.uniform(0,5,1)
    acc = (np.exp(-(prop-1.0)**2/2.0) + np.exp(-(prop-4.0)**2/2.0)) \
    / ( (np.exp(-(old-1.0)**2/2.0) + np.exp(-(old-4.0)**2/2.0)) )
    if np.random.uniform(size=1) < acc:
        sim[i] = prop
    else:
        sim[i] = old