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

时间 2020-07-11 标签分布参数估计方法实现

幂律分布

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

若是随机变量X的密度函数为：web

p (x) \propto x^{- α}

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

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

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

对于连续型随机变量：svg

p (x) = Pr (x \leq X < x + d x) = C x^{- α}

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

p (x) = \frac{α - 1}{x_{min}} {(\frac{x}{x_{min}})}^{- α}

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

F (x) = \int_{- \infty}^{x} p (x) d x = 1 - {(\frac{x}{x_{min}})}^{- α + 1}

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

对于离散型随机变量：
函数

p (x) = Pr (X = x) = C x^{- α}

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

p (x) = \frac{x^{- α}}{ς (α, x_{min})}

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

F (x) = 1 - \frac{ς (α, x)}{ς (α, x_{min})} = 1 - \frac{\sum_{n = 0}^{\infty} {(n + x)}^{- α}}{\sum_{n = 0}^{\infty} {(n + 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 }}}$ atom

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

如何估计参数 $\alpha$ ?

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

\hat{α} = 1 + n {[\sum_{i = 1}^{n} I n \frac{x_{i}}{x_{min}}]}^{- 1}

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

σ = \frac{\hat{α} - 1}{\sqrt{n}} + O (\frac{1}{n})

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

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

\hat{α} ≃ 1 + n {[\sum_{i = 1}^{n} I n \frac{x_{i}}{x_{min} - \frac{1}{2}}]}^{- 1}

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

σ = \frac{\hat{α} - 1}{\sqrt{n [\frac{ς^{″} (\hat{α}, x_{min})}{ς (\hat{α}, x_{min})} - {(\frac{ς^{'} (\hat{α}, x_{min})}{ς (\hat{α}, x_{min})})}^{2}]}}

$\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^{*} = max_{x \geq x_{min}} \frac{| S (x) - F (x) |}{\sqrt{F (x) (1 - F (x))}}

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

幂律分布数据分析指南：

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

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

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 - α}}

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

离散型幂律分布：

x = r o u n d (x_{min} (1 - u)^{- \frac{1}{1 - α}})

$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()