A while back I began working on a project tasked with producing large loss expectations and ranges for our insured lines of business. The targeted expectations of interest were counts (i.e., how many claims will a line of business incur in excess 1M USD in the next 12 months) and aggregate losses (the sum of the loss dollars associated with the large count projection).
For modeling large count expectations, we would employ one of the following distributions based on the data’s variance to mean ratio:

- If \(Var[N]/E[N] \gt 1\), model large counts using the Negative Binomial Distribution

- If \(Var[N]/E[N] \lt 1\), model large counts using the Binomial Distribution

- If \(Var[N]/E[N] = 1\), model large counts using the Poisson Distribution

Once the details of the modeling approach were agreed upon, the next step was to create a prototype in Excel prior to developing a more robust implementation intended for production. It was at this point that we learned Excel does not have a function to compute the cumulative distribution function for the negative binomial distribution. Digging deeper, we found that the negative binomial CDF can be obtained using the *regularized incomplete beta function*. At this time, I was familiar with the beta function, but not the incomplete beta function, or what *regualrized* meant in this context.

This post is an attempt to highlight variants of the beta function, and reinforce that understanding by implementing a fuctiion that reproduces CDF output from `scipy.stats.nbinom`

for a given parameterization.

(Note: The Wikipedia page for the beta function contains a section titled Software implementation that contains an example of how calculate the complete beta value in VBA using Excel’s builtin `GammaLn`

).

## The Gamma function

No discussion of the Beta function would be complete without first introducing the *gamma function*.
The *gamma function* is an extension of the factorial function, with its argument shifted down by 1, to real
and complex numbers:

If \(n\) is a positive integer, the gamma function reduces to:

We make a point to mention the gamma function since the Beta function can be represented in terms of it, which we demonstrate next.

## The Beta function

The *beta function* is given by:

for real \((a,b) \gt 0\).

We can represent the beta function in terms of the gamma function as follows:

When \(a\) and \(b\) are positive integers, the expression becomes:

Notice that in the first representation of the beta function the limits of integration were \((0,1)\). The *incomplete beta function* is a generalization of the beta function which allows the upper limit of integration to take on values in the range \([0,1]\). Symbolically, the incomplete beta function can be represented as:

When \(x=1\), **the incomplete beta function and the beta function are identical**. Put another way, the beta function is the incomplete beta function evaluated at \(x=1\).

Having described the beta function and the incomplete beta function, we are now in a position to introduce the *regularized incomplete beta function*, also referred to as the *regularized beta function*. It is defined as the ratio of the incomplete beta function to the beta function, evaluated at \(x\):

The regularized beta function is the cumulative distribution function for the beta distribution, and as already mentioned, can be used to calculate the CDF for both the negative binomial and binomial distributions.

For a binomial random variable \(X\), to determine the probability of \(k\) or fewer successes in \(n\) independent trials, \(k \leq n\), the CDF is given by:

To demonstrate, assume \(p=.25\), \(n=5\) and we’re interested in determining the probability of 3 or fewer successes.
Expansion of the binomial summation results in:

Alternatively, calculation with the regularized incomplete beta function yields:

For the negative binomial distribution, assuming the common \(r,k\) parameterization in which the PDF is given by

the expression for the negative binomial CDF can be reduced to:

## Implementation

We now have all the necessary pieces to implement a function to compute the CDF of the negative binomial distribution using functions available in `scipy.special`

. We will compare our results to `scipy.stats.nbinom.cdf`

for a set of inputs to assess the correctness of the implementation.

```
"""
Calculation of negative binomial CDF using regularized incomplete
beta function.
"""
from scipy.special import gamma, beta, betainc
from scipy.stats import nbinom
import numpy as np
# The call signature for nbinom.cdf is:
# nbinom.cdf(<nbr_failures>, <nbr_successes>, <prob_success>)
nbr_successes = 10
nbr_failures = 10
prob_success = .5
# Vectorize nbinom.pmf.
nb_pdfs = np.vectorize(nbinom.pmf)
# Determine sum of nbinom PDF from 0-10.
sum_pdfs = nbinom_pdf(np.arange(nbr_failures+1),nbr_successes,prob_success).sum() # 0.58809852600097756
# Determine nbinom CDF using nbinom.cdf.
actual_cdf = nbinom.cdf(nbr_failures,nbr_successes,prob_success) # 0.58809852600097678
def nb_cdf(k, n, p)
"""
Implementation of negative binomial CDF using the
regularized incomplete beta function.
"""
I = betainc(n+1, k, p)
return(1-I)
nb_cdf(nbr_failures, nbr_successes, prob_success)
# returns: 0.58809852600097678
```

As expected, we find that all three approaches arrive at the same negative binomial CDF for a given set of inputs.