# [NumPy-Tickets] [NumPy] #2030: random Power Law distribution a global approach

NumPy Trac numpy-tickets@scipy....
Wed Jan 25 04:04:54 CST 2012

```#2030: random Power Law distribution a global approach
--------------------------+-------------------------------------------------
Reporter:  tristan       |       Owner:  somebody
Type:  enhancement   |      Status:  new
Priority:  normal        |   Milestone:  Unscheduled
Component:  numpy.random  |     Version:  1.6.1
Keywords:  power law     |
--------------------------+-------------------------------------------------
Hi,

It's my first post in Numpy so if I'm not clear or if you need more
informations, let me know.

I search for a Numpy function which permit to draw a random distribution
from a power law: N(x) = x**alpha.
I'm surprised to not find an easy way to make it with Numpy.

Of course there are numpy.random.power (only for alpha = a-1 > -1) and
numpy.random.pareto (only for alpha = -a-1 < -1) respectively in
[http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.power.html]
and
[http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.pareto.html].
In the Pareto distribution, there is no way (at least as explain in the
documentation) to change the shape 'm'.

Finally, there is no solution to have a power law with alpha = -1!

So my proposition is quiet easy.
I propose to add a new numpy.random function (called for example
'genpower' for GENeral POWER law) that permit us to draw a random
distribution from a power law where alpha is a real.
In this aim, I wrote a simple function based on the IDL randomP.pro
([http://astro.uni-tuebingen.de/software/idl/astrolib/math/randomp.html]).

{{{
import numpy as np

def genpower(alpha, min, max, N=1000, seed=None):
""" Random Power Law distribution generator

PDF(x,alpha) = x^alpha

Input:
arg:
- alpha      : Power Law index
- min        : Minimal value for X
- max        : Maximal value for X
*arg:
- N          : Number of value to draw (default N=1000)
- sedd       : Value for the seed generator (default None)
"""
## Convert min, max a, and alpha to float
min = float(min)
max = float(max)
alpha = float(alpha)

## Fixe seed
if seed != None: np.random.seed(seed)

## Get uniform vector
rand = np.random.uniform(0.0, 1.0, N)

## Test alpha != 1
if alpha != -1.0:
## General case
pow = alpha + 1.0
X = ( min**pow + rand*(max**pow - min**pow) )**(1.0/pow)
# Normalisation
integral = 1.0/pow * (max**pow - min**pow)
else:
## Particular case alpha = -1
X = min * (max/min)**rand
# Normalisation
integral = np.log(max/min)
#endIF

## Return
#return X, integral
return X
#endDEF

if __name__ == "__main__":
alpha = -1.5
min = 3.
max = 1e3
N = 1e5
X = genpower(alpha, min, max, N=N)

## Show distribution
import pylab as py
py.figure()
ax = py.subplot(2, 1, 1)
n, b, p = py.hist(X, 100, histtype='step', color='k')
bin = np.array([ (b[i]+b[i+1])/2. for i in range(0,len(b)-1) ]) # Get
centered bin value
powerlaw = bin**alpha * np.sum(n)/np.sum(bin**alpha)
py.plot(bin, powerlaw)
ax = py.subplot(2, 1, 2)
n, b, p = py.hist(X, 100, histtype='step', color='k')
bin = np.array([ (b[i]+b[i+1])/2. for i in range(0,len(b)-1) ]) # Get
centered bin value
powerlaw = bin**alpha * np.sum(n)/np.sum(bin**alpha)
py.loglog(bin, powerlaw)
py.show()
#endMAIN
}}}

I the attached a file 'randomPL.eps', which shows some distribution for
different value of alpha.

I hope to help Numpy community.

Tristan

--
Ticket URL: <http://projects.scipy.org/numpy/ticket/2030>
NumPy <http://projects.scipy.org/numpy>
My example project
```