PDiracDelta

optimized_result = fmin_l_bfgs_b(log_likelihood_binned,
                 x0=initial_params,
                 # fprime=log_likelihood_deriv,
                 args=(X,),
                 approx_grad=True,
                 bounds=bounds)

MLE_params = {'size': optimized_result[0][0], 'prob': optimized_result[0][1]}

PDiracDelta

from scipy.optimize import fmin_l_bfgs_b

""" Get reasonable initial values """
# estimates from fitdistr function in R
m = np.mean(X)
v = np.var(X)
size = (m ** 2) / (v - m) if v > m else 1
# convert mu/size parametrization to prob/size
p0 = size / ((size + m) if size + m != 0 else 1)
r0 = size

PDiracDelta

def log_likelihood_binned(params, *args):
    r, p = params
    occurrences_binned = args[0]
    observed_values = np.arange(len(X_binned))
    N = np.sum(occurrences_binned)

    result = np.sum(np.multiply(occurrences_binned, gammaln(observed_values + r))) \
             - np.sum(np.multiply(occurrences_binned, np.log(factorial(observed_values)))) \
             - N * (gammaln(r)) \
             + np.sum(np.multiply(occurrences_binned, observed_values * np.log(1 - (p if p < 1 else 1 - infinitesimal)))) \

PDiracDelta

# inspired by https://github.com/gokceneraslan/fit_nbinom/blob/master/fit_nbinom.py
import numpy as np
from scipy.special import gammaln
from scipy.misc import factorial

infinitesimal = np.finfo(np.float).eps

def log_likelihood(params, *args):
    r, p = params
    X = args[0]