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)

为什么需要> 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

Answer 1

在其域上集成的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可测量的距离,则是,您可以将表单的PDF定义C * x**a为否定a.但是你必须说出你想要的域名,而且我认为还没有可用的PDF文件scipy.stats.

Answer 2

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

返回