Python : generating random numbers from a power law distribution

Question

I want to draw a random variable between 2 to 15, from a power law distribution with negative exponent (a = -2). I found the following :

r = scipy.stats.powerlaw.rvs(a, loc = 2, scale = 13, size = 1000)

But it does not take negative numbers for a.

Anyone knows a way out?

Someone answered it right here: http://stackoverflow.com/questions/17882907/python-scipy-stats-powerlaw-negative-exponent - if you understand something from it. — Faflok, Jun 29 '15 at 11:27
@RobFoley, I need to check it for different negative exponents.. -1.5, -2, -3 etc.. hence i was hoping if i could get a more general advice :( @Faflok I read that.. but it didn't help me get rid of the `a>0` condition. — Panchi, Jun 29 '15 at 11:37

score 20 · Accepted Answer · edited Sep 09 '17 at 09:56

Power law distribution as defined in numpy.random and scipy.stats are not defined for negative a in the mathematical sense as explained in the answer to this question: they are not normalizable because of the singularity at zero. So, sadly, the math says 'no'.

You can define a distribution with pdf proportional to x^{g-1} with g < 0 on an interval which does not contain zero, if that's what you are after.

For pdf(x) = const * x**(g-1) for a <= x <= b, the transformation from a uniform variate (np.random.random) is:

In [3]: def rndm(a, b, g, size=1):
    """Power-law gen for pdf(x)\propto x^{g-1} for a<=x<=b"""
   ...:     r = np.random.random(size=size)
   ...:     ag, bg = a**g, b**g
   ...:     return (ag + (bg - ag)*r)**(1./g)

Then you can do, for example,

In [4]: xx = rndm(1, 2, g=-2, size=10000)

and so on.

For completeness, here is pdf:

In [5]: def pdf(x, a, b, g):
    ag, bg = a**g, b**g
   ....:     return g * x**(g-1) / (bg - ag)

This all assumes that a < b and g != 0. These formulas should agree with numpy.power and scipy.stats.powerlaw for a=0, b=1 and g > 0.

Python : generating random numbers from a power law distribution

1 Answers1