我正在尝试将逆伽马分布拟合到直方图。我可以通过玩弄它的两个参数（参见图 1 中的橙色线）轻松挑选出一组合适的参数。
但是，当我尝试使用 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]）。更改步数和/或烧录样本的数量无济于事。
如果我增加步行者的数量，它只会更好地采样输入位置 - 他们似乎永远不会去任何地方。
有这个帖子链接，但我认为它不适用于这里 - 逆伽马分布可以说有点烦人，但曲线拟合问题应该得到很好的定义..？
- - - - - - - - - -编辑 - - - - - - - - - -
我知道这个问题：删除flux为零的数据点会导致步行者开始移动，我得到了一个结果。我觉得这种行为有点奇怪，所以更新一下问题：
为什么两个单独的异常值会阻止司仪移动？