MCMC Выборка кривой Максвелла с использованием ведущего Python

Question

Я пытаюсь познакомиться с MCMC с помощью emcee. Я хочу просто взять пример из дистрибутива Максвелла Больцмана, используя набор примеров кода на github, https://github.com/dfm/emcee/blob/master/examples/quickstart.py.

Пример кода действительно отличный, но когда я меняю распределение с гауссовского на максвелловское, я получаю ошибку, TypeError: lnprob () принимает ровно 2 аргумента (дано 3)

Однако он не вызывается нигде, где ему не даны соответствующие параметры? Необходимы некоторые указания относительно того, как определить кривую Максвелла и как она вписывается в этот пример кода.

Вот что у меня есть;

    from __future__ import print_function
import numpy as np
import emcee

try:
    xrange
except NameError:
    xrange = range
def lnprob(x, a, icov):
    pi = np.pi
    return np.sqrt(2/pi)*x**2*np.exp(-x**2/(2.*a**2))/a**3

ndim = 2
means = np.random.rand(ndim)

cov  = 0.5-np.random.rand(ndim**2).reshape((ndim, ndim))
cov  = np.triu(cov)
cov += cov.T - np.diag(cov.diagonal())
cov  = np.dot(cov,cov)


icov = np.linalg.inv(cov)


nwalkers = 50


p0 = [np.random.rand(ndim) for i in xrange(nwalkers)]


sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[means, icov])

pos, prob, state = sampler.run_mcmc(p0, 5000)

sampler.reset()

sampler.run_mcmc(pos, 100000, rstate0=state)

Спасибо

Answer 1

Я думаю, что есть пара проблем, которые я вижу. Главный из них заключается в том, что ведущий хочет, чтобы вы задали натуральный логарифм функции распределения вероятностей, которую вы хотите выбрать. Итак, вместо того, чтобы:

def lnprob(x, a, icov):
    pi = np.pi
    return np.sqrt(2/pi)*x**2*np.exp(-x**2/(2.*a**2))/a**3

вы бы вместо этого хотели, например

def lnprob(x, a):
    pi = np.pi
    if x < 0:
        return -np.inf
    else:
        return 0.5*np.log(2./pi) + 2.*np.log(x) - (x**2/(2.*a**2)) - 3.*np.log(a)

где оператор if... else... должен явно сказать, что отрицательные значения x имеют нулевую вероятность (или -infinity в лог-пространстве).

Вы также не должны рассчитывать icov и передать его lnprob поскольку это необходимо только для случая Гаусса в примере, на который вы ссылаетесь.

Когда вы звоните:

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[means, icov])

args значение должно быть просто дополнительные аргументы, которые ваш lnprob функция требует, так что в вашем случае это будет значение a что вы хотите установить ваше распределение Максвелла-Больцмана с помощью. Это должно быть одно значение, а не два случайно инициализированных значения, которые вы установили при создании mean,

В целом, следующее должно работать для вас:

from __future__ import print_function

import emcee
import numpy as np
from numpy import pi as pi

# define the natural log of the Maxwell-Boltzmann distribution
def lnprob(x, a):               
    if x < 0:                                                                
        return -np.inf
    else:
        return 0.5*np.log(2./pi) + 2.*np.log(x) - (x**2/(2.*a**2)) - 3.*np.log(a)

# choose a value of 'a' for the distributions
a = 5. # I'm choosing 5!

# choose the number of walkers
nwalkers = 50

# set some initial points at which to calculate the lnprob
p0 = [np.random.rand(1) for i in xrange(nwalkers)]

# initialise the sampler
sampler = emcee.EnsembleSampler(nwalkers, 1, lnprob, args=[a])

# Run 5000 steps as a burn-in.
pos, prob, state = sampler.run_mcmc(p0, 5000)

# Reset the chain to remove the burn-in samples.
sampler.reset()

# Starting from the final position in the burn-in chain, sample for 100000 steps.
sampler.run_mcmc(pos, 100000, rstate0=state)

# lets check the samples look right
mbmean = 2.*a*np.sqrt(2./pi) # mean of Maxwell-Boltzmann distribution
print("Sample mean = {}, analytical mean = {}".format(np.mean(sampler.flatchain[:,0]), mbmean))
mbstd = np.sqrt(a**2*(3*np.pi-8.)/np.pi) # std. dev. of M-B distribution
print("Sample standard deviation = {}, analytical = {}".format(np.std(sampler.flatchain[:,0]), mbstd))