Random negative binomial variates.

def rnegative_binomial(mu, alpha, size=None):
    """
    Random negative binomial variates.
    """
    # 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)