Реализация алгоритма мегаполиса-Хастинга - PullRequest
Новая
       10

Реализация алгоритма мегаполиса-Хастинга

0 голосов
/ 05 апреля 2019

В течение последних нескольких недель я пытался понять MCMC и Metropolis-Hastings, но у меня не получалось каждый раз, когда я пытался это реализовать.

Поэтому я пытаюсь использовать алгоритм Метрополиса-Гастингса, чтобы получить распределение Больцмана из равномерного распределения, но оно не работает.

Вот краткое изложение того, что я делаю:

  1. Я рисую случайное число из равномерного распределения m.
  2. Я рисую другое случайное число из равномерногораспределение n.
  3. Я устанавливаю dU = нм.
  4. Если dU <0, я принимаю dU, устанавливаю m = n и повторяю. </li>
  5. Если dU> 0, явычислите w = exp (-b * dU), где b равно 1 / kT, и нарисуйте случайное число из равномерного распределения r.
  6. Если w> r, я принимаю dU, установите m = n иповторение.7 Если w

Я новичок в этом поле и в 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

1 Ответ

0 голосов
/ 06 апреля 2019

Ваша проблема двоякая.Во-первых, вы выбираете новую энергию ('n = random ()') как безразмерную величину, которая противоречит всему, что вы делаете (ваша температура в Кельвинах, кБ в Дж / К и т. Д.).И, во-вторых, использование таких значений, как 10 23 и обратное, не подходит для физического моделирования - вам лучше оказаться где-то в пределах 0 ... 1 диапазона и изменить его масштаб позже.Ниже я сделал код, который работает в электрон-вольтах, сэмплировал также новую энергию в эВ и произвел нечто, напоминающее правду. 

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

enter image description here

...