我目前正在尝试使用PyMC来确定适合给定数据的幂律的参数。我正在使用取自的pdf公式：
A. Clauset，C。R. Shalizi和M.E. J. Newman，“权力法”
经验数据的分布”，《暹罗修订版》，第51卷，第4卷，第4页。
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
最佳答案
我认为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。）
关于python - 用PyMC拟合幂律功能，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/24402834/