Approximating the Sum of Independent Non-Identical Binomial Random Variables

The distribution of sum of independent non-identical binomial random variables is frequently encountered in areas such as genomics, healthcare, and operations research. Analytical solutions to the density and distribution are usually cumbersome to find and difficult to compute. Several methods have been developed to approximate the distribution, and among these is the saddlepoint approximation. However, implementation of the saddlepoint approximation is non-trivial and, to our knowledge, an R package is still lacking. In this paper, we implemented the saddlepoint approximation in the \textbf{sinib} package. We provide two examples to illustrate its usage. One example uses simulated data while the other uses real-world healthcare data. The \textbf{sinib} package addresses the gap between the theory and the implementation of approximating the sum of independent non-identical binomials.


Introduction
Convolution of independent non-identical binomial random variables appears in a variety of applications, such as analysis of variant-region overlap in genomics [1], calculation of bundle compliance statistics in healthcare organizations [2], and reliability analysis in operations research [3].
Computating the exact distribution of the sum of non-identical independent binomial random variables requires enumeration of all possible combinations of binomial outcomes that satisfy the totality constraint. However, analytical solutions are often difficult to find for sums of greater than two binomial random variables. Several studies have proposed approximate solutions [4,5]. In particular, Eisinga et al. examined the saddlepoint approximation, and compared them to exact solutions [6]. They note that, in practice, these approximations are often as good as the exact solution and are easy to implement in most statistical software.
Despite the theoretical development of aforementioned approximate solutions, a software implementation in R is still lacking. The stats package includes functions for the most frequently used distribution such as dbinom and dnorm. In addition, it also include less frequent distributions such as pbirthday. However, it does not contain functions for distribution of sum of independent non-identical binomials. Another package, EQL, provides a saddlepoint function to approximate the mean of i.i.d random variables. However, this function does not apply to the case where the random variables are not identical. In this paper, we address this deficiency by implementing a saddlepoint approximation in the open source package sinib (sum of independent non-identical binomial random variables). The package provides the standard suite of probability (psinib, distribution (dsinib, quantile (qsinib, and random number generator (rsinib functions. The package is accompanied by a detailed documentation, and can be easily integrated into existing applications.
The remainder of this paper is organized as follows, section 2 formulates the distribution of sum of independent non-identical binomial random variables. Section 3 gives an overview of the saddlepoint approximation. Section 4 describes the design and implementation of the sinib package. Section 5 uses two examples to illustrate its usage. As a bonus, we assess the accuracty of saddlepoint approximation in both examples. Section 6 draws final conclusion and discusses possible future development of the package.
2 Overview of the distribution Suppose X 1 ,...,X m are independent non-identical binomial random variables, and S m = m i=1 X i . We are interested in finding the distribution of S m .
In the special case of m = 2, the probability simplifies to Computation of the exact distribution often involves enumerating all possible combinations of each variable that sums to a given value, which becomes infeasible when n is large. A fast recursion method to compute the exact distribution has been proposed [7,8]. The algorithm is as follows: 1. Compute the exact distribution of each X i .
2. Calculate the distribution of S 2 = X 1 + X 2 using equation 2 and cache the result.
Although the recursion speeds up the calculation, studies have shown that the result may be numerically unstable due to round-off error in computing P (S r = 0) if r is large [9,6]. Therefore, approximation methods are still widely used in literature.

Saddlepoint approximation
The saddlepoint approximation, first proposed by [10] and later extended by [11], provides highly accurate approximations for the probability and density of any distribution. In brief, let M (u) be the moment generating function, and K(u) = log(M (u)) be the cumulant generating function. The saddlepoint approximation to the PDF of the distribution is given as:P whereû is the unique value that satisfies K (û) = s.
[6] applied the saddlepoint approximation to sum of independent non-identical binomial random variables. Suppose that X i ∼ Binomial(n i , p i ) for i = 1, 2, ..., m. The cumulant generating function of S m = m i=1 X i is: The first-and second-order derivative of K(u) are: where The saddlepoint ofû can be obtained by solving K (û) = s. A unique root can always be found because K(u) is strictly convex and therefore K (u) is monotonically increasing on the real line.
The above shows the first-order approximation of the distribution. The approximation can be improved by adding a second-order correction term [12]. where Although the saddlepoint equation cannot be solved at boundaries s = 0 and s = m i=1 n i , their exact probabilities can be computed easily: Incorporation of boundary solutions into the approximation gives: We implemented equation 10 as the final approxmation of the probability density function. For the cumulative density, [12] gave the following approximator: The letters Φ and φ denotes the probability and density of the standard normal distribution.
The accuracy can be improved by adding a second-order continuity correction: We implemented equation 12 in the package to approximate the cumulative distribution.

The sinib package
The package used only functions in base R and the stats package to minimize compatibility issues. The arguments for the functions in the sinib package are designed to have similar meaning to those in the stats package, thereby minimizing the learning required. To illustrate, we compare the arguments of the *binom and the *sinib functions. From the help page of the binomial distribution: • x, q: vector of quantiles.
• p: vector of probabilities.
• n: number of observations.
• size: number of trials.
• prob: probability of success on each trial.
• log, log.p: logical; if TRUE, probabilities p are given as log(p).
• lower.tail: Since the distribution of sum of independent non-identical binomials is defined by a vector of trial and probability pairs (each pair for one constituent binomial), it was neccessary to redefine these arguments in the *sinib functions. Therefore, the following two arguments need to be redefined: • size: integer vector of number of trials.
• prob: numeric vector of success probabilities.
All other arguments remain the same. It is worth noting that when size and prob arguments are given as vectors of length 1, the *sinib functions redues to *binom functions: 5 Two examples

Sum of two binomials
We use two examples to illustrate the use of this package, starting from the simplest case of two binomial random variables with the same mean but different sizes, X ∼ Bin(n, p) and Y ∼ Bin(m, p). The distribution of S = X + Y has an analytical solution, S ∼ Bin(m + n, p). We can therefore use different combinations of (m, n, p) to assess the accuracy of the saddlepoint approximation of the CDF. We use m, n = {10, 100, 1000} and p = {0.1, 0.5, 0.9} to assess the approximation. The ranges of m and n are chosen to be large and the value of p are chosen to represent both boundaries.

Healthcare monitoring
In the second example, we used a health system monitoring dataset by [2]. Suppose n i and p i take the following values.
size=as.integer(c (12,14,4,2,20,17,11,1,8,11)) prob=c(0.074, 0.039, 0.095, 0.039, 0.053, 0.043, 0.067, 0.018, 0.099, 0.045) Since it is difficult to find an analytical solution to the density, we estimated the density with simulation (1e8 trials) and treated it as the ground truth. We then compared simulations with 1e3, 1e4, 1e5, and 1e6 trials, as well as the saddlepoint approximation, to the ground truth.    ggplot(merged,aes(s,diff))+ geom_point()+ facet_grid(~type)+ theme_bw()+ xlab('Quantile')+ylab('Truth-Approx') Figure 6 shows that that the saddlepoint method and the simulation with 1e6 both provide good approximations, while simulations of smaller sizes show clear deviations. Note that the saddlepoint approximation is 5 times faster than the simulation of 1e6 trials.  In this paper, we presented an implementation of the saddlepoint method to approximate the distribution of sum of independent and non-identical binomials. We assessed the accuracy of the method by, first, comparing it with the analytical solution on a simple case of two binomials, and second, with the simulated ground truth on a real-world dataset in healthcare monitoring. These assessments suggest that, while the saddlepoint approximation deviates from the ground truth around the means, it generally provides an approximation superior to simulation in terms of both speed and accuracy. Overall, the sinib package addresses the gap between the theory and implementation on the approximation of sum of independent and non-identical binomial random variables.
In the future, we aim to explore other approximation methods such as the Kolmogorov approximation and the Pearson curve approximation described by [7].