emcee样本初始位置

Question

我试着把逆伽马分布拟合成直方图。我可以很容易地选择一组好的拟合参数，通过使用它的两个参数(c.f )。(Fig1中的橙色行)。

然而，当我试图用emcee估计后验时，我只是简单地抽样了链的初始起点。

import numpy as np 
import corner 
import emcee 
from scipy.special import gamma as gamma_function
import matplotlib.pyplot as plt

def _inverse_gamma_distribution(x0, alpha, beta):
    return beta**alpha/gamma_function(alpha)* (1/x0)**(alpha + 1)* np.exp(-beta/x0)


flux =          np.array([0.,            0.,         0.08753462, 0.48757902, 0.59385076, 0.32146012, 0.20280991, 0.06542532, 0.01888836, 0.00369042, 0.00133481,  0.,            0.,         0.10504154, 0.44777665])
bin_center =    np.array([0.06898463,    0.12137053, 0.21353749, 0.37569469, 0.66099163, 1.16293883, 2.04605725, 3.59980263, 6.33343911, 11.1429583, 19.60475495, 0.06898463,    0.12137053, 0.21353749, 0.37569469])
error =         np.array([0.,            0.,         0.03914667, 0.06965415, 0.0579539,  0.03214601, 0.01924986, 0.00824282, 0.00333902, 0.0011127,  0.0005045,   0.,            0.,         0.04288303, 0.0667506 ])


def lnprior(theta):
    alpha, beta         = theta
    if 1.0 < alpha < 2.0 and 1.0 < beta < 2.0:
        return 0.0
    return -np.inf

def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)

def lnlike(theta, x, y, yerr):
    alpha, beta         = theta
    model               = np.array(_inverse_gamma_distribution(x, alpha, beta))
    return -0.5*np.sum((y - model)**2/yerr**2)

p0                      = [1.5, 1.5]

ndim, nwalkers          = 2, 100
pos                     = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)]

sampler                 = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(bin_center, flux, error))
sampler.run_mcmc(pos, 1000)
samples                 = sampler.chain[:, 500:, :].reshape((-1, ndim))

print("mean posterior: ", samples.T[0].mean(), samples.T[1].mean()) 
print("std posterior: ",  samples.T[0].std(),  samples.T[1].std())

fig                     = corner.corner(samples)

fig, ax                 = plt.subplots()
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip')
ax.set_ylim([0.0001, 10])
x                       = np.linspace(0.1, 20, 1000)

ax.errorbar(bin_center, flux, yerr=error, fmt=".")
ax.plot(x, _inverse_gamma_distribution(x, *p0))
for alpha, beta in samples[np.random.randint(len(samples), size=100)]:
    ax.plot(x, _inverse_gamma_distribution(x, alpha, beta), "-", color="grey", alpha=0.1)
plt.show()

后验均值和标准差对应于我在pos = [np.array(p0) + 5e-1*np.random.randn(ndim) for i in range(nwalkers)]中设定的值(即p0 +- 0.5)。改变步骤的数量或燃烧样本的数量都没有帮助。

​

​

如果我增加步行者的数量，它只是更好地采样输入位置-他们似乎永远不会去任何地方。

​

​

有一个帖子link，但我不认为它适用于这里-逆伽马分布可以说是有点烦人，但曲线拟合问题应该是最好的定义。

-------------------EDIT--------------------

我知道问题所在:删除flux为零的数据点会导致步行者开始移动，我得到了结果。我觉得这个行为有点奇怪，所以要更新这个问题：为什么两个孤立点阻止emcee移动？

Answer 1

步行者从未开始移动的原因在于错误数组：

error = np.array([0., 0., 0.03914667, 0.06965415, ... ])

error的前两个值是0.。因此，似然函数-0.5*np.sum((y - model)**2/yerr**2)对于(alpha, beta)的任何值都是-np.inf，而步行者永远不会开始移动。移除这些值或将错误设置为任何非零值，让我们放松一下，fit就会像它应该的那样迅速收敛！