random Power Law distribution a global approach (migrated from Trac #2030)

[thouis]

Original ticket http://projects.scipy.org/numpy/ticket/2030 Reported 2012-01-25 by trac user tristan, assigned to unknown.

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

[thouis]

Attachment in Trac by trac user tristan, 2012-01-25: randomPL.eps

[thouis]

Attachment in Trac by trac user tristan, 2012-01-25: genpower.py