# Power Law Probability Distribution From Observations

This is probably a very simple question: I am trying to understand a power-law fitting technique by Aaron Clauset (http://tuvalu.santafe.edu/~aaronc/powerlaws/), but to do this I need to understand how a probability distribution is calculated from a series of observations, for example if my observations are:

x = [1, 40, 100, 200, 2000, 4, 10, 50, 75];

How do I calculate the probability distribution P(x) ?

I would like to do this so that I can fit a power law to the probability distribution and determine if a power law fit is acceptable and accurate.

Thanks for any help.

Regards,

Fido

-
On second thought, am I misinterpreting the issue. It's not that I calculate $$P(x)$$ and then fit a power law but rather attempt to fit a power law distribution $$P(x) \propto x^\alpha$$ and then determine the accuracy of the fitted distribution. Is that the correct procedure? – Fido Jun 4 '13 at 12:44

The first step is to fit the power-law model to your data. The model has two parameters alpha and xmin whose values need to be chosen. The usual definition of the power-law distribution is $p(x) \propto x^{-\alpha}$ where $\alpha>1$ and for $x\geq x_{\rm min}$.

As Clauset, Shalizi and Newman point out, choosing the correct value for xmin is slightly tricky, so for now just set xmin=1, the smallest value in your data set. Given this choice, and the fact that your data are integers (rather than values with decimals), the maximum likelihood estimate for $\alpha$, which we denote as $\hat{\alpha}$ to indicate that it is an estimated derived from the data rather than the true value, is roughly

$\hat{\alpha} = 1 + \frac{n}{\sum_{i=1}^{n} \ln x_i/(x_{\rm min}-1/2)}$

This expression (Eq. 3.7 in the arxiv article linked above) is an approximation. The exact estimator is given by the ratio of two incomplete Hurwitz zeta functions (Eq. 3.4), or by directly maximizing the log-likelihood function for the discrete power-law distribution (Eq. 3.5). But, if you just want a rough estimate, the above equation should be fine.

In your case, $n=9$ since your vector $x$ has that many integers in it.

Once you've fitted the model, then you would need to evaluate whether the fit is statistically significant (or rather, whether the model is terrible, given the data). The results is a value $p$ that tells you whether you can rule out the fitted power-law model. When $p<0.1$, you reject it. Otherwise, you say it's statistically plausible. Clauset, Shalizi and Newman describe how to do this using a semi-parametric bootstrap, which is a little complicated to write out here.

For the data you give, however, there is no point in doing it because a data set with $n=9$ is so small that basically any distribution, including the power law, will be statistically plausible.

Also, I should point out that Clauset and Shalizi published Matlab and R implementations of their methods, which you could use to do this analysis.

-