7

我想为 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

谢谢！

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

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

4

5 回答 5

5

在其域上集成的 PDF 必须等于 1。换句话说，概率密度函数曲线下的面积必须等于 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

幂律密度函数有一个从 0 <= x <= 1 开始的域。在这个域上，对于任何> -1 ， 的积分x**b都是有限的。越小，b在附近炸得太快。所以当 时，它不是一个有效的概率密度函数。bx**bx = 0b <= -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 的可测量距离，那么是的，您可以定义一个 PDF 格式C * x**a为负数a。但是您必须说明您想要的域，而且我认为（目前）还没有可用的 PDF scipy.stats。

于 2013-07-26T14:06:19.663 回答
5

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

返回

幂律

于 2018-03-04T16:56:12.913 回答
3

如果 r 是均匀随机偏差 U(0,1)，则以下表达式中的 x 是幂律分布随机偏差：

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

其中 xmin 是幂律分布保持的最小（正）值，alpha 是分布的指数。

于 2017-01-03T05:15:47.763 回答
0

我的回答几乎和上面的 Virgil 一样，关键的区别在于 alpha 实际上是幂律分布的负指数

因此，如果 r 是均匀随机偏差 U(0,1)，则以下表达式中的 x 是幂律分布随机偏差：

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

其中 xmin 是幂律分布保持的最小（正）值，alpha 是分布的负指数，即 P(x) = [constant] * x**-alpha

于 2019-12-17T19:24:48.683 回答
0

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

p(k) = k^(-gamma)

y是介于 0和1 之间的变量 uniform。

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 和平均连通性。

于 2017-09-06T00:29:33.877 回答