Python : generating random numbers from a power law distribution [duplicate]

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?

Answer 1

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.