To date, Bayesian inferences for the negative binomial distribution (NBD) have relied on computationally intensive numerical methods (e.g., Markov chain Monte Carlo) as it is thought that the posterior densities of interest are not amenable to closed-form integration. In this article, we present a “closed-form” solution to the Bayesian inference problem for the NBD that can be written as a sum of polynomial terms. The key insight is to approximate the ratio of two gamma functions using a polynomial expansion, which then allows for the use of a conjugate prior. Given this approximation, we arrive at closed-form expressions for the moments of both the marginal posterior densities and the predictive distribution by integrating the terms of the polynomial expansion in turn (now feasible due to conjugacy). We demonstrate via a large-scale simulation that this approach is very accurate and that the corresponding gains in computing time are quite substantial. Furthermore, even in cases where the computing gains are more modest our approach provides a method for obtaining starting values for other algorithms, and a method for data exploration
