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

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

Answer 2

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.

Answer 3

If you want to generate power-law distribution, you can use a random deviation. You just have to generate a random number between [0,1] and apply the inverse method (Wolfran). In this case, the probability density function is:

p(k) = k^(-gamma)

and y is the variable uniform between 0 and 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)

Now to generate a distribution, you just have to create an array

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)

This will work to generate a power-law distribution with gamma=3.0, if you want to fixe the average of distribution, you have to study Complex Networks cause the k_min depends of k_max and the average conectivity.

Answer 4

If r is a uniform random deviate U(0,1), then x in the following expression is a power-law distributed random deviate:

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

where xmin is the smallest (positive) value above which the power-law distribution holds, and alpha is the exponent of the distribution.