2.2.2. Probability-Generating Function

The probability-generating function of the random variable *X* ∼ BNDL(*p*) can be obtained using its moment-generating function which is equivalent to calculating *E tX* ; therefore, the probability-generating function of the random variable *X* is

$$G\_X(t) = E\left(t^X\right) = M\_X(\log(t)) = \frac{1 - p(t - 2) + p^2(t - 1)}{\left(2 - t + p(t - 1)\right)^2 (p + 1)}$$

.

Since,

$$G\_X^{(k)}(t) = \frac{d^k G\_X(t)}{dt^k} = E\left\{X(X-1)(X-2)\dots(X-k+1)t^{X-k}\right\}.$$

Therefore, at *t* = 1, we can obatin

$$G\_X^{(k)}(1) = \left. \frac{d^k G\_X(t)}{dt^k} \right|\_{t=1} = E\{X(X-1)(X-2)\dots(X-k+1)\},$$

where *μ*(*k*) = *E*{*X*(*X* − 1)(*X* − 2)...(*X* − *k* + 1)} is the *k*th factorial moment of *X*.

**Figure 2.** Histograms of the BNDL model for simulated data.

2.2.3. Non-Central Moments and Variance

If *X* ∼ BNDL(*p*) distribution, then the *k*th moment about zero of *X* is given by

$$\mu\_k' = E(X^r) = \sum\_{x=0}^{\infty} x^k p\_x = \sum\_{x=0}^{\infty} x^k \frac{(1-p)^x \left(1+x+2p-p^2\right)}{\left(p+1\right) \left(2-p\right)^{x+2}}.$$

The first four raw moments can be obtained as follows

$$
\mu\_1' = E(X) = \frac{(p+2)(1-p)}{p+1},
$$

$$
\mu\_2' = E\left(X^2\right) = \frac{(1-p)\left(8-3p-2p^2\right)}{p+1},
$$

$$
\mu\_3' = E\left(X^3\right) = \frac{(1-p)\left(44-53p+6p^2+6p^3\right)}{p+1},
$$

$$
\mu\_4 = \frac{1}{2}
$$

and

$$
\mu\_4' = E\left(X^4\right) = \frac{(1-p)\left(308 - 516p + 346p^2 - 12p^3 - 24p^4\right)}{p+1}.
$$

The variance in the random variable *X* is

$$Var(X) = E\left(X^2\right) - \left[E(X)\right]^2 = \frac{(1-p)\left(4+5p-2p^2-p^3\right)}{\left(p+1\right)^2}.$$

#### 2.2.4. Central Moments

The *k*th moment about the mean of *X* is

$$\mu\_r = E\left[\left(X - \mu\_1'\right)^k\right] = \sum\_{x=0}^{\infty} \left(\mathbf{x} - \mu\_1'\right)^k p\_x(\mathbf{x}) = \sum\_{x=0}^{\infty} \left(\mathbf{x} - \mu\_1'\right)^k \frac{(1-p)^x \left(1+\mathbf{x}+2p-p^2\right)}{\left(p+1\right)\left(2-p\right)^{x+2}}.$$

Therefore, the second, third and fourth central moments of the random variable *X* are

$$\mu\_2 = \frac{(1-p)\left(4+5p-2p^2-p^3\right)}{\left(p+1\right)^2},$$

$$\mu\_3 = \frac{\left(1-p\right)\left(12+21p-7p^2-21p^3+5p^4+2p^5\right)}{\left(p+1\right)^3},$$

and

$$\mu\_4 = \frac{(1-p)\left(100+181p-132p^2-285p^3+50p^4+137p^5-27p^6-9p^7\right)}{\left(p+1\right)^4}$$

#### 2.2.5. Skewness and Kurtosis

The coefficient of skewness and the coefficient of kurtosis of the of BNDL distribution are, respectively,

$$\beta\_1 = \frac{\mu\_3}{\sqrt{\mu\_2^3}} = \frac{(1-p)\left(12+21p-7p^2-21p^3+5p^4+2p^5\right)}{\left(4+p-7p^2+p^3+p^4\right)^{3/2}}.$$

$$\beta\_2 = \frac{\mu\_4}{\mu\_2^2} = \frac{100+181p-132p^2-285p^3+50p^4+137p^5-27p^6-9p^7}{\left(1-p\right)\left(4+5p-2p^2-p^3\right)^2}.$$

#### 2.2.6. Index of Dispersion

The index of dispersion (ID) indicates whether a certain distribution is suitable for under- or over-dispersed datasets. For example, ID = 1 for the Poisson distribution where the variance is equal to the mean, for the geometric distribution and the negative binomial distribution ID > 1, while the binomial distribution has ID < 1.

**Theorem 1.** *If X* ∼ *BNDL*(*p*)*, then Var*(*X*) > *E*(*X*) *for all p* ∈ (0, 1).

**Proof.** We have

$$\text{ID}(X) = \frac{\text{Var}(X)}{E(X)} = \frac{4 + 5p - 2p^2 - p^3}{p^2 + 3p + 2}.$$

This function is a monotonic decreasing function as *p* ∈ (0, 1) increases. It converges to 2 when *p* → 0, while it tends to 1 as *p* → 1; therefore, ID(*X*) ∈ (1, 2), which means that ID(*X*) > 1, and hence, *Var*(*X*) > *E*(*X*).

From Theorem 1, BNDL distribution should only be used in the count data analysis with over-dispersion. In Table 1, some of the empirical findings of these measured are due for considerations.


**Table 1.** Mean, Variance, Skewness, kurtosis and ID of the BNDL distribution for different values of the parameter *p*.

#### 2.2.7. Log-Concavity

A necessary and sufficient condition that *px* be strongly unimodal is that it has to be log-concave, i.e., *p*<sup>2</sup> *<sup>x</sup>*+<sup>1</sup> ≥ *px px*+<sup>2</sup> for all *x* (see Keilson and Gerber [20])).

**Theorem 2.** *The pmf of the BNDL distribution in (3) is log-concave.*

**Proof.** From (3), we can directly reach

$$p\_{x+1}^2 = \frac{(1-p)^{2x+2} \left(2+x+2p-p^2\right)^2}{\left(p+1\right)^2 \left(2-p\right)^{2x+6}}.$$

and

$$p\_x p\_{x+2} = \frac{(1-p)^{2x+2} \left(1+x+2p-p^2\right) \left(3+x+2p-p^2\right)}{\left(p+1\right)^2 \left(2-p\right)^{2x+6}}.$$

After some algebraic operations, we find that

$$p\_{x+1}^2 - p\_x p\_{x+2} = \frac{(1-p)^{2x+2}}{(p+1)^2 (2-p)^{2x+6}} > 0,$$

for all *x* and for all choices *p* ∈ (0, 1).

Theorem 2 confirms that the BNDL distribution is strongly unimodal.

*2.3. Reliability Properties of the BNDL Distribution*

## 2.3.1. Survival Function

If *X* ∼ BNDL(*p*) distribution, then from (4), the survival function of *X* is

$$S(\mathbf{x}; p) = P(X \ge \mathbf{x}) = \frac{(1 - p)^{\mathbf{x} + 1} \left( \mathbf{3} + \mathbf{x} + p - p^2 \right)}{\left( p + 1 \right) \left( 2 - p \right)^{\mathbf{x} + 2}}.$$

2.3.2. Hazard Rate and Mean Residual Life Functions

The hazard (failure) rate function is the probability that an item has survived time *x*, given that it has survived to at least time *x*. If *X* ∼ BNDL(*p*) distribution, then its hazard rate (failure rate) function is given as

$$r(\mathbf{x};p) = P(X=\mathbf{x}|X>\mathbf{x}) = \frac{p\_{\mathbf{x}}(\mathbf{x};p)}{S(\mathbf{x};p)} = \frac{1+\mathbf{x}+2p-p^2}{(1-p)(3+\mathbf{x}+p-p^2)}.$$

Obviously, the upper limit of the failure rate function is <sup>1</sup> <sup>1</sup>−*<sup>p</sup>* , i.e., lim*x*→∞*r*(*x*; *<sup>p</sup>*) <sup>=</sup> <sup>1</sup> <sup>1</sup>−*<sup>p</sup>* . Graphical illustrations of hazard rate function are presented in Figure 3 while descriptive measures are presented in Figure 4.

**Figure 3.** Plots of hazard rate of BNDL distribution for some choices of *p*.

The mean residual life function of *X* is given by

$$m(\mathbf{x}; p) = P(X - \mathbf{x} | X > \mathbf{x}) = \frac{\sum\_{t=x+1}^{\infty} S(t; p)}{S(\mathbf{x}; p)} = \frac{(p - 1)\left(p^2 - \mathbf{x} - \mathbf{5}\right)}{3 + p - p^2 + \mathbf{x}}.$$

**Corollary 1.** *If X* ∼ *BNDL*(*p*) *distribution, then it has an increasing failure rate and decreasing mean residual life.*

As we explained through Theorem 2, the BNDL distribution has a property of logconcavity; therefore, according to Gupta et al. [21], the BNDL distribution has an IFR property. According to Kemp [22], the next chain is verified

IFR ⇒ IFRA ⇒ NBU ⇒ NBUE ⇒ DMRL.

So, the BNDL distribution is


**Figure 4.** Plots of the BNDL model for (**a**) Mean, (**b**) Variance, (**c**) Skewness, (**d**) Kurtosis and (**e**) ID.

#### *2.4. Stochastic Orderings*

Stochastic orders are important measures to judge comparative behaviors of random variables. Shaked and Shanthikumar [8] showed that many stochastic orders exist and have various applications. Given two random variables *X* and *Y*, we say that *X* is smaller than *Y* in the


$$X \leq\_{lr} Y \Rightarrow X \leq\_{hr} Y \Rightarrow X \leq\_{st} Y \; \prime$$

and

$$X \leq\_{lr} Y \Rightarrow X \leq\_{rh} Y \Rightarrow X \leq\_{st} Y$$

also,

$$X \leq\_{lr} Y \Rightarrow X \leq\_{surl} Y.$$

**Theorem 3.** *Let X* ∼ *BNDL*(*p*1) *and Y* ∼ *BNDL*(*p*2)*; then, X* ≤*lr Y for all p*<sup>1</sup> > *p*2.

**Proof.** Let

$$L(\mathbf{x}; p\_1, p\_2) = \frac{p\_X(\mathbf{x}; p\_1)}{p\_Y(\mathbf{x}; p\_2)}.$$

Now,

$$L(\mathbf{x}; p\_1, p\_2) = \frac{(p\_2 + 1)(2 - p\_2)^{\mathbf{x} + 2}(1 - p\_1)^{\mathbf{x}}\left(1 + \mathbf{x} + 2p\_1 - p\_1^2\right)}{(p\_1 + 1)(2 - p\_1)^{\mathbf{x} + 2}(1 - p\_2)^{\mathbf{x}}\left(1 + \mathbf{x} + 2p\_2 - p\_2^2\right)}.$$

and

$$L(\mathbf{x}+1;\ p\_1,p\_2) = \frac{(p\_2+1)(2-p\_2)^{\mathbf{x}+3}(1-p\_1)^{\mathbf{x}+1}\left(2+\mathbf{x}+2p\_1-p\_1^2\right)}{\left(p\_1+1\right)\left(2-p\_1\right)^{\mathbf{x}+3}\left(1-p\_2\right)^{\mathbf{x}+1}\left(2+\mathbf{x}+2p\_2-p\_2^2\right)}.$$

Therefore,

$$\frac{L(\mathbf{x}+1;\ p\_1,p\_2)}{L(\mathbf{x};p\_1,p\_2)} = \frac{(2-p\_2)(1-p\_1)\left(2+\mathbf{x}+2p\_1-p\_1^2\right)\left(1+\mathbf{x}+2p\_2-p\_2^2\right)}{(2-p\_1)(1-p\_2)\left(2+\mathbf{x}+2p\_2-p\_2^2\right)\left(1+\mathbf{x}+2p\_1-p\_1^2\right)}\tag{5}$$

,

Let *p*<sup>1</sup> = 1 − *δ* and *p*<sup>2</sup> = 1 − *δ* − *ε*, where 0 < *δ* < 1 and 0 < *ε* < 1 − *δ*. After substitution of the values *p*<sup>1</sup> and *p*<sup>2</sup> in (5), we obtain

$$\frac{L(\mathbf{x} + \mathbf{1}; \, p\_1, p\_2)}{L(\mathbf{x}; p\_1, p\_2)} = \frac{\eta\_1 \left(\delta + \delta^2 + \delta \varepsilon\right)}{\eta\_2 \left(\delta + \delta \varepsilon + \delta^2 + \varepsilon\right)}$$

where

$$
\eta\_1 = \left(\mathfrak{Z} + \mathfrak{x} - \delta^2\right) \left(\mathfrak{Z} + \mathfrak{x} - \left(\delta + \varepsilon\right)^2\right),
$$

and

$$
\eta\_2 = \left( 3 + \mathfrak{x} - \left( \delta + \mathfrak{e} \right)^2 \right) \left( 2 + \mathfrak{x} - \left( \delta \right)^2 \right) .
$$

After some algebraic operations, we find that

$$
\eta\_1 - \eta\_2 = -\varepsilon (2\delta + \varepsilon) < 0 \Rightarrow \eta\_1 < \eta\_2.
$$

Therefore,

$$
\eta\_1 \left( \delta + \delta^2 + \delta \varepsilon \right) < \eta\_2 \left( \delta + \delta \varepsilon + \delta^2 + \varepsilon \right).
$$

This implies that

$$\frac{L(\mathbf{x}+1;\ p\_1,p\_2)}{L(\mathbf{x};p\_1,p\_2)} < 1 \Rightarrow L(\mathbf{x}+1;\ p\_1,p\_2) < L(\mathbf{x};p\_1,p\_2).$$
