docs: def nb_fit(data, P_init=None, R_init=None, epsilon=1e-8, max_iters=100):
      """
      Fits the NB distribution to data using method of moments.
  
      Args:
          data (array): genes x cells
          P_init (array, optional): NB success prob param - genes x 1
          R_init (array, optional): NB stopping param - genes x 1
  
      Returns:
          P, R - fit to data
      """
      means = data.mean(1)
      variances = data.var(1)
      if (means > variances).any():
          raise ValueError("For NB fit, means must be less than variances")
      genes, cells = data.shape
      # method of moments
      P = 1.0 - means/variances
      R = means*(1-P)/P
      for i in range(genes):
          result = minimize(nb_ll_row, [P[i], R[i]], args=(data[i,:],),
                  bounds = [(0, 1), (eps, None)])
          params = result.x
          P[i] = params[0]
          R[i] = params[1]
          #R[i] = fsolve(nb_r_deriv, R[i], args = (data[i,:],))
          #P[i] = data[i,:].mean()/(data[i,:].mean() + R[i])
      return P,R