Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have several query frequencies, and I need to estimate the coefficient of Zipf's law. These are the top frequencies:

share|improve this question
according to wikipedia page Zipf's law has two parameters. Number of the elements $N$ and $s$ the exponent. What is $N$ in your case, 10? And frequencies can be calculated by dividing your supplied values by the sum of all the supplied values? – mpiktas Feb 1 '11 at 13:44
let it's ten, and frequencies can be calculated by dividing your supplied values by the sum of all the supplied values.. how can I estimate? – Diegolo Feb 1 '11 at 13:51
up vote 19 down vote accepted

Update I've updated the code with maximum likelihood estimator as per @whuber suggestion. Minimizing sum of squares of differences between log theoretical probabilities and log frequencies though gives an answer would be a statistical procedure if it could be shown that it is some kind of M-estimator. Unfortunately I could not think of any which could give the same results.

Here is my attempt. I calculate logarithms of the frequencies and try to fit them to logarithms of theoretical probabilities given by this formula. The final result seems reasonable. Here is my code in R.

fr <- c(26486, 12053, 5052, 3033, 2536, 2391, 1444, 1220, 1152, 1039)

p <- fr/sum(fr)

lzipf <- function(s,N) -s*log(1:N)-log(sum(1/(1:N)^s))

opt.f <- function(s) sum((log(p)-lzipf(s,length(p)))^2)

opt <- optimize(opt.f,c(0.5,10))

> opt
[1] 1.463946

[1] 0.1346248

The best quadratic fit then is $s=1.47$.

The maximum likelihood in R can be performed with mle function (from stats4 package), which helpfully calculates standard errors (if correct negative maximum likelihood function is supplied):

ll <- function(s) sum(fr*(s*log(1:10)+log(sum(1/(1:10)^s))))

fit <- mle(ll,start=list(s=1))

> summary(fit)
Maximum likelihood estimation

mle(minuslogl = ll, start = list(s = 1))

  Estimate  Std. Error
s 1.451385 0.005715046

-2 log L: 188093.4 

Here is the graph of the fit in log-log scale (again as @whuber suggested):

s.sq <- opt$minimum
s.ll <- coef(fit)


Red line is sum of squares fit, green line is maximum-likelihood fit.

Log-log graph of fits

share|improve this answer
There's also an R package zipfR I haven't tried it though. – onestop Feb 1 '11 at 14:51
@onestop, thanks for the link. It would be nice if someone would answer this question using this package. My solution definitely lacks depth, though it gives some kind of answer. – mpiktas Feb 1 '11 at 15:09
(+1) You're really impressive. So many good contributions in so many different statistical fields! – chl Feb 1 '11 at 22:22
@chl, thanks! I certainly feel that I am not the only one with such characteristics in this site though ;) – mpiktas Feb 2 '11 at 7:37

There are several issues before us in any estimation problem:

  1. Estimate the parameter.

  2. Assess the quality of that estimate.

  3. Explore the data.

  4. Evaluate the fit.

For those who would use statistical methods for understanding and communication, the first should never be done without the others.

For estimation it is convenient to use maximimum likelihood. The frequencies are so large we can expect the well-known asymptotic properties to hold. The log likelihood for data $\{f_i\},\ i=1,\ldots,n$ is

$$\Lambda = -s \sum_{j=1}^n{f_j \log(j)} - \left(\sum_{j=1}^n{f_j}\right) \log\left(\sum_{i=1}^n{i^{-s}}\right).$$

Numerical minimization with these data yields $\hat{s} = 1.45041$ and $\Lambda(\hat{s}) = -94046.7$. This is significantly better (but just barely so) than the least squares solution (based on log frequencies) of $\hat{s}_{ls} = 1.463946$ with $\Lambda(\hat{s}_{ls}) = -94049.5$. (The optimization can be done with a minor change to the elegant, clear R code provided by mpiktas.)

ML will also estimate confidence limits for $s$ in the usual ways. The chi-square approximation gives $[1.43922, 1.46162]$ (if I did the calculations correctly :-).

Given the nature of Zipf's law, the right way to graph this fit is on a log-log plot, where the fit will be linear (by definition):

enter image description here

To evaluate the goodness of fit and explore the data, look at the residuals (data/fit, log-log axes again):

enter image description here

This is not too great: although there's no evident serial correlation or heteroscedasticity in the residuals, they typically are around 10% (away from 1.0). With frequencies in the thousands, we wouldn't expect deviations by more than a few percent. The goodness of fit is readily tested with chi square. We obtain $\chi^2 = 656.476$ with 10 - 1 = 9 degrees of freedom; this is highly significant evidence of departures from Zipf's Law.

Because the residuals appear random, in some applications we might be content to accept Zipf's Law (and our estimate of the parameter) as an acceptable albeit rough description of the frequencies. This analysis shows, though, that it would be a mistake to suppose this estimate has any explanatory or predictive value for the dataset examined here.

share|improve this answer
@whuber, I might humbly suggest a little caution with the formulation given above. Zipf's law is usually stated as a relative-frequency result. It is not (normally considered) the distribution from which an iid sample is drawn. An iid framework is probably not the best idea for these data. Perhaps I'll post more on this later. – cardinal Feb 1 '11 at 18:09
@cardinal I look forward to what you have to say. If you don't have time for a thorough response, even a sketch of what you think might be the "best idea for these data" would be most welcome. I can guess where you're going with this: the data have been ranked, a process which creates dependencies and should require me to defend a likelihood derived without recognizing the potential effects of the ranking. It would be nice to see an estimation procedure with sounder justification. I am hopeful, though, that my analysis can be rescued by the sheer size of the dataset. – whuber Feb 1 '11 at 19:20
@whuber +1 for this detailed explanation of the modeling approach, and ensuing comments. – chl Feb 1 '11 at 22:26
@cardinal Easily. For example, you collect frequencies and identify and rank the ten highest. You hypothesize Zipf's Law. You collect a new set of frequencies and report them according to the previous ranking. That's the iid situation, to which my analysis is perfectly suited, contingent upon the new ranks agreeing with the old ones. – whuber Feb 2 '11 at 16:27
@cardinal Not pedantry; correctness. My analysis blithely overlooked the important issues that you raise. I agree that the last argument is a bit of hand-waving, too. It would still be interesting (hint) to see an approach to the original question that is directly suited to how such frequency data are typically collected. – whuber Feb 2 '11 at 19:48

Here is my attempt to fit the data, evaluate and explore the results using VGAM:


freq <- dzipf(1:100, N = 100, s = 1)*1000 #randomizing values
freq <- freq  + abs(rnorm(n=1,m=0, sd=100)) #adding noize

zdata <- data.frame(y = rank(-freq, ties.method = "first") , ofreq = freq)
fit = vglm(y ~ 1, zipf, zdata, trace = TRUE,weight = ofreq,crit = "coef")

s <- (shat <- Coef(fit)) # the coefficient we've found
probs <- dzipf(zdata$y, N = length(freq), s = s) # expected values
chisq.test(zdata$ofreq, p = probs) 
plot(zdata$y,(zdata$ofreq),log="xy") #log log graph
lines(zdata$y, (probs)*sum(zdata$ofreq),  col="red") # red line, num of predicted frequency

enter image description here

    Chi-squared test for given probabilities

data:  zdata$ofreq
X-squared = 99.756, df = 99, p-value = 0.4598

In our case Chi square's null hypotheses is that the data is distributed according to zipf's law, hence larger p-values support the claim that the data is distributed according to it. Note that even very large p-values are not a proof, just a an indicator.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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