3
$\begingroup$

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.

$\endgroup$
0

1 Answer 1

1
$\begingroup$

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');

enter image description here

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');

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ Did I get the PDF for powerlaw wrong then? Because in your answer the two PDFs are the same, but for me the powerlaw PDF has an opposite sign. Also, in your example, beta = -4, and I was given the impression in the papers powerlaw refers to that it is greater than 1. Conversely, if I'm working with a distribution with beta > 1, then alpha under scipy.stats.powerlaw would be negative, and they explicitly say it must be positive. The one last part I forgot to mention is that in scipy.stats.powerlaw 0 <= x <= 1, whereas in powerlaw (and what I'm trying to simulate) x >= x_min. $\endgroup$
    – Massif
    Mar 29, 2021 at 12:59
  • $\begingroup$ I'm not confident enough with powerlaw package to be able to help you in your specific case. I can only suggest you to take a look at the code implementation github.com/jeffalstott/powerlaw/blob/… and find directly there the PDF that better suits your needs. $\endgroup$ Mar 29, 2021 at 13:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.