幂律分布的参数估计方法及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
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 }}} }}
其中:ς(α,xmin)=∑n=0∞(n+xmin)−ας(α,xmin)=∑n=0∞(n+xmin)−α\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]} }}
其中:ς(α,xmin)=∑n=0∞(n+xmin)−ας(α,xmin)=∑n=0∞(n+xmin)−α\varsigma (\alpha ,{x_{\min }}) = \sum\limits_{n = 0}^\infty {{{(n + {x_{\min }})}^{ - \alpha }}}
*连续与离散公式不可混用
实际上,由于样本量有限(尤其对于尾部的数据),分布函数CDF比密度函数PDF更稳健
如何估计XminXminX_{min}?
计算KS统计量,取得XminXminX_{min}使D∗D∗D^*最小
{D^*} = \mathop {\max }\limits_{x \ge {x_{\min }}} \frac{{\left| {S(x) - F(x)} \right|}}{{\sqrt {F(x)(1 - F(x))} }}
幂律分布数据分析指南:
- 估计XminXminX_{min}和αα\alpha的值
- 计算数据与幂律分布之间的拟合优度,如果该拟合优度>0.1则认为数据服从幂律分布
- 进行幂律分布与备择假设分布的似然比检验,如果似然比显著不为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)F−1(u){F^{ - 1}}(u)得到
应用以上定律,连续型幂律分布:
x = {x_{\min }}{(1 - u)^{ - \frac{1}{{1 - \alpha }}}}
离散型幂律分布:
x =round( {x_{\min }}{(1 - u)^{ - \frac{1}{{1 - \alpha }}}})
该方法较为简便,且在xminxmin {x_{\min }}>5时,引入的误差小于1%
幂律分布图像的Python代码
import scipy.stats as st
import numpy as np
import matplotlib.pyplot as pltx=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()
幂律分布的参数估计方法及R实现相关推荐
- 幂律分布参数估计幂律分布公式计算
如果你看到本帖,哦那恭喜你,其他帖子就不用看了,我都替你筛选过了,都用处不大,没有一个能实际求出参数的,都是生成随机分布的点,然后做拟合,然后就没然后了,把人能气死,别提了就,要么就是给你介绍一堆什么 ...
- 从幂律分布到特征数据概率分布——12个常用概率分布
在机器学习领域,概率分布对于数据的认识有着非常重要的作用.不管是有效数据还是噪声数据,如果知道了数据的分布,那么在数据建模过程中会得到很大的启示. 首先,如下图所示8个特征数据概率分布情况(已经做归一 ...
- Power law and Power law distribution(幂律和幂律分布)
原文:<Power-law distribution in empirical data> 1. Introduction 有些分布可以很好的描述,比如成年男性的身高,某物体的重量等,它们 ...
- 幂律分布图matlab代码,关于幂律分布,你还应该知道如何用代码实现!| 集智百科...
今天我们继续学习幂律分布的基本概念--幂律概率分布,以及如何用代码实现幂律分布.内容来自集智百科,集智百科是复杂系统领域的百科全书,涵盖复杂系统领域的基本概念(持续完善中). 我们正在组织撰写翻译相应 ...
- 数理统计10.15 | 幂律分布
数理统计10.15 | 幂律分布 定义 示例 幂律与"长尾" 克莱伯定律(Kleiber's Law) Zipf定律:书籍中单词频率的分布 Pareto定律(帕累托定律) 性质 标 ...
- 可推导出幂律分布的模型的文献小结
可推导出幂律分布的模型的文献小结 看了一些paper后,一直想写一个关于幂律分布的文献综述,但近年来研究复杂系统,特别是复杂网络的文献增长迅速,而只要是涉及复杂网络的,基本上都谈到了类幂律分布,因此, ...
- Python数据可视化:幂律分布
1.公式推导 对幂律分布公式: 对公式两边同时取以10为底的对数: 令,且为常数,所以公式变为: 所以对于幂律公式,对X,Y取对数后,在坐标轴上为线性方程. 2.可视化 从图形上来说 ...
- 【转载】关于幂律分布的一个笔记
关于幂律分布的一个笔记 原文转自:http://blog.sina.com.cn/s/blog_55954cfb0100ps89.html 0:题外话或补记 最早知道二八法则,还是一本介绍犹太民族杰出 ...
- 关于幂律分布的一个笔记_哈克_新浪博客
关于幂律分布的一个笔记_哈克_新浪博客 关于幂律分布的一个笔记 (2011-03-02 18:12:27) 转载▼ 标签: 幂律 二八法则 杂谈 ...
最新文章
- java struts技术_java技术框架之:struts
- 春泥棒(偷春人) — ヨルシカ(MV + 歌词、汉译、罗马音)
- Linux下使用expect实现跳板机自动跳转/免密登录/自动登录(转)
- python的concat用法_Pandas串联操作concat()用法介绍
- python日历模块_Python日历模块| firstweekday()方法与示例
- vs2010 sp1 安装 Silverlight4_Tools 提示 错误 解决办法
- 【ACL2020】最新效果显著的关系抽取框架了解一下?
- 通讯(transport)
- Oracle 19c 安装教程
- c语言20%3c=10,C语言 练习题(2)
- 高清晰卫星图片:东京、法兰克福机场、华盛顿机场、金字塔、凯旋门
- 3dmax 保存慢 卡死
- 密歇根州立大学联合京东提出深度强化学习算法DeepPage用于分页推荐
- 大学计算机基础四大专业课,《大学计算机基础》课程教学大纲.doc
- 小麦苗博客用到的图片
- Python使用matplotlib可视化哑铃图、强调从一个点到另一个点的变化、数量的变化、客户满意度的变化等(Dumbbell Plot)
- Netflix如何在上万台机器中管理微服务?(史上最全)
- TC118S/TC118H单通道直流马达驱动IC
- 数据分析与数据挖掘实战案例本地房价预测(716):
- ios swift MVVM实例(Model-View-ViewModel)