使用PyMC拟合幂律函数
我现在正在尝试使用PyMC来确定给定数据的幂律拟合参数。我使用的概率密度函数（pdf）公式来自以下文献：
A. Clauset, C. R. Shalizi, 和 M. E. J. Newman, "经验数据中的幂律分布," Siam rev., vol. 51, iss. 4, pp. 661-703, 2009.
为了生成具有特定参数的样本数据来测试我的代码，我使用了以下Python的powerlaw包，它实现了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
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)