Negative binomial distribution with Python scipy.stats

Question

I have a dataset of counts on which I tried to fit a Poisson distribution, but my variance is larger than the average so I decided to use a negative binomial distribution.

I use these formulas

to estimate r and p based on the mean and variance of my dataset. However, the nbinom.pmf function requires n and p as parameters. How can I estimate n based on r? The plot is not right if I use r as n.

Answer 1

def convert_params(mu, alpha):
    """ 
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    Parameters
    ----------
    mu : float 
       Mean of NB distribution.
    alpha : float
       Overdispersion parameter used for variance calculation.

    See https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations
    """
    var = mu + alpha * mu ** 2
    p = (var - mu) / var
    r = mu ** 2 / (var - mu)
    return r, p

Answer 2

scipy.stats.nbinom uses a different n, p parameterisation than the one given. See the definition here https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.nbinom.html

You can use the following code to fit the parameters used by nbinom to your sample:

# Estimate parameters
mu = np.mean(sample)  # Mean
sigma_sqr = np.var(sample)  # Variance

# Convert mean and variance to n, p parameterisation
n = mu**2 / (sigma_sqr - mu)
p = mu / sigma_sqr

If you want to test that the estimates actually work, compare them to true parameter values from simulated data:

# Generate sample data
n = 1000
p = 0.009
sample = nbinom.rvs(n=n, p=p, size=10000)

# Estimate parameters
mu = np.mean(sample)  # Mean
sigma_sqr = np.var(sample)  # Variance

# Convert mean and variance to n, p parameterisation
n_est = mu**2 / (sigma_sqr - mu)
p_est = mu / sigma_sqr

# Print results
print("""
                   {:<3}   {:<3}
True parameters:  {:<3}   {:<3}
Estimates      :  {:<3}   {:<3}""".format('n', 'p', n, p,
                                  np.round(n_est, 0), np.round(p_est, 4)))

```