1
$\begingroup$

Background

I have data that roughly follows a power law with a negative exponent (up to a point; also, the parameters of the "fit" were just guesstimated by eye as a demonstration):

enter image description here

Now I tried using the histogram to fit a power law in the shape of

$$f(x) = k \cdot x^a,$$

but I quickly realized that this strongly depends not only on the bin size, but also on whether or not I normalize the histogram (to drop $k$ in the process). So I already inquired about how to approach this here, where I was (I guess unsurprisingly) told that it's not a good idea to fit a histogram at all, and to use maximum likelihood estimators (MLEs) instead.

Question

Now I've started to read up about MLEs, and I'm starting to understand what they are all about, but mind you, I'm not a mathematician. First, I tried to derive the MLEs on my own, but ultimately got stuck. If

$$f(x_i; k, a) = k \cdot x_i^a,$$

then I think

$$L(x_i; k, a) = k^n \cdot \left(\prod_i^n x_i\right)^a,$$

which means

$$\frac{\partial L}{\partial k} = n \cdot k^{n-1} \cdot \left(\prod_i^n x_i\right)^a, \hspace{1cm} \frac{\partial L}{\partial a} = k^n \cdot \left(\prod_i^n x_i\right)^a \cdot \sum_i^n \ln{x_i}.$$

When I set them equal to zero, I can solve for $k$:

$$k = \frac{n}{\sum_i^n \ln{x_i}},$$

but that's as far as I got, and I'm not even sure if that's correct. I didn't manage to sove for $a$, but then I came across this post, which is very related to my question, but unfortunately, I didn't really undestand the answers. I did look into the referenced paper by Clauset et al. though. But it seems that their MLEs are not concerned with my problem:

We assume in these calculations that α > 1, since distributions with α ≤ 1 are not normalizable and hence cannot occur in nature. It is possible for a probability distribution to go as x−α with α ≤ 1 if the range of x is bounded above by some cutoff, but different maximum likelihood estimators are needed to fit such a distribution.

Then I found this paper by White et al., which I think is concerned with my problem, but it's all going over my head at the moment.

So in conclusion, I'm looking for MLEs for $k$ and $a$, though I guess I would also need to address the cutoff for small $x$ somehow.

$\endgroup$
16
  • $\begingroup$ How committed are you to fitting a power law? It won't work well and will be terrible for about half the particles. What is the objective of fitting a distribution? $\endgroup$
    – whuber
    Commented Feb 10, 2022 at 21:49
  • $\begingroup$ I guess rather quite commited? I know it sounds silly, but that's what people in my field did (don't ask me how; from what I've learned now I would guess they actually just fitted histograms like I did), and I would like to compare their results to mine to see if I find similar exponents. $\endgroup$
    – mapf
    Commented Feb 10, 2022 at 22:27
  • $\begingroup$ The power law distribution as typically written has a negative exponent: $kx^{-\alpha}$, and your eyeball estimate is at $\alpha = 3$ in that formulation, so the Clauset et al. paper would apply. $\endgroup$
    – jbowman
    Commented Feb 10, 2022 at 23:14
  • 2
    $\begingroup$ It is probably not to often that a sample from a power law distribution with $\alpha \approx 3$ will have the two leftmost bins with smaller frequencies than the third bin when the bin width is around 0.02 as in your example data. The point is that it seems that a power law would not be a good fit for your data - unless maybe you're only interested in the righthand tail which (I think) is part of why @whuber is asking you about your objective. $\endgroup$
    – JimB
    Commented Feb 11, 2022 at 6:03
  • 1
    $\begingroup$ One way to approach this problem is to fit a likely distribution that has this unimodal shape and exhibits a power tail to the right. It's hard to nominate any distributions without understanding the origin of the data, though. Many other approaches are possible, but they require more guidance to be successful: exactly how would you trade off errors in fitting the tail against bad fits at the left? What has been done in your literature could provide guidance. In many cases what has been published is statistically poor, so following the literature blindly is never a good idea. $\endgroup$
    – whuber
    Commented Feb 11, 2022 at 16:00

1 Answer 1

2
$\begingroup$

Ok, after getting statistical consultation from my university and reading Clauset et al. and White et al. more carefully, I think I finally understand what is the right approach. So adapting Clauset et al.'s terminology, assume I am trying to fit a power law of the following form:

$$f(x) = x^{-\alpha}$$

Previous Issues

The first issue I had was that I assumed that power laws would generally have positive exponents, but apparently this is not so; it seems like in literature, the term "power law" typically referrs to power laws with negative exponents. In the beginning however, this confused me, so because of the positive $\alpha$ values, I thought Clauset at al. did not apply to my case.

The second issue was that I misunderstood what Clauset at al.'s $x_\text{min}$ was used for, as I was still thinking of power laws with a positive exponent, so the lower limit didn't make much sense to me. Now I get that for negative exponents, $f(x)$ diverges for $x \rightarrow 0$, and in order to have a probability density function $p(x)$, we need a lower limit for $x$ such that $\int_{x_\text{min}}^\infty p(x) dx = 1$.

The third issue was that I thought in my original formulation of $f(x)=k \cdot x^{-a}$, I should be able to fit $k$ and $\alpha$ independently. But now I understand that $k$ is not an independent parameter, but the normalization constant which depends on $x_\text{min}$ and $a$. So that's why I was never able to determine the MLEs of both $k$ and $a$.

Finally, the last issue was that I misunderstood what Clauset et al. referred to as the continuous and discrete cases. I thought this was referring to whether the data I was trying to fit was continous or discrete--in which case it was discrete of course--, when it was actually referring to whether the underlying distribution the data were drawn from is continuous or discrete--in which case it is continuous.

Solution

So following Clauset et al.'s definition,

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

where $C$ is the normalization constant. For the continuous case this means that

$$C = (\alpha - 1) \cdot x_\text{min}^{\alpha - 1}.$$

More importantly however, instead of using the approximation suggested for the discrete case (which I originally thought I had to use but which doesn't work for my data), for the continuous case, the MLE for $\alpha$ is easy to compute via:

$$\hat{\alpha} = 1 + n \left[\sum_{i=1}^n \text{ln} \frac{x_i}{x_\text{min}}\right]^{-1},$$

"where $x_i$, $i = 1, ..., n$, are the observed values of $x$ such that $x_i \ge x_\text{min}$".

With this, it is only a matter of estimating the correct $x_\text{min}$. Clauset et al. for example suggests to fit the data using each $x$ value as $x_\text{min}$, and choosing the one that produced the smallest Kolmogorov-Smirnov distance.

Luckily, there already exists a Python library that can do all this: powerlaw (accompanying paper). (The same and similar packages are also awailable for several other programming languages, see e.g. here.)

All that being said however, I guess one should be careful if a power law is even the right model in the first place, or if instead for example an exponential or lognormal distribution is the better choice. Next to the papers, this is for example also discussed by Aaron Clauset here, and the powerlaw package offers likelihood tests to check this. Though in my case, I think a power law still makes the most sense when considering the underlying physicals laws that generate the particle distribution:

enter image description here

Lastly for those interested, there is also another more recent paper on the power law fitting by Virkar and Clauset.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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