MCMC使用Python的emcee對麥克斯韋曲線進行採樣

Question

我想向自我介紹MCMC與司儀採樣。我想簡單地使用github上的一組示例代碼，從Maxwell Boltzmann分佈中抽取一個樣本，https://github.com/dfm/emcee/blob/master/examples/quickstart.py。

的示例代碼是真的優良，但是當我從高斯改變分配到一個麥克斯韋，收到錯誤，類型錯誤：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功能需要，所以在你的情況下，這是您要設置您的麥克斯韋Boltxmann的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))