python scipy.stats.powerlaw negative exponent

Question

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

Answer 1

A PDF, integrated over its domain, must equal one. In other words, the area under a probability density function's curve must equal one.

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

The powerlaw density function has a domain from 0 <= x <= 1. On this domain, the integral of x**b is finite for any b > -1. When b is smaller, x**b blows up too rapidly near x = 0. So it is not a valid probability density function when 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

Thus for x**(a-1), a must satisfy a-1 > -1 or equivalently, a > 0.

The first constant a in a * x**(a-1) is the normalizing constant which makes the integral of a * x**(a-1) over the domain [0,1] equal to 1. So you don't get to choose this constant independent of a.

Now if you change the domain to be a measurable distance away from 0, then yes, you could define a PDF of the form C * x**a for negative a. But you'd have to state what domain you want, and I don't think there is (yet) a PDF available in scipy.stats for this.