Characterizing/Fitting Word Count Data into Zipf / Power Law / LogNormal Using NLTK and Pandas, I was able to process some text files and generate word count data for them, and finally create a histogram describing word frequency.
However, I'm wondering what kind of analysis should I do in order to characterize this distribution. I'm not sure how should I proceed in order to characterize it. I know for a fact that it wouldn't be possible to fit it into a Poisson distribution as the mean is different from the variance.
Any pointers on how can I find out which type of distribution could this data be fit into? Since we're talking about discrete data and looking at the histogram, my initial guesses were Poisson or Negative Binomial. However, mean is different from variance so that would leave me with negative binomial, or binomial.
I tend to think that Negative Binomial is more likely, however I still have to figure out a way to test this assumption.
 A: The distribution of word frequencies is often characterized by Zipf's law, which states that it has Pareto distribution $p(k) \sim k^{-s}$, so-called power law.
This power law can be well seen as a straight line on the log-log plot of word counts:
import nltk.corpus
import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
# may need nltk.download() to use Brown corpus
counter_of_words = Counter(nltk.corpus.brown.words())
counter_of_counts = Counter(counter_of_words.values())
word_counts = np.array(list(counter_of_counts.keys()))
freq_of_word_counts = np.array(list(counter_of_counts.values()))
plt.scatter(np.log(word_counts), np.log(freq_of_word_counts))
plt.xlabel('Log of word frequency')
plt.ylabel('Log of number of such words')
plt.title('Power law for word frequencies')
plt.show();


The negated slope of this line (roughly 0.5) corresponds to the parameter $s$ of the Zipf law. You can estimate this value with maximizing likelihood:
def neg_zipf_likelihood(s):
    n = sum(freq_of_word_counts)
    # for each word count, find the probability that a random word has such word count
    probas = word_counts ** (-s) / np.sum(np.arange(1, n+1) **(-s))
    log_likelihood = sum(np.log(probas) * word_counts)
    return -log_likelihood

from scipy.optimize import minimize_scalar
s_best = minimize_scalar(neg_zipf_likelihood, [0.1, 3.0] )
print(s_best.x)

which gives you the value of 0.5366.
If you are still not sure whether you need Zipf's distribution or any other distribution, you can compare log likelihood of your data under different distribution, or choose one using Kolmogorov-Smirnov test.
A: I have used the example of David Dale for implementing a MLE power law discrete function.
import nltk.corpus
import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
# may need nltk.download() to use Brown corpus
counter_of_words = Counter(nltk.corpus.brown.words())
counter_of_counts = Counter(counter_of_words.values())

# We sort data
counter_of_counts = sorted(counter_of_counts.items(), key=lambda pair: pair[1], reverse=True)
word_counts = np.asarray(counter_of_counts)[:,0]
freq_of_word_counts = np.asarray(counter_of_counts)[:,1]


f,ax = plt.subplots()
ax.scatter(word_counts, freq_of_word_counts, label = "data")
ax.set_xlabel('Word frequency')
ax.set_ylabel('Number of such words')
ax.set_xscale("log")
ax.set_yscale("log")



def loglik(b):  
    # Power law function
    Probabilities = word_counts**(-b)

    # Normalized
    Probabilities = Probabilities/Probabilities.sum()

    # Log Likelihoood
    Lvector = np.log(Probabilities)

    # Multiply the vector by frequencies
    Lvector = np.log(Probabilities) * freq_of_word_counts

    # LL is the sum
    L = Lvector.sum()

    # We want to maximize LogLikelihood or minimize (-1)*LogLikelihood
    return(-L)

s_best = minimize(loglik, [2])
print(s_best)
ax.plot(word_counts[0:2*10**2], 4*10**4*word_counts[0:2*10**2]**-s_best.x, '--', color="orange", lw=3, label = "fitted MLE")
ax.legend()

The result give us a slope of 1.62 which visually fits the data very well.

