实施metropolis-hasting算法

Question

在过去的几周里，我一直在努力了解MCMC和Metropolis-Hastings，但每次尝试实施它都失败了。

所以我试图使用metropolis-Hastings算法从均匀分布中得到Boltzmann分布，但它不起作用。

以下是我正在做的事情的摘要：

1. 我从均匀分布m中绘制一个随机数。
2. 我从均匀分布n中绘制另一个随机数。
3. 我设置dU = n-m。
4. 如果dU <0，我接受dU，设置m = n，然后重复。
5. 如果dU> 0，我计算w = exp（-b * dU），其中b是1 / kT，并从均匀分布r中绘制一个随机数。
6. 如果w> r，我接受dU，设置m = n，然后重复。 7如果w <r，我拒绝dU，设置m = m，然后重复。

我是这个领域和python的初学者，所以我不确定代码是错还是算法错了（可能两者都有。）

我的代码附在下面。谢谢。

import random
%matplotlib inline
import numpy as np
import scipy
import matplotlib.pyplot as plt
from scipy import stats
k = 1.38064852 * 10**(-23)
t = 298
b = 1/(t*k)   U = []
m = np.random.uniform(0, 1)
for j in range(100000):
    n = np.random.uniform(0, 1)
    du = n-m
    if du<0:
            U.append(du)
            m = n 
    elif du > 0:
            w = np.exp(-b*du)
            r = np.random.uniform(0, 1)
            if w > r:
                    U.append(du)
                    m = n
            else:
                    U.append(du)
                    m = m

Answer 1

你的问题是双重的。一个是你将新能量（'n = random（）'）作为无量纲量进行采样，这与你正在做的其他事情相矛盾（你的温度是开尔文，kB是J / K等）。第二，使用像1023和反向值这样的值在物理模拟中并不好 - 你最好在0 ... 1范围内的某个地方，然后重新调整。下面我制作的代码在电子伏特中起作用，也在eV中采样新的能量，并产生类似于真相的东西。

import numpy as np
import matplotlib.pyplot as plt

kB = 1.0/11600. # eV/K
T  = 300 # K
b  = 1.0/(kB * T) # inverse temperature, eV^-1

np.random.seed(76543217) # for reproducibility

N = 100000
EE = np.empty(N+1) # energy
DE = np.empty(N+1) # delta energy

Ei = 1.0 # initial energy, 1eV
EE[0] = Ei
DE[0] = 0.0
for k in range(N):
    E  = np.random.random()/b # sample energy, in eV
    dE = E - Ei
    if dE < 0.0:
        Ei = E
        EE[k+1] = Ei
        DE[k+1] = dE
    elif dE > 0.0:
        w = np.exp(-b*dE)
        r = np.random.random()
        if w > r:
            Ei = E
            EE[k+1] = Ei
            DE[k+1] = dE
        else:
            EE[k+1] = Ei
            DE[k+1] = 0.0

x = np.linspace(0, N+1, num=N+1)

print(EE[N-30:])
print(np.mean(EE[N-1000:]))
print(np.mean(EE[N-2000:]))

fig, ax = plt.subplots(1, 1)

ax.plot(x[N-1000:], EE[N-1000:], 'r-', lw=5, alpha=0.6, label='Energy')
ax.plot(x[N-1000:], DE[N-1000:], 'go', lw=5, alpha=0.6, label='Delta Energy')
plt.show()

打印是1000和2000样品的两个最后平均值，看起来热化给我

0.010188070423940562
0.010666101150488673

和E / dE的图表