关于幂律分布，你还应该知道如何用代码实现！| 集智百科
今天我们继续学习幂律分布的基本概念——幂律概率分布，以及如何用代码实现幂律分布。内容来自集智百科，集智百科是复杂系统领域的百科全书，涵盖复杂系统领域的基本概念（持续完善中）。
我们正在组织撰写翻译相应的维基词条，并附上代码实现。想要自己创建词条，一起贡献知识的小伙伴们可以通过链接报名哦。点击「编辑」，做些改变，按下「保存」，你将影响世界！
幂律概率分布
（Power-law probability distributions）
广义上，幂律概率分布是一个密度函数（或离散情况下的概率质量函数）具有以下形式的分布： 对于较大的
，
，其中
且
是一个慢变函数(Slowly varying function)，对于任何正因子
，它都满足
的这个属性来自于
渐进的标度不变性。因此，
仅控制左尾的形状和有限范围。如果
是常量因子函数，并且我们有一个幂律适用于所有的
值，在许多情况下，可以很容易地依据幂律假设出一个下限。结合这两种情况，当
是一个连续变量，幂律有以下形式：
其中，
是标准化常量因子
下面我们来讨论这个分布的性质。 首先，它的矩可表示为：
当
，定义是完备的； 当
，发散： 当
，均值与高阶矩都是无穷大； 当
，均值存在，但方差和高阶矩都是无穷大。 如果从这种分布中抽取有限样本，意味着中心矩估计永远不会收敛——并且随着数据的增多，他们还有增大的趋势。这种幂律概率分布又被称为帕累托型分布，具有帕累托尾部特征的分布，或是具有规则变化的分布。
帕累托型分布 |来源：https://www.google.com
一种不满足上面的一般形式的修改，即指数截止幂律分布。
在这种分布中，指数衰减项
最终会在较大的处超过正常的幂律分布。这种分布无法成比例缩放，因此并不是幂律;不过，它会在截止前的有限区域内近似地缩放。（注意，一般的幂律分布是这种分布的简单形式，即
的指数截止幂律分布。）这种分布是渐近幂律分布的常见替代方法，因为它考虑了有限大小的影响。
Tweedie分布是一族统计模型，其特征是基于可加（additive）与可再生（reproductive）卷积以及标度变换（scale transformation）的闭包（closure）。因此，这些模型都表达了方差和均值之间的幂律关系。这些模型作为数学收敛的焦点，类似于正态分布在中心极限定理中所扮演的角色。这种收敛效应解释了为什么在自然过程中, 方差-平均幂律表现得如此广泛, 就像泰勒在生态学中的定律和在物理学中的涨落标度。还可以证明，使用扩展箱（ expanding bins）的方法时，这种方差- 均值幂律分布（variance-to-mean power law）意味着存在1 /f噪声，而1/ f噪声可能是由于Tweedie收敛效应（Tweedie convergence effect）而产生的。
图形验证法
检验幂律概率分布
虽然人们已经提出了更成熟更稳健的方法，但通过随机样本检验幂律概率分布的最常用的图形方法还是帕累托双分位图（Pareto quantile-quantile plots ）（或帕累托Q-Q图），平均剩余寿命图（mean residual life plot）和双对数图（Paretoquantile-quantile plots）。
（log-log图）
在双对数图上呈现直线是必要的，但并没有足够的证据证明他们的幂律关系，直线的斜率就对应于幂律指数。
另一种更强大的图形检验法是利用bundles of residual quantile functions 残余分位函数束。（注意，幂律分布也称为帕累托分布。）这里假设从概率分布中获得随机样本，并且我们想知道分布的尾部是否遵循幂律（换句话说，我们想知道分布是否有“帕累托尾”）。此处随机样本也被称为“数据”。
帕累托Q-Q图是这样绘制的：它将取对数后（样本）数据的分位数与取均值为1的指数分布对应的分位数（或标准帕累托分布的位数）进行比较。如果得到的散点图表现是“渐近收敛”为直线，就应该怀疑其服从幂律分布。帕累托 Q-Q图的局限是它在尾部指数（也称为帕累托指数）接近于0时表现不佳，因为帕累托Q-Q图难以检验尾部是缓慢变化的分布。
帕累托分布 |来源：https://www.google.com
另一种检验幂律概率分布的方法是平均剩余寿命图，它包含以下步骤：首先对数据取对数，然后将高于第i阶统计量的数据平均值与第i阶统计量进行比较绘制，从i=1,…,n,其中n是随机样本容量。如果绘制出的散点图走势呈现为一条“稳定”的水平直线，那么应该考虑其服从幂律分布。但由于平均剩余寿命图对异常值非常敏感（它并不稳健），所以它通常会产生一些难以解释的图形;而这些图形通常被称为 Hill horror plots 。
Hill horror plots |来源：https://www.google.com
双对数图是使用随机样本以图形方式检验尾部分布的另一种方式。使用这个方法必须要谨慎，因为双对数图中呈现直线对幂律概率分布是必要不充分条件，许多非幂律分布在双对数图上也显示为直线。这个方法是将特定数在该分布中的概率估计量的对数 | 对比这个数的对数 | 进行绘图。通常，此估计量是该数据在数据集中出现的次数的比例。如果图中的点在x较大时倾向于“收敛”为直线，则可得出结论，该分布具有“幂律尾”（power-law tail）。目前这些类型的绘图的应用示例已经发表。但这种方法的局限是，需要大量的数据才能使结果可靠。此外，它仅适用于离散（或分组）数据。
log-log Graph 双对数曲线图 | 来源：https://www.google.com
不过，目前已经提出了使用随机样本检验幂律概率分布的另一种图形方法。该方法包括绘制对数变换样本的束，是最早提出使用随机样本探索矩的存在和矩生成函数的工具，基于残差分位函数（RQF）（也称为残差百分位函数），它提供了许多众所周知的概率分布的尾部行为的完整表征，包括幂律分布与其他类型的重尾，甚至非重尾分布的分布。这种方法绘制的图形没有上面提到的平均剩余寿命图、双对数图和帕累托 Q-Q图的缺点，它们对异常值很敏感，能够直观地检验具有小值的幂律，并且不适用于分析大量数据。此外，其他分布类型的尾部也可以用这个方法观察检验。
绘制幂律分布
（Plotting power-law distributions）
一般来说，幂律分布是在双对数坐标轴上绘制的，强调右尾部分。最简便直观的方法是通过（互补）累积分布函数（cumulativedistribution function，缩写为cdf）说明：
累积分布函数 |来源：https://www.google.com
注意，cdf也是幂律函数，只是它的标度指数较小。从数据处理角度，cdf的等价形式是rank-frequency 分布，即先按升序排列的观察值，再将它们与矢量
对应。
尽管便于记录数据，抑或是便于拟合平滑概率密度（质量）函数，但这些方法在数据表示中引入了隐式偏差，因此应该避免。另一方面，所述的cdf法对处理这些隐式偏差更稳健（但并非没有偏误）并且保留了在双对数图形上的线性特征。
虽然在同时用线性最小二乘法拟合幂律时，使用cdf绘制优于pdf（概率密度函数），但其不可避免地在数学上有不准确性。因此，在估计幂律分布的指数时，建议使用最大似然估计。
从经验数据估计指数
有许多方法可以估算幂律尾部的标度指数值，但并非所有方法都能产生无偏且一致的结果。一些最可靠的技术通常基于最大似然估计。替代方法通常基于双对数概率，双对数累积分布函数或对数分组数据进行线性回归，但是，应该避免这些方法，因为它们都可能导致对标度系数的具有显著偏误的估计。
关于这些方法，以及能够使用它们的条件，可以进一步发现，这篇文章全面而详细地提供了可用的代码（Matlab、Python、R和C++）来评估和测试幂律分布的过程。
详细代码如下：
# coding: utf-8
# # 用numpy生成0，1之间的幂律分布
#
# ### 概率密度函数为
# f(x) = a*x^(a-1)
a = 0.4
# 采样数量
samples = 10000
s = np.random.power(a, samples)
# 绘图展示结果
import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(s, bins=50)
x = np.linspace(0, 1, 100)
y = a*x**(a-1.)
normed_y = samples*np.diff(bins)[0]*y
plt.plot(x, normed_y)
plt.show()
# # 使用原生方法生成0，1之间的幂律分布
import math
# 分布函数的反函数
def rev(x,a):
return math.exp(math.log(x) / a)
# 生成分布
s1 = []
for i in range(samples):
s1.append(rev(np.random.uniform(0,1),a))
# 绘图
count, bins, ignored = plt.hist(s1,bins=50)
x = np.linspace(0, 1, 100)
y = a*x**(a-1.)
normed_y = samples*np.diff(bins)[0]*y
plt.plot(x, normed_y)
plt.show()
# # 线性拟合生成结果
# 统计不同区间的数据数量
divide_num = 100
ys = np.zeros(divide_num)
xs = np.linspace(0,1,divide_num)
for i in range(len(s)):
ys[int(s[i] * 100)] += 1
# 使用sklearn包中的回归工具
from sklearn import linear_model
# 回归
x_log = np.log(xs)
y_log = np.log(ys)
#线性拟合数据准备
X_para=[]
Y_para=[]
for x ,y in zip(x_log[1:],y_log[1:]):
X_para.append([float(x)])
Y_para.append(float(y))
# 使用sklearn的线性拟合函数进行拟合
regr = linear_model.LinearRegression()
regr.fit(X_para, Y_para)
#
plt.title(“fit the log data”)
plt.scatter(x_log,y_log,color =”black”)
plt.plot(X_para, regr.predict(X_para),color=’blue’,linewidth=3)
plt.show()
通过R函数估计指数, 并绘制双对数数据拟合线：
详细代码如下：
pwrdist <- function(u,…) {
# u is vector of event counts, e.g. how many
# crimes was a given perpetrator charged for by the police
fx <- table(u)
i <- as.numeric(names(fx))
y <- rep(0,max(i))
y[i] <- fx
m0 <– glm(y~log(1:max(i)),family=quasipoisson())
print(summary(m0))
sub <-
paste(“s=”,round(m0$coef[2],2),”lambda=”,sum(u),”/”,length(u))
plot(i,fx,log=”xy”,xlab=”x”,sub=sub,ylab=”counts”,…)
grid()
lines(1:max(i),(fitted(m0)),type=”b”)
return(m0)
}
如果你有更好的代码，欢迎贡献到集智百科上，或者联系我们，成为我们的志愿者！
来源：集智百科
地址：
http://wiki.swarma.net/index.php?title=%E5%B9%82%E5%BE%8B%E5%88%86%E5%B8%83&variant=zh-hans
编辑：孟婕
推荐阅读
推荐课程
商务合作及投稿转载｜swarma@swarma.org
◆◆◆
搜索公众号：集智俱乐部
加入“没有围墙的研究所”
让苹果砸得更猛烈些吧！
原文始发于微信公众号（ 集智俱乐部 ）：集智
微信扫一扫,分享到朋友圈
