python scipy.stats.powerlaw 负指数 [英] python scipy.stats.powerlaw negative exponent

问题描述

我想为 scipy.stats.powerlaw 例程提供负指数，例如a=-1.5，为了抽取随机样本:

<代码>"""powerlaw.pdf(x, a) = a * x**(a-1)"""从 scipy.stats 导入幂律R = powerlaw.rvs(a, size=100)

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

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

文档在这里

I want to supply a negative exponent for the scipy.stats.powerlaw routine, e.g. a=-1.5, in order to draw random samples:

"""
powerlaw.pdf(x, a) = a * x**(a-1)
"""

from scipy.stats import powerlaw
R = powerlaw.rvs(a, size=100)

Why is a > 0 required, how can I supply a negative a in order to generate the random samples, and how can I supply a normalization coefficient/transform, i.e.

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

The documentation is here

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

Thanks!

EDIT: I should add that I'm trying to replicate IDL's RANDOMP function:

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

解决方案

The Python package powerlaw can do this. Consider for a>1 a power law distribution with probability density function

f(x) = c * x^(-a) 

for x > x_min and f(x) = 0 otherwise. Here c is a normalization factor and is determined as

c = (a-1) * x_min^(a-1).

In the example below it is a = 1.5 and x_min = 1.0 and comparing the probability density function estimated from the random sample with the PDF from the expression above gives the expected result.

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

returns

这篇关于python scipy.stats.powerlaw 负指数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！