Fit a distribution to a Counter in scipy

Question

I have a collections.Counter object with a count of the occurrences of different values like this:

1:193260
2:51794
3:19112
4:9250
5:6486

How can I fit a probability distribution to this data in scipy? scipy.stats.expon.fit() seems to want a list of numbers. It seems wasteful to create a list with 193260 [1]s, 51794 [2]s, etc. Is there a more elegant or efficient way?

Answer 1

It looks like scipy.stats.expon.fit is basically a small wrapper over scipy.optimize.minimize, where it first creates a function to compute neg-log-likelihood, and then uses scipy.optimize.minimize to fit the pdf parameters.

So, I think what you need to do here is write your own function that computes the neg-log-likelihood of the counter object, and then call scipy.optimize.minimize yourself.

More specifically, scipy defines the expon 'scale' parameter here http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html

So, the pdf is:

pdf(x) = 1 / scale * exp ( - x / scale)

So, taking the logarithm of both sides we get:

log_pdf(x) = - log(scale) - x / scale

Therefore, the negative-log-likelikehood of everything in your counter object would be:

def neg_log_likelihood(scale):
    total = 0.0
    for x, count in counter.iteritems():
       total += (math.log(scale) + x / scale) * count
    return total

Here is a program to try this out.

import scipy.stats
import scipy.optimize
import math
import collections

def fit1(counter):
    def neg_log_likelihood(scale):
        total = 0.0
        for x, count in counter.iteritems():
           total += (math.log(scale) + x / scale) * count
        return total

    optimize_result = scipy.optimize.minimize(neg_log_likelihood, [1.0])
    if not optimize_result.success:
        raise Exception(optimize_result.message)
    return optimize_result.x[0]

def fit2(counter):
    data = []
    # Create an array where each key is repeated as many times
    # as the value of the counter.
    for x, count in counter.iteritems():
        data += [x] * count
    fit_result = scipy.stats.expon.fit(data, floc = 0)
    return fit_result[-1]    

def test(): 
    c = collections.Counter()
    c[1] = 193260
    c[2] = 51794
    c[3] = 19112
    c[4] = 9250
    c[5] = 6486

    print "fit1 'scale' is %f " % fit1(c)
    print "fit2 'scale' is %f " % fit2(c)

test()

Here is the output:

fit1 'scale' is 1.513437 
fit2 'scale' is 1.513438