python scipy.stats.powerlaw否定指数
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了python scipy.stats.powerlaw否定指数相关的知识,希望对你有一定的参考价值。
我想为scipy.stats.powerlaw例程提供一个负指数,例如a = -1.5,为了绘制随机样本:
"""
powerlaw.pdf(x, a) = a * x**(a-1)
"""
from scipy.stats import powerlaw
R = powerlaw.rvs(a, size=100)
为什么需要> 0,如何提供负a以生成随机样本,以及如何提供归一化系数/变换,即
PDF(x,C,a) = C * x**a
文档在这里
http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.powerlaw.html
谢谢!
编辑:我应该补充一点,我正在尝试复制IDL的RANDOMP函数:
在其域上集成的PDF必须等于一个。换句话说,概率密度函数曲线下的面积必须等于1。
In [36]: import scipy.integrate as integrate
In [40]: y, err = integrate.quad(lambda x: 0.5*x**(-0.5), 0, 1)
In [41]: y
Out[41]: 0.9999999999999998 # The integral is close to 1
powerlaw密度函数具有0 <= x <= 1的域。在该域上,x**b
的积分对于任何b
> -1都是有限的。当b
较小时,x**b
在x = 0
附近爆炸太快。所以当b <= -1
时它不是一个有效的概率密度函数。
In [38]: integrate.quad(lambda x: x**(-1), 0, 1)
UserWarning: The maximum number of subdivisions (50) has been achieved...
# The integral blows up
因此,对于x**(a-1)
,a
必须满足a-1 > -1
或等同地,a > 0
。
a
中的第一个常数a * x**(a-1)
是归一化常数,它使得a * x**(a-1)
在域[0,1]上的积分等于1.所以你不能选择这个常数独立于a
。
现在,如果您将域更改为距离0可测量的距离,则是,您可以为负C * x**a
定义a
形式的PDF。但是你必须说出你想要的域名,我认为scipy.stats
中还没有(可用)PDF。
如果r是均匀随机偏差U(0,1),则下面表达式中的x是幂律分布随机偏差:
x = xmin * (1-r) ** (-1/(alpha-1))
其中xmin是幂律分布所在的最小(正)值,而alpha是分布的指数。
Python包powerlaw可以做到这一点。考虑a>1
具有概率密度函数的幂律分布
f(x) = c * x^(-a)
对于x > x_min
和f(x) = 0
否则。这里c
是一个归一化因子,并被确定为
c = (a-1) * x_min^(a-1).
在下面的示例中,它是a = 1.5
和x_min = 1.0
,并将从随机样本估计的概率密度函数与上面表达式中的PDF进行比较,得出预期结果。
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as pl
import numpy as np
import powerlaw
a, xmin = 1.5, 1.0
N = 10000
# generates random variates of power law distribution
vrs = powerlaw.Power_Law(xmin=xmin, parameters=[a]).generate_random(N)
# plotting the PDF estimated from variates
bin_min, bin_max = np.min(vrs), np.max(vrs)
bins = 10**(np.linspace(np.log10(bin_min), np.log10(bin_max), 100))
counts, edges = np.histogram(vrs, bins, density=True)
centers = (edges[1:] + edges[:-1])/2.
# plotting the expected PDF
xs = np.linspace(bin_min, bin_max, 100000)
pl.plot(xs, [(a-1)*xmin**(a-1)*x**(-a) for x in xs], color='red')
pl.plot(centers, counts, '.')
pl.xscale('log')
pl.yscale('log')
pl.savefig('powerlaw_variates.png')
回报
如果要生成幂律分布,可以使用随机偏差。你只需要在[0,1]之间生成一个随机数并应用逆方法(Wolfran)。在这种情况下,概率密度函数是:
p(k)= k ^( - gamma)
y是0和1之间的变量均匀。
y~U(0,1)
import numpy as np
def power_law(k_min, k_max, y, gamma):
return ((k_max**(-gamma+1) - k_min**(-gamma+1))*y + k_min**(-gamma+1.0))**1.0/(-gamma + 1.0)
现在要生成分发,您只需创建一个数组
nodes = 1000
scale_free_distribution = np.zeros(nodes, float)
k_min = 1.0
k_max = 100*k_min
gamma = 3.0
for n in range(nodes):
scale_free_distribution[n] = power_law(k_min, k_max,np.random.uniform(0,1), gamma)
这将产生伽马= 3.0的幂律分布,如果你想要确定分布的平均值,你必须研究Complex Networks导致k_min取决于k_max和平均锥度。
以上是关于python scipy.stats.powerlaw否定指数的主要内容,如果未能解决你的问题,请参考以下文章