Plotting confidence intervals for Maximum Likelihood Estimate

numpy, python, scipy, statistics, statsmodels

Solution

Looks like you're ok on the first part, so I'll tackle your second and third points.

There are plenty of ways to fit smooth curves, with scipy.interpolate and splines, or with scipy.optimize.curve_fit. Personally, I prefer `curve_fit`, because you can supply your own function and let it fit the parameters for you.

Alternatively, if you don't want to learn a parametric function, you could do simple rolling-window smoothing with numpy.convolve.

As for code quality: you're not taking advantage of numpy's speed, because you're doing things in pure python. I would write your (existing) code like this:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

# N is the true number of books.
# t is the number of weeks.
# unk is the true number of repeats found 
t = 30
unk = 3
def numberrepeats(N, t, iters):
    rand = np.random.randint(0, N, size=(t, iters))
    return t - np.array([len(set(r)) for r in rand])

iters = 1000
ydata = np.empty(500-10)
for N in xrange(10,500):
    sampledunk = np.count_nonzero(numberrepeats(N,t,iters) == unk)
    ydata[N-10] = sampledunk/iters

print "MLE is", np.argmax(ydata)
xdata = range(10, 500)
print len(xdata), len(ydata)
plt.plot(xdata,ydata)
plt.show()

It's probably possible to optimize this even more, but this change brings your code's runtime from ~30 seconds to ~2 seconds on my machine.

Problem

I am trying to write code to produce confidence intervals for the number of different books in a library (as well as produce an informative plot). My cousin is at elementary school and every week is given a book by his teacher. He then reads it and returns it in time to get another one the next week. After a while we started noticing that he was getting books he had read before and this became gradually more common over time. Say the true number of books in the library is N and the teacher picks one uniformly at random (with replacement) to give to you each week. If at week t the number of occasions on which you have received a book you have read is x, then I can produce a maximum likelihood estimate for the number of books in the library following https://math.stackexchange.com/questions/615464/how-many-books-are-in-a-library . Example: Consider a library with five books A, B, C, D, and E. If you receive books [A, B, A, C, B, B, D] in seven successive weeks, then the value for x (the number of duplicates) will be [0, 0, 1, 1, 2, 3, 3] after each of those weeks, meaning after seven weeks, you have received a book you have already read on three occasions. To visualise the likelihood function (assuming I have understood what one is correctly) I have written the following code which I believe plots the likelihood function. The maximum is around 135 which is indeed the maximum likelihood estimate according to the MSE link above. ``` from __future__ import division import random import matplotlib.pyplot as plt import numpy as np #N is the true number of books. t is the number of weeks.unk is the true number of repeats found t = 30 unk = 3 def numberrepeats(N, t): return t - len(set([random.randint(0,N) for i in xrange(t)])) iters = 1000 ydata = [] for N in xrange(10,500): sampledunk = [numberrepeats(N,t) for i in xrange(iters)].count(unk) ydata.append(sampledunk/iters) print "MLE is", np.argmax(ydata) xdata = range(10, 500) print len(xdata), len(ydata) plt.plot(xdata,ydata) plt.show() ``` The output looks like My questions are these: - Is there an easy way to get a 95% confidence interval and plot it on the diagram? - How can you superimpose a smoothed curve over the plot? - Is there a better way my code should have been written? It isn't very elegant and is also quite slow. Finding the 95% confidence interval means finding the range of the x axis so that 95% of the time the empirical maximum likelihood estimate we get by sampling (which should theoretically be 135 in this example) will fall within it. The answer @mbatchkarov has given does not currently do this correctly. There is now a mathematical answer at https://math.stackexchange.com/questions/656101/how-to-find-a-confidence-interval-for-a-maximum-likelihood-estimate .

Original source