# 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] N = X.size result = np.sum(gammaln(X + r)) \ - np.sum(np.log(factorial(X))) \ - N * (gammaln(r)) \ + np.sum(X * np.log(1 - (p if p < 1 else 1 - infinitesimal))) \ + N * r * np.log(p) return -result