Zipf’s Distribution: Formula, Worked Example, and How to Fit It

Zipf’s distribution is a discrete probability distribution in which the probability of the item at rank k falls as a power of k: the most common item is twice as likely as the second only when the exponent is 1, and in general p(k) is proportional to 1/k^s. This page shows the exact formula in both of its forms (finite and infinite), works one example through by hand, and covers the step most introductions skip: fitting the exponent without the bias that log-log regression introduces.

Definition: what Zipf’s distribution is

Zipf’s distribution assigns probability to the ranks 1, 2, 3, … of a set of items so that probability declines as a power of rank. Rank 1 is the most frequent item, rank 2 the next, and so on. The shape is governed by a single exponent, written s here (SciPy calls it a). The larger the exponent, the more the top few ranks dominate. The pattern is named for George Kingsley Zipf, who described it for word frequencies in Human Behavior and the Principle of Least Effort (1949).

Zipf’s distribution vs Zipf’s law vs the zeta distribution

These three terms overlap, and software uses them inconsistently. Zipf’s law is the empirical claim that a ranked frequency list follows 1/k^s, usually with s near 1. Zipf’s distribution is the probability model built on that claim. The zeta distribution is the version with no upper limit on rank. SciPy’s documentation states that the Zipf distribution is also known as the zeta distribution, which is a special case of the Zipfian distribution. In that library scipy.stats.zipf is the infinite (zeta) form and scipy.stats.zipfian is the finite form with N items.

Zipf’s law is also closely tied to the power law and Pareto distribution. Newman (2005) notes that a power law is “also known variously as Zipf’s law or the Pareto distribution”, and lists city sizes, earthquakes and personal fortunes among the examples.

Formula, with every symbol explained

Finite form (N items): p(k) = 1 / (k^s × H_N,s), for k = 1, 2, …, N.

• k: the rank, a whole number from 1 to N.
• s: the exponent. It can be any value from 0 upward in the finite form. At s = 0 every rank is equally likely.
• N: the number of distinct items (for example, the size of a vocabulary).
• H_N,s: the generalised harmonic number, 1/1^s + 1/2^s + … + 1/N^s. It is the normalising constant that makes the N probabilities sum to 1.

Infinite form (zeta distribution): p(k) = 1 / (ζ(a) × k^a), for k = 1, 2, 3, …, with a > 1 required.

• ζ(a): the Riemann zeta function, the sum of 1/k^a over all k. It is finite only when a > 1, which is why the infinite form is undefined at a = 1 and below.

The SciPy documentation gives exactly these two probability mass functions: f(k, a) = 1/(ζ(a) k^a) with k ≥ 1 and a > 1 for zipf, and f(k, a, n) = 1/(H_n,a k^a) with a ≥ 0 for zipfian. The finite form converges to the infinite form as N grows, which is why SciPy’s own example shows zipfian with a very large n matching zipf.

Worked example with the arithmetic shown

Take N = 5 items and s = 1.

1. Compute the normaliser: H_5,1 = 1 + 1/2 + 1/3 + 1/4 + 1/5 = 2.2833.
2. Divide each term by it. p(1) = 1 / 2.2833 = 0.4380. p(2) = 0.5 / 2.2833 = 0.2190. p(3) = 0.3333 / 2.2833 = 0.1460. p(4) = 0.25 / 2.2833 = 0.1095. p(5) = 0.2 / 2.2833 = 0.0876.
3. Check: 0.4380 + 0.2190 + 0.1460 + 0.1095 + 0.0876 = 1.0001, which is 1 within rounding.

Rank k	p(k)	Cumulative
1	0.4380	0.4380
2	0.2190	0.6569
3	0.1460	0.8029
4	0.1095	0.9124
5	0.0876	1.0000

With s = 1 the ratio p(1)/p(2) is exactly 2 and p(1)/p(3) is exactly 3: each item’s probability is the top probability divided by its rank. For the infinite form, take a = 2: p(1) = 1/ζ(2) = 6/π² = 0.6079, then p(2) = 0.1520 and p(3) = 0.0675.

Interpretation: what the exponent tells you

The exponent controls how concentrated the probability is in the first ranks. The table uses the same N = 5 as above.

Exponent s	H5,s	p(1)	p(5)	Reading
0.5	3.2317	0.3094	0.1384	Gentle slope; top item about 2.2 times the fifth
1	2.2833	0.4380	0.0876	Classic Zipf; top item exactly 5 times the fifth
2	1.4636	0.6832	0.0273	Steep; the top item holds over two thirds of the mass

Three consequences of the infinite form are worth knowing before you interpret a fitted exponent. First, it needs a > 1 to be a probability distribution at all. Second, the mean exists only when a > 2, because the sum of k × k^-a diverges otherwise. Third, the variance exists only when a > 3. A sample mean or standard deviation from data with an exponent near 2 therefore does not settle toward a stable value as the sample grows, so describe such data with the exponent and rank-frequency plots, not with mean and standard deviation (see descriptive versus inferential statistics).

Assumptions and common mistakes

• Fitting by a straight line on a log-log plot. It is the most common way to estimate the exponent, and it is biased. Clauset, Shalizi and Newman (2009) report that least-squares fitting can produce substantially inaccurate estimates for power-law distributions and gives no indication of whether the data follow a power law at all. They present a framework that combines maximum-likelihood fitting with goodness-of-fit tests based on the Kolmogorov-Smirnov statistic (see the Kolmogorov-Smirnov test guide, including the trap when parameters are estimated from the same data).
• Assuming power-law behaviour holds across the whole range. The same authors point to the difficulty of identifying the range over which power-law behaviour holds, and to large fluctuations in the tail. Report which ranks or values you fitted.
• Mixing up the two SciPy functions. zipf rejects a ≤ 1; zipfian accepts a = 1 but needs N. Using the wrong one changes the normaliser.
• Treating a good-looking plot as a test. A near-straight log-log line does not rule out other heavy-tailed shapes. Check for extreme values separately (see winsorizing versus trimming), and test fit formally rather than by eye (compare normality testing, the same logic with a different reference distribution).

How to report it (APA-style example)

State the model, the estimation method, the data used and the fitted exponent. An example sentence for a results section, using the simulated check described in the next section:

A discrete power-law (zeta) distribution was fitted by maximum likelihood to 5,000 simulated counts generated with a = 2. The estimated exponent was â = 2.01.

For real data, add the range of values fitted, the sample size, the goodness-of-fit result and an interval for the exponent if your software provides one. Put the formula in the methods section so the reader knows whether you used the finite or the infinite form.

How to calculate it in Excel, R or Python

Excel

Put the ranks 1 to 5 in A2:A6 and the exponent in B1. In B2 enter =1/A2^$B$1 and fill down to B6. In C2 enter =B2/SUM($B$2:$B$6) and fill down. With B1 = 1, C2 returns 0.4380.

R

k <- 1:5
s <- 1
w <- 1 / k^s
p <- w / sum(w)
round(p, 4)   # 0.4380 0.2190 0.1460 0.1095 0.0876

Python (SciPy)

import numpy as np
from scipy import stats, special, optimize

# finite form: N = 5, s = 1
print(np.round(stats.zipfian.pmf(np.arange(1, 6), a=1, n=5), 4))
# infinite form: a = 2, p(1) = 6 / pi^2
print(stats.zipf.pmf(1, 2))

# maximum-likelihood exponent for the zeta form
x = stats.zipf.rvs(2.0, size=5000, random_state=1)
nll = lambda a: a * np.log(x).sum() + len(x) * np.log(special.zeta(a))
print(optimize.minimize_scalar(nll, bounds=(1.01, 6), method='bounded').x)  # about 2.006

The first two lines reproduce the worked values above. The last block minimises the negative log-likelihood of the zeta form, which is the maximum-likelihood approach of Clauset, Shalizi and Newman (2009).

FAQ

Is Zipf’s distribution the same as Zipf’s law?

Not exactly. Zipf’s law is the empirical observation that frequency falls roughly as 1/rank^s; Zipf’s distribution is the probability model that formalises it, with a normalising constant so probabilities sum to 1.

What is the difference between zipf and zipfian in SciPy?

zipf is the infinite (zeta) form and needs a > 1. zipfian is the finite form with N items and accepts a ≥ 0.

Why must the exponent exceed 1 in the infinite form?

The normaliser ζ(a) is a sum over all ranks and is finite only when a > 1. At a = 1 it is the harmonic series, which grows without bound.

Can I estimate the exponent from a log-log plot?

You can use the plot to look at the data, but not as your estimator. Use maximum likelihood and a goodness-of-fit test, as in Clauset, Shalizi and Newman (2009).

Sources

• SciPy documentation, scipy.stats.zipf and scipy.stats.zipfian (checked 2026-10-29).
• Newman, M. E. J. (2005). Power laws, Pareto distributions and Zipf’s law. Contemporary Physics, 46, 323-351. arXiv:cond-mat/0412004 (abstract checked 2026-10-29).
• Clauset, A., Shalizi, C. R., and Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Review, 51(4), 661-703. https://doi.org/10.1137/070710111; preprint https://arxiv.org/abs/0706.1062 (abstract checked 2026-10-29).
• Zipf, G. K. (1949). Human behavior and the principle of least effort. Addison-Wesley.