# 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