MCMC 使用 Python 的司仪对麦克斯韦曲线进行采样 [英] MCMC Sampling a Maxwellian Curve Using Python's emcee

问题描述

我正在尝试向司仪介绍自己的 MCMC 抽样.我想简单地使用 github 上的一组示例代码从 Maxwell Boltzmann 分布中抽取一个样本，https://github.com/dfm/emcee/blob/master/examples/quickstart.py.

示例代码非常出色，但是当我将分布从高斯分布更改为麦克斯韦分布时，我收到错误消息，TypeError: lnprob() 只需要 2 个参数(给出 3 个)

然而，在没有给定适当参数的任何地方都不会调用它?需要一些关于如何定义麦克斯韦曲线并使其适合此示例代码的指导.

这是我所拥有的；

 from __future__ import print_function将 numpy 导入为 np进口司仪尝试:范围除了名称错误:xrange = 范围def lnprob(x, a, icov):pi = np.pi返回 np.sqrt(2/pi)*x**2*np.exp(-x**2/(2.*a**2))/a**3ndim = 2意味着 = 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)步行者 = 50p0 = [np.random.rand(ndim) for i in xrange(nwalkers)]采样器 = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[means, icov])pos, prob, state = sampler.run_mcmc(p0, 5000)采样器.reset()sampler.run_mcmc(pos, 100000, rstate0=state)

谢谢

解决方案

我认为我看到了一些问题.主要的一个是司仪希望你给它你想要采样的概率分布函数的自然对数.所以，而不是:

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

你会想要，例如

def lnprob(x, a):pi = np.pi如果 x <0:返回-np.inf别的:返回 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 的值你想用它来设置你的 Maxwell-Boltxmann 分布.这应该是单个值，而不是您在创建 mean 时设置的两个随机初始化的值.

总的来说，以下内容应该适合您:

from __future__ import print_function进口司仪将 numpy 导入为 np从 numpy 导入 pi 作为 pi# 定义麦克斯韦-玻尔兹曼分布的自然对数def lnprob(x, a):如果 x <0:返回-np.inf别的:返回 0.5*np.log(2./pi) + 2.*np.log(x) - (x**2/(2.*a**2)) - 3.*np.log(a)# 为分布选择一个值a"a = 5. # 我选择 5！# 选择步行者的数量步行者 = 50# 设置一些初始点来计算lnprobp0 = [np.random.rand(1) for i in xrange(nwalkers)]# 初始化采样器采样器 = emcee.EnsembleSampler(nwalkers, 1, lnprob, args=[a])# 运行 5000 步作为老化.pos, prob, state = sampler.run_mcmc(p0, 5000)# 重置链以移除老化样本.采样器.reset()# 从老化链的最后一个位置开始，采样 100000 步.sampler.run_mcmc(pos, 100000, rstate0=state)# 让我们检查样本是否正确mbmean = 2.*a*np.sqrt(2./pi) # Maxwell-Boltzmann 分布的均值打印(样本均值 = {}，分析均值 = {}".格式(np.mean(sampler.flatchain[:,0])，mbmean))mbstd = np.sqrt(a**2*(3*np.pi-8.)/np.pi) # std.开发M-B 分布打印(样本标准偏差 = {}，分析 = {}".格式(np.std(sampler.flatchain[:,0])，mbstd))

I am trying to introduce myself to MCMC sampling with emcee. I want to simply take a sample from a Maxwell Boltzmann distribution using a set of example code on github, https://github.com/dfm/emcee/blob/master/examples/quickstart.py.

The example code is really excellent, but when I change the distribution from a Gaussian to a Maxwellian, I receive the error, TypeError: lnprob() takes exactly 2 arguments (3 given)

However it is not called anywhere where it is not given the appropriate parameters? In need of some guidance as to how to define a Maxwellian Curve and have it fit into this example code.

Here is what I have;

    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)

Thanks

解决方案

I think there are a couple of problems that I see. The main one is that emcee wants you to give it the natural logarithm of the probability distribution function that you want to sample. So, rather than having:

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

you would instead want, e.g.

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)

where the if...else... statement is to explicitly say that negative values of x have zero probability (or -infinity in log-space).

You also shouldn't have to calculate icov and pass it to lnprob as that's only needed for the Gaussian case in the example you link to.

When you call:

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

the args value should just be any additional arguments that your lnprob function requires, so in your case this would be the value of a that you want to set your Maxwell-Boltxmann distribution up with. This should be a single value rather than the two randomly initialised values you have set when creating mean.

Overall, the following should work for you:

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

这篇关于MCMC 使用 Python 的司仪对麦克斯韦曲线进行采样的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！