用 PyMC 拟合幂律函数

Question

我目前正在尝试使用 PyMC 来确定适合给定数据的幂律参数。我使用的 pdf 公式取自：

A. Clauset、C. R. Shalizi 和 M. E. J. Newman，“幂律 经验数据中的分布，”Siam rev., vol. 51, iss. 4, pp. 661-703，2009 年。

为了生成具有特定给定参数的示例数据以测试我的代码，我使用了以下 Python 幂律包，该包实现了 Clauset 等人的方法：

https://pypi.python.org/pypi/powerlaw

如果我使用固定的 xmin 值（即幂律函数的下限），我的代码运行良好。但是，一旦我尝试确定 xmin 值，输出就会产生过高的 xmin 值。我已将相应的 xmin 部分注释掉：

test = powerlaw.Power_Law(xmin = 1., parameters = [1.5])
simulated = test.generate_random(1000)
fit = powerlaw.Fit(simulated, xmin=1.)
print fit.alpha
print fit.xmin

xmin = 1.

#alpha = mc.Uniform('alpha', 0,6, value=1.5) 
alpha = mc.Exponential('alpha', 1.5)
#xmin = mc.Uniform('xmin', min(simulated), max(simulated), value=min(simulated))
#xmin = mc.Exponential('xmin', 1.)
#print xmin.value

@mc.stochastic(observed=True)
def power_law(value=simulated, alpha=alpha, xmin=xmin):
    #value = value[value >= xmin]
    return np.sum(np.log((alpha-1) * xmin**(alpha-1) * value**-alpha))

model = mc.MCMC([alpha,xmin,power_law])
model.sample(iter=5000)

print(model.stats()['alpha']['mean'])
#print(model.stats()['xmin']['mean'])

alpha_samples = model.trace('alpha')[:]
#xmin_samples = model.trace('xmin')[:]

figsize(12.5,10)

ax = plt.subplot(311)
ax.set_autoscaley_on(False)

plt.hist(alpha_samples, histtype='stepfilled', bins=20, label="posterior of alpha", color="#A60628", normed=True)
plt.legend(loc="upper left")
plt.xlim([0,2])
plt.xlabel("alpha value")

#plt.subplot(312)

#plt.hist(xmin_samples, histtype='stepfilled', bins=20, 
#         label="posterior of xmin", color="#A60628", normed=True)
#plt.legend(loc="upper left")
#plt.xlim([0,500])
#plt.xlabel("xmin value")

我认为一个问题是我应该始终只考虑我的 power_law 函数中的数据 >= xmin。如果这样做，当我还确定 xmin 时，我会得到“正确”的 alpha 值，但 xmin 仍然太高。我也觉得这是一个不公平的比较，因为要在 MCMC 过程中查看的数据样本大小不同，因此可能性比较也存在偏差。

也许有人知道我该如何处理这个问题。

更新：我当前的代码在这里可用： http://www.philippsinger.info/notebooks/pl_pymc.html

Answer 1

我认为当xmin 小于某些数据值时，您的可能性存在问题是正确的。我的解决方案是明确禁止这种情况，当它出现时返回 -np.inf 的对数似然：

@mc.stochastic(observed=True)
def power_law(value=simulated, alpha=alpha, xmin=xmin):
    if value.min() < xmin:
        return -np.inf
    return np.sum(np.log((alpha-1) * xmin**(alpha-1) * value**-alpha))

我还建议使用总样本一半的老化期，并以图形方式检查收敛性，如下所示：

model.sample(iter=5000, burn=2500)
pm.Matplot.plot(model)

（见this SO answer for a PyMC3 example of how I like to use burn-in and graphical convergence checks。）