如何在Python中正确地将数据拟合到幂律分布？

Question

我正在考虑《白鲸记》小说中独特单词出现次数，并使用powerlaw Python包将单词频率拟合到幂律分布上。

我不确定为什么我无法复制Clauset等人之前的工作结果，因为p值和KS得分都很“差”。

这个想法是将独特单词的频率拟合成幂律分布。然而，由scipy.stats.kstest计算的Kolmogorov-Smirnov拟合优度检验看起来很糟糕。

我有以下函数将数据拟合到幂律分布上：

import numpy as np
import powerlaw
import scipy
from scipy import stats

def fit_x(x):
    fit = powerlaw.Fit(x, discrete=True)
    alpha = fit.power_law.alpha
    xmin  = fit.power_law.xmin
    print('powerlaw', scipy.stats.kstest(x, "powerlaw", args=(alpha, xmin), N=len(x)))
    print('lognorm', scipy.stats.kstest(x, "lognorm", args=(np.mean(x), np.std(x)), N=len(x)))

---

下载赫尔曼·梅尔维尔(Herman Melville)小说《白鲸》中独特单词的频率（根据Aaron Clauset等人的说法，应该遵循幂律分布）:

wget http://tuvalu.santafe.edu/~aaronc/powerlaws/data/words.txt

Python脚本：

x =  np.loadtxt('./words.txt')
fit_x(x)

结果：

结果：

('powerlaw', KstestResult(statistic=0.862264651286131, pvalue=0.0))
('log norm', KstestResult(statistic=0.9910368602492707, pvalue=0.0))

---

当我比较期望结果并在相同的Moby Dick数据集上遵循这个 R教程时，我得到了一个不错的p值和KS检验值：

library("poweRlaw")
data("moby", package="poweRlaw")
m_pl = displ$new(moby)
est = estimate_xmin(m_pl)
m_pl$setXmin(est)
bs_p = bootstrap_p(m_pl)
bs_p$p
## [1] 0.6738

---

当计算KS检验值并通过powerlaw Python库进行后处理时，我可能会错过什么？PDF和CDF看起来没问题，但KS检验结果似乎有误。

- 0x90

---

空手俱乐部网络规模相对较小。我认为你在其中没有足够的样本来做出非常强有力的统计陈述。 - Paul Brodersen

@PaulBrodersen 根据作者提供的网站链接：https://arxiv.org/abs/0706.1062 《实证数据中的幂律分布》一文讨论了“常用于分析幂律数据的方法，如最小二乘拟合，可能会产生相当不准确的幂律分布参数估计，即使在这些方法返回准确答案的情况下，它们仍然不令人满意，因为它们没有表明数据是否符合幂律分布。” 然而，我不确定这里的结果是怎么回事。 - Scott Staniewicz

这是一篇来自一家声誉良好的期刊的最新文章，阐述了幂律采样的结果。然而，这个结果本身已经很老了（出于某种原因，它仍然不太为人所知）。 - Paul Brodersen

我并不是在试图解决您拟合问题的指向，只是在引用定理（pdf第2页）：“尽管它们被称为幂律或无标度分布，但它们在子采样下并不具备不变性。换句话说，如果P(s)遵循幂律分布，则P_{sub}(s)不是一个幂律分布[...]”我认为这个定理对于许多发表在文献中的结果的解释有重要意义（通常相当毁灭性），应该警告您不要过度解释自己数据中可能出现的任何无标度性质。 - Paul Brodersen

为适应您的标量和指数，我会在双对数轴上绘制度数，并从图像的直线部分估算参数。 - Paul Brodersen

Answer 1

我认为你应该注意数据是连续还是离散的，然后选择适当的测试方法；此外，正如前面所说，数据的大小将对结果产生一定影响，希望这能帮到你。

- 鲍业文

Answer 2

我仍然不清楚如何使用scipy.stats.kstest和powerlaw库来确定显著性和拟合度。

虽然，powerlaw实现了自己的distribution_compare功能，该功能返回似然比R和R的p-val（请参见Aaron Clauset在此处的一些内容）：

R：浮点数 两个分布对数据的拟合的对数似然比。如果大于0，则更喜欢第一个分布。如果小于0，则更喜欢第二个分布。 p：浮点数 R的显著性

from numpy import genfromtxt
import urllib
import powerlaw

urllib.urlretrieve('https://raw.github.com/jeffalstott/powerlaw/master/manuscript/words.txt', 'words.txt')
words = genfromtxt('words.txt')

fit = powerlaw.Fit(words, discrete=True)

print(fit.distribution_compare('power_law', 'exponential', normalized_ratio=True))
(9.135914718776998, 6.485614241379581e-20)
print(fit.distribution_compare('power_law', 'truncated_power_law'))
(-0.917123083373983, 0.1756268316869548)
print(fit.distribution_compare('power_law', 'truncated_power_law'))
(-0.917123083373983, 0.1756268316869548)
print(fit.distribution_compare('power_law', 'lognormal'))
(0.008785246720842022, 0.9492243713193919)

- 0x90