Goodness-of-Fit Test for the Bivariate Negative Binomial Distribution

Abstract

When modeling real-world data, we face the challenge of determining which probability distribution best represents the data. To address this intricate problem, we rely on goodness-of-fit tests. However, when the data come from a bivariate negative binomial distribution, the literature reveals no existing goodness-of-fit test for this distribution. For this reason, in this article, we propose and study a computationally convenient goodness-of-fit test for the bivariate negative binomial distribution. This test is based on a bootstrap approximation and a parallelization strategy. To this end, we use a reparameterization technique based on the probability generating function and a Cramér-von Mises-type statistic. From the simulation studies, we conclude that the results converge to the established nominal levels as the sample size increases, and in all cases considered, the parametric bootstrap method provides an accurate approximation of the null distribution of the statistic we propose. Additionally, we verify the power of the proposed test, as well as its application to five real datasets. To accelerate the massive computational work, we employ the parallelization strategy that, according to Novoa-Muñoz (2024), was the most efficient among the techniques he analyzed.

1. Introduction

In today’s world, where large volumes of data proliferate, it is common to encounter count data, with the Univariate Negative Binomial Distribution (UNBD) being a useful discrete method for representing such data. This distribution has been used to model insurance policy allocation in actuarial practice (see Shi and Valdez [1]), football betting (see Van Gemert and Van Ophem [2]), and, interestingly, the population dynamics of invasive species (see Krkošek et al. [3]).

By extending the idea of modeling, the Bivariate Negative Binomial Distribution (BNBD) was developed, which is used to describe phenomena involving workplace accidents (see Arbus and Kerrich [4]; Edwards and Gurland [5]), ecological data (see Holgate [6]), migration data related to traffic accidents (see Maher [7]), and bus accident data within London’s vehicle fleet (see Kopocińska [8]).

But in a real-life situation, when bivariate count data are available, how can one determine if the BNBD is appropriate for modeling such data? According to González-Albornoz and Novoa-Muñoz [9], a crucial aspect of data analysis is testing the goodness of fit (gof) of given observations with a probabilistic model. In this regard, and to the best of our knowledge, the only goodness-of-fit tests related to the Negative Binomial distribution that we have found in the literature are those developed by (a) Heller [10], which was applied to a set of microbiological data; (b) Beltrán-Beltrán and O’Reilly [11], who applied it to two datasets (deaths caused by horse kicks and laboratory mice); and (c) Meintanis [12], who applied it to accident data in the Netherlands over the period 1997–2004. However, cases (a) and (b) are tests for univariate data, while case (c) is bivariate but corresponds to a Poisson–Negative Binomial distribution model. More recently, Hudecová et al. [13] and Wang et al. [14] published goodness-of-fit tests for bivariate count time series based on a bivariate Poisson model. Moreover, as far as we know, there is no literature on gof tests for the BNBD. For this reason, the objective of this article is to propose and study a gof test for the BNBD that is consistent against any fixed alternative.

The test we propose employs the probability-generating function (pgf) of the BNBD and its counterpart, the empirical probability-generating function (epgf). The foundation of this proposal lies in the fact that the pgf characterizes the distribution of a random vector and can be consistently estimated by the epgf (see for example Novoa-Muñoz [15]; Novoa-Muñoz and Jiménez-Gamero [16]). This statistical test compares the epgf of the data with an estimator of the pgf of the BNBD.

To decide when to reject the null hypothesis, we need to know the null distribution of the statistic or at least an approximation of it. Since the proposed statistic is of the Cramér-von Mises type, obtaining the null distribution for finite sample sizes is not feasible. The asymptotic alternative we used was to approximate the null distribution using a parametric bootstrap estimator, whose properties can be found in Novoa-Muñoz and Jiménez-Gamero [17].

Since infinite sample sizes do not exist in the real world, we conducted a simulation study to evaluate the performance of the proposed test for finite sample sizes. Due to the computational intensity of the parametric bootstrap method, we used the parallel algorithm parRapply to reduce the simulation runtime. This algorithm proved to be the most efficient among the parallel versions analyzed in Novoa-Muñoz [18].

The present study follows the following sequence: In Section 2, we present the pgf of the DBNB and a reparameterization that makes it computationally more efficient. In Section 3, we develop the test statistic we propose, which, since it does not have a known distribution, we show can be approximated by its null asymptotic distribution. However, as this approximation still depends on the true parameter, in Section 4, we propose a parametric bootstrap estimator. Section 5 is dedicated to presenting the three techniques we used to estimate the parameters of the BNBD. In Section 6, we present and discuss the results of a simulation study and the application of the proposed test to two real datasets. The simulation study was conducted to evaluate the performance of the proposed test and analyze its power.

The notation we use below is detailed in the abbreviations listed at the end of the manuscript.

2. Reparametrization of the Probability Generating Function of the BNBD

From now on, we will consider that $\gamma =({\gamma }_0,{\gamma }_1,{\gamma }_2)\in \mathrm{\Theta }$, where

$$\mathrm{\Theta }=\left\{\left({\gamma }_0,{\gamma }_1,{\gamma }_2\right)\in {\mathbb{R}}^3:{\gamma }_0>{\gamma }_2,{\gamma }_1>{\gamma }_2,{\gamma }_2>0\right\},$$

is the respective parameter space.

On the other hand, according to Kocherlakota and Kocherlakota [19], the genesis of the BNBD depends on the underlying random mechanism. One such mechanism, composition, gives rise to the following probability generating function (pgf):

$$g(t_1,t_2)={(1-p_1-p_2-p_3)}^v{\left(1-p_1{t}_1-p_2{t}_2-p_3{t}_1{t}_2\right)}^{-v},$$

(1)

where $p_1,p_2,p_3\in [0,1]$, $v\in \mathbb{N}$ represents the number of trials before the first success or failure.

From (1), Edwards and Gurland [5] and Subrahmaniam [20] derived the probability mass function of the BNBD given below (see Min et al. [21]):

$$P(X_1=r,X_2=s)={(1-p_1-p_2-p_3)}^v\sum _{i=0}^{min(r,s)}{\displaystyle \frac{\mathrm{\Gamma }(v+r+s-i)}{\mathrm{\Gamma }\left(v\right)i!(r-i)!(s-i)!}}{p}_1^{r-i}{p}_2^{s-i}{p}_3^i,$$

where $\mathrm{\Gamma }$ is the Gamma function.

For computational purposes and to facilitate the estimation of the distribution parameters, it is highly convenient to incorporate the following reparameterization:

$$p_1={\displaystyle \frac{{\gamma }_0-{\gamma }_2}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}},p_2={\displaystyle \frac{{\gamma }_1-{\gamma }_2}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}},p_3={\displaystyle \frac{{\gamma }_2}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}}.$$

Thus, the random vector $\left(X_1,X_2\right)$ follows a Bivariate Negative Binomial distribution with parameter $\gamma =\left({\gamma }_0,{\gamma }_1,{\gamma }_2\right)$, which we will denote as $\left(X_1,X_2\right)\sim BNB\left(\gamma \right)$, and its probability mass function (pmf), with which we will work in this article, is given by

$$P_{\gamma }(X_1=r,X_2=s)={\displaystyle \frac{{({\gamma }_0-{\gamma }_2)}^r{({\gamma }_1-{\gamma }_2)}^s\mathrm{\Gamma }(v+r+s)}{r!s!\mathrm{\Gamma }\left(v\right){(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}^{v+r+s}}}S(r,s),$$

(2)

where

$$S(r,s)=\sum _{i=0}^{min(r,s)}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{r}{i}\right)\left(\genfrac{}{}{0pt}{}{s}{i}\right)}{\left(\genfrac{}{}{0pt}{}{v+r+s-1}{i}\right)}}{\tau }^i\mathrm{with}\tau ={\displaystyle \frac{{\gamma }_2(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}{({\gamma }_0-{\gamma }_2)({\gamma }_1-{\gamma }_2)}}.$$

With this reparameterization and considering the parameter $\gamma$, the pgf (1) of the BNBD can be written as

$$g(t;\gamma ):=E\left(t_1^{X_1}{t}_2^{X_2}\right)={\left(1+(1-t_1){\gamma }_0+(1-t_2){\gamma }_1-(t_1-1)(t_2-1){\gamma }_2\right)}^{-v},$$

where $t=(t_1,t_2)\in {\mathbb{R}}^2$ and $\gamma \in \mathrm{\Theta }$.

3. Cramér-Von Mises-Type Statistic and Its Asymptotic Null Distribution

In this section, we introduce the test statistic, which is of the Cramér-von Mises type. The essence of the Cramér-von Mises statistic lies in measuring the distance between two continuous distribution functions. In particular, we use it to measure the distance between the pgf, $g(t;\gamma )$, and its empirical counterpart, the epgf $g_n\left(t\right)$, which is defined below. For a detailed explanation of the Cramér-von Mises distance, see Baringhaus and Henze [22].

In order to test the hypothesis

$$H_0:(X_1,X_2)\sim BNB\left(\gamma \right),\mathrm{for}\mathrm{some}\gamma \in \mathrm{\Theta },$$

against the alternative

$$H_1:(X_1,X_2)\nsim BNB\left(\gamma \right),\forall \gamma \in \mathrm{\Theta },$$

let ${\mathit{X}}_1=(X_{11},X_{12}),\dots ,{\mathit{X}}_n=(X_{n1},X_{n2})$ be independent and identically distributed (iid) random vectors defined on a probability space $(\mathrm{\Omega },\mathcal{A},P)$ and taking values in ${\mathbb{N}}_0^2$. Moreover, let

$${\displaystyle g_n\left(t\right)=\frac{1}{n}\sum _{i=1}^n{t}_1^{X_{i1}}{t}_2^{X_{i2}},}$$

be the epgf of ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ for some appropriate $W\subseteq {\mathbb{R}}^2$.

To develop our statistical test, the following result is fundamental, and its proof, for the d-dimensional case, can be reviewed in Novoa-Muñoz and Jiménez-Gamero [17]:

Proposition 1.

Let ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ be iid from a random vector $\mathit{X}=(X_1,X_2)\in {\mathbb{N}}_0^2$. Let $g\left(t\right)=E\left(t_1^{X_1}{t}_2^{X_2}\right)$ be the pgf of $\mathit{X}$, which is defined on $W\subseteq {\mathbb{R}}^2$. Let $0\le b_j\le c_j<\infty ,j=1,2$ such that $R=[b_1,c_1]\times [b_2,c_2]\subseteq W$; then,

$$\underset{t\in R}{sup}|g_n\left(t\right)-g\left(t\right)|\stackrel{a.s.}{⟶}0.$$

According to Proposition 1, the empirical probability-generating function (epgf) consistently estimates the pgf. Moreover, if ${\widehat{\gamma }}_n$ is a consistent estimator of $\gamma$ and $H_0$ is true, then the pgf is consistently estimated by $g(t;{\widehat{\gamma }}_n)$, since g is a continuous function. On the other hand, because the distribution of $\mathit{X}$ is determined solely by its pgf, $g\left(t\right),t\in {[0,1]}^2$, it is reasonable to think that $H_0$ should reject the null hypothesis for large values of $B_{n,w}\left({\widehat{\gamma }}_n\right)$, which are defined by

$$B_{n,w}\left({\widehat{\gamma }}_n\right)={\int }_{{[0,1]}^2}{B}_n^2(t;{\widehat{\gamma }}_n)w\left(t\right)dt,$$

(3)

where

$$B_n(t;\gamma )=\sqrt{n}\left\{g_n\left(t\right)-g(t;\gamma )\right\},$$

${\widehat{\gamma }}_n={\widehat{\gamma }}_n({\mathit{X}}_1,{\mathit{X}}_2,\dots ,{\mathit{X}}_n)$ is a consistent estimator of $\gamma$. Additionally, $w\left(t\right)$ is a measurable and non-negative weight function to which we will impose the condition of being finite to ensure that the double integral in (3) is finite for each fixed n. That is,

$$0\le {\int }_{{[0,1]}^2}w\left(t\right)dt<\infty .$$

(4)

The next step is to determine what constitutes large samples of $B_{n,w}\left({\widehat{\gamma }}_n\right)$. To do this, the first alternative is to understand its null distribution. However, since $B_{n,w}\left({\widehat{\gamma }}_n\right)$ is a Cramér-von Mises-type statistic, it does not have a known distribution. Therefore, we will estimate it through its asymptotic null distribution. To achieve this, we will assume that the estimator ${\widehat{\gamma }}_n$ satisfies the following regularity condition:

Assumption 1.

Under $H_0$, if $\gamma =({\gamma }_0,{\gamma }_1,{\gamma }_2)\in \mathrm{\Theta }$ denotes the true parameter value, then

$$\sqrt{n}\left({\widehat{\gamma }}_n-\gamma \right)={\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^n\ell \left({\mathit{X}}_i;\gamma \right)+{\mathbf{o}}_P\left(1\right),$$

where $\ell :{\mathbb{N}}_0^2\times \mathrm{\Theta }⟶{\mathbb{R}}^3$ is such that $E_{\gamma }\left\{\ell ({\mathit{X}}_1;\gamma )\right\}=\mathbf{0}$ and $J\left(\gamma \right)=E_{\gamma }\left\{\ell {\left({\mathit{X}}_1;\gamma \right)}^{\top }\ell \left({\mathit{X}}_1;\gamma \right)\right\}<\infty$. Here, ${\mathbf{o}}_P\left(1\right)$ is a vector consisting of three $o_P\left(1\right)$ elements. 

Assumption 1 is fulfilled by most commonly used estimators; see [19,23].

The regularity conditions for a parameter estimator have important practical implications, as they ensure the validity and good behavior of the estimators obtained. These conditions allow for the derivation of asymptotic properties such as consistency, asymptotic normality, and efficiency. Some key practical implications are the following:

• Ensuring the consistency of the estimator: This guarantees that the estimator converges to the true value of the parameter as the sample size increases. This is crucial for ensuring that the model’s results are representative of the underlying phenomenon.
• Asymptotic normality: This allows the distribution of the estimator to approximate a normal distribution as the sample size grows. This facilitates inference, such as hypothesis testing and confidence interval construction.
• Efficiency of the estimator: This can enable the estimator to reach the Cramér–Rao bound, meaning it has the lowest possible variance among unbiased estimators.
• Mathematical interpretability and proper modeling: They ensure that the probability or likelihood functions have appropriate properties, such as being differentiable, non-degenerate, and well-defined. This enables reliable analytical derivations and numerical computations.
• Robustness to minor violations: Although some regularity conditions may be difficult to verify in practice, this provides a solid theoretical framework that allows models to adapt to small violations or approximations.

In summary, regularity conditions ensure that estimators are reliable, accurate, and suitable for making valid inferences, which is essential in practical statistical applications, such as scientific research, economics, biology, and other fields.

The next result gives the asymptotic null distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$.

Theorem 1.

Let ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ be iid from $\mathit{X}=(X_1,X_2)\sim BNB\left(\gamma \right)$. Suppose that Assumption 1 holds. Then,

$$B_{n,w}\left({\widehat{\gamma }}_n\right)=\left|\right|V_n{\left|\right|}_{\mathcal{H}}^2+o_P\left(1\right),$$

where $V_n\left(t\right)={\displaystyle \frac{1}{\sqrt{n}}}{\sum }_{i=1}^n{B}^0({\mathit{X}}_i;t;\gamma )$, with

$$B^0({\mathit{X}}_i;t;\gamma )=t_1^{X_{i1}}{t}_2^{X_{i2}}-g(t;\gamma )-v{\left\{g(t;\gamma )\right\}}^{{\displaystyle \frac{v+1}{v}}}\left(t_1-1,t_2-1,(t_1-1)(t_2-1)\right)\ell {\left({\mathit{X}}_i;\gamma \right)}^{\top },$$

$i=1,\dots ,n$. Moreover,

$$B_{n,w}\left({\widehat{\gamma }}_n\right)\stackrel{L}{⟶}\sum _{j\ge 1}{\lambda }_j{\chi }_{1j}^2,$$

(5)

where ${\chi }_{11}^2,{\chi }_{12}^2,\dots$ are independent ${\chi }^2$ variates with one degree of freedom, and the set $\left\{{\lambda }_j\right\}$ is the non-null eigenvalues of the operator $C\left(\gamma \right)$ defined on the function space $\{\xi :{\mathbb{N}}_0^2\to \mathbb{R}, \mathit{such}\mathit{that}{E}_{\gamma }\left\{{\xi }^2\left(\mathit{X}\right)\right\}<\infty ,\forall \gamma \in \mathrm{\Theta }\}$ as follows:

$$C\left(\gamma \right)\xi \left(\mathit{x}\right)=E_{\gamma }\left\{h(\mathit{x},\mathit{Y};\gamma )\xi \left(\mathit{Y}\right)\right\},$$

where

$$h(\mathit{x},\mathit{y};\gamma )={\int }_{{[0,1]}^2}{B}^0(\mathit{x};t;\gamma )B^0(\mathit{y};t;\gamma )w\left(t\right)dt.$$

(6)

Proof. 

By definition, $B_{n,w}\left({\widehat{\gamma }}_n\right)={\parallel B_n(t;{\widehat{\gamma }}_n)\parallel }_{\mathcal{H}}^2$. Note that

$$B_n(t;{\widehat{\gamma }}_n)={\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^{n}B({\mathit{X}}_i;t;{\widehat{\gamma }}_n),\mathrm{with}B({\mathit{X}}_i;t;\gamma )=t_1^{X_{i1}}{t}_2^{X_{i2}}-g(t;\gamma ).$$

(7)

By Taylor expansion of $B({\mathit{X}}_i;t;{\widehat{\gamma }}_n)$ around $\gamma$,

$$\begin{array}{cc}\hfill B_n(t;{\widehat{\gamma }}_n)={\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^{n}B({\mathit{X}}_i;t;\gamma )& +{\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}^{\left(1\right)}({\mathit{X}}_i;t;\gamma )\sqrt{n}{({\widehat{\gamma }}_n-\gamma )}^{\top }\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\sqrt{n}({\widehat{\gamma }}_n-\gamma ){\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}^{\left(2\right)}({\mathit{X}}_i;t;\gamma )\sqrt{n}{({\widehat{\gamma }}_n-\gamma )}^{\top }+r_n,\hfill \end{array}$$

(8)

where

$$r_n={\left.{\displaystyle \frac{1}{3!}}{\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^n\sum _{k_1+k_2+k_3=3}{\displaystyle \frac{{\partial }^{3}B({\mathit{X}}_i;t;\gamma )}{\partial {\gamma }_0^{k_1}\partial {\gamma }_1^{k_2}\partial {\gamma }_2^{k_3}}}\right|}_{\gamma =\tilde{\gamma }}\left({\widehat{\gamma }}_{0n}-{\gamma }_0\right)\left({\widehat{\gamma }}_{1n}-{\gamma }_1\right)\left({\widehat{\gamma }}_{2n}-{\gamma }_2\right),$$

and $\tilde{\gamma }=\alpha {\widehat{\gamma }}_n+(1-\alpha )\gamma$, for some $0<\alpha <1$, $B^{\left(1\right)}(\mathit{x};t;\tau )$, is the vector of the first derivatives, and $B^{\left(2\right)}(\mathit{x};t;\tau )$ is the matrix of the second derivatives of $B(\mathit{x};t;\tau )$ with respect to $\tau =({\tau }_0,{\tau }_1,{\tau }_2)$. That is,

$$B^{\left(1\right)}(x;t;\tau )=\left(B_1^{\left(1\right)}(x;t;\tau ),B_2^{\left(1\right)}(x;t;\tau ),B_3^{\left(1\right)}(x;t;\tau )\right),B_j^{\left(1\right)}(x;t;\tau )={\displaystyle \frac{\partial B(x;t;\tau )}{\partial {\tau }_{j-1}}},j=1,2,3.$$

As the first derivatives are

$$B_1^{\left(1\right)}(x;t;\tau )=v(1-t_1)q,B_2^{\left(1\right)}(x;t;\tau )=v(1-t_2)q,B_3^{\left(1\right)}(x;t;\tau )=-v(t_1-1)(t_2-1)q,$$

(9)

where $q={\left\{g(t;\tau )\right\}}^{{\displaystyle \frac{v+1}{v}}}$, then $\left|B_j^{\left(1\right)}\left({\mathit{X}}_1;t;\tau \right)\right|\le v<\infty$ for $j=1,2,3$$\forall t\in {[0,1]}^2$.

Thus, considering (4), it results that

$$E_{\gamma }\left\{{∥B_j^{\left(1\right)}\left({\mathit{X}}_1;t;\gamma \right)∥}_{\mathcal{H}}^2\right\}<\infty ,j=1,2,3.$$

(10)

Since, for $j=1,2,3$, we have that

$$E_{\gamma }\left[{\left\{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}_j^{\left(1\right)}({\mathit{X}}_i;t;\gamma )\right\}}^2\right]={\displaystyle \frac{1}{n}}{E}_{\gamma }\left[{\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}}^2\right]+{\displaystyle \frac{n-1}{n}}{\left[E_{\gamma }\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}\right]}^2,$$

then

$$\begin{array}{cc}\hfill E_{\gamma }\left[{\left\{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}_j^{\left(1\right)}({\mathit{X}}_i;t;\gamma )-E_{\gamma }\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}\right\}}^2\right]& \hfill \\ \hfill & ={\displaystyle \frac{1}{n}}{E}_{\gamma }\left[{\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}}^2\right]-{\displaystyle \frac{1}{n}}{\left[E_{\gamma }\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}\right]}^2\hfill \\ \hfill & \le {\displaystyle \frac{1}{n}}{E}_{\gamma }\left[{\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}}^2\right].\hfill \end{array}$$

Using the Markov inequality and (10) for $j=1,2,3$, we obtain

$$P_{\gamma }\left[{∥{\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}_j^{\left(1\right)}({\mathit{X}}_i;t;\gamma )-E_{\gamma }\left\{B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}∥}_{\mathcal{H}}>\epsilon \right]\le {\displaystyle \frac{1}{n{\epsilon }^2}}{E}_{\gamma }\left[{∥B_j^{\left(1\right)}({\mathit{X}}_1;t;\gamma )∥}_{\mathcal{H}}^2\right]\to 0.$$

Thus, in the space $\mathcal{H}$, we obtain that

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}^{\left(1\right)}({\mathit{X}}_i;t;\gamma )\stackrel{P}{⟶}{E}_{\gamma }\left\{B^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}=-v{\left\{g(t;\gamma )\right\}}^{{\displaystyle \frac{v+1}{v}}}(t_1-1,t_2-1,(t_1-1)(t_2-1)).$$

For the second derivatives, it follows from (9) that

$$B^{\left(2\right)}(x;t;\tau )=-v(v+1){\left\{g(t;\tau )\right\}}^{{\displaystyle \frac{v+2}{v}}}{A}_t,$$

where $A_t=\left(\begin{array}{ccc}{(t_1-1)}^2& (t_1-1)(t_2-1)& {(t_1-1)}^2(t_2-1)\\ (t_1-1)(t_2-1)& {(t_2-1)}^2& (t_1-1){(t_2-1)}^2\\ {(t_1-1)}^2(t_2-1)& (t_1-1){(t_2-1)}^2& {\left\{(t_1-1)(t_2-1)\right\}}^2\end{array}\right)$.

Therefore, $\left|B_{jk}^{\left(2\right)}\left({\mathit{X}}_1;t;\tau \right)\right|\le v(v+1)<\infty$ for $j,k\in \{1,2,3\}$$\forall t\in {[0,1]}^2$.

Thus, considering (4), it results that

$$E_{\gamma }\left\{{∥B_{jk}^{\left(2\right)}\left({\mathit{X}}_1;t;\gamma \right)∥}_{\mathcal{H}}^2\right\}<\infty ,j,k\in \{1,2,3\}.$$

Following similar steps as those for the vector $B^{\left(1\right)}({\mathit{X}}_i;t;\gamma )$ for the matrix of second derivatives, we obtain that in the space $\mathcal{H}$ it results in

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}^{\left(2\right)}({\mathit{X}}_i;t;\gamma )\stackrel{P}{⟶}{E}_{\gamma }\left\{B^{\left(2\right)}({\mathit{X}}_1;t;\gamma )\right\}=-v(v+1){\left\{g(t;\gamma )\right\}}^{{\displaystyle \frac{v+2}{v}}}{A}_t.$$

Additionally, using Assumption 1, we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}\sqrt{n}({\widehat{\gamma }}_n-\gamma ){\displaystyle \frac{1}{n}}\sum _{i=1}^n{B}^{\left(2\right)}({\mathit{X}}_i;t;\gamma )& \sqrt{n}{({\widehat{\gamma }}_n-\gamma )}^{\top }=\hfill \\ \hfill & {\displaystyle \frac{1}{2}}{E}_{\gamma }\left\{\ell ({\mathit{X}}_1;\gamma )\right\}{E}_{\gamma }\left\{B^{\left(2\right)}({\mathit{X}}_1;t;\gamma )\right\}{\displaystyle \frac{1}{\sqrt{n}}}\sum _{k=1}^n\ell {\left({\mathit{X}}_k;\gamma \right)}^{\top }\hfill \\ \hfill & +E_{\gamma }\left\{\ell ({\mathit{X}}_1;\gamma )\right\}{E}_{\gamma }\left\{B^{\left(2\right)}({\mathit{X}}_1;t;\gamma )\right\}{\mathbf{o}}_P{\left(1\right)}^{\top }+o_P\left(1\right)\hfill \end{array}$$

On the other hand, using similar arguments to the previous ones, we obtain $\parallel r_n{\parallel }_{\mathcal{H}}=o_P\left(1\right)$; then, using Assumption 1, (8) can be written as

$$B_n(t;{\widehat{\gamma }}_n)=B_n(t;\gamma )+b_n,$$

where $\parallel b_n{\parallel }_{\mathcal{H}}=o_P\left(1\right)$, and 

$$B_n(t;\gamma )={\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^n\left[B({\mathit{X}}_i;t;\gamma )+E_{\gamma }\left\{B^{\left(1\right)}({\mathit{X}}_1;t;\gamma )\right\}\ell {\left({\mathit{X}}_i;\gamma \right)}^{\top }\right].$$

On the other hand, observe that

$$\parallel B_n{\left(\gamma \right)\parallel }_{\mathcal{H}}^2={\displaystyle \frac{1}{n}}\sum _{i=1}^n\sum _{j=1}^{n}h({\mathit{X}}_i,{\mathit{X}}_j;\gamma ),$$

where $h(\mathit{x},\mathit{y};\gamma )$ is defined in (6) and satisfies $h(\mathit{x},\mathit{y};\gamma )=h(\mathit{y},\mathit{x};\gamma )$, $E_{\gamma }\left\{h^2({\mathit{X}}_1,{\mathit{X}}_2;\gamma )\right\}<\infty$, $E_{\gamma }\left\{\left|h\right({\mathit{X}}_1,{\mathit{X}}_1;\gamma \left)\right|\right\}<\infty$ and $E_{\gamma }\left\{h({\mathit{X}}_1,{\mathit{X}}_2;\gamma )\right\}=0$. Thus, from Theorem 6.4.1.B in Serfling [24],

$$\parallel B_n{\left(\gamma \right)\parallel }_{\mathcal{H}}^2\stackrel{L}{⟶}\sum _{j\ge 1}{\lambda }_j{\chi }_{1j}^2$$

where ${\chi }_{11}^2,{\chi }_{12}^2,\dots$ and the set $\left\{{\lambda }_j\right\}$ is as defined in the statement of the Theorem. In particular, $\parallel B_n{\left(\gamma \right)\parallel }_{\mathcal{H}}^2=O_P\left(1\right)$, which implies (5).   □

As seen in Theorem 1, the null asymptotic distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$ still depends on the true value of the parameter $\gamma$, which, being unknown, does not provide a useful solution. This issue is resolved by replacing $\gamma$ with $\widehat{\gamma }$.

However, a greater challenge lies in determining the set ${\left\{{\lambda }_j\right\}}_{j\ge 1}$, since calculating the eigenvalues of an operator is not straightforward. Moreover, it is also necessary to obtain the expression for $h(\mathit{x},\mathit{y};\gamma )$, which is difficult to derive, because it depends on the function ℓ, which generally does not have a simple expression.

Thus, in the following section, we employ the parametric bootstrap method to approximate the null distribution of the test statistic.

4. Bootstrap Approximation of the Null Distribution

The parametric bootstrap method provides an alternative way to estimate the null distribution of the proposed statistic. For this, let ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ be iid random vectors taking values in ${\mathbb{N}}_0^2$, and suppose that ${\widehat{\gamma }}_n={\widehat{\gamma }}_n({\mathit{X}}_1,\dots ,{\mathit{X}}_n)\in \mathrm{\Theta }$.

Furthermore, to use the bootstrap method, we must consider that ${\mathit{X}}_1^{\ast },\dots ,{\mathit{X}}_n^{\ast }$ are iid random vectors from a population distributed according to the $BNB\left({\widehat{\gamma }}_n\right)$ law given that ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$. Also, let $B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)$ be the bootstrap version of $B_{n,w}\left({\widehat{\gamma }}_n\right)$, which is obtained by replacing ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ and ${\widehat{\gamma }}_n={\widehat{\gamma }}_n({\mathit{X}}_1,\dots ,{\mathit{X}}_n)$ with ${\mathit{X}}_1^{\ast },\dots ,{\mathit{X}}_n^{\ast }$ and ${\widehat{\gamma }}_n^{\ast }={\widehat{\gamma }}_n({\mathit{X}}_1^{\ast },\dots ,{\mathit{X}}_n^{\ast })$, respectively, in the expression for $B_{n,w}\left({\widehat{\gamma }}_n\right)$.

Theorem 2, presented later, demonstrates that the parametric bootstrap method consistently approximates the null distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$. To achieve this, we must assume the following hypotheses, which are satisfied by commonly used estimators, as mentioned after Assumption 1.

Assumption 2.

Assumption 1 holds and the functions ℓ and J satisfy 

(1) 

${sup}_{\varrho \in {\mathrm{\Theta }}_0}{E}_{\varrho }\left[{\parallel \ell (\mathit{X};\varrho )\parallel }^{2}I\left\{\parallel \ell (\mathit{X};\varrho )\parallel >\tau \right\}\right]⟶0$ as $\tau \to \infty$, where ${\mathrm{\Theta }}_0\subseteq \mathrm{\Theta }$ is an open neighborhood of γ. 

(2) 

$\ell (\mathit{X};\varrho )$ and $J\left(\varrho \right)$ are continuous as functions of ϱ at $\varrho =\gamma$, and $J\left(\varrho \right)$ is finite $\forall \varrho \in {\mathrm{\Theta }}_0$. 

Theorem 2.

Let ${\mathit{X}}_1,\dots ,{\mathit{X}}_n$ be iid from a random vector $\mathit{X}=(X_1,X_2)\in {\mathbb{N}}_0^2$. Suppose that Assumption 2 holds and that ${\widehat{\gamma }}_n=\gamma +o\left(1\right)$ for some $\gamma \in \mathrm{\Theta }$. Then,

$$\underset{x\in \mathbb{R}}{sup}\left|P_{\ast }\left\{B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)\le x\right\}-P_{\gamma }\left\{B_{n,w}\left({\widehat{\gamma }}_n\right)\le x\right\}\right|\stackrel{a.s.}{⟶}0.$$

Proof. 

By definition, $B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)={\parallel B_n^{\ast }(t;{\widehat{\gamma }}_n^{\ast })\parallel }_{\mathcal{H}}^2$, with 

$$B_n^{\ast }(t;{\widehat{\gamma }}_n^{\ast })={\displaystyle \frac{1}{\sqrt{n}}}\sum _{i=1}^{n}B({\mathit{X}}_i^{\ast };t;{\widehat{\gamma }}_n^{\ast })$$

and $B(\mathit{X};t;\gamma )$ defined in (7).

Following similar steps to those given in the proof of Theorem 1, it can be seen that $B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)={\parallel V_n^{\ast }\parallel }_{\mathcal{H}}^2+o_{P_{\ast }}\left(1\right)$, where $V_n^{\ast }\left(t\right)$ is defined as $V_n\left(t\right)$, with ${\mathit{X}}_i$ and $\gamma$ replaced by ${\mathit{X}}_i^{\ast }$ and ${\widehat{\gamma }}_n$, respectively.

Now, the result follows by applying Theorem 1.1 in Kundu et al. [25] to $V_n^{\ast }$ and the continuous mapping theorem.    □

From Theorem 2, the test function

$${\mathrm{\Psi }}^{\ast }=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}{B}_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)\ge b_{n,w,\alpha }^{\ast },\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

or, equivalently, the test that rejects $H_0$ when $p^{\ast }=P_{\ast }\{B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)\ge B_{obs}\}\le \alpha ,$ is asymptotically correct in the sense that, when $H_0$ is true, $lim{P}_{\gamma }({\mathrm{\Psi }}^{\ast }=1)=\alpha$, where $b_{n,w,\alpha }^{\ast }=inf\{x:P_{\ast }(B_{n,w}^{\ast }\left({\widehat{\gamma }}_n^{\ast }\right)\ge x)\le \alpha \}$ is the $\alpha$ upper percentile of the bootstrap distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$, and $B_{obs}$ is the observed value of the test statistic.

5. Parameter Estimation

In this research, three techniques were analyzed to estimate the parameters established in parameterization (2), namely, the (a) method of moments, (b) zero–zero cell frequency, and (c) maximum likelihood. In what follows, we will use different notation to denote the estimator of $\gamma$ depending on the method employed: $\tilde{\gamma }$ for the method of moments, $\tilde{\tilde{\gamma }}$ for the zero–zero cell frequency method, and $\widetilde{\tilde{\tilde{\gamma }}}$ for the maximum likelihood. For this, let us consider $(X_1,X_2)\sim BNB\left(\gamma \right)$.

5.1. Method of Moments

The usual procedure for obtaining the moment estimators is to equate the sample moments to the population moments and solve the resulting equations. According to [26], we shall use the two first marginal moments and the first mixed central moment $m_{1,1}$.

The moments of the distribution in (2) are

$$\begin{array}{cc}\hfill {\mu }_{1,0}^{\prime }& =v{\gamma }_0,{\mu }_{0,1}^{\prime }=v{\gamma }_1,{\mu }_{1,1}=v({\gamma }_2+{\gamma }_0{\gamma }_1),{\mu }_{2,0}=v{\gamma }_0(1+{\gamma }_0),\hfill \\ \hfill {\mu }_{0,2}& =v{\gamma }_1(1+{\gamma }_1),{\mu }_{2,1}=v(2{\gamma }_0+1)({\gamma }_2+{\gamma }_0{\gamma }_1),{\mu }_{1,2}=v(2{\gamma }_0+1)({\gamma }_1+{\gamma }_0{\gamma }_1),\hfill \\ \hfill {\mu }_{2,2}& =v\left\{2\left(v+1\right){\gamma }_2^2+2\left(v+2\right){\gamma }_0{\gamma }_1{\gamma }_2+3\left(v+2\right){\gamma }_0^2{\gamma }_1^2+{\gamma }_2\left(2{\gamma }_0+2{\gamma }_1+1\right)\right.\hfill \\ \hfill & +\left.\left(v+1\right){\gamma }_0{\gamma }_1\left(1+{\gamma }_0+{\gamma }_1\right)+{\gamma }_0{\gamma }_1\left({\gamma }_0+{\gamma }_1\right)\right\}.\hfill \end{array}$$

5.1.1. Estimators of $\tilde{\gamma }$

Using the sample marginal first moments ${\overline{X}}_1$ and ${\overline{X}}_2$, as well as the mixed moment $m_{1,1}$, we have the moment estimators

$${\tilde{\gamma }}_0={\displaystyle \frac{{\overline{X}}_1}{v}},{\tilde{\gamma }}_1={\displaystyle \frac{{\overline{X}}_2}{v}},{\tilde{\gamma }}_2={\displaystyle \frac{m_{1,1}}{v}}-{\displaystyle \frac{{\overline{X}}_1{\overline{X}}_2}{v^2}}.$$

5.1.2. Asymptotic Variance Matrix of $\tilde{\gamma }$

The asymptotic variance–covariance matrix of the estimators obtained through the method of moments, ${\mathrm{\Sigma }}_{MM}$, is given by

$${\mathrm{\Sigma }}_{MM}={\displaystyle \frac{1}{nv}}\left(\begin{array}{ccc}{\gamma }_0\left(1+{\gamma }_0\right)& {\gamma }_2+{\gamma }_0{\gamma }_1& \left(2{\gamma }_0+1\right)\left({\gamma }_2+{\gamma }_0{\gamma }_1\right)\\ {\gamma }_2+{\gamma }_0{\gamma }_1& {\gamma }_1\left(1+{\gamma }_1\right)& \left(2{\gamma }_0+1\right)\left({\gamma }_1+{\gamma }_0{\gamma }_1\right)\\ \left(2{\gamma }_0+1\right)\left({\gamma }_2+{\gamma }_0{\gamma }_1\right)& \left(2{\gamma }_0+1\right)\left({\gamma }_1+{\gamma }_0{\gamma }_1\right)& \varsigma \end{array}\right),$$

where $\varsigma =(v+2){\gamma }_2^2+{\gamma }_0{\gamma }_1\left\{4{\gamma }_2+2(v+3){\gamma }_0{\gamma }_1+(v+2)({\gamma }_0+{\gamma }_1)+v+1\right\}$.

5.1.3. Asymptotic Distribution of $\tilde{\gamma }$

Since $\tilde{\gamma }$ is a function of sample moments, applying Theorem 11.2.14 (Delta Method) given in Lehmann and Romano [27], it follows that

$$\sqrt{n}\left(\tilde{\gamma }-\gamma \right)\stackrel{L}{⟶}{N}_3(\mathbf{0},{\mathrm{\Sigma }}_{MM}).$$

5.1.4. Asymptotic Properties of $\tilde{\gamma }$

According to the results of Section 5.1.3, $\tilde{\gamma }$ is a consistent and asymptotically normal estimator.

5.2. Zero–Zero Cell Frequency

5.2.1. Estimators of $\tilde{\tilde{\gamma }}$

From (2),

$$P_{\gamma }(0,0)={(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}^{-v}.$$

A combination of ${\overline{X}}_1$, ${\overline{X}}_2$ and the observed zero–zero cell relative frequency ${\displaystyle \frac{n_{00}}{n}}$ can be used to estimate the parameters ${\gamma }_0$, ${\gamma }_1$ and ${\gamma }_2$:

$${\tilde{\tilde{\gamma }}}_0={\displaystyle \frac{{\overline{X}}_1}{v}},{\tilde{\tilde{\gamma }}}_1={\displaystyle \frac{{\overline{X}}_2}{v}},{\tilde{\tilde{\gamma }}}_2=1+{\displaystyle \frac{{\overline{X}}_1}{v}}+{\displaystyle \frac{{\overline{X}}_2}{v}}-{\left\{{\displaystyle \frac{n_{00}}{n}}\right\}}^{-{\displaystyle \frac{1}{v}}}.$$

5.2.2. Asymptotic Variance Matrix of $\tilde{\tilde{\gamma }}$

The asymptotic variance–covariance matrix of the estimators obtained through the zero–zero cell frequency method, ${\mathrm{\Sigma }}_{ZZ}$, is given by

$${\mathrm{\Sigma }}_{ZZ}={\displaystyle \frac{1}{vn}}\left(\begin{array}{ccc}{\gamma }_0(1+{\gamma }_0)& {\gamma }_2+{\gamma }_0{\gamma }_1& {\gamma }_2(1+{\gamma }_0)\\ {\gamma }_2+{\gamma }_0{\gamma }_1& {\gamma }_1(1+{\gamma }_1)& {\gamma }_2(1+{\gamma }_1)\\ {\gamma }_2(1+{\gamma }_0)& {\gamma }_2(1+{\gamma }_1)& \left(1+{\gamma }_0+{\gamma }_1\right)\left(2{\gamma }_2-{\gamma }_0-{\gamma }_1\right)+\zeta \end{array}\right),$$

where $\zeta ={\displaystyle \frac{1}{v}}\left({\displaystyle \frac{1}{P_{\gamma }(0,0)}}-1\right){(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}^2$.

5.2.3. Asymptotic Distribution of $\tilde{\tilde{\gamma }}$

Since $\tilde{\tilde{\gamma }}$ is a function of sample moments, applying Theorem 11.2.14 (Delta Method) given in Lehmann and Romano [27], it follows that

$$\sqrt{n}\left(\tilde{\tilde{\gamma }}-\gamma \right)\stackrel{L}{⟶}{N}_3(\mathbf{0},{\mathrm{\Sigma }}_{ZZ}).$$

5.2.4. Asymptotic Properties of $\tilde{\tilde{\gamma }}$

According to the results of Section 5.2.3 $\tilde{\tilde{\gamma }}$ is a consistent and asymptotically normal estimator.

5.3. Maximum Likelihood Estimators

5.3.1. Estimators of $\widetilde{\tilde{\tilde{\gamma }}}$

Differentiating the pmf given in (2) with respect to ${\gamma }_0$, ${\gamma }_1$, and ${\gamma }_2$ successively, we have

$${\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_0}}={\displaystyle \frac{r}{{\gamma }_0-{\gamma }_2}}-{\displaystyle \frac{v+r+s}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}}-{\displaystyle \frac{(1+{\gamma }_1)R(r,s)}{({\gamma }_0-{\gamma }_2)(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}},$$

(11)

where $R(r,s)={\displaystyle \frac{S^{\prime }(r,s)}{S(r,s)}}$, and

$$S^{\prime }(r,s)=\sum _{i=0}^{min(r,s)}{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{r}{i}\right)\left(\genfrac{}{}{0pt}{}{s}{i}\right)}{\left(\genfrac{}{}{0pt}{}{v+r+s-1}{i}\right)}}i{\tau }^i.$$

A similar procedure yields

$${\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_1}}={\displaystyle \frac{s}{{\gamma }_1-{\gamma }_2}}-{\displaystyle \frac{v+r+s}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}}-{\displaystyle \frac{(1+{\gamma }_0)R(r,s)}{({\gamma }_1-{\gamma }_2)(1+{\gamma }_0+{\gamma }_1-{\gamma }_2)}},$$

(12)

and

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_2}}& =-{\displaystyle \frac{r}{{\gamma }_0-{\gamma }_2}}-{\displaystyle \frac{s}{{\gamma }_1-{\gamma }_2}}-{\displaystyle \frac{v+r+s}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}}\hfill \\ & +\left({\displaystyle \frac{1}{{\gamma }_2}}+{\displaystyle \frac{1}{{\gamma }_0-{\gamma }_2}}+{\displaystyle \frac{1}{{\gamma }_1-{\gamma }_2}}-{\displaystyle \frac{1}{1+{\gamma }_0+{\gamma }_1-{\gamma }_2}}\right)R(r,s).\hfill \end{array}$$

(13)

The estimating equations are

$$\sum _{r,s}{\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_0}}=0,\sum _{r,s}{\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_1}}=0,\sum _{r,s}{\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_1}}=0,$$

which yield the maximum likelihood estimators $\widetilde{\tilde{\tilde{\gamma }}}$ as

$${\displaystyle \frac{{\overline{X}}_1}{{\widetilde{\tilde{\tilde{\gamma }}}}_0-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}-{\displaystyle \frac{v+{\overline{X}}_1+{\overline{X}}_2}{\xi }}-{\displaystyle \frac{\left(1+{\widetilde{\tilde{\tilde{\gamma }}}}_1\right)\overline{R}}{\xi \left({\widetilde{\tilde{\tilde{\gamma }}}}_0-{\widetilde{\tilde{\tilde{\gamma }}}}_2\right)}}=0,$$

$$-{\displaystyle \frac{v+{\overline{X}}_1+{\overline{X}}_2}{\xi }}+{\displaystyle \frac{{\overline{X}}_2}{{\widetilde{\tilde{\tilde{\gamma }}}}_1-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}-{\displaystyle \frac{\left(1+{\widetilde{\tilde{\tilde{\gamma }}}}_0\right)\overline{R}}{\xi \left({\widetilde{\tilde{\tilde{\gamma }}}}_1-{\widetilde{\tilde{\tilde{\gamma }}}}_2\right)}}=0,$$

$$-{\displaystyle \frac{{\overline{X}}_1}{{\widetilde{\tilde{\tilde{\gamma }}}}_0-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}-{\displaystyle \frac{{\overline{X}}_2}{{\widetilde{\tilde{\tilde{\gamma }}}}_1-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}+{\displaystyle \frac{v+{\overline{X}}_1+{\overline{X}}_2}{\xi }}+\left\{{\displaystyle \frac{1}{{\widetilde{\tilde{\tilde{\gamma }}}}_2}}+{\displaystyle \frac{1}{{\widetilde{\tilde{\tilde{\gamma }}}}_0-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}+{\displaystyle \frac{1}{{\widetilde{\tilde{\tilde{\gamma }}}}_1-{\widetilde{\tilde{\tilde{\gamma }}}}_2}}-{\displaystyle \frac{1}{\xi }}\right\}\overline{R}=0,$$

where   $\xi =1+{\widetilde{\tilde{\tilde{\gamma }}}}_0+{\widetilde{\tilde{\tilde{\gamma }}}}_1-{\widetilde{\tilde{\tilde{\gamma }}}}_2$.

Here, $\overline{R}$ is the sample mean of $R(x,y)$. Solving these equations leads to

$${\widetilde{\tilde{\tilde{\gamma }}}}_0={\displaystyle \frac{{\overline{X}}_1}{v}},{\widetilde{\tilde{\tilde{\gamma }}}}_1={\displaystyle \frac{{\overline{X}}_2}{v}},$$

while ${\widetilde{\tilde{\tilde{\gamma }}}}_2$ is the solution to the equation

$$v{\widetilde{\tilde{\tilde{\gamma }}}}_2=\overline{R}.$$

This equation for ${\widetilde{\tilde{\tilde{\gamma }}}}_2$ can be solved iteratively using either the moment or zero–zero cell frequency estimate as an initial value.

5.3.2. Asymptotic Variance Matrix of $\widetilde{\tilde{\tilde{\gamma }}}$

By definition, The third-order information matrix, $F_3$, is given by (for large n)

$$F_3=n\left(\sum _{r,s}{\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_i}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_j}}\right),\mathrm{with}i,j\in \{0,1,2\}.$$

Returning to Equations (11)–(13) and recalling that the asymptotic variance–covariance matrix of the estimators obtained by the maximum likelihood method, ${\mathrm{\Sigma }}_{ML}$, is the inverse of the information matrix, then

$${\mathrm{\Sigma }}_{ML}={\displaystyle \frac{1}{n}}{\left(\sum _{r,s}{Q}_i(r,s)Q_j(r,s)P_{\gamma }(r,s)\right)}^{-1},i,j\in \{0,1,2\},$$

where

$$Q_i(r,s)={\displaystyle \frac{1}{P_{\gamma }(r,s)}}{\displaystyle \frac{\partial P_{\gamma }(r,s)}{\partial {\gamma }_i}},i=0,1,2.$$

5.3.3. Asymptotic Distribution of $\widetilde{\tilde{\tilde{\gamma }}}$

Since $\widetilde{\tilde{\tilde{\gamma }}}$ is a function of sample moments, applying Theorem 11.2.14 (Delta Method) given in Lehmann and Romano [27], it follows that

$$\sqrt{n}\left(\widetilde{\tilde{\tilde{\gamma }}}-\gamma \right)\stackrel{L}{⟶}{N}_3(\mathbf{0},{\mathrm{\Sigma }}_{ML}).$$

5.3.4. Asymptotic Properties of $\widetilde{\tilde{\tilde{\gamma }}}$

According to the results of Section 5.3.3, $\widetilde{\tilde{\tilde{\gamma }}}$ is a consistent and asymptotically normal estimator.

6. Numerical Results and Discussion

The properties of the statistic $B_{n,w}\left({\widehat{\gamma }}_n\right)$ are asymptotic (see, for example, [17]). However, we need to apply the test statistic in the real world, i.e., for finite sample sizes. This section is dedicated to investigating the behavior of the bootstrap approximation of the proposed statistic in a finite scenario. To this end, we conducted a simulation study and present the results obtained.

It is important to note that, as we did not find another consistent test statistic for the BNBD, the numerical study presented has no point of comparison. This is the reason why the study was carried out using three estimation techniques.

Given the massive volume of calculations performed, we minimized computation time by using the parallelization algorithm called parRapply, which has been suggested in [18]. All computational work was done through codes written in the R programming language [28].

It is important to point out that parRapply is a function from the parallel package, which is available in [18]. parRapply is a computational technique based on task simultaneity, allowing for faster computation times by partitioning the work to be executed across the different cores available on a computer. Its main goal is to reduce the execution time of a program. For more detailed information about parRapply, you can refer to the documentation available in [18].

To calculate $B_{n,w}\left({\widehat{\gamma }}_n\right)$, we used the weight function

$$w(t;a_1,a_2)=t_1^{a_1}{t}_2^{a_2},$$

(14)

with $a_1>-1$ and $a_2>-1$, resulting in the following explicit form of the test statistic:

$$\begin{array}{cc}\hfill & B_{n,w}\left({\widehat{\gamma }}_n\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n\sum _{j=1}^n\left({\displaystyle \frac{1}{X_{i1}+X_{j1}+a_1+1}}\right)\left({\displaystyle \frac{1}{X_{i2}+X_{j2}+a_2+1}}\right)\hfill \\ \hfill & +q^v\sum _{k=0}^{\infty }\sum _{l=0}^k\sum _{m=0}^{k-l}{\displaystyle \frac{{({\widehat{\gamma }}_0-{\widehat{\gamma }}_1)}^{k-l-m}{({\widehat{\gamma }}_1-{\widehat{\gamma }}_2)}^l{\widehat{\gamma }}_2^m}{l!m!(k-m-l)!q^k}}\left({\displaystyle \frac{n{q}^v(2v+k-1)!}{(2v-1)!A_{kl}{C}_{lm}}}-{\displaystyle \frac{2(v+k-1)!}{(v-1)!}}\sum _{i=1}^n{\displaystyle \frac{1}{D_{iklm}}}\right),\hfill \end{array}$$

where   $q=1+{\widehat{\gamma }}_0+{\widehat{\gamma }}_1-{\widehat{\gamma }}_2$,   $A_{kl}=a_1+k-l+1$,   $C_{lm}=a_2+l+m+1$,  and  $D_{iklm}=(X_{i1}+A_{kl})(X_{i2}+C_{lm})$.

Remark 1.

For practical purposes, truncating the infinite series to the first 20 terms yielded sufficiently accurate values for the test statistic and good performance of the corresponding subroutine. 

6.1. Simulated Data for Type I Error

In practice, the exact bootstrap estimator of the null distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$ cannot be calculated. As usual, we approximate it by simulation as follows:

1.

Estimate $\gamma$ through $\widehat{\gamma }$ and compute the observed value of the test statistic $B_{obs}$.

2.

For some large integer K, repeat for every $k\in \{1,\dots ,K\}$:

(a)

Generate ${\mathit{X}}^{\ast k}=({\mathit{X}}_1^{\ast k},\dots ,{\mathit{X}}_n^{\ast k})$, where ${\mathit{X}}_1^{\ast k},\dots ,{\mathit{X}}_n^{\ast k}$ are iid from a $BNB\left(\widehat{\gamma }\right)$

(b)

Calculate the test statistic evaluated at ${\mathit{X}}^{\ast k}$, obtaining $B_n^{\ast k}\left({\widehat{\gamma }}^{\ast k}\right)$.

3.

Approximate the p-value by $\widehat{p}={\displaystyle \frac{1}{K}}{\sum }_{k=1}^{K}I\{B_n^{\ast k}\left({\widehat{\gamma }}^{\ast k}\right)>B_{obs}\}$.

It is important to note that in Appendix A, we present the computational implementation of the parametric bootstrap algorithm written in R code.

To study the goodness of the bootstrap approximation of $B_{n,w}\left({\widehat{\gamma }}_n\right)$ for finite sample sizes, we generated samples of size $n=30\left(20\right)70$ from $BNB({\gamma }_0,{\gamma }_1,{\gamma }_2)$, with ${\gamma }_0={\gamma }_1$, where ${\gamma }_0,{\gamma }_1\in \{0.29,0.30,0.31,0.32\}$ and the parameter ${\gamma }_2$ were chosen such that the correlation coefficient, $\rho ={\displaystyle \frac{{\gamma }_0{\gamma }_1+{\gamma }_2}{\sqrt{{\gamma }_0{\gamma }_1(1+{\gamma }_0)(1+{\gamma }_1)}}}$, equaled to 0.25, 0.50, and 0.75 in order to evaluate the quality of the approximations for data with low, moderate, and high correlations, respectively.

To estimate the parameter $\gamma$, we employed the three techniques described in Section 5. Then, we approximated the bootstrap p-values, $p^{\ast }$, for the proposed test. For this purpose, we used the weighting function given in (14) for $(a_1,a_2)\in \{(0,0),(1,0),(0,1)\}$ and generated $K=500$ bootstrap samples.

The above procedure was repeated 1000 times, and we calculated the fraction of the estimated p-values that were less than or equal to 0.05 and 0.10. These fractions correspond to the estimated probabilities of Type I error for $\alpha =0.05,0.10$, respectively.

Next, we repeat the previous experiment for ${\gamma }_0\ne {\gamma }_1$, where ${\gamma }_0,{\gamma }_1\in \{0.29,0.30,0.31,0.32\}$ and ${\gamma }_2$ were chosen such that the correlation coefficient, $\rho ={\displaystyle \frac{{\gamma }_0{\gamma }_1+{\gamma }_2}{\sqrt{{\gamma }_0{\gamma }_1(1+{\gamma }_0)(1+{\gamma }_1)}}}$, was approximately equal to 0.25, 0.50, and 0.75.

In both scenarios, ${\gamma }_0={\gamma }_1$ and ${\gamma }_0\ne {\gamma }_1$, we experimented with different numbers of trials before the first event occurred, v. We obtained similar results for the different simulated cases; here, we only report the results for $v=5$.

Table 1. Simulation results for the probability of type I error for $\rho =0.25$  and $E\left(X_1\right)=E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.29,0.009425)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.033	0.080	0.036	0.087	0.061	0.110
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.034	0.081	0.035	0.088	0.060	0.109
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.032	0.080	0.036	0.087	0.058	0.111
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.083	0.042	0.089	0.056	0.107
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.083	0.043	0.090	0.057	0.108
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.040	0.082	0.044	0.090	0.057	0.109
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.088	0.045	0.095	0.051	0.102
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.090	0.047	0.097	0.050	0.101
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.087	0.046	0.097	0.051	0.102
$\gamma =(0.30,0.30,0.007500)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.030	0.084	0.035	0.085	0.061	0.112
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.034	0.083	0.035	0.087	0.062	0.113
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.080	0.038	0.089	0.059	0.112
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.086	0.042	0.091	0.056	0.107
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.085	0.041	0.093	0.057	0.106
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.085	0.044	0.093	0.055	0.105
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.090	0.047	0.095	0.052	0.102
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.087	0.048	0.094	0.051	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.092	0.049	0.095	0.050	0.101
$\gamma =(0.32,0.32,0.003200)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.030	0.083	0.037	0.086	0.060	0.108
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.031	0.082	0.035	0.089	0.061	0.109
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.032	0.081	0.036	0.087	0.060	0.107
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.087	0.041	0.091	0.058	0.107
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.087	0.043	0.093	0.057	0.105
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.036	0.086	0.043	0.093	0.055	0.104
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.089	0.045	0.096	0.053	0.102
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.091	0.046	0.095	0.053	0.103
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.093	0.045	0.096	0.052	0.101

,

Table 2. Simulation results for the probability of type I error for $\rho \approx$ 0.25  and $E\left(X_1\right)\ne E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.30,0.008492)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.031	0.079	0.039	0.085	0.058	0.113
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.034	0.081	0.039	0.087	0.056	0.111
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.075	0.035	0.084	0.057	0.114
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.086	0.042	0.089	0.056	0.110
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.085	0.041	0.089	0.056	0.109
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.085	0.042	0.088	0.055	0.110
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.092	0.046	0.096	0.053	0.101
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.093	0.045	0.095	0.052	0.103
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.090	0.044	0.096	0.052	0.100
$\gamma =(0.29,0.31,0.007543)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.032	0.079	0.035	0.083	0.059	0.115
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.033	0.083	0.035	0.085	0.058	0.112
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.031	0.077	0.033	0.081	0.059	0.116
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.081	0.040	0.087	0.057	0.110
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.036	0.087	0.043	0.089	0.055	0.109
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.085	0.044	0.090	0.056	0.108
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.093	0.044	0.095	0.053	0.103
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.090	0.045	0.094	0.053	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.089	0.046	0.092	0.050	0.100
$\gamma =(0.30,0.31,0.006492)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.031	0.081	0.037	0.084	0.060	0.112
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.033	0.078	0.038	0.084	0.059	0.112
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.076	0.037	0.086	0.060	0.110
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.035	0.081	0.041	0.086	0.055	0.109
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.083	0.043	0.087	0.055	0.107
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.035	0.082	0.041	0.084	0.056	0.111
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.086	0.044	0.092	0.054	0.102
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.085	0.045	0.095	0.051	0.101
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.085	0.045	0.095	0.052	0.101

,

Table 3. Simulation results for the probability of type I error for $\rho =0.50$  and $E\left(X_1\right)=E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.29,0.102950)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.033	0.075	0.035	0.083	0.060	0.115
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.031	0.083	0.037	0.084	0.059	0.112
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.082	0.038	0.086	0.058	0.113
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.082	0.043	0.084	0.055	0.112
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.083	0.044	0.087	0.055	0.111
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.085	0.042	0.089	0.055	0.111
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.088	0.045	0.095	0.052	0.104
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.090	0.047	0.096	0.051	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.093	0.045	0.097	0.052	0.100
$\gamma =(0.30,0.30,0.105000)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.034	0.074	0.038	0.084	0.059	0.114
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.033	0.075	0.037	0.085	0.060	0.112
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.034	0.076	0.041	0.084	0.058	0.114
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.083	0.044	0.086	0.054	0.109
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.085	0.042	0.087	0.056	0.111
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.082	0.042	0.086	0.055	0.110
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.089	0.046	0.095	0.053	0.103
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.091	0.044	0.095	0.053	0.103
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.093	0.045	0.096	0.050	0.102
$\gamma =(0.32,0.32,0.108800)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.032	0.080	0.039	0.084	0.060	0.113
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.035	0.076	0.036	0.083	0.061	0.114
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.030	0.078	0.038	0.084	0.060	0.113
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.081	0.042	0.087	0.057	0.109
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.084	0.042	0.086	0.056	0.109
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.082	0.041	0.085	0.057	0.108
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.087	0.044	0.095	0.051	0.100
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.091	0.045	0.094	0.053	0.104
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.093	0.046	0.095	0.051	0.103

,

Table 4. Simulation results for the probability of type I error for $\rho \approx$ 0.50  and $E\left(X_1\right)\ne E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.31,0.104990)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.032	0.078	0.038	0.083	0.060	0.114
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.032	0.076	0.037	0.084	0.060	0.113
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.034	0.081	0.037	0.083	0.059	0.113
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.080	0.041	0.084	0.057	0.111
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.081	0.043	0.085	0.055	0.110
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.082	0.042	0.087	0.056	0.109
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.086	0.046	0.096	0.053	0.102
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.090	0.044	0.095	0.053	0103
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.091	0.047	0.094	0.051	0.103
$\gamma =(0.30,0.31,0.105984)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.031	0.072	0.037	0.083	0.059	0.115
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.033	0.079	0.035	0.083	0.060	0.115
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.074	0.038	0.085	0.058	0.114
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.080	0.042	0.084	0.056	0.109
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.081	0.043	0.085	0.055	0.111
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.082	0.041	0.087	0.056	0.109
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.087	0.045	0.094	0.052	0.103
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.090	0.044	0.093	0.053	0.103
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.040	0.087	0.047	0.095	0.051	0.101
$\gamma =(0.30,0.32,0.106939)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.031	0.079	0.034	0.081	0.061	0.115
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.032	0.076	0.035	0.083	0.060	0.112
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.032	0.077	0.035	0.082	0.061	0.114
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.082	0.042	0.087	0.054	0.111
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.082	0.041	0.088	0.053	0.106
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.081	0.043	0.088	0.054	0.107
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.088	0.046	0.095	0.051	0.104
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.086	0.044	0.096	0.050	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.088	0.044	0.095	0.053	0.103

,

Table 5. Simulation results for the probability of type I error for $\rho =0.75$  and $E\left(X_1\right)=E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=0.05$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.29,0.196475)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.034	0.085	0.038	0.089	0.059	0.110
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.033	0.084	0.036	0.088	0.060	0.111
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.032	0.084	0.037	0.087	0.058	0.111
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.085	0.042	0.089	0.057	0.107
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.086	0.045	0.090	0.054	0.104
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.084	0.043	0.089	0.056	0.109
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.044	0.088	0.047	0.095	0.051	0.104
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.090	0.044	0.096	0.052	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.093	0.045	0.097	0.052	0.100
$\gamma =(0.30,0.30,0.202500)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.035	0.076	0.039	0.086	0.058	0.110
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.034	0.075	0.039	0.085	0.059	0.111
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.034	0.075	0.040	0.084	0.057	0.112
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.036	0.083	0.043	0.087	0.054	0.108
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.085	0.044	0.088	0.055	0.109
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.086	0.044	0.089	0.055	0.107
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.090	0.046	0.095	0.052	0.101
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.091	0.046	0.095	0.052	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.092	0.045	0.096	0.051	0.101
$\gamma =(0.31,0.31,0.208475)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.033	0.078	0.039	0.085	0.059	0.112
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.035	0.079	0.039	0.086	0.060	0.113
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.033	0.078	0.038	0.084	0.050	0.112
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.084	0.043	0.088	0.054	0.107
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.040	0.085	0.045	0.089	0.053	0.108
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.040	0.086	0.046	0.090	0.052	0.106
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.089	0.046	0.095	0.051	0.100
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.091	0.046	0.094	0.052	0.101
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.041	0.092	0.044	0.096	0.050	0.100

and

Table 6. Simulation results for the probability of type I error for $\rho \approx 0.75$  and $E\left(X_1\right)\ne E\left(X_2\right)$.

	$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=70$
	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$	$\mathbf{\alpha }=\mathbf{0.05}$	$\mathbf{\alpha }=\mathbf{0.10}$
$\gamma =(0.29,0.30,0.199475)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.031	0.084	0.037	0.089	0.060	0.109
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.030	0.083	0.036	0.088	0.061	0.108
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.030	0.084	0.036	0.089	0.059	0.109
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.036	0.084	0.044	0.089	0.057	0.106
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.086	0.045	0.091	0.056	0.105
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.085	0.045	0.089	0.055	0.107
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.089	0.047	0.094	0.052	0.104
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.090	0.045	0.096	0.053	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.092	0.045	0.096	0.053	0.101
$\gamma =(0.30,0.31,0.205476)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.034	0.076	0.039	0.085	0.058	0.111
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.034	0.075	0.038	0.086	0.057	0.110
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.035	0.075	0.039	0.086	0.058	0.111
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.084	0.043	0.087	0.054	0.109
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.086	0.044	0.089	0.055	0.108
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.039	0.086	0.044	0.089	0.054	0.108
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.044	0.091	0.046	0.096	0.052	0.101
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.091	0.045	0.095	0.053	0.102
$B_{n,(0,1)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.044	0.092	0.045	0.096	0.052	0.101
$\gamma =(0.30,0.32,0.208408)$
$B_{n,(0,0)}\left(\tilde{\gamma }\right)$	0.035	0.079	0.039	0.085	0.058	0.112
$B_{n,(1,0)}\left(\tilde{\gamma }\right)$	0.035	0.080	0.039	0.087	0.059	0.110
$B_{n,(0,1)}\left(\tilde{\gamma }\right)$	0.034	0.079	0.039	0.085	0.056	0.112
$B_{n,(0,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.084	0.044	0.088	0.054	0.106
$B_{n,(1,0)}\left(\tilde{\tilde{\gamma }}\right)$	0.037	0.084	0.044	0.088	0.054	0.109
$B_{n,(0,1)}\left(\tilde{\tilde{\gamma }}\right)$	0.038	0.085	0.045	0.089	0.053	0.107
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.090	0.046	0.094	0.051	0.100
$B_{n,(1,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.042	0.091	0.046	0.094	0.052	0.101
$B_{n,(0,0)}\left(\widetilde{\tilde{\tilde{\gamma }}}\right)$	0.043	0.092	0.047	0.096	0.050	0.101

summarize the obtained results. We denote the test statistic by $B_{n,a}\left(\widehat{\gamma }\right)$ when the weight function w takes the form given by (14) for some $a=(a_1,a_2)$. Additionally, note that we used the notation detailed in Section 5 to represent the technique employed to estimate the parameter $\gamma$.

In turn, as documented in the literature, the tables show that the maximum likelihood method provided better results than the other two estimation methods analyzed, followed by the method of moments. Furthermore, it is observed that the results converged to the established nominal level as the sample size increased.

These two conclusions can be clearly visualized in the graphs presented in Appendix A. On the one hand, it can be observed that as the sample size increased, the results approached the established nominal level: in the left area of each graph, the convergence toward $\alpha =0.05$ is evident, while in the right area, the convergence toward $\alpha =0.10$ is observable. Furthermore, it can be detected that using the maximum likelihood method yielded results closer to the established nominal levels ($\alpha =0.05$ on the left side and $\alpha =0.10$ on the right side) compared to the other two methods, with the method of moments following in accuracy.

Observing the values presented in the tables, we conclude that the parametric bootstrap method provides an accurate approximation to the null distribution of $B_{n,w}\left({\widehat{\gamma }}_n\right)$ in all the cases considered. This is also evident in the graphs in Appendix B.

6.2. Power of the Proposed Test Statistic

To analyze the power, we repeated the previous experiment for samples of size $n=50$ across four alternatives (see, for example, Kocherlakota and Kocherlakota [19] for a description of these distributions):

• Bivariate Binomial distribution $BB(m;p_1,p_2,p_3)$, where $p_1+p_2-p_3\le 1$, $p_1\ge p_3$, $p_2\ge p_3$, and $p_3>0$;
• Bivariate Hermite distribution $BH(\mu ,{\sigma }^2;{\lambda }_1,{\lambda }_2,{\lambda }_3)$, where $\mu >{\sigma }^2({\lambda }_1+{\lambda }_2+{\lambda }_3)$;
• Bivariate Neyman type A distribution $BNTA(\lambda ;{\lambda }_1,{\lambda }_2,{\lambda }_3)$, where $0<{\lambda }_1+{\lambda }_2+{\lambda }_3\le 1$;
• Bivariate Poisson distribution $BP({\lambda }_1,{\lambda }_2,{\lambda }_3)$, where ${\lambda }_1>{\lambda }_3$, and ${\lambda }_2>{\lambda }_3>0$.

The bivariate Binomial distribution and the bivariate Neyman type A distribution have been used as alternatives in other related articles (see, for example, Loukas and Kemp [29] and Rayner and Best [30]). The other alternative distributions employed are analogous to those used by Gürtler and Henze [31] in the univariate case.

To generate random samples from the bivariate distributions used as alternatives, we implemented computational algorithms following the simulation procedures provided in Kocherlakota and Kocherlakota [19].

Table 7. Simulation results for the power. The values are in the form of percentages rounded to the nearest integer.

Alternative	${\mathit{B}}_{(0,0)}$	${\mathit{B}}_{(1,0)}$	${\mathit{B}}_{(0,1)}$
$BB(1;0.41,0.02,0.01)$	95	91	88
$BB(1;0.55,0.03,0.02)$	97	92	89
$BB(2;0.42,0.02,0.01)$	98	94	90
$BB(2;0.51,0.01,0.01)$	94	89	88
$BB(2;0.71,0.04,0.03)$	92	86	85
$BH(0.99;1,0.66,0.10,0.10)$	96	87	92
$BH(1.00;1,0.30,0.30,0.25)$	95	86	90
$BH(1.00;1,0.31,0.33,0.24)$	92	86	86
$BH(1.30;1,0.61,0.03,0.02)$	93	84	85
$BH(1.40;1,1.00,0.26,0.12)$	95	86	88
$BNTA(0.15;0.35,0.33,0.32)$	94	87	93
$BNTA(0.42;0.32,0.35,0.33)$	92	89	91
$BNTA(0.50;0.31,0.34,0.35)$	93	89	92
$BNTA(0.70;0.35,0.30,0.33)$	91	88	92
$BNTA(0.75;0.32,0.34,0.34)$	94	89	93
$BP(1.00,1.00,0.25)$	97	87	90
$BP(1.00,1.00,0.50)$	94	83	92
$BP(1.00,1.00,0.75)$	96	82	90
$BP(1.50,1.00,0.31)$	95	86	89
$BP(1.50,1.00,0.92)$	96	86	91

displays the alternatives considered and the estimated power for nominal significance level $\alpha =0.05$. Additionally, to estimate the parameter $\gamma$, we only used the maximum likelihood method, as it yielded better results in the type I error study. Analyzing this table, we can conclude that all the considered tests, denoted by $B_{(a_1,a_2)}$, were able to detect the alternatives studied and with good power results, especially when $B_{(0,0)}$ was used.

All calculations performed in this research were programmed using the open-source software R [28]. To leverage the potential of parallel programming (see Novoa-Muñoz [18]), the codes were executed utilizing the supercomputing infrastructure of the NLHPC (ECM–02)² belonging to the Mathematical Modeling Center of the University of Chile [32]. This infrastructure includes, among other resources, 132 HP computing nodes (128 HP SL230 nodes and 4 HP SL250 nodes), each equipped with two 10-core Intel Xeon Ivy Bridge E5-2660V2 processors.The manufacturer of the NLHPC supercomputer was Atos, located in the city of Bezons, France.

6.3. Real Data Set

The test statistic was applied to two real datasets:

First Case: Biometric Data. This case was analyzed by Dunn [33] to illustrate and represent a real case of the BNBD. The sample was obtained from Arbous and Kerrich [4] in a study on the treatment and analysis of the propensity for workplace accidents in the railway sector among a group of experienced workers. The study involved 122 workers, with the first variable representing a six-year period (1937–1942) and the second variable representing a five-year period (1943–1947). The results correspond to the number of workers who had accidents in the two non-overlapping time periods.

Second Case: Agronomic–Botanical Data. This dataset is based on the number of plants of the species Lacistema aggregatum and Protium guianense in each of 100 contiguous quadrants. The data were first presented by Holgate [6] and later analyzed by various researchers up to 1995, including Gillings [34], Crockett [35], Loukas and Kemp [29], Kocherlakota and Kocherlakota [19], and Rayner and Best [30]. They were also studied by Novoa-Muñoz and Jiménez-Gamero [17]. Researchers Crockett [35], Loukas and Kemp [29], Rayner and Best [30], and Novoa-Muñoz and Jiménez-Gamero [17] examined the data to determine whether they corresponded to a bivariate Poisson model and concluded that the bivariate Poisson distribution did not correctly model the mentioned data.

The following three cases are part of the kaggle.com [36] repository, where companies from around the world can upload their databases for a community of data scientists to generate predictive models. It is a challenge of skills and competencies among developers who seek to identify patterns and trends in the data. These cases correspond to current data recorded in recent times.

kaggle.com [36] also offers modeling and evaluation tools to develop and compare models. It is an online community for learning, sharing, and competing in data science.

It is important to note that, based on the information available to date, we have no evidence that any of these three datasets have previously been analyzed to determine whether they follow a probability distribution.

Third Case: Black Friday Sales Data. This is a public dataset containing information about purchases made during Black Friday. It includes buyer variables such as age, gender, purchase amount, category, and quantity of products purchased. The dataset contains over 550,000 rows and 12 columns. Two columns containing the quantity of products purchased were randomly selected. The dataset can be used to identify patterns and trends in Black Friday purchases.

Fourth Case: Jamboree Education Data. This dataset contains information on students’ academic performance in admission exams. It includes variables such as scores in different sections of the exam, the student’s age, and gender. The objective is to predict students’ academic performance. This dataset is useful for analyzing and improving educational outcomes.

Fifth Case: Insurance Data. This database contains information about insurance customers. It includes variables such as age, gender, marital status, and income level. It also contains information about insurance policies, such as premium amount and insurance type. The dataset aims to predict the likelihood of a customer renewing their insurance policy. It is useful for analyzing and improving customer retention in the insurance industry.

Table 8. p-values for the real datasets rounded to three decimal places.

Test Statistics	Biometric	Agronomic	Sales	Education	Insurance
$B_{n\left(0,0\right)}$	0.434	0.277	0.003	0.892	0.493
$B_{n\left(1,0\right)}$	0.443	0.327	0.037	0.836	0.665
$B_{n\left(0,1\right)}$	0.440	0.303	0.023	0.878	0.258

shows the p-values for testing the goodness of fit to a BNBD for the five real datasets using the statistical test we have proposed. For this, we used $(a_1,a_2)\in \{(0,0),(1,0),(0,1)\}$. Additionally, for the five real cases, we applied a nominal significance level of $\alpha =0.05$ and considered $v=5$, which was employing the same number of trials before the first event occurred.

The results displayed in this table allow us to conclude the following: For the Black Friday Sales data, the null hypothesis is rejected, meaning that Black Friday sales are not well modeled by a DBNB. On the other hand, for the other four datasets, it is not possible to reject the null hypothesis. Therefore, we conclude that the biometric, agronomic, education, and insurance data can be well modeled by a DBNB.

6.4. Comparisons with Other Goodness-of-Fit Tests

The only goodness-of-fit tests we have found in the statistical literature for bivariate count data are for the bivariate Poisson distribution and the bivariate Hermite distribution. The similarity is that these tests can be defined as follows:

• Consistent: As the sample size increases, the probability of correctly rejecting a false null hypothesis approaches 1. In other words, a consistent test has the ability to detect any deviation from the null hypothesis when sufficient information is available (i.e., with a sufficiently large sample size).
• Based on the probability generating function.
• Of the Cramér-von Mises type.
• With null asymptotic distribution approximated by the parametric bootstrap method.
• Specific tests for bivariate count data.

However, their implementation complexity may be a disadvantage compared to more generic tests like Kolmogorov–Smirnov, Anderson–Darling, or chi-squared.