How to Find the MLEs of a Negative Binomial Distribution

The maximum-likelihood estimates depend on which negative-binomial convention you use and whether the size parameter is known. This article uses the convention that each observation counts failures before the $r$th success. It gives the closed-form estimate when $r$ is fixed and the equation to solve when both parameters are unknown.

Model and likelihood

Let $X_1,\ldots,X_n$ be independent observations, where each $X_i$ is the number of failures before $r$ successes. Write $p$ for the success probability. The probability mass function is

$$P(X=x)=\binom{x+r-1}{x}p^r(1-p)^x,\qquad x=0,1,2,\ldots$$

For positive integer $r$, the binomial coefficient has its usual meaning. In estimation, it is also common to allow a positive real-valued size parameter $r$, using the gamma-function extension. Under this parameterization,

$$\mathbb{E}[X]=\frac{r(1-p)}{p}.$$

For observations $x_1,\ldots,x_n$, let $S=\sum_{i=1}^n x_i$ and let $\bar{x}=S/n$. The likelihood is

$$L(r,p)=\prod_{i=1}^n \binom{x_i+r-1}{x_i}p^{nr}(1-p)^S.$$

The combinatorial factors do not depend on $p$, which makes its estimate straightforward when $r$ is fixed.

MLE of $p$ when $r$ is known

The part of the log-likelihood involving $p$ is

$$\ell(p)=nr\log p+S\log(1-p)+\text{constant}.$$

Differentiate and set the result to zero:

$$\frac{d\ell}{dp}=\frac{nr}{p}-\frac{S}{1-p}=0.$$

Solving gives

$$\hat{p}=\frac{nr}{nr+S}=\frac{r}{r+\bar{x}}.$$

Thus, for known $r$, the MLE of the success probability is the size divided by the size plus the sample mean. The fitted mean is $r(1-\hat{p})/\hat{p}=\bar{x}$. If $S=0$, this formula gives the boundary estimate $\hat{p}=1$. If the parameter space is restricted to $0<p<1$, the likelihood instead has a supremum as $p$ approaches 1, but no interior MLE.

MLE when both $r$ and $p$ are unknown

If $S=0$, all observations are zero. The likelihood is maximized at the boundary $p=1$ if that value is allowed, and every positive $r$ gives the same likelihood, so $r$ is not identified. If $0<p<1$ is required, the likelihood supremum is approached as $p$ tends to 1.

For $S>0$, for each proposed value of $r$, the maximizing value of $p$ is still $r/(r+\bar{x})$. Substitute this into the log-likelihood, then solve the score equation for $r$:

$$\sum_{i=1}^n\left[\psi(x_i+r)-\psi(r)\right] +n\log\left(\frac{r}{r+\bar{x}}\right)=0,$$

where $\psi$ is the digamma function, the derivative of $\log\Gamma$. This equation generally has no closed-form solution, so use a numerical root finder or maximize the profile log-likelihood directly. A finite root need not exist: compare the profile likelihood with its limit as $r$ tends to infinity, which corresponds to a Poisson model with mean $\bar{x}$. If that limit is favored, there is no finite-$r$ MLE under the negative-binomial model alone. After finding a finite $\hat{r}$, calculate

$$\hat{p}=\frac{\hat{r}}{\hat{r}+\bar{x}}.$$

When the observations are nonnegative integers, the digamma differences can be written as finite sums:

$$\psi(x_i+r)-\psi(r)=\sum_{j=0}^{x_i-1}\frac{1}{r+j}.$$

For an observation equal to zero, this sum is empty and contributes zero. This form can make numerical evaluation easier.

Worked example

Suppose the observed failure counts are $0,0,2,6$, and both parameters are to be estimated. There are $n=4$ observations and $\bar{x}=2$. The score equation becomes

$$\frac{2}{r}+\frac{2}{r+1}+\frac{1}{r+2}+\frac{1}{r+3}+\frac{1}{r+4}+\frac{1}{r+5} +4\log\left(\frac{r}{r+2}\right)=0.$$

Numerically, its solution is approximately $\hat{r}=0.495$. Substituting this into the formula for $p$ gives

$$\hat{p}=\frac{0.495}{0.495+2}\approx 0.198.$$

The fitted mean is approximately $0.495(1-0.198)/0.198=2$, matching the sample mean, as expected. In practice, a software optimizer should also be checked for convergence and for whether it found a maximum rather than another stationary point.

If your convention is different

Some sources define a negative-binomial variable as the number of successes before a fixed number of failures, or write the model in terms of its mean $\mu$ and size $r$. Under the mean-size form,

$$P(X=x)=\frac{\Gamma(x+r)}{\Gamma(r)x!} \left(\frac{r}{r+\mu}\right)^r \left(\frac{\mu}{r+\mu}\right)^x.$$

For fixed $r$, the MLE of $\mu$ is $\bar{x}$. This is the same result expressed in different parameters: $\mu=r(1-p)/p$. Always check what the variable counts and which probability is called $p$ before applying a formula.

Common mistakes

• Swapping $p$ and $1-p$. In the convention used here, $p$ is the probability of success, while the observed count is failures.
• Treating the size estimate as closed form. The estimate of $p$ is explicit for fixed $r$, but estimating both parameters usually requires numerical optimization.
• Using the wrong mean. The sample mean estimates $r(1-p)/p$, not $p$ itself.
• Ignoring the Poisson limit. As $r$ grows without bound while the mean stays fixed, the negative binomial approaches a Poisson distribution. If the profile likelihood is maximized at this limit, there is no finite-size MLE under the negative-binomial model alone. Inspect the profile likelihood instead of forcing a finite estimate.