code: def rnegative_binomial(mu, alpha, size=None):
      
      # Using gamma-poisson mixture rather than numpy directly
      # because numpy apparently rounds
      mu = np.asarray(mu, dtype=float)
      pois_mu = np.random.gamma(alpha, mu / alpha, size)
      return np.random.poisson(pois_mu, size)
docs: Random negative binomial variates.