alternative parametrization of the negative binomial in scipy
from scipy.stats import nbinomdef convert_params(mu, theta): """ Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports See https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations """ r = theta var = mu + 1 / r * mu ** 2 p = (var - mu) / var return r, 1 - pdef pmf(counts, mu, theta): """ >>> import numpy as np >>> from scipy.stats import poisson >>> np.isclose(pmf(10, 10, 10000), poisson.pmf(10, 10), atol=1e-3) True """ return nbinom.pmf(counts, *convert_params(mu, theta))def logpmf(counts, mu, theta): return nbinom.logpmf(counts, *convert_params(mu, theta))def cdf(counts, mu, theta): return nbinom.cdf(counts, *convert_params(mu, theta))def sf(counts, mu, theta): return nbinom.sf(counts, *convert_params(mu, theta))
The Wikipedia page you linked given a precise formula for p and r in terms of mu and sigma, see the very last bullet item in the Alternative parametrization section,https://en.m.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations