Parameter estimation for power-law distributions by maximum likelihood
  methods
Parameter estimation for power-law distributions by maximum likelihood methods
Theoretical Physics
Preprint
Parameter estimation for power-law distributions
by maximum likelihood methods
Heiko Bauke
Rudolf Peierls Centre for Theoretical Physics, University of
Oxford, 1 Keble Road, Oxford, OX1 3NP, United Kingdom∗
Distributions following a power-law are an ubiquitous phenomenon. Methods for determining the exponent
of a power-law tail by graphical means are often used in practice but are intrinsically unreliable. Maximum
likelihood estimators for the exponent are a mathematically sound alternative to graphical methods.
http://arXiv.org/abs/0704.1867v2
Parameter estimation for power-law distributions
by maximum likelihood methods
Heiko Bauke
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, 1 Keble Road, Oxford, OX1 3NP, United Kingdom∗
(Dated: August 11, 2007)
Distributions following a power-law are an ubiquitous phenomenon. Methods for determining the exponent
of a power-law tail by graphical means are often used in practice but are intrinsically unreliable. Maximum
likelihood estimators for the exponent are a mathematically sound alternative to graphical methods.
PACS numbers: 02.50.Tt, 89.75.-k
1. Introduction
The distribution of a discrete random variable is referred to as
a distribution with a power-law tail if it falls as
p(k) ∼ k−γ (1)
for k ∈ N and k ≥ kmin. Power-laws are ubiquitous distribu-
tions that can be found in many systems from different disci-
plines, see [1, 2] and references therein for some examples.
Experimental data of quantities that follow a power-law are
usually very noisy; and therefore obtaining reliable estimates
for the exponent γ is notoriously difficult. Estimates that are
based on graphical methods are certainly used most often in
practice. But simple graphical methods are intrinsically unre-
liable and not able to establish a reliable estimate of the expo-
nent γ .
For that reason, the authors of [3] introduced an alternative
approach based on a maximum likelihood estimator for the
exponent γ . Unfortunately the authors concentrate on a rather
idealized type of power-law distributions, namely
p1(k;γ) =
ζ (γ,1)
with k ∈ N, where the normalization constant ζ (γ,1) is given
by the Hurwitz-ζ -function which is defined for γ > 1 and a >
ζ (γ,a) =
(i+ a)γ
. (3)
The distribution (2) is characterized by one parameter only,
and therefore all properties of this distribution (e. g. its mean)
are determined solely by the exponent γ . In many applications
the power-law (2) is too restrictive.
If one states that a quantity follows a power-law, then this
means usually that the tail (k ≥ kmin) of the distribution p(k)
falls proportionally to k−γ . Probabilities p(k) for k < kmin may
differ from the power-law and admit the possibility to tune
the mean or other characteristics independently of γ . In some
∗E-mail:heiko.bauke@physics.ox.ac.uk
situations probabilities p(k) may differ from a power-law for
k ≥ kmax as well, e. g. the distribution may have an exponential
cut-off.
Therefore, I will generalize the maximum likelihood ap-
proach introduced in [3] to distributions that follow a power-
law within a certain range kmin ≤ k < kmax but differ from a
power-law outside this range in an arbitrary way. Furthermore,
I will give some statements about the large sample properties
of the estimate of the power-law exponent and present a nu-
merical procedure to identify the power-law regime of the dis-
tribution p(k). But first, let us see what is wrong with popular
graphical methods.
2. Trouble with graphical methods
All graphical methods for estimating power-law exponents are
based on a linear least squares fit of some empirical data points
(x1,y1), (x2,y2),. . . , (xM,yM) to the function
y(x) = a0 + a1x . (4)
The linear least squares fit minimizes the residual
(yi −a0 −a1xi)2 . (5)
Estimates â0 and â1 of the parameters a0 and a1 are given by
â0 =
∑Mi=1 yi
∑Mi=1 x
∑Mi=1 xi
∑Mi=1 yixi
∑Mi=1 x
∑Mi=1 xi
â1 =
∑Mi=1 yixi
∑Mi=1 xi
∑Mi=1 yi
∑Mi=1 x
∑Mi=1 xi
. (7)
The ansatz for the residual (5) and derivation of (6) and
(7) are based on several assumptions regarding the data points
(xi,yi). It is assumed that there are no statistical uncertain-
ties in xi, but yi may contain some statistical error. The errors
in different yi are independent identically distributed random
mailto:heiko.bauke@physics.ox.ac.uk
variables with mean zero. In particular the standard devia-
tion of the error is independent of xi. For various graphical
methods for the estimation of the exponent of a power-law
distribution these conditions are not met, leading to the poor
performance of these methods.
To illustrate the failure of graphical methods by a computer
experiment N = 10000 random numbers mi had been drawn
from distribution (2) with γ = 2.5 and an estimate γ̂ for the ex-
ponent γ was determined by various graphical methods. The
estimator is a random variable and its distribution depends on
the method that has been used to obtain the estimate. Impor-
tant measures of the quality of an estimator are its mean and
its standard deviation. If the mean of the estimator equals the
true exponent γ then the estimator is unbiased and estimators
with a distribution that is concentrated around γ are desirable.
For each graphical method a histogram of the distribution of
the estimator was calculated to rate the quality of the estimator
by repeating the numerical experiment 500 times.
The most straight forward (and most unreliable) graphical
approach is based on a plot of the empirical probability dis-
tribution p̂(k) on a double-logarithmic scale. Introducing the
indicator function I [·], which is one if the statement in the
brackets is true and else zero, the empirical probability distri-
bution is given by
p̂(k) =
I [mi = k] . (8)
An estimate γ̂ for the power-law exponent γ is established by
a least squares fit to
(xi,yi) = (lnk, ln p̂(k)) for all k ∈ N with p̂(k) > 0, (9)
γ̂ equals the estimate (7) for the slope, see Figure 1 a. Because
the lack of data points in the tail of the empirical distribution
this procedure underestimates systematically the exponent γ ,
see Table I.
There are two ways to deal with the sparseness in the tail of
the empirical distribution, logarithmic binning and consider-
ing the empirical cumulative distribution P̂(k) instead of p̂(k).
The cumulative probability distribution of (2) is defined by
P(k) =
ζ (γ,1)
. (10)
If p(k) has a power-law tail with exponent γ then P(k) follows
approximately a power-law with exponent γ − 1 because for
k ≫ 1 the distribution P(k) can be approximated by
P(k) ≈
ζ (γ,1)
(γ −1)ζ (γ,1)
. (11)
The empirical cumulative probability distribution is given by
P̂(k) =
I [mi ≥ k] . (12)
It is less sensitive to the noise in the tail of the distribution and
therefore a fit of
(xi,yi) = (lnk, ln P̂(k)) for all k ∈ N with P̂(k) > 0 (13)
Table I: Mean and standard deviation of the distribution of the esti-
mate for the power-law exponent γ for various methods. All methods
have been applied to the same data sets of random numbers from
distribution (2) with γ = 2.5. See text for details.
mean standard deviation
method estimate of estimate
fit on empirical distribution 1.597 0.167
fit on cumulative empirical distri-
bution
2.395 0.304
fit on empirical distribution with
logarithmic binninga
2.397 0.080
fit on cumulative empirical distri-
bution with logarithmic binning
2.544 0.127
maximum likelihood 2.500 0.016
aIn [3] a similar experiment is reported. For a fit of the logarithmically
binned empirical probability distribution the authors find a systematical bias
of 29 %. I cannot reconstruct such a strong bias, instead I get a bias of 5 %
only. Probably the quality of this method depends on the details of the binning
procedure.
to a straight line gives much better estimates for the exponent,
see Figure 1 b. But there is still a small bias to too small values
and the distribution of this estimate is rather broad, see Table I.
Logarithmic binning reduces the noise in the tail of the em-
pirical distributions p̂(k) and P̂(k) by merging data points into
groups. By introducing the logarithmically scaled boundaries
bi = roundc
i with some c > 1 (14)
(The function roundx rounds x to the nearest integer.) a linear
least squares fit is performed to
(xi,yi) =
bi + bi+1 −1
bi+1−1
p̂(k)
bi+1 −bi
(xi,yi) =
bi + bi+1 −1
bi+1−1
P̂(k)
bi+1 −bi
, (16)
respectively. As a consequence of the binning the width of
the distribution of the estimate γ̂ of the power-law exponent
γ is reduced, see Figure 1 c, 1 d and Table I. According to
the numerical experiments a fit of the logarithmically binned
cumulative distribution gives the best results among graphical
methods. It shows the smallest systematic bias.
All the methods that have been considered so far have a
common weakness. In the deviation of (6) and (7) it was as-
sumed that the standard deviation of the distribution of the
error in yi is the same for all data points (xi,yi). But this is
obviously not the case. For fixed k the empirical distribution
p̂(k) is a random variable with mean p1(k;γ) and standard
deviation
p1(k;γ)(1− p1(k;γ))/N. For the corresponding
data on a logarithmic scale the standard deviation is approxi-
mately given by the quotient
p1(k;γ)(1− p1(k;γ))/N
p1(k;γ)
1− p1(k;γ)
N p1(k;γ)
. (17)
fit on the distribution p̂(k) fit on the cumulative distribution P̂(k) = ∑i≥k p̂(i)
a) b)
c) d)
Figure 1: Comparison of various methods for estimating the exponent of a power-law. Each figure shows data for a single data set of N = 10000
samples drawn from distribution (2) with γ = 2.5. Insets present histograms of estimates for γ for 500 different data sets.
A power-law distribution p1(k;γ) is a monotonically decreas-
ing function and therefore (17) is an increasing function of k.
Because the variation of the statistical error is not taken into
account, the distribution of the estimate γ̂ is very broad.
Methods that deal with the cumulative distribution have an
additional weakness. Cumulation has the side-effect that the
statistical errors in yi are not independent any more, which
violates another assumption of the deviation of (6) and (7).
To sum up, estimates of exponents of power-law distribu-
tions based on a linear least squares fit are intrinsically inac-
curate and lack a sound mathematical justification.
3. Maximum likelihood estimators
Maximum likelihood estimators offer a solid alternative to
graphical methods. Let p(k;θ ) denote a single parameter prob-
ability distribution. The maximum likelihood estimator θ̂N for
the unknown parameter based on a sample m1,m2, . . . ,mN of
size N is given by
θ̂N = argmax
[L(θ )] = argmax
[lnL(θ )] , (18)
where
L(θ ) =
p(mi;θ ) (19)
denotes the likelihood function. In the limit of asymptotically
large samples and under some regularity conditions maximum
likelihood estimators share some desirable features [5, 6].
• The estimator θ̂N exists and is unique.
• The estimator θ̂N is consistent, that means for every ε > 0
|θ̂N −θ |< ε
= 1 , (20)
where P
|θ̂N −θ |< ε
denotes the probability that the dif-
ference |θ̂N −θ | is less than ε .
• The estimator θ̂N is asymptotically normal with mean θ
and variance
(∆θ̂N)2
ln p(k;θ )
, (21)
where E [·] indicates the expectation value of the quantity
in the brackets.
• Maximum likelihood estimators have asymptotically mini-
mal variance among all asymptotically unbiased estimators.
One says, they are asymptotically efficient.
4. Maximum likelihood estimators for
genuine power-laws
The most general discrete genuine power-law distribution has
a lower as well as an upper bound and is given by
pkmin,kmax(k;γ) =
ζ (γ,kmin,kmax)
for k ∈ N with kmin ≤ k < kmax. Where the non-standard nota-
ζ (γ,kmin,kmax) := ζ (γ,kmin)− ζ (γ,kmax) (23)
has been introduced. If the upper bound is missing the distri-
bution
pkmin(k;γ) =
ζ (γ,kmin)
has to be considered for k ∈N with k ≥ kmin. The distributions
(22) and (24) are generalizations of (2) and will be useful for
the analysis of more general distributions that show a power-
law behavior only in a certain range but have an arbitrary pro-
file outside the power-law regime. This kind of distributions
will be considered in section 5.
The maximum likelihood estimator γ̂N for the parameter γ
of the distribution (22) follows from (18) and is given by
γ̂N = argmax
−N lnζ (γ,kmin,kmax)
or equivalently by the implicit equation
ζ ′(γ̂N ,kmin,kmax)
ζ (γ̂N ,kmin,kmax)
lnmi = 0 , (26)
which has to be solved numerically. The prime denotes the
derivative with respect to γ . The asymptotic variance of this
estimator γ̂N follows from (21) and equals
(∆γ̂N)2
ζ (γ,kmin,kmax)2
ζ ′′(γ,kmin,kmax)ζ (γ,kmin,kmax)− ζ ′(γ,kmin,kmax)2
1 2 3 4 5 6 7 8
Figure 2: Asymptotic standard deviation for maximum likelihood esti-
mators for the exponent of a power-law distribution (24).
Figure 3: Empirical distribution of the maximum likelihood estima-
tor (histogram) versus its theoretical asymptotic distribution, which is
given by a normal distribution with mean γ = 2.5 and variance (27).
The histogram has been obtained from the same data as in Figure 1.
In the limit kmax → ∞ equations (25), (26), and (27) give the
maximum likelihood estimator and the asymptotic variance
of this estimator for power-law distributions lacking an upper
cut-off (24). A graphical representation of the standard devia-
tion (27) in the limit kmax → ∞ is given in Figure 2. For each
fixed kmin the quantity ∆γ̂N
N grows faster than linear with γ .
Therefore the larger the exponent γ the larger the sample size
that is necessary to get an estimate within a given error bound.
If the maximum likelihood method is applied (assuming
a distribution (24)) to the same data as in section 2, numeri-
cal experiments show that the estimates for the exponent are
much more precise. The estimate has no identifiable system-
atic bias, the standard deviation of the distribution of the esti-
mate is smaller by an order of magnitude compared to graphi-
cal methods, see Table I and Figure 3.
5. Maximum likelihood method for
general power-law distributions
The maximum likelihood procedure outlined in section 4 can
be generalized further to distributions p(k) that are no pure
power-laws (22) or (24) but follow a power-law within a cer-
tain finite range or follow a power-law in the whole tail of the
distribution and have an arbitrary profile outside the power-
law regime. The popurse of this section is to establish meth-
ods for identifying the power-law regime and for estimating
the exponent of the power-law regime without making special
assumptions about the profile of the probability distribution
beyond the power-law regime.
The main problem for a generalization of the maximum
likelihood approach is that there might be no good hypothe-
sis for the profile of the probability distribution beyond the
power-law regime. To overcome this difficulty the empirical
data set is restricted to a window kcmin ≤ mi < kcmax. (The
following discussion covers the case of a power-law tail distri-
butions as well by setting kcmax = ∞.) Assuming that p(k) has
a power-law profile for kcmin ≤ k < kcmax then the probabil-
ity distribution of the restricted data set is given by (22) with
kmin = kcmin and kmax = kcmax (or by (24) with kmin = kcmin)
and some unknown exponent γ . This allows to estimate the
power-law exponent by the application of the maximum like-
lihood method on the restricted data set of size N′ as presented
in section 4 without making a hypothesis about the profile of
the probability distribution beyond the power-law regime.
In order to apply the maximum likelihood method one has
to determine the cut-off points kcmin and kcmax first. Here it
has to be taken into account that if the window kcmin ≤ mi <
kcmax is chosen too large the estimate γ̂ is systematically bi-
ased, but on the other hand if it is too small the statistical error
is larger than necessary. In some cases one can make conser-
vative estimates for kcmin and kcmax by plotting the empirical
probability distribution (8) on a double-logarithmic scale. An
appropriate window can also be found by determining esti-
mates γ̂N′(kcmin) as a function of the window and a χ2-test.
Assuming the empirical data is drawn from a distribution
with a power-law tail (no upper cut-off) the lower cut-off point
kcmin can be determined in the following systematic way. By
varying the parameter kcmin the maximum likelihood approach
gives a sequence of estimates γ̂N′(kcmin). If kcmin is very large
the estimate will be quite inaccurate because only a tiny frac-
tion of the experimental data is taken into account; but the
smaller the cut-off kcmin the more accurate the estimate of the
exponent. If kcmin approaches the point from above (but is
still above) where the probability distribution starts do differ
from a power-law γ̂N′(kcmin) will give a very precise estimate
for the exponent γ . On the other hand, if kcmin is too small the
hypothesis that the (restricted) empirical data is drawn from a
power-law distribution is violated which causes a significant
change of the estimate of the power-law exponent.
If the empirical data is drawn from a distribution having
both a lower crossover point as well as an upper crossover
point a sequence of estimates γ̂N′(kcmin) is determined by re-
stricting the data to a sliding window kcmin ≤ mi < wkcmin =
kcmax with w > 1. As long as the window lies completely
within the power-law regime the maximum likelihood esti-
mate obtained from the restricted data set will give a reliable
estimate of the power-law exponent. If the window lies at least
partly outside the power-law regime the estimate is systemati-
cally biased.
To illustrate the procedures outlined above I generated two
data sets from two distributions having a power-law regime.
The first data set of N = 10000 samples was drawn from a
distribution with a power-law tail which is given by
p(k) ∼
5−2.5 for 1 ≤ k ≤ 5
k−2.5 for k > 5 .
Plotting the sequence of estimates γ̂N′(kcmin) against the pa-
rameter kcmin reveals the exponent γ = 2.5 as wells as the
crossover point k = 5 very clearly, see Figure 4. The second
data set of N = 100000 samples was drawn from a distribution
with two crossover points, viz.
p(k) ∼
5−2.25 for 1 ≤ k ≤ 5
k−2.25 for 5 ≤ k ≤ 100
100−2.25e−0.05(k−100) for k > 100 .
Figure 5 shows the sequence of estimates γ̂N′(kcmin) that had
been determined from restricted data sets of samples within
the sliding window kcmin ≤ mi < 5kcmin. This sequence ex-
hibits a broad plateau that corresponds to the power-law expo-
nent γ = 2.25. If the window does not lie completely inside
c min
Figure 4: Sequence of estimates γ̂N ′(kcmin) as a function of the cut-
off kcmin for a data set of N = 10000 samples from distribution (28).
Filled symbols mark where the χ2-test has rejected the hypothesis
that the restricted data follows a power-law (24) with exponent γ =
γ̂N ′(kcmin). An error probability of α = 0.001 was chosen.
c min
2kc max
Figure 5: Sequence of estimates γ̂N ′(kcmin) as a function of the lower
cut-off kcmin for a data set of N = 100000 samples from distribution
(29). Filled symbols mark where the χ2-test has rejected the hypoth-
esis that the restricted data follows a power-law (22) with exponent
γ = γ̂N ′(kcmin). An error probability of α = 0.001 and the window
width w = 5 had been chosen.
the power-law regime the estimate γ̂N′(kcmin) deviates system-
atically from the known exponent.
Apart from a visual inspection of the γ̂N′(kcmin) plot the
crossover point(s) to the power-law regime can be determined
by means of a χ2-test. To apply a χ2-test the data set has to be
divided into some bins and the following binning turned out
to be appropriate: The data is partitioned into a small number
b, say b = 6, of bins. In the case of a distribution with a power-
law tail this means each bin j collects n j items such that
n1 = N
′ p̂(kcmin) q1 = pkcmin(kcmin; γ̂N′(kcmin)) (30)
n2 = N
′ p̂(kcmin + 1) q2 = pkcmin(kcmin + 1; γ̂N′(kcmin))
and finally
nb = N
k=kcmin+b−1
p̂(k) qb =
k=kcmin+b−1
pkcmin(k; γ̂N′ (kcmin)) ,
where q j denotes the probability that a data point falls into
bin j under the assumption that the (restricted) data follows
the power-law (24) with kmin = kcmin and the exponent γ =
γ̂N′(kcmin). For distributions with a finite power-law regime
the binning procedure can be carried out in a similar way. In
this case the summation index in (32) is bounded by kcmin +
b−1≤ k < kcmax and the probability pkcmin,kcmax(k; γ̂N′(kcmin))
has to be considered instead of pkcmin(k; γ̂N′(kcmin)).
The test statistic of the χ2-test is given by
(n j −N′q j)2
N′q j
. (33)
If the to kcmin ≤ mi < kcmax restricted data is given by the
power-law (22) or (24) with kmin = kcmin, kmax = kcmax,
and the exponent γ = γ̂N′(kcmin) then the statistic c2 follows
asymptotically a χ2-distribution with ν = (b− 1) degrees of
freedom, which is given by
pχ2(x,ν) =
xν/2−1e−x/2
2−ν/2 Γ(ν/2)
. (34)
Let χ2α be the (1−α)-quantile of the distribution (34). The hy-
pothesis that the restricted data is given by the power-law (22)
or (24), respectively, with kmin = kcmin, kmax = kcmax, and the
exponent γ = γ̂N′(kcmin) is accepted with the error probability
α if c2 ≤ χ2α . If the window kcmin ≤ mi < kcmax lies not com-
pletely within the power-law regime this hypothesis will be
rejected by the χ2-test and one can detect the upper crossover
point as well as the lower crossover point (where the power-
law loses its validity) in a reliable way, see Figure 4 and Fig-
ure 5.
6. Computational remarks
The normalizing factors of the probability distributions (22)
and (24) are given by the Hurwitz-ζ -function. This func-
tion is less common than other special functions and may not
be available in the reader’s favorite statistical software pack-
age but the GNU Scientific Library [7] offers an open source
implementation of this function. A direct calculation of the
Hurwitz-ζ -function by truncating the sum (3) gives unsatis-
factory results.
The maximum likelihood estimator of the exponent can be
computed numerically either by solving (25) or (26). Equation
(25) has the advantage that it can be solved without calculating
derivatives of the Hurwitz-ζ -function [8], whereas the solu-
tion of (26) involves its first derivative (e. g. bisection method)
or even higher derivatives (e. g. Newton-Raphson method).
An explicit implementation of these derivatives is often not
available but may be calculated numerically.
7. Conclusion
Methods based on a least squares fit are not suited to establish
estimates for power-law distribution exponents because least
squares fits rely on assumptions about the data set that are
not fulfilled by empirical data from power-law distributions.
In this paper maximum likelihood estimators have been intro-
duced as a reliable alternative to graphical methods. These es-
timators are asymptotically efficient and can be applied to data
from a wide class of distributions having a power-law regime.
The crossover points that separate the power-law regime from
the rest of the distribution can be determined by a procedure
based on a χ2-test.
Finally I would like to mention that the idea to plot a se-
quence of estimates γ̂N′(kcmin) as shown in Figure 4 is re-
lated to so-called Hill plots [9, 10]. The Hill estimator is
a maximum likelihood estimator for the inverse of the ex-
ponent of the continuous Pareto distribution p(k) = (γ −
1)(k/kmin)
−γ/kmin, see [10] for a detailed discussion.
Acknowledgments
Work sponsored by the European Community’s FP6 Infor-
mation Society Technologies programme under contract IST-
001935, EVERGROW.
[1] M.E.J. Newman, Contemporary Physics 46(5), 323 (2005)
[2] E.F. Keller, BioEssays 27(10), 1060 (2005)
[3] M.L. Goldstein, S.A. Morris, G.G. Yen, The European Physical
Journal B 41(2), 255 (2004)
[4] F. Ramsey, D. Schafer, The Statistical Sleuth: A Course in Meth-
ods of Data Analysis, 2nd edn. (Duxbury Press, Pacific Grove,
CA, 2002)
[5] L.J. Bain, M. Engelhardt, Introduction to Probability and Math-
ematical Statistics, Duxbury Classic Series, 2nd edn. (Duxbury
Press, 2000)
[6] Y. Pawitan, In all likelihood: statistical modelling and inference
using likelihood, Oxford science publications (Oxford Univer-
sity Press, 2001)
[7] GNU Scientific Library, http://www.gnu.org/software/gsl/
[8] R.P. Brent, Algorithms for Minimization without Derivatives
(Prentice-Hall, Englewood Cliffs, New Jersey, 1973)
[9] B.M. Hill, Annals of Statistics 3(5), 1163 (1975)
[10] H. Drees, L. de Haan, S. Resnick, Annals of Statistics 28(25),
254 (2000)
http://www.gnu.org/software/gsl/
