[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


More information about the NumPy-Tickets mailing list