2.1. Group testing setup

We now state the group testing problem in slightly more formal language and introduce some notation similar to that of [Reference Aldridge, Johnson and Scarlett5]. We will write $n$ for the total population of individuals, which we will refer to as “items.” Instead of referring to infected and healthy individuals, we will follow standard group testing terminology by calling them “defective” and “nondefective,” respectively. We write ${{\mathcal {K}}}$ for the set of defective items (or defective set for short) and $k = | {{\mathcal {K}}} |$ for the total number of defective items, and $T$ for the number of tests.

As is standard, we can represent a nonadaptive testing strategy by a binary $T \times n$ test matrix $X$, with rows corresponding to tests and columns corresponding to items. Here, the entry $X_{ti} = 1$ means that item $i$ appears in test $t$. In this paper, we focus on Bernoulli testing, where the $(X_{ti})$ are independent Bernoulli random variables with parameter $p$. We will refer to this as a “Bernoulli test design with parameter $p$,” and this design will apply throughout. Of particular interest will be the case $p = 1/k$.

The outcome of test $t$ is represented as a binary value $Y_t$ (where $Y_t = 1$ means a positive test) which can be calculated for this test matrix as

(1)\begin{equation} Y_t = \bigvee_{i \in {{\mathcal{K}}}} X_{ti}, \end{equation}

where $\bigvee$ represents a standard binary OR, capturing the fact that each test is positive if and only if it contains at least one of the items in ${{\mathcal {K}}}$.

As in [Reference Aldridge, Baldassini and Johnson4] we write $q_0 = (1-p)^{k}$, noting that each test is positive independently with probability $1-q_0$ (since it is negative, if and only if it contains none of the $k$ defectives).

Again, as in [Reference Aldridge, Baldassini and Johnson4,Reference Coja-Oghlan, Gebhard, Hahn-Klimroth and Loick10] and other papers, we will often study what is referred to in [Reference Aldridge, Johnson and Scarlett5] as the sparse regime. In this setting, the number of items $n$ tends to infinity and the number of defectives $k = k(n) = n^{\theta }$ for some explicit parameter $\theta \in (0,1)$. In this context, it is natural to consider an asymptotic regime where the number of tests $T= (c/q_0) k \log (n)$ for some constant $c$, noting that if $p=1/k$, then asymptotically $q_0$ converges to $e^{-1}$. Here and throughout our work, we write $\log$ for the natural logarithm.

Another asymptotic setting of interest (see, e.g., [Reference Aldridge1] and [Reference Aldridge, Johnson and Scarlett5] Sect. 5.5) is referred to as the linear regime. Here, again $n$ tends to infinity, and the number of defectives $k = k(n) = \beta n$ for some explicit parameter $\beta \in (0,1)$. In this context, we consider an asymptotic regime where the number of tests $T=c n$ for some constant $c$. Although, in this regime, Aldridge [Reference Aldridge1] proved that no nonadaptive algorithm can outperform individual testing, there is still interest in understanding the performance of two-stage or adaptive algorithms, and analysis of $G$ of the kind presented in this paper can help with this.

However, in many practical group testing contexts, we are also interested in what we refer to as the finite blocklength setting, following terminology popularized for example by the work of Polyanskiy et al. [Reference Polyanskiy, Poor and Verdú27]. In this context, we wish to understand the performance of algorithms in solving concrete problems such as $n=500$, $k=10$ (see [Reference Aldridge, Baldassini and Johnson4]), where asymptotic bounds may not necessarily give the best guide to actual performance. Interest in finite blocklength problems was particularly prompted by the COVID pandemic, where for example the use of 96-well PCR plates [Reference Erlich, Gilbert, Ngo, Rudra, Thierry-Mieg, Wootters, Zielinski and Zuk16] means that the number of tests may be bounded by (a multiple of) 96. This setting has typically been less well explored than asymptotic settings such as the sparse or linear regimes described above, however we provide some bounds in this context.

At various points, in our analysis, we will find it useful to work in terms of the falling moments of various random variables, and we will write $M_{({s})}( {Y} ) := {{\mathbb {E}}} (Y)_{(s)} = {{\mathbb {E}}}Y(Y-1) \ldots (Y-s+1) = {{\mathbb {E}}}Y!/(Y-s)!$ for the $s$th falling moment of random variable $Y$. We will also write $w({\mathbf {u}})$ for the Hamming weight of the binary vector ${\mathbf {u}}$.

2.2. COMP algorithm

The COMP algorithm uses the test outputs $Y$ and matrix $X$ to produce an estimate of the defective set ${{\mathcal {K}}}$ which we will write as $\widehat {{{\mathcal {K}}}}_{\rm COMP}$. In fact, it is easier to consider the complement of $\widehat {{{\mathcal {K}}}}_{\rm COMP}$: an item will appear in this complement if it appears in some negative test. Formally speaking

(2)\begin{equation} \widehat{{{\mathcal{K}}}}_{\rm COMP}^{c} = \{ i: Y_s = 0 \text{ for some } s \text{ with } X_{si} = 1\}. \end{equation}

Notice that if item $j$ really is defective (i.e., $j \in {{\mathcal {K}}}$), then for every test $s$ with $X_{sj} = 1$ then $Y_s = 1$ (by the definition of the group testing action in (1)), so that (2) implies that $j$ is not in $\widehat {{{\mathcal {K}}}}_{\rm COMP}$. In other words, ${{\mathcal {K}}} \subseteq \widehat {{{\mathcal {K}}}}_{\rm COMP}$ (see [Reference Aldridge, Johnson and Scarlett5] Lem. 2.3).

COMP is an attractive algorithm because it is simple to perform and interpret, and because of this performance guarantee in one direction. However, in practice if we do not perform enough tests, then $\widehat {{{\mathcal {K}}}}_{\rm COMP}$ can be significantly larger than ${{\mathcal {K}}}$, meaning that many nondefective items would be misclassified, potentially causing problems in healthcare-related situations where unnecessary quarantine could result.

For this reason, Aldridge [Reference Aldridge2] proposed what he refers to as a conservative two-stage algorithm. Here, given a total budget of $T$ tests, we should use $T_1$ of them to perform the COMP algorithm in the usual way, and then the remaining $T_2 := T-T_1$ tests to perform individual testing of each of the items in $\widehat {{{\mathcal {K}}}}_{\rm COMP}$, which have not been classified as nondefective. While this algorithm may be inferior in performance to a two-stage algorithm which uses the $T_1$ tests to perform the DD algorithm [Reference Aldridge, Baldassini and Johnson4] followed by individual testing, it has the advantage of being transparent to perform for healthcare professionals without a mathematical background.

Clearly, the conservative two-stage algorithm of Aldridge [Reference Aldridge2] will succeed in finding all the defective items if the number of second stage tests is greater than or equal to the number of items to be tested. That is, it will succeed when $T_2 = T-T_1 \geq |\widehat {{{\mathcal {K}}}}_{\rm COMP}|$, meaning that we would like to find the size of $\widehat {{{\mathcal {K}}}}_{\rm COMP}$. Furthermore, for a given budget of $T$ tests, since larger values of $T_1$ give smaller $\widehat {{{\mathcal {K}}}}_{\rm COMP}$ (more stage one tests allow more nondefective items to be screened out) but leave fewer tests available in the second stage, we would like to find a sensible choice of $T_1$ that manages this tradeoff.

2.3. Basic properties of intruding items

We now define the key property we will study in this paper:

More formally, if nondefective item $\ell$ is intruding, then $Y_s = 1$ for every $s$ with $X_{s\ell }=1$, so that (by (2)) we know $\ell \notin \widehat {{{\mathcal {K}}}}_{\rm COMP}^{c}$, so $\ell \in \widehat {{{\mathcal {K}}}}_{\rm COMP}$. In other words, if a nondefective item is intruding, then COMP will mistakenly declare it to be defective. Since nonintruding items are declared to be nondefective by COMP, we know that the size of $\widehat {{{\mathcal {K}}}}_{\rm COMP}$ (which determines the success of the two-stage algorithm of Aldridge [Reference Aldridge2]) is exactly $|\widehat {{{\mathcal {K}}}}_{\rm COMP}| = k + G$.

For a given nondefective item $i$, we can work out the marginal distribution of $G_i$ relatively easily:

3.1. Moments and associated random variables

If the $G_i$ were independent, then the analysis of the distribution of $G$ would be easy: using Lemma 2.2, then $G$ would be binomial with parameters $n-k$ and $(1-p q_0)^{T}$. However, there is a dependence between the $G_i$. In fact, if we learn that a given item is intruding, that suggests there might be more positive tests than average, which would imply that other items are more likely to be intruding.

We formalize this intuition by showing that the $G_i$ are more likely be equal to 1 together than independence would imply. That means that if COMP fails, it is more likely to fail badly (with a large total $G$) than a naive analysis based on Lemma 2.2 might suggest. We prove this in two ways, first by showing in Corollary 3.3 that the $G_i$ are pairwise positively correlated, and second by proving a stronger result (Proposition 3.5) which shows that the $G_i$ have the property of association (see Definition 3.4), which is stronger than positive correlation (see the discussion in [Reference Esary, Proschan and Walkup17] for example).

In fact, we will deduce the positive correlation result (Corollary 3.3) from an expression for the falling moments of $G$ (Proposition 3.1), which may be of independent interest, extending the result for ${{\mathbb {E}}}G$ implicit in [Reference Aldridge2] Thm. 1). We describe the distribution of $G$ as a binomial mixture of binomial distributions as follows. As in [Reference Aldridge, Baldassini and Johnson4], if we let $M_{0}$ be the number of negative tests, then we know that:

(3)\begin{align} M_{0} & \sim {{\rm Bin}}(T,(1-p)^{k}), \end{align}

(4)\begin{align} G \mid {M_{0}=m} & \sim {{\rm Bin}} (n-k,(1-p)^{m}). \end{align}

The first result follows because a test is negative if and only if it contains no defective items, and for Bernoulli testing this occurs independently across tests with probability $q_0 = (1-p)^{k}$. Moreover, the second result follows because a nondefective item is intruding if and only if it does not appear in any of the $M_0$ negative tests, and each nondefective item is present in a given test independently with probability $p$. We can use this to prove the following result:

Proposition 3.1. Under a Bernoulli test design with parameter $p$, the falling moments of $G$ are given by

(5)\begin{equation} M_{({s})}( {G} ) = {{\mathbb{E}}}G(G-1) \ldots (G-s+1) = \binom{n-k}{s} s! ( 1 - q_0 ( 1- (1-p)^{s} ))^{T}, \end{equation}

for any integer $s \geq 0$.

We prove this result using the following intermediate lemma:

Lemma 3.2 (Falling Moments of Binomial Distribution)

Suppose $X \sim {{\rm Bin}}(L,t)$. Then, the $s$th falling moment of $X$ is given by:

(6)\begin{equation} M_{({s})}( {X} ) = \binom{L}{s} s! \cdot t^{s}. \end{equation}

Proposition 3.1 follows using Lemma 3.2 when we recall the distributions of $M_0$ and $G \mid {M_0 = m}$ from (3) and (4), and apply the law of iterated expectation to express $M_{({s})}( {G} ) = {{\mathbb {E}}} [ {{\mathbb {E}}} ( (G)_{s} \,|\, {M_0} ) ]$. We omit the details for brevity.

Note that we can give an alternative proof of Proposition 3.1, using [Reference Harremoës, Johnson and Kontoyiannis20] Lem. 2.2, which gives a multinomial-type expansion for falling factorials based on the Vandermonde identity. This expansion simplifies in the case of binary random variables to give the fact that the falling moment can be expressed as a sum over sets:

(7)\begin{equation} M_{({s})}( {G} ) = s! \sum_{S: |S| = s} {{\mathbb{E}}} \left( \prod_{i \in S} G_{i} \right). \end{equation}

The summation over sets contributes $s! \binom {n-k}{s}$ equal expectation terms, each one of which equals the probability that all elements of a specified set are intruding, which corresponds to the event that none of them ever appear in a negative test, so the result follows by independence of all test items.

Using Proposition 3.1, we can deduce the following result that shows that the $G_i$ have positive pairwise correlation:

Corollary 3.3. For nondefective items $i \neq j$:

(8)\begin{equation} {{\rm Cov}}(G_i, G_j) = ( 1 - q_0 ( 2p - p^{2} ))^{T} - ( 1 - q_0 p)^{2T} \geq 0.\end{equation}

Proof. We consider the variance of $G$ in two different ways:

(9)\begin{equation} {{\rm Var}}(G) = M_{({2})}( {G} ) + {{\mathbb{E}}}G - ({{\mathbb{E}}}G)^{2}, \end{equation}

and (using the symmetry between pairs of $G_i$ implied by the Bernoulli matrix design)

(10)\begin{align} {{\rm Var}}(G) & = \sum_{i} {{\rm Var}}(G_i) + \sum_{i \neq j} {{\rm Cov}}(G_i,G_j) \nonumber\\ & = (n-k) ( {{\mathbb{E}}}G_i - ({{\mathbb{E}}}G_i)^{2} ) + (n-k)(n-k-1) {{\rm Cov}}(G_i, G_j) \nonumber\\ & = {{\mathbb{E}}}G - \frac{({{\mathbb{E}}}G)^{2}}{n-k} + (n-k)(n-k-1) {{\rm Cov}}(G_i, G_j), \end{align}

using the fact that for a binary random variable $Y$ we have ${{\rm Var}}(Y) = {{\mathbb {E}}}Y^{2} - ({{\mathbb {E}}}Y)^{2} = {{\mathbb {E}}}Y - ({{\mathbb {E}}}Y)^{2}$ and that $(n-k) {{\mathbb {E}}}G_i = {{\mathbb {E}}}G$. Now, equating (9) and (10), we obtain that

\begin{align*} & (n-k)(n-k-1) {{\rm Cov}}(G_i, G_j) = M_{({2})}( {G} ) - ({{\mathbb{E}}}G)^{2} \left(1- \frac{1}{n-k} \right)\\ & \quad = (n-k)(n-k-1) [ ( 1 - q_0 ( 1- (1-p)^{2} ))^{T} - ( 1 - q_0 ( 1- (1-p) ) )^{2T}], \end{align*}

using the expressions for $M_{({2})}( {G} )$ and $M_{({1})}( {G} )$ from Proposition 3.1, and the result follows on cancellation.

However, we can prove a stronger property than just positive correlation. Recall the following definition:

Definition 3.4. [17]

Random variables ${\mathbf {X}} = (X_1, X_2, \ldots, X_n)$ are (positively) associated if, for all increasing functions $f$ and $g$,

(11)\begin{equation} {{\mathbb{E}}} ( f({\mathbf{X}}) g({\mathbf{X}}) ) \geq {{\mathbb{E}}} ( f({\mathbf{X}}) ) {{\mathbb{E}}} ( g({\mathbf{X}}) ). \end{equation}

We will prove the following proposition:

Combining Proposition 3.5 with Theorem 3.1 of [Reference Denuit, Dhaene and Ribas12], we find that $G$ is larger, in a convex sense, than a binomial random variable $H$ with parameters $n-k$ and $(1-pq_0)^{T}$. That is, we have that ${{\mathbb {E}}}g(G)\geq {{\mathbb {E}}}g(H)$ for all real-valued functions $g$ with $g(x+1)-2g(x)+g(x-1)\geq 0$ for all positive integers $x$, and for which the expectations exist (where we note from Property 3.4 of [Reference Denuit, Lefèvre and Utev13] that convex ordering on the integers is equivalent to convex ordering on the real line). This formalizes the notion that $G$ is more variable than it would be if the $G_i$ were independent.

The function $g(x)=x!/(x-s)!$ satisfies this convexity condition, since direct calculation gives that $g(x+1)-2g(x)+g(x-1) = s(s-1) (x-1) \ldots (x-s+2)$ in this case. We thus deduce that (see Lemma 3.2 for the value of $M_{({s})}( {H} ))$

(12)\begin{equation} M_{({s})}( {G} )\geq M_{({s})}( {H} )=\binom{n-k}{s}s!(1-pq_0)^{sT}, \end{equation}

that is, the falling moments of $H$ act as lower bounds for those of $G$. Indeed, direct calculation gives that the ratio

(13)\begin{equation} \frac{M_{({s})}( {G} )}{M_{({s})}( {H} )} = \left( \frac{ 1- q_0 (1- (1-p)^{s})}{(1-p q_0)^{s}} \right)^{T} =: R(s)^{T}, \end{equation}

where $R(s+1) - R(s) = (p q_0 (1-q_0)) ( 1 - (1-p)^{s} ) (1-p q_0)^{-s-1} \geq 0$, so that the ratio between successive falling moments of $G$ and $H$ is increasing in $s$.

Some numerical illustration of this is given in Table 1, along with comparison of falling moments of $G$ with those of other distributions which are more suitable than $H$ as approximations of $G$ and which we now discuss in more detail.

Table 1. First four falling moments of approximating distributions in the case $n=500$, $k=10$, $p=0.1$, $T=100$.

True falling moments $M_{({s})}( {G} )$are given by (5). Negative binomial falling moments $M_{({s})}( {Z} )$ are given by (15) with parameters are given by (18). Poisson falling moments are given by $M_{({s})}( {Y} ) = \lambda ^{s}$, where $\lambda = M_{({1})}( {G} )$ is chosen to match the first moment of $G$. Geometric falling moments are given by $M_{({s})}( {X} ) = s! (1/\alpha - 1)^{s}$, where $\alpha = 1/(1+M_{({1})}( {G} ))$ is chosen to match the first moment of $G$. The binomial random variable $H$ and its falling moments are as given in the discussion following Proposition 3.5.

3.2. Negative binomial approximation

In this subsection, we start to show that the distribution of $G$ can be approximated by $Z$, where $Z \sim \text {NB}(r,q)$ follows the negative binomial distribution. For concreteness, we use the parameterization where the probability mass function for the negative binomial distribution is

(14)\begin{equation} {{\mathbb{P}}}(Z=z) = f(z; r, q) := \frac{\Gamma(z+r)}{\Gamma(r) z !} q^{r}(1-q)^{z} \quad \text{for } z=0,1,2 \ldots. \end{equation}

Note that in the case of integer $r$, the normalization constant ${\Gamma (z+r)}/{(\Gamma (r) z !)} = \binom {z+r-1}{z}$, and $Z$ can be interpreted in terms of the number of failures to see the $r$th success in a sequence of Bernoulli trials. However, we do not require $r$ to be an integer here.

The parameter $r$ is sometimes referred to as the dispersion. In this sense, it is worth noting that the case $r=1$ corresponds to a geometric random variable (as mentioned above) and the limiting regime $r \rightarrow \infty$ and $q = r/(r+\lambda )$ (which corresponds to a mean-preserving limit—see Lemma 3.6) gives the mass function of a Poisson random variable with mean $\lambda$. However, in the finite blocklength setting, we will see that $G$ is often well approximated by a distribution with $1 \ll r \ll \infty$, meaning that neither the geometric nor Poisson approximation are valuable.

One natural question is that of which negative binomial distribution (which choice of parameters) to use. We argue that one natural choice is based on a standard moment matching argument—since we need two parameters, we need to set ${{\mathbb {E}}}G = {{\mathbb {E}}}Z$ and ${{\mathbb {E}}}G^{2} = {{\mathbb {E}}}Z^{2}$. In fact, equivalently, since Proposition 3.1 and Lemma 3.6 give closed form expressions for the falling moments of $G$ and $Z$, it is actually easier to solve $M_{({s})}( {G} ) = M_{({s})}( {Z} )$ for $s=1,2$. It is a straightforward exercise involving the Gamma function to prove the following:

Lemma 3.6. (Falling moments of negative binomial distribution)

The falling moments of $Z \sim \text {NB}(r,q)$ are given by:

(15)\begin{equation} M_{({s})}( {Z} ) = \frac{\Gamma(s+r)}{\Gamma(r)}\left(\frac{1-q}{q}\right)^{s}. \end{equation}

Hence, combining Proposition 3.1 and Lemma 3.6, we have successfully matched the moments if

(16)\begin{align} \frac{r (1-q)}{q} & = (n-k) ( 1 - q_0 p )^{T} = M_{({1})}( {G} ), \end{align}

(17)\begin{align} \frac{r (r+1) (1-q)^{2}}{q^{2}} & = (n-k)(n-k-1) ( 1 - q_0 ( 2p-p^{2} ) )^{T} = M_{({2})}( {G} ), \end{align}

or, equivalently, if

(18)\begin{equation} r= \frac{ M_{({1})}( {G} )^{2}}{ M_{({2})}( {G} ) - M_{({1})}( {G} )^{2}} \quad \text{and}\quad q = \frac{ M_{({1})}( {G} )}{ M_{({2})}( {G} ) + M_{({1})}( {G} ) - M_{({1})}( {G} )^{2}}. \end{equation}

In Figure 1, we illustrate the quality of this negative binomial approximation for $G$ in the case $n=500$, $k=10$, and $p=0.1$. We plot the mass function of the negative binomial random variable $Z$ (with parameter choices as in (18)) and the corresponding estimated mass function for $G$ obtained by simulation, for several values of the number of tests $T$ in a nonadaptive algorithm. These experiments show that this negative binomial distribution gives a good approximation to the distribution of $G$ for various different values of $T$.

Figure 1. Probability mass functions of $G$ (blue, obtained by simulation) and approximating negative binomial random variable $Z$ (black). In this case, the group testing parameters are $n=500$, $k=10$, and $p=0.1$. The four plots correspond to $T=60, 80, 100, 120$.

We also find that for many finite blocklength examples the extra flexibility offered by a two parameter approximation means that the negative binomial distribution with parameters given by (18) approximates the low-order falling moments of $G$ (and hence the variance, skewness, and kurtosis) better than either the Poisson or geometric distributions with a single parameter chosen by matching means. This is illustrated in Table 1, again in the case $n=500$, $k=10$, $p=0.1$, $T=100$. In this setting, the dispersion parameter $r$ of $Z$ is $3.66$, which is well separated from both $r=\infty$, corresponding to the Poisson, and $r=1$, corresponding to the geometric.

In general, direct calculation shows that for any $n,p,q,T$ such that $r \geq 1$, then writing $X$, $Z$, $Y$, and $H$ for the geometric, negative binomial, Poisson, and binomial, respectively, defined in Table 1, we have

(19)\begin{equation} M_{({s})}( {X} ) \geq M_{({s})}( {Z} ) \geq M_{({s})}( {Y} ) \geq M_{({s})}( {H} )\quad \text{for all } s \geq 1. \end{equation}

This follows by rewriting the expressions in (19) in the form

$$s! \left( \frac{r(1-q)}{q} \right)^{s} \geq \frac{\Gamma(s+r)}{\Gamma(r)} \left( \frac{(1-q)}{q} \right)^{s} \geq \left( \frac{r(1-q)}{q} \right)^{s} \geq \left( \frac{r(1-q)}{q} \right)^{s} \frac{\binom{n-k}{s} s!}{(n-k)^{s}}$$

which simplifies to

$$s! r^{s} \geq (s+r-1)(s+r-2) \ldots r \geq r^{s} \geq r^{s} \frac{ (n-k) \ldots (n-k-s+1)}{(n-k)^{s}}.$$

Recall from (12) above that $M_{({s})}( {G} ) \geq M_{({s})}( {H} )$ for all $s$. However, it is not the case that $M_{({s})}( {G} ) \geq M_{({s})}( {Y} )$ for all $s$, since $M_{({s})}( {G} ) = 0$ for $s > n-k$, whereas $M_{({s})}( {Y} ) > 0$ for all $s$, although from Table 1 it appears that $M_{({s})}( {G} ) \geq M_{({s})}( {Y} )$ for small $s$ at least.

In the sparse regime described above (in which $k = n^{\theta }$, $p=1/k$ and $T = (c/q_0) k \log (n)$), we note that the dispersion parameter $r$ of (18) tends to infinity as $n\to \infty$. This is equivalent to the fact that $M_{(2)}(G)/M_{(1)}(G)^{2}$ tends to $1$. We can deduce that this holds from (18) by writing

(20)\begin{align} \frac{M_{(2)}(G)}{M_{(1)}(G)^{2}} & = \frac{(n-k)(n-k-1) (1- q_0 (2p-p^{2}) )^{T}}{(n-k)^{2} (1- q_0 p )^{2T}} \nonumber\\ & \simeq \left( 1 + \frac{p^{2} q_0(1-q_0)}{(1-q_0 p)^{2}} \right)^{T} \nonumber\\ & \simeq \exp \left( \frac{p^{2} T q_0(1-q_0)}{(1-q_0 p)^{2}} \right)\nonumber\\ & \simeq \exp \left( \frac{c(1-q_0) \log(n)}{k} \right) \rightarrow 1, \end{align}

since $q_0 p^{2} T = c \log (n)/k$.

Similarly, in the linear regime ($k=\beta n$, $T= c n$) using (20), we obtain

$$\frac{M_{(2)}(G)}{M_{(1)}(G)^{2}} \simeq \exp \left( \frac{ c q_0(1-q_0)}{\beta^{2} n} \right) \rightarrow 1,$$

since $p^{2} T = c/(\beta ^{2} n)$. A Poisson approximation to the distribution of $G$ may thus be appropriate in the large-$n$ limit for the sparse and linear regimes, but the numerical results of this section make it clear that a negative binomial approximation is a more natural choice for finite blocklength applications.

3.3. Convergence of falling moments

Next, we show that, in this framework, matching the first two moments of $G$ and $Z$ ensures that all the falling moments converge.

We can see this informally by simplifying (5) and (15), respectively. The former gives (to leading order)

(21)\begin{align} M_{({s})}( {G} ) & = \binom{n-k}{s} s! ( 1 - q_0 ( 1- (1-p)^{s}))^{T} \nonumber\\ & \simeq ( n-k )^{s} ( 1 - q_0 ( p s ) )^{T} \nonumber\\ & \simeq ( n-k )^{s} \exp ( - q_0 p s T ) = ( (n-k) \exp({-}q_0 p T))^{s}, \end{align}

using the fact that $1-(1-p)^{s} = 1 - (1-p s + O(p^{2})) = ps + O(p^{2})$. The latter gives

(22)\begin{equation} M_{({s})}( {Z} ) = \frac{\Gamma(s+r)}{\Gamma(r)}\left(\frac{1-q}{q}\right)^{s} \simeq \left( \frac{r(1-q)}{q} \right)^{s}. \end{equation}

Note that the moment matching condition of (16) gives that $r(1-q)/q = (n-k) (1-q_0 p)^{T} \simeq (n-k) \exp (-q_0 p T)$, meaning that (21) and (22) agree. A more formal comparison of falling moments is given in the following theorem, where we recall that with the usual choice $p=1/k$ we have $q_0=(1-p)^{k}\approx e^{-1}$:

Theorem 3.7. Consider the number of intruding defectives $G$ and negative binomial $Z$ with parameters given by moment matching, satisfying (18). Under a Bernoulli test design, for any integer $s \geq 1$, if we write $C = q_0(1-q_0)/(1-q_0 p)^{2}$, then the moment ratio satisfies

(23)\begin{equation} \frac{M_{({s})}( {G} )}{M_{({s})}( {Z} )} \geq \left( \frac{ n-k-s}{(n-k)(1 + (s-1)/(2r))} \right)^{s} \left( 1 + \frac{1}{2 }s(s-1) C p^{2} \left( 1 - \frac{(s-2) (1-2q_0) p}{3 (1-q_0 p)} \right) \right)^{T}, \end{equation}

and

(24)\begin{equation} \frac{M_{({s})}( {G} )}{M_{({s})}( {Z} )} \leq \exp ( s(s-1) C p^{2} T (1-q_0 p)^{2-s} ). \end{equation}

Hence, for any fixed $s$, the ratio $M_{({s})}( {G} )/M_{({s})}( {Z} ) \rightarrow 1$, for both

1. 1. the sparse regime with $k = n^{\theta }$ for some $\theta \in (0,1)$ and $T= c e k \log n$, and

2. 2. the linear regime with $k = \beta n$ and $T = c n$.

Here, we control the upper bound in (24) using the fact that (see Section 3.2) the $p^{2} T$ is $c \log n/(q_0 k)$ or $c/(\beta ^{2} n)$, respectively. Similarly, we control the lower bound (23) using the fact that in both regimes the $n-k$ and dispersion parameter $r$ tend to infinity (again see Section 3.2).

Additionally, the bounds (23) and (24) allow us to control the ratio $M_{({s})}( {G} )/M_{({s})}( {Z} )$ in the finite blocklength regime, deducing bounds that can be compared with the concrete values given in Table 1 for example.

3.4. Stein–Chen method

Having seen in Sections 3.2 and 3.3 that a negative binomial distribution seems to be a reasonable approximation for the distribution of $G$, in this section, we adapt the Stein–Chen method to prove explicit error bounds in the approximation of $G$ by a negative binomial distribution. We emphasize that these bounds apply for any finite blocklength application. The error in our approximation of $G$ by $Z$ will be measured in total variation distance, defined by

(25)\begin{equation} d_{{\rm TV}}(G,Z)=\sup_{A\subseteq{{\mathbb{Z}}}^{+}}|{{\mathbb{P}}}(G\in A)-{{\mathbb{P}}}(Z\in A)|=\inf_{(G,Z)}{{\mathbb{P}}}(G\not=Z), \end{equation}

where ${{\mathbb {Z}}}^{+}=\{0,1,\ldots \}$ and the infimum is taken over all couplings of $G$ and $Z$.

First, we briefly review the Stein–Chen method in the context of negative binomial approximation. For a more detailed introduction to this technique more generally, we refer the reader to [Reference Ross28]. Recall from [Reference Brown and Phillips8] that $Z \sim \text {NB}(r,q)$ if and only if

(26)\begin{equation} {{\mathbb{E}}} [ (1-q) (r+Z) g(Z+1) - Z g(Z) ] = 0, \end{equation}

for all test functions $g:\mathbb {Z}^{+}\to \mathbb {R}$ for which the expectation exists. This follows as a consequence of the fact that the negative binomial probability mass function of (14) satisfies $z {{\mathbb {P}}}(Z = z) = (1-q) (r+z-1) {{\mathbb {P}}}(Z = z-1)$.

The key to the analysis is that for each set $A \subseteq {{\mathbb {Z}}}^{+}$ we can define a function $f_A:\mathbb {Z}^{+}\to \mathbb {R}$ which satisfies $f_A(0)=0$ and the Stein–Chen equation

(27)\begin{equation} (1-q)(r+z) f_A(z+1) - z f_A(z) = {{\mathbb{I}}}(z \in A) - {{\mathbb{P}}}(Z \in A), \end{equation}

for all $z\in \mathbb {Z}^{+}$. Then, for any random variable $Y$, we can take the expectation of (27) over $Y$ to obtain

(28)\begin{equation} {{\mathbb{E}}} ( (1-q)(r+Y) f_A(Y+1) - Y f_A(Y) ) = {{\mathbb{P}}}(Y \in A) - {{\mathbb{P}}}(Z \in A). \end{equation}

Note that (as expected) if $Y$ were negative binomial, then both RHS and LHS of (28) would be zero (the latter due to the characterization in (26)). However, we can deduce that if the LHS of (28) is small, then so is the RHS. Indeed, if the LHS of (28) is small uniformly over choices of $A$, then we can deduce a bound in total variation distance since, combining (25) and (27), we have

(29)\begin{equation} d_{{\rm TV}}(Y, W) = \sup_{A \subseteq{{\mathbb{Z}}}^{+}} | {{\mathbb{E}}} ( (1-q)(r+Y) f_A(Y+1) - Y f_A(Y) )|. \end{equation}

Having reviewed the Stein–Chen method in general, we will now describe how we bound (28) in this specific case. For our approximating negative binomial distribution for $G$, we will make the following choices of the parameters $q$ and $r$:

(30)\begin{equation} q=\frac{\mu}{\sigma^{2}}\quad\text{and}\quad r=\frac{\mu^{2}}{\sigma^{2}-\mu}=\frac{1}{e^{Tp^{2} q_0}-1}, \end{equation}

where

$$\mu=(n-k)e^{{-}Tp q_0}\quad\text{and}\quad\sigma^{2}=(n-k)^{2}(e^{{-}Tp(2-p)q_0}-e^{{-}2Tp q_0} )+(n-k)e^{{-}Tp q_0},$$

and where, as before, we write $q_0 = (1-p)^{k}$.

We remark that these parameter choices do not match the first two moments of $Z$ with those of $G$, unlike those of (18). Instead, the parameters $q$ and $r$ are chosen to match the first two moments of $Z$ with those of $G^{\prime \prime }$, to be defined precisely below, in which the binomial mixture which defines $G$ is replaced by a particular Poisson mixture. This may seem a little unnatural at first, but makes sense in the setting of the proof in Appendix C, and the ultimate effect should be negligible since the distributions of $G$ and $G^{\prime \prime }$ are close, as our proof demonstrates.

Our main result is the following.

Theorem 3.8. Let $G$ be as above, and let $Z$ have a negative binomial distribution with parameters $q$ and $r$ given by (30). Then, defining $K = e^{T p q_0}$, we have

(31)\begin{align} d_{{\rm TV}}(G,Z)& \leq 2\min\left\{\frac{q_0}{4\sqrt{1-q_0}},\frac{1}{\sqrt{T}}\alpha(q_0)+\frac{1}{\sqrt{2\pi e}}\log\left(\frac{1}{\sqrt{1-q_0}}\right)\right\}+\frac{1}{K}\nonumber\\ & \quad + \frac{(2-q)(n-k)}{1-q} \left(e^{r+1} K^{r} \exp({-}K r) +\int_0^{1}\left|\widehat{\Gamma}\left(\left\lceil\frac{\log(x)}{\log(1-p)}\right\rceil,Tq_0 \right)-\widehat{\Gamma}(r,K r x)\right|dx\right), \end{align}

where

$$\alpha(q_0)=\frac{0.4748[\sqrt{1-q_0}(1+2q _0^{2} e^{{-}q_0})+ q_0^{2}+(1-q_0)^{2}]}{\sqrt{q_0(1-q_0)}},$$

and $\widehat {\Gamma }(\cdot,\cdot )$ is the normalized upper incomplete gamma function, defined by

$$\widehat{\Gamma}(s,y)= \frac{1}{\Gamma(s)} \int_y^{\infty} t^{{-}s-1}e^{{-}t}\,dt,$$

where $\Gamma (\cdot )$ is the gamma function.

We conclude this section with some discussion and numerical illustration of the bound of Theorem 3.8. We will discuss further aspects of the convergence of this bound in Remark 3.9, but first use numerical illustrations to gain some initial understanding of its behavior. Firstly, we note that although our bound applies for any finite blocklength, there are examples in which it is worse than the trivial bound $d_{{\rm TV}}(G,Z)\leq 1$, or performs poorly compared with the simulation results we observed in Section 3.2. For example, with $n=500$, $k=10$, $p=0.1$, and $T=100$, the bound of Theorem 3.8 gives $d_{{\rm TV}}(G,Z)\leq 1.80$, which is clearly uninformative. We give some further examples of our upper bound in Table 2, all in the case $n=2500$ and for various values of $k$, $p$, and $T$, from which it is clear that there are some cases in which the bound of Theorem 3.8 performs very well, and demonstrates proximity of the distribution of $G$ to negative binomial.

Table 2. The upper bound of Theorem 3.8 in the case $n=2500$ for various values of $k$, $p$, and $T$.

The symbol “—” indicates that the upper bound is larger than 1, and therefore uninformative.

It is interesting to note that in many cases the largest share of the error estimate in Theorem 3.8 comes from the first term of the upper bound, which arises from the approximation of the binomial random variable $M_0$ by a Poisson random variable of the same mean (see the proof in Appendix C for details). Nevertheless, the bound we would obtain without making this approximation generally performs worse than our Theorem 3.8. We conjecture that this is because the integral in the final term of the upper bound is made smaller by this approximation of $M_0$ (compared with the corresponding term without this approximation), resulting in a smaller upper bound overall because of the relatively large factors multiplying this integral.

Remark 3.9. While the incomplete gamma functions in (31) make the integral in that bound not straightforward to interpret directly, we can provide an upper bound on this quantity using concentration of measure inequalities.

We split the region of integration in three parts, with breaks at $(1 \pm \epsilon )/K$. On the region $( (1-\epsilon )/K, (1+\epsilon )/K )$, we simply bound the integrand by 1, to give $2 \epsilon /K$. Using the fact that for $0 \leq u,v \leq 1$ we can bound $|u-v| \leq \max (u,v)$ and $|u-v| \leq \max (1-u,1-v)$, we can hence bound the integral in (31) by

\begin{align*} & \frac{2 \epsilon}{K} + \frac{1}{K} \max \left( {{\mathbb{P}}} \left( \xi^{\prime} < \frac{(1-\epsilon)}{K} \right), {{\mathbb{P}}} \left(\eta^{\prime} < \frac{(1-\epsilon)}{K} \right) \right)\\ & \quad + \max \left( {{\mathbb{P}}} \left( \xi^{\prime} > \frac{(1+\epsilon)}{K} \right), {{\mathbb{P}}} \left(\eta^{\prime} > \frac{(1+\epsilon)}{K} \right) \right), \end{align*}

where $\eta ^{\prime } \sim \Gamma (r,rK)$ and $\xi ^{\prime }=(1-p)^{M^{\prime }}$, with $M^{\prime }\sim \text {Po}(Tq_0)$; see also (C.9) in the proof of Theorem 3.8.

The inequalities (C.5) and (C.6) below give us that we may upper bound both ${{\mathbb {P}}}(\eta '>z)$ and ${{\mathbb {P}}}(\eta '< z)$ by

$$(K z e)^{r} \exp({-}r K z) = \exp ( r (1 - K z + \log(K z)) ).$$

Hence, for example, we know by writing $m(s) = s - \log (1+s) \simeq s^{2}/2$ that

$${{\mathbb{P}}} \left(\eta^{\prime} > \frac{(1+\epsilon)}{K} \right) \leq \exp ( - r m(\epsilon) )\quad \text{and}\quad {{\mathbb{P}}} \left(\eta^{\prime} < \frac{(1-\epsilon)}{K} \right) \leq \exp ({-}r m(-\epsilon) ).$$

A similar standard Chernoff bounding argument, as used for (C.5) and (C.6), gives that, for $Y \sim {{\rm Po}}(\lambda )$, both ${{\mathbb {P}}}( Y > y)$ and ${{\mathbb {P}}}(Y < y)$ are bounded above by $\exp ( - \lambda h(y/\lambda - 1) )$, where $h(s) = (1+s) \log (1+s) - s \simeq s^{2}/2$. Hence, for example, we can bound

\begin{align*} {{\mathbb{P}}}( \xi^{\prime} >x) & = {{\mathbb{P}}} \left( M' < \frac{-\log x}{ -\log(1-p)} \right) \leq {{\mathbb{P}}} \left( M' < \frac{-\log x}{ p} \right)\\ & \leq \exp \left(- T q_0 h \left( \frac{-\log x}{T q_0 p} - 1 \right) \right), \end{align*}

if $-\log x/p \leq {{\mathbb {E}}}M' = T q_0$, since $M' \sim {{\rm Po}}(T q_0)$. Hence, taking $x= (1+\epsilon )/K = (1+\epsilon ) e^{-T p q_0}$, we obtain

$${{\mathbb{P}}} \left( \xi^{\prime} > \frac{(1+\epsilon)}{K} \right) \leq \exp \left(- T q_0 h \left( - \frac{\log (1+\epsilon)}{T q_0 p} \right) \right).$$

Similarly, we can bound

$${{\mathbb{P}}}( \xi^{\prime} < x) = {{\mathbb{P}}} \left( M' > \frac{-\log x}{ -\log(1-p)} \right) \leq \exp \left(- T q_0 h \left( \frac{-\log x}{-T q_0 \log(1-p)} - 1 \right) \right),$$

if $-\log x/-\log (1-p) \leq {{\mathbb {E}}}M' = T q_0$, since $M' \sim {{\rm Po}}(T q_0)$. Hence, taking $x= (1-\epsilon )/K = (1-\epsilon ) e^{-T p q_0}$ and writing $m(-p) = - p - \log (1-p) \geq 0$ as above, we obtain that

$${{\mathbb{P}}} \left( \xi^{\prime} < \frac{(1-\epsilon)}{K} \right) \leq \exp \left( - T q_0 h \left( - \frac{-\log(1-\epsilon) - T q_0 m({-}p)}{-T q_0 \log(1-p)} \right) \right),$$

assuming that the numerator is positive, which holds, for example, if $\epsilon > T q_0 m(-p)$. This condition is satisfied, for example, for large $n$ in the linear regime where $T=c n$, $k = \beta n$, and $p = 1/k$, since $m(p) \simeq p^{2}/2$ so $T q_0 m(-p) \sim {\rm const.}/n$ in this regime.