Difference Between SciPy.Stats.powerlaw and powerlaw package

Question

I am trying to simulate random variables that are power law distributed based on my understanding of the definition in this Wikipedia article and several other resources where the consensus is that a "power law distributed random variable" has the probability density function (PDF) of the form formula1 and in particular, I'm interested in the case where x_min=1, which reduces to formula2.

However, I've noticed that while the powerlaw python package seems to use this definition, scipy.stats.powerlaw uses a slightly different definition for powerlaw distributions and that is formula3 where alpha is positive and x is between 0 and 1 (inclusive), which we can rearrange if we let formula4 to the form formula5, which would match the form of formula2 (above) but with a negative sign in front of it, with the caveat that alpha would then be negative in the powerlaw version.

The reason I was interested in using the SciPy version is that the Scipy package defines a "percent point function" (ppf) which can be used to generate a set of random values from uniformly distributed probabilities (i.e. the PPF is the inverse of the CDF).

So I'm trying to understand the difference between how the two packages define powerlaw distributions (are they even talking about the same thing?) so that I can translate from my distribution parameters (x_min and alpha) which seem to work under the Powerlaw model to the SciPy.Stats.powerlaw model in order to be able to use the PPF function.

Alternatively, if there is a way to use the Powerlaw package to generate random values from a set of uniformly distributed probabilities, I would greatly appreciate any pointers.

Answer 1

I think you're right.

scipy.stats.powerlaw defines

$$ p(x, \alpha) = \alpha x^{\alpha - 1} $$

powerlaw is much more complex and I don't know it very well but (as I can understand) when you generate random variates from a continuous distribution with $x_{min}=1$, it defines a PDF

$$ p(x, \beta) = -(\beta - 1) x^{-\beta} $$

so that $\beta = 1 - \alpha$.

You can verify this.

import numpy as np
import scipy.stats as sps
import powerlaw
import matplotlib.pyplot as plt

# scipy.stats.powerlaw
alpha = 5
# powerlaw
beta = 1 - alpha

Let's define scipy.stats.powerlaw distribution with $\alpha$

# define scipy distribution
dist = sps.powerlaw(a=alpha)
x = np.linspace(dist.ppf(.0001), dist.ppf(.9999), 1000)
pdf = dist.pdf(x)
plt.plot(x, pdf)
plt.xlabel('x')
plt.ylabel('PDF');

and take 1000 random variates

# take 1000 random variates
np.random.seed(0)
rvs = dist.rvs(1000)

then, let's define powerlaw distribution with $\beta$ and take 1000 random variates

# define powerlaw distribution
theoretical_distribution = powerlaw.Power_Law(
    xmin=1., parameters=[beta], discrete=False
)
# take 1000 random variates
np.random.seed(0)
simulated_data = theoretical_distribution.generate_random(1000)

and plot both

plt.hist(simulated_data, histtype='step', label='powerlaw')
plt.hist(rvs, histtype='step', label='scipy.stats')
plt.title('Random variates')
plt.legend(loc='upper left');

As you can see, they come from the same distribution, and we can check fitting the random variates obtained with powerlaw to scipy.stats.powerlaw

# fit powerlaw random variates with scipy.stats
fit_simulated_data = sps.powerlaw.fit(simulated_data, loc=0, scale=1)
print('alpha:', fit_simulated_data[0])

that gives

alpha: 4.948952195656542

which is the $\alpha$ we defined for scipy.stats.powerlaw.