2

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

A. Clauset、CR Shalizi 和 MEJ Newman，“经验数据中的幂律分布”，Siam rev.，第一卷。51 号 4，第 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

4

1 回答 1

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)

（有关我喜欢如何使用老化和图形收敛检查的 PyMC3 示例，请参阅此 SO 答案。）

于 2014-06-25T15:53:16.067 回答