Solved – Exemplar MLE for negative binomial

Question

I recently compared MLE estimates for a negative binomial fit using two different pieces of software, and got different results. I'd like to determine which (if either) is correct.

To do that, I'd like an example of a simple set of data along with the corresponding MLE solution. I'd really like an analytic proof of the optimality of the MLE solution, since determining it with some third software library seems a bit circular. However, we are free to choose any data that might make the analysis work out nicely.

For the record, I noticed the discrepancy through these snippets of code (first library in R, second library in Python). The data was {2,4} (i.e., I only fit two data points).

library(MASS)
x=fitdistr(c(2,4),"Negative Binomial")
p = x$estimate['size'][[1]] / (x$estimate['size'][[1]] + x$estimate['mu'][[1]])
# Results:
# > p
# [1] 0.9708738
# > x
#      size         mu    
#   100.00020     3.00000 
#  (505.70845) (  1.24298)

and

import fit_nbinom
x = fit_nbinom.fit_nbinom(np.array([2,4]))
mu = ((1.0-x['prob'])/x['prob'])  * x['size']
# Results:
# In []: x
# Out[]: {'prob': 0.98892868281051893, 'size': 270.37321834245245}
#
# In []: mu
# Out[]: 3.026899423427535

Answer 1

The negative binomial density is:

$$f(x) = \Gamma(x + n )/(\Gamma(n) x!) p^n (1-p)^x$$

So for your sample:

$$L([2,4]) = f(2)f(4) \propto \Gamma(2 + n )\Gamma(4 + n )/2\Gamma(n)( p^{2n} (1-p)^6 ) $$

which we are to maximize

The mean is 3 by default as others have noted. Define a few useful R wrappers:

loglik <- function(size) sum(dnbinom(c(2,4), size=size, mu=3, log=T))

and just visualize what happens over a range of "size"

plot(1:1000, sapply(1:1000, loglik), type='l', xlab='Size parameter', ylab='Trace of log likelihood')

You can see the likelihood approaches an asymptote and has no unique, finite maximum.