MCMC Pythonのemceeを使ったMaxwellian Curveのサンプリング

Question

私はemceeでMCMCサンプリングを紹介しようとしています。 githubのサンプルコードを使って、Maxwell Boltzmannのディストリビューションからサンプルを取ってみたい、https://github.com/dfm/emcee/blob/master/examples/quickstart.py。

例コードが本当に優れているが、私はマクスウェルにガウスの分布を変更すると、私はエラーを受け取り、TypeError例外：lnprob（）は、正確に2つの引数（3所与）

を要するが、それは適切なパラメータが与えられていない場所では呼び出されませんか？ Maxwellian Curveを定義し、このサンプルコードに適合させる方法についていくつかの指針が必要です。

これは私が持っているものです。私は私が見る問題のカップルがあると思います

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

。主なものは、あなたがサンプリングしたい確率分布関数の自然対数を与えることです。たとえば、以下のように指定するのではなく、

のようにすることができます。

...他にあれば...文は、明示的に xの負の値は、（対数空間でまたは-Infinity）ゼロの確率を持っていると言うことです

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) 

。

icovを計算して、lnprobに渡す必要はありません。これは、リンク先の例のガウス関数の場合にのみ必要です。

お問い合わせの際：

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

をお使いの場合には、これはあなたがあなたのマックスウェル・Boltxmannを設定したいaの値になりますのでargs値だけで、あなたのlnprob機能が必要であることを任意の追加の引数でなければなりませんとの分配。これは、meanを作成するときに設定した2つのランダムに初期化された値ではなく、1つの値にする必要があります。

全体的に、次はあなたのために働く必要があります。

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))