How to calculate Zipf's law coefficient from a set of top frequencies?

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

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039

-
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

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
$minimum [1] 1.463946$objective
[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

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

Coefficients:
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) plot(1:10,p,log="xy") lines(1:10,exp(lzipf(s.sq,10)),col=2) lines(1:10,exp(lzipf(s.ll,10)),col=3)  Red line is sum of squares fit, green line is maximum-likelihood fit. - There's also an R package zipfR cran.r-project.org/web/packages/zipfR/index.html 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): To evaluate the goodness of fit and explore the data, look at the residuals (data/fit, log-log axes again): 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. - @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: require("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") summary(fit) 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   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.

-