Python中的scipy.stats.powerlaw负指数

Question

我希望为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)

为什么需要 a > 0，如何提供负的 a 来生成随机样本，以及如何提供归一化系数/转换，即。

PDF(x,C,a) = C * x**a

文档在这里

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.powerlaw.html

Thanks！

编辑：我应该补充一下，我正在尝试复制IDL的RANDOMP函数：

http://idlastro.gsfc.nasa.gov/ftp/pro/math/randomp.pro

- jtlz2

---

这里的'a'是什么，它是如何计算的？从scipy文档中没有理解@unutbu。 - Sharvari Gc

Answer 1

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')

返回值

- Felix Hoffmann

Answer 2

PDF（概率密度函数）在其定义域内的积分必须等于一。换句话说，概率密度函数曲线下的面积必须等于一。

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

幂律密度函数的定义域为0 <= x <= 1。在此定义域内，对于任何 b > -1，x**b 的积分都是有限的。当 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 * x**(a-1) 中的第一个常数 a 是使 a * x**(a-1) 在区间 [0,1] 上的积分等于 1 的标准化常数。因此您不能独立选择这个常数而不考虑 a。

现在如果您将区域更改为距离 0 可测量的距离，则可以定义形式为负 a 的 C * x**a 的概率密度函数。但是您必须说明您想要的区域，并且我认为（目前）在 scipy.stats 中没有可用的 PDF。

- unutbu

---

到最后一部分：使用正位置“loc”我们可以移动分布。从这个解释中可以得出，对于位置“loc”，对“a”的限制可以放宽为一个函数。这值得一些测试，并且在scipy.stats中扩展应该是可能的。 - Josef

虽然您可以使用 loc 来移动分布，但由于移动是在生成基础分布之后执行的，因此限制条件 a > 0 仍将存在。 - unutbu

你说得对，我没有正确考虑这个问题。它需要一个额外的形状参数来移动和扩大分布的支持。 - Josef

我认为在scipy中这还不可能，但我真的很想知道为什么这么难。这个pdf文件具有解析累积分布函数，因此可以轻松地使用其他人指出的方法，在教科书或这里中概述的方法中解决。 - Rho Phi

Answer 3

如果r是一个均匀分布的随机变量U(0,1)，那么以下表达式中的x就是符合幂律分布的随机变量：

x = xmin * (1-r) ** (-1/(alpha-1))

xmin是大于该值的最小（正）值，超过该值后符合幂律分布，alpha是分布的指数。

- Virgil

---

为什么不直接写成 x = xmin * (r) ** (-alpha) 呢？ - theQman

不知道。这只是我从Aaron Clauset那里得到的公式。 - Virgil

Answer 4

如果您想生成幂律分布，可以使用随机偏差。您只需要在[0,1]之间生成一个随机数，并应用反方法（Wolfram）。在这种情况下，概率密度函数为：

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)

如果你想生成一个gamma=3.0的幂律分布，这将有效。如果你想固定分布的平均值，你必须学习复杂网络，因为k_min取决于k_max和平均连通性。

- Emanuel Fontelles

Answer 5

我的答案与Virgil的几乎相同，但有一个关键区别，即alpha实际上是powerlaw分布的负指数

因此，如果r是均匀随机变量U（0,1），则以下表达式中的x是power-law分布的随机变量：

x = xmin * (1-r) ** (-1/(alpha-1))

xmin是幂律分布成立的最小（正）值，alpha是分布的负指数，即P(x) = [常数] * x**-alpha

- Michael Bietenholz