# Estimating exponent of Zipf distribution using MLE vs fitting linear regression on log-transformed rank and frequency data

I'm having trouble understanding why I get radically different results if I try to find the parameter of a Zipf distribution when I use the methods proposed by Clauset et al. (2009) as opposed to using the log-transformed rank and frequency data and fitting a linear regression.

This is the code I'm using:

# Creating 500 random samples from Zipf distribution with parameter 1.5
frequencies = np.random.zipf(1.5, 100)
n = len(frequencies)

# Continuous MLE
alpha_hat_cont = 1 + n / sum(np.log(frequencies))

# Approximation of discrete MLE
alpha_hat_discr = 1 + n / sum(np.log(frequencies/0.5))

# Fitting linear regression rank vs frequency to log-transformed data
ranks = np.arange(1, n+1)
slope=-np.polyfit(np.log(ranks), np.log(np.array(sorted(frequencies, reverse=True))), 1)[0]

print(alpha_hat_cont, alpha_hat_discr, slope)


One can also use the powerlaw package in Python, as such:

fit = powerlaw.Fit(frequencies, discrete=True, xmin=1)
print(fit.alpha)


which gives the exact same result as alpha_hat_discr above (if the argument xmin is specified to be equal to 1). I know results don't have to be the same (Clauset et al. suggest using MLE because the OLS on the log-transformed data is a bad approximation), but these are radically different.

For context, I'm trying to find the Zipf exponent of the rank-frequency distribution of a corpus.

Thank you very much for your help!

1. This is a Zeta distribution, not a Zipf distribution, despite Numpy's naming.

2. The issue appears to be that the frequencies variable does not contain the frequencies of the values but the values themselves. I am not a Python programmer, so you'll have to put up with a little R:

x <- rzeta(500, s=1.5)
[1] 1 1 1 1 1 1 2 1 1 1

xf <- table(x)
x
1   2   3   4   5   6   7  11  12  13
395  53  19  11   6   3   4   1   1   2

frequencies <- as.numeric(xf)   # The values in xf:  395, 53, ...
values <- as.numeric(names(xf)) # The labels in xf:  1, 2, 3, ...

summary(lm(log(frequencies/500)~log(values)))

*** stuff ***

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.5042     0.4509  -3.336  0.00594 **
log(values)  -1.5713     0.1977  -7.947 4.03e-06 ***


... and you can see the coefficient estimate for log(values) is pretty close to -s, the parameter we are trying to estimate.

• Thank you jbowman! However, I'm still confused. I thought a power law relation exists also between the frequency and the rank of the word, while what you're doing is the relationship between the amount of words appearing n times and the number n. I guess this is also a power law relationship, but it's not what I'm trying to do. Or am I missing something major here? Sorry for the confusion. Apr 9, 2023 at 8:10
• I'm not sure you quite understand what the random number generator produces. In this case, the word "1" appears 395 times, the word "2" appears 53 times, etc. "1" denotes the most common word (in terms of the underlying probability), "2" is the second most common word, etc. xf above is the table of frequencies of words, not of the number of words appearing n times. Apr 9, 2023 at 14:57