If \(p \sim \mathrm{Beta}(\alpha, \beta)\) and
\(X \sim \mathrm{NegBinomial}(r, p)\), then
\(X \sim \mathrm{BetaNegBinomial}(r, \alpha, \beta)\).
Probability mass function
$$
f(x) = \frac{\Gamma(r+x)}{x! \,\Gamma(r)}
\frac{\mathrm{B}(\alpha+r, \beta+x)}{\mathrm{B}(\alpha, \beta)}
$$
Cumulative distribution function is calculated using recursive algorithm that employs the fact that
\(\Gamma(x) = (x - 1)!\) and
\(
\mathrm{B}(x, y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
\). This enables re-writing probability mass function as
$$
f(x) = \frac{(r+x-1)!}{x! \, \Gamma(r)} \frac{\frac{(\alpha+r-1)!\,(\beta+x-1)!}{(\alpha+\beta+r+x-1)!}}{\mathrm{B}(\alpha,\beta)}
$$
what makes recursive updating from \(x\) to \(x+1\) easy using the properties of factorials
$$
f(x+1) = \frac{(r+x-1)!\,(r+x)}{x!\,(x+1) \, \Gamma(r)} \frac{\frac{(\alpha+r-1)!\,(\beta+x-1)!\,(\beta+x)}{(\alpha+\beta+r+x-1)!\,(\alpha+\beta+r+x)}}{\mathrm{B}(\alpha,\beta)}
$$
and let's us efficiently calculate cumulative distribution function as a sum of probability mass functions
$$F(x) = \sum_{k=0}^x f(k)$$