docs: def nb_ll(data, P, R):
      """
      Returns the negative binomial log-likelihood of the data.
  
      Args:
          data (array): genes x cells
          P (array): NB success probability param - genes x clusters
          R (array): NB stopping param - genes x clusters
  
      Returns:
          cells x clusters array of log-likelihoods
      """
      # TODO: include factorial...
      #data = data + eps
      genes, cells = data.shape
      clusters = P.shape[1]
      lls = np.zeros((cells, clusters))
      for c in range(clusters):
          P_c = P[:,c].reshape((genes, 1))
          R_c = R[:,c].reshape((genes, 1))
          # don't need constant factors...
          ll = gammaln(R_c + data) - gammaln(R_c) #- gammaln(data + 1)
          ll += data*np.log(P_c) + xlog1py(R_c, -P_c)
          #new_ll = np.sum(nbinom.logpmf(data, R_c, P_c), 0)
          lls[:,c] = ll.sum(0)
      return lls