How to calculate the maximum likelihood of the sum of independent negative binomial variables? In the paper Small-sample estimation of negative binomial dispersion, with applications to SAGE data, section 4.1 Conditional maximum likelihood:
For RNA sequencing data, assume the counts for a single tag across $n$ libraries is a negative binomial random variable. Consider
$Y_1, \cdots , Y_n$ as independent and $NB(\mu_i = m_i \lambda, \phi)$, where $m_i$ is the library size (i.e. total number of tags sequenced for library $i$) and $\lambda$ represents the proportion of the library that is a particular tag. The probability mass function is:
$$f(y_i;\mu,\phi)=P(Y=y_i)=\frac{\Gamma(\phi^{-1}+y_i)}{\Gamma(\phi^{-1})\Gamma(y_i+1)}\left(\frac{1}{\phi^{-1}\mu^{-1}+1}\right)^{y_i}\left(\frac{1}{\phi\mu+1}\right)^{\phi^{-1}} \tag 1$$
$$\text{E}(Y)=\mu$$ and $$\text{Var}(Y)=\mu+\phi\mu^2$$
If all libraries are the same size (i.e. $m_i ≡ m$), the sum
$$Z = Y_1 + \cdots + Y_n \sim NB(nm\lambda, \phi n^{−1})$$
How to derive conditional maximum likelihood for $\phi$ that is independent of $\lambda$?
$$l_{Y|Z=z}(\phi)=\left[\sum_{i=1}^{n}\log\Gamma(y_i+\phi^{-1})\right]+\log\Gamma(n\phi^{-1})-\log\Gamma(z+n\phi^{-1})-n\log\Gamma(\phi^{-1})$$

If we write out the log likelihood, dropping terms that don't involve $\phi$, we get:
$$\sum\log \Gamma(y_i+\phi^{-1}) - n \log\Gamma(\phi^{-1}) + \sum y_i\log\left({\phi\mu \over 1 + \phi\mu}\right) + n\phi^{-1}\log\left({1 \over 1 + \phi\mu}\right)$$
where we have multiplied the numerator and denominator of $1/(\phi^{-1}\mu^{-1}+1)$ by $\phi\mu$.  Rewriting the $\log$ terms and rearranging gives us:
$$ \sum\log \Gamma(y_i+\phi^{-1}) - n \log\Gamma(\phi^{-1}) + \sum y_i\log\phi\mu - \left(\sum y_i+n\phi^{-1}\right)\log(1+\phi\mu)$$
Some more rearranging:
$$\sum\log \Gamma(y_i+\phi^{-1}) - n \log\Gamma(\phi^{-1}) + z\log\phi\mu -(z+n\phi^{-1})\log(1+\phi\mu)$$
I do not believe the equation as presented in the paper is correct, and referring back to the original formula for the probability distribution, you can see there are two terms on the r.h.s. that a) involve $\phi$ and b) are not arguments to the Gamma function, that have disappeared in their final expression, to be replaced by two terms involving $\log \Gamma$ that appear nowhere in the probability distribution as arguments to the Gamma function.
Substituting $\bar{y} = (1/n)\sum y_i$ for $\mu$, as it is the MLE for $\mu$ and is what enables us to get rid of $\lambda$ (which is hidden inside $\mu$) gives us the one-dimensional optimization problem:
$$\max_{\phi} \left[\sum\log \Gamma(y_i+\phi^{-1}) - n \log\Gamma(\phi^{-1}) + n\bar{y}\log\phi\bar{y} -n(\bar{y}+\phi^{-1})\log(1+\phi\bar{y})\right]$$
There is no analytic solution to this problem, only numerical ones.
