幂律分布的参数估计方法及R实现

幂律分布

幂律分布出现在许多自然以及人为的现象中，如城市的人口、地震的强度以及停电的影响范围等。但其检验及特征描述可能由于长尾部分的波动以及幂律分布适用的范围而变得复杂，常用的方法，如最小二乘拟合，在这方面往往无能为力（他既不能判定数据是否服从幂律分布，又可能给出不准确的参数估计）。Clauset、Shalizi和Newman给出了一个用于识别与测度幂律现象的新框架：该方法基于Kolmogorov-Smirnov(KS)统计量与最大似然比，结合了极大似然估计方法与拟合优度检验。

如果随机变量X的密度函数为：

$$p(x) \propto {x^{ - \alpha }}$$

则称X服从幂律分布，其中$\alpha$一般在(2,3)范围内

在我们一般所见的现象中，X不会在其整个取值范围内服从幂律分布，更可能在大于某个数Xmin的范围内服从幂律分布，我们称X尾部的分布服从幂律分布。

对于连续型随机变量：

$$p(x) = \Pr (x \le X < x + dx) = C{x^{ - \alpha }}$$

$$p(x) = \frac{{\alpha - 1}}{{{x_{\min }}}}{\left( {\frac{x}{{{x_{\min }}}}} \right)^{ - \alpha }}$$

$$F(x) = \int_{ - \infty }^x {p(x)} dx = 1 - {\left( {\frac{x}{{{x_{\min }}}}} \right)^{ - \alpha + 1}}$$

对于离散型随机变量：

$$p(x) = \Pr (X = x) = C{x^{ - \alpha }}$$

$$p(x) = \frac{{{x^{ - \alpha }}}}{{\varsigma (\alpha ,{x_{\min }})}}$$

$$F(x) = 1 - \frac{{\varsigma (\alpha ,x)}}{{\varsigma (\alpha ,{x_{\min }})}} = 1 - \frac{{\sum\limits_{n = 0}^\infty {{{(n + x)}^{ - \alpha }}} }}{{\sum\limits_{n = 0}^\infty {{{(n + {x_{\min }})}^{ - \alpha }}} }}$$

其中：$\varsigma (\alpha ,{x_{\min }}) = \sum\limits_{n = 0}^\infty {{{(n + {x_{\min }})}^{ - \alpha }}}$

*一种近似方法是把离散分布看作是连续分布抽样的取整

---

如何估计参数$\alpha$?

在连续情况下，$\alpha$的极大似然估计与标准误为：

$$\hat \alpha {\rm{ = }}1{\rm{ + }}n{\left[ {\sum\limits_{i = 1}^n {In\frac{{{x_i}}}{{{x_{\min }}}}} } \right]^{ - 1}}$$

$$\sigma = \frac{{\hat \alpha - 1}}{{\sqrt n }} + O\left( {\frac{1}{n}} \right)$$

在离散情况下，$\alpha$的极大似然估计与标准误为：

$$\hat \alpha \simeq 1{\rm{ + }}n{\left[ {\sum\limits_{i = 1}^n {In\frac{{{x_i}}}{{{x_{\min }} - \frac{1}{2}}}} } \right]^{ - 1}}$$

$$\sigma = \frac{{\hat \alpha - 1}}{{\sqrt {n\left[ {\frac{{\varsigma ''(\hat \alpha ,{x_{\min }})}}{{\varsigma (\hat \alpha ,{x_{\min }})}} - {{\left( {\frac{{\varsigma '(\hat \alpha ,{x_{\min }})}}{{\varsigma (\hat \alpha ,{x_{\min }})}}} \right)}^2}} \right]} }}$$

其中：$\varsigma (\alpha ,{x_{\min }}) = \sum\limits_{n = 0}^\infty {{{(n + {x_{\min }})}^{ - \alpha }}}$

*连续与离散公式不可混用

实际上，由于样本量有限(尤其对于尾部的数据)，分布函数CDF比密度函数PDF更稳健

---

如何估计$X_{min}$?

计算KS统计量，取得$X_{min}$使$D^*$最小

$${D^*} = \mathop {\max }\limits_{x \ge {x_{\min }}} \frac{{\left| {S(x) - F(x)} \right|}}{{\sqrt {F(x)(1 - F(x))} }}$$

---

幂律分布数据分析指南：

1. 估计$X_{min}$和$\alpha$的值
2. 计算数据与幂律分布之间的拟合优度，如果该拟合优度>0.1则认为数据服从幂律分布
3. 进行幂律分布与备择假设分布的似然比检验，如果似然比显著不为0，则似然比符号决定取哪个分布（似然比检验可由其他模型比较方法代替：如贝叶斯方法、交叉验证方法等）

---

幂律分布参数估计R代码（仅包括第一、二步）

library(poweRlaw)
data('moby')
#首先建立幂律对象，xmin被设为1，尺度参数被设为空
m_m=displ$new(moby)
m_m$getXmin()#返回预设的Xmin,其值为1

#Xmin与alpha参数的调节方法
#m_m$setXmin(5)
#m_m$setPars(2)

#通过实际分布函数与理论分布之间的距离最小化，求出Xmin
est = estimate_xmin(m_m)
m_m$setXmin(est)
est = estimate_pars(m_m)

plot(m_m)
lines(m_m,col=2,lw=3)

#如何得到图像点的数据
dd = plot(m_m) #返回数据框
head(dd, 3)

##################################################################
#解决Xmin的不确定性：bootstrap方法
#N <- 数据集中的样本数
#for(i in 1:B){
#  sample(,N)
#  estx = estimate_xmin(m_m)
#  m_m$setXmin(estx)
#  estpar = estimate_pars(m_m)
#}
parallel::detectCores()#查看有几个线程
bs = bootstrap(m_m, no_of_sims=200, threads=4)
plot(jitter(bs$bootstraps[,2], factor=1.2), bs$bootstraps[,3])
###################################################################

###################################################################
#是否真的为幂律分布：bootstrap方法
#先计算Xmin与Pars
#ksd=原始数据集的ks距离
#ntail=大于xmin的样本个数
#for(i in 1:B){
#  从二项分布B(n,ntail/n)中抽一个样n1
#  从小于xmin的数中抽n-n1个样
#  从指数为pars的幂律分布中抽n1个样
#  计算ks统计量
#  if ks>ksd P=P+1
#}
#p=P/B
#p值<0.05则拒绝幂律分布
bs_p = bootstrap_p(m_m, no_of_sims=100, threads=4)
###################################################################

---

如何生成一个幂律分布？

随机模拟第一定律：有随机变量u服从[0,1]上的均匀分布，任意随机分布都可由 ${F^{ - 1}}(u)$得到

应用以上定律，连续型幂律分布：

$$x = {x_{\min }}{(1 - u)^{ - \frac{1}{{1 - \alpha }}}}$$

离散型幂律分布：

$$x =round( {x_{\min }}{(1 - u)^{ - \frac{1}{{1 - \alpha }}}})$$

该方法较为简便，且在${x_{\min }}$>5时，引入的误差小于1%

---

幂律分布图像的Python代码

import scipy.stats as st
import numpy as np
import matplotlib.pyplot as plt

x=np.linspace(1,8,800)

fig=plt.figure(figsize=(18,6))
ax1=fig.add_subplot(1,2,1)
ax2=fig.add_subplot(1,2,2)

xmin=1
alphas=[2,2.5,3]
for alpha in alphas:
    y1=((alpha-1)/xmin)*(np.power((x/xmin),-alpha))
    ax1.plot(x,y1,label='alpha=%s'%(alpha))
    ax1.legend(loc="best")
    ax1.set_title('Power Law Distribution xmin=1')
    ax1.set_xlabel("x")
    ax1.set_ylabel('p(x)')

xmins=[1,2,3]
alpha=2.5
for xmin in xmins:
    y2=((alpha-1)/xmin)*(np.power((x/xmin),-alpha))
    ax2.plot(x,y2,label='xmin= %s'%(xmin))
    ax2.legend(loc="best")
    ax2.set_title('Power Law Distribution alpha=2.5')
    ax2.set_xlabel("x")
    ax2.set_ylabel('p(x)')

plt.show()