On some aspects of a bivariate alternative zero-inflated logarithmic series distribution

Abstract

In this paper, we discuss some important aspects of the bivariate alternative zero-inflated logarithmic series distribution (BAZILSD) of which the marginals are the alternative zero-inflated logarithmic series distributions of Kumar and Riyaz (2015. An alternative version of zero-inflated logarithmic series distribution and some of its applications. Journal of Statistical Computation and Simulation, 85(6), 1117–1127). We study some important properties of the distribution by deriving expressions for its probability mass function, factorial moments, conditional probability generating functions, and recursion formulae for its probabilities, raw moments and factorial moments. The parameters of the BAZILSD are estimated by the method of maximum likelihood and certain test procedures are also considered. Further certain real-life data applications are cited for illustrating the usefulness of the model. A simulation study is conducted for assessing the performance of the maximum likelihood estimators of the parameters of the BAZILSD.

KEYWORDS:

• Factorial moments
• generalized likelihood ratio test
• probability generating function
• Rao’s score test
• zero-inflated distributions

1. Introduction

Bivariate discrete distributions have received much attention in the literature. For example, see Ghosh and Balakrishnan (2015), Hassan and El-Bassiouni (2013), Kemp (2013), Kumar (2008), Kocherlakota and Kocherlakota (1992) and references therein. Due to the extensive applications of logarithmic series distribution in various areas of scientific research especially in biology, ecology, meteorology, etc., the bivariate logarithmic series distribution (BLSD) is of particular interest. Chapter 7 of Kocherlakota and Kocherlakota (1992) is fully devoted to the BLSD. Subrahmaniam (1966) defined the BLSD through the following probability generating function (pgf) (1.1) $A(t_1,t_2)=\frac{-\log (1-{\theta }_1{t}_1-{\theta }_2{t}_2-{\theta }_3{t}_1{t}_2)}{-\log (1-{\theta }_1-{\theta }_2-{\theta }_3)},$(1.1) in which ${\theta }_1>0$, ${\theta }_2>0$ and ${\theta }_3\ge 0$ such that ${\theta }_1+{\theta }_2+{\theta }_3<1$. An important drawback of the BLSD in practical point of view is that it excludes the (0, 0)-th observation from its support. To overcome this difficulty, Kumar and Riyaz (2014) considered a class of bivariate distribution namely the ‘bivariate zero-inflated logarithmic series distribution (BZILSD)’ through the following probability mass function (pmf), for any non-negative integers $m$ and $n$, ${\theta }_1>0$, ${\theta }_2>0$ and ${\theta }_3\ge 0$ such that ${\theta }_1+{\theta }_2+{\theta }_3\le 1$. (1.2) $f(m,n)=\delta {\theta }_1^m{\theta }_2^n\sum _{r=0}^{min(m,n)}\frac{D_r^{\ast }}{r!(m-r)!(n-r)!}{\left(\frac{{\theta }_3}{{\theta }_1{\theta }_2}\right)}^r,$(1.2) where $D_r^{\ast }=\underset{j=0}{\overset{m+n-r-1}{\mathrm{\Pi }}}\frac{{(j+1)}^2}{(j+2)}$ and $\delta =\left[F_{2,1}\right(1,1;2;{\theta }_1+{\theta }_2+{\theta }_3){]}^{-1}$ in which $F_{2,1}(a,b;c;z)$ is the Gauss hypergeometric function (cf. Mathai & Haubold, 2008).

Kumar and Riyaz (2013) considered the zero-inflated logarithmic series distribution (ZILSD) through the following pgf, in which $A=\theta [-\ln (1-\theta ){]}^{-1}$ with $\theta \in (0,1)$. (1.3) $G_1\left(t\right)=A{F}_{2,1}(1,1;2;\theta t),$(1.3) or equivalently, (1.4) $G_1\left(t\right)=-A{t}^{-1}\ln (1-\theta t).$(1.4) Kumar and Riyaz (2015) considered another zero-inflated logarithmic series distribution, which they termed as ‘the alternative zero-inflated logarithmic series distribution (AZILSD)’, through the following pmf, for $x=0,1,2,\dots$, (1.5) $g\left(x\right)=B\frac{F_{2,1}(1+x,1+x;2+x;\alpha ){\theta }^x}{1+x},$(1.5) in which $B=[-\ln (1-\theta -\alpha ){]}^{-1}(\theta +\alpha )$, $\theta >0$, $\alpha \ge -1$ and $|\theta +\alpha |<1$ such that $\theta \ne -\alpha$. The pgf of the AZILSD with pmf (1.5) is (1.6) $G_2\left(t\right)=B{F}_{2,1}(1,1;2;\theta t+\alpha )$(1.6) or equivalently, (1.7) $G_2\left(t\right)=B\frac{-\ln (1-\theta t-\alpha )}{(\theta t+\alpha )}.$(1.7) Kumar and Riyaz (2017) studied an extended version of AZILSD and its important properties. Kumar and Riyaz (2016) considered an order $k$ version of AZILSD and studied its important applications.

Through this paper, we consider a bivariate version of the AZILSD through the name ‘the bivariate alternative zero-inflated logarithmic series distribution’ or, in short ‘the BAZILSD’, and discuss some of its important aspects. In Section 2, we derive the BAZILSD as a bivariate random sum distribution of independent and identically distributed bivariate Bernoulli random variables and show that the marginal distributions of the BAZILSD are AZILSD. We obtain expressions for its pmf, mean, covariance, factorial moments and conditional pgfs which are included in Section 2. In Section 3, we derive certain recursion formulae for probabilities, raw moments and factorial moments of the BAZILSD. In Section 4, we describe the estimation of the parameters of the BAZILSD by method of maximum likelihood and certain test procedures are suggested. And in Section 5, we illustrate the usefulness of the BAZILSD through fitting the distribution to certain real-life data sets. In Section 6, a brief simulation study is conducted for examining the performance of the maximum likelihood estimators of the parameters of the BAZILSD.

It is important to note that the BAZILSD possesses a bivariate random sum structure as shown in Section 2. Certain bivariate random sum distributions are studied in the literature. For example, see Kumar (2007, 2013). The random sum structure arises in several areas of scientific research particularly in actuarial science, agricultural science, biological science and physical science. Chapter 9 of Johnson et al. (2005) fully devoted to univariate random sum distributions.

For simplicity in the notations, we adopt the following notations throughout in the manuscript. (1.8) $\begin{array}{rl}{R}_j\left(\theta \right)& =F_{2,1}(1+j,1+j;2+j;\theta +\alpha ),\end{array}$(1.8) (1.9) $\begin{array}{rl}\mathrm{\Lambda }& =R_0^{-1}\left(\theta \right),\end{array}$(1.9) (1.10) $\begin{array}{rl}{h}_j& =R_j({\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha ),\end{array}$(1.10) (1.11) $\begin{array}{rl}{\psi }_j& =R_j\left(\theta \right),\end{array}$(1.11) (1.12) $\begin{array}{rl}{\text{β}}_j& =R_j\left(0\right),\end{array}$(1.12) (1.13) $\begin{array}{rl}\text{and}H(t_1,t_2)& =\mathrm{\Lambda }{R}_0({\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha ).\end{array}$(1.13)

2. A genesis and some properties of the BAZILSD

First, we derive the BAZILSD in the following and discuss some of its properties.

Consider the sequence $\left\{Y_n=(Y_{1n},Y_{2n});n\ge 1\right\}$ of independent and identically distributed bivariate Bernoulli random vectors, each with pgf $P(t_1,t_2)={\lambda }_1{t}_1+{\lambda }_2{t}_2+{\lambda }_3{t}_1{t}_2,$in which ${\lambda }_j=\frac{{\theta }_j}{\theta }$, $j=1,2,3$ with $\theta ={\theta }_1+{\theta }_2+{\theta }_3$ such that ${\theta }_1>0$, ${\theta }_2>0$ and ${\theta }_3\ge 0.$ Let $X$ be a non-negative integer valued random variable having AZILSD with pgf (1.6), in which $\theta ={\theta }_1+{\theta }_2+{\theta }_3$. Assume that $\left\{Y_n:n\ge 1\right\}$ and $X$’s are independent. Define $S_n=(S_{1n},S_{2n})$, for each $n\ge 0$ in which $(S_{10},S_{20})=(0,0)$ and $S_{rm}=\sum _{j=1}^m{Y}_{rj}$, for $r=1,2$ and $m\ge 1$. Set $S_X=\sum _{n=0}^{\mathrm{\infty }}{S}_n{I}_{\left[X=n\right]}$ where $I_{\left[X=n\right]}$ denotes the indicator function of an event $[X=n]$. Then the pgf of $S_X$ is (2.1) $\begin{array}{rl}H(t_1,t_2)& =G_2\left\{P\right(t_1,t_2\left)\right\}\\ & =\mathrm{\Lambda }{F}_{2,1}(1,1;2;{\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha ),\end{array}$(2.1) where $\mathrm{\Lambda }$ is defined in (1.9).

We call a distribution with pgf (2.1) ‘the bivariate alternative zero-inflated logarithmic series distribution’ or, in short ‘the BAZILSD’. Clearly when $\alpha =0$, the pgf given in (2.1) reduces to the following pgf of the BZILSD with pmf (1.2). (2.2) $B(t_1,t_2)=\frac{F_{2,1}(1,1;2;{\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2)}{F_{2,1}(1,1;2;{\theta }_1+{\theta }_2+{\theta }_3)},$(2.2) which shows that the proposed bivariate model of the AZILSD can be considered as a more flexible model in practical point of view compared to the BZILSD. Further, it can be noted that the marginals of the BAZILSD are AZILSD whereas the marginals of the BZILSD are not ZILSD.

Proposition 2.1.

If $V=(V_1,V_2)$ follows the BAZILSD, then the marginal distribution of $V_j$ for $j=1,2$ is AZILSD with pgf given below. $H_{V_1}\left(t\right)=\mathrm{\Lambda }{F}_{2,1}[1,1;2;({\theta }_1+{\theta }_3)t+{\theta }_2+\alpha ]$and $H_{V_2}\left(t\right)=\mathrm{\Lambda }{F}_{2,1}[1,1;2;({\theta }_2+{\theta }_3)t+{\theta }_1+\alpha ].$

The proof follows from the fact that $H_{V_1}\left(t\right)=H(t,1)$ and $H_{V_2}\left(t\right)=H(1,t)$.

Proposition 2.2.

The pgf of the conditional distribution of $V_1$ given $V_2=v$ is the following: for any non-negative integer $v$, (2.3) $H_{V_1|V_2}\left(t\right)={\left(\frac{{\theta }_2+{\theta }_{3}t}{{\theta }_2+{\theta }_3}\right)}^v\frac{F_{2,1}(1+v,1+v;2+v;{\theta }_{1}t+\alpha )}{F_{2,1}(1+v,1+v;2+v;{\theta }_1+\alpha )}.$(2.3)

Proof:

For any non-negative integer $v,$ assume that $P(V_2=v)>0$. Now, we have the following partial derivatives of order $(0,v)$ of $H(t_1,t_2)$ with respect to $t_2$ evaluated at $(t_1,t_2)=(t,0)$. (2.4) $H^{(0,v)}(t,0)=\mathrm{\Lambda }({\theta }_2+{\theta }_{3}t{)}^v(\prod _{j=0}^{v-1}{D}_j)R_v\left({\theta }_{1}t\right),$(2.4) where for $j=0,1,2,\dots ,$ (2.5) $D_j=\frac{{(j+1)}^2}{(j+2)}$(2.5) and $R_j\left(t\right)$ is defined in (1.8).

Now, applying the formula for the conditional pgf in terms of partial derivatives of the joint pgf developed by Subrahmaniam (1966), we obtain the conditional pgf of $V_1$ given $V_2=v$ as $\begin{array}{rl}{H}_{V_1|V_2=v}\left(t\right)& =\frac{H^{(0,v)}(t,0)}{H^{(0,v)}(1,0)}\\ & =\frac{{({\theta }_2+{\theta }_{3}t)}^v}{{({\theta }_2+{\theta }_3)}^v}\frac{R_v\left({\theta }_{1}t\right)}{R_v\left({\theta }_1\right)},\end{array}$which implies (2.3) in the light of (1.8).

Remark 2.1.

The conditional distribution of $V_1$ given $V_2=v$ as given in (2.3) can be written as $H_{V_1|V_2}\left(t\right)=H_{Z_1}\left(t\right)H_{Z_2}\left(t\right)$ where $H_{Z_1}\left(t\right)$ is the pgf of a binomial random variable with parameters $z_1$ and $p=\frac{{\theta }_3}{{\theta }_2+{\theta }_3}$ and $H_{Z_2}\left(t\right)$ is the pgf of a random variable following the AZILSD with parameters $v$, ${\theta }_1$ and $\alpha$. Thus clearly, the conditional distribution $V_1$ given $V_2=v$ is the distribution of the sum of two independent random variables $Z_1$ and $Z_2$.

By using Remark 2.1, we obtain the following proposition.

Proposition 2.3:

Let $V=(V_1,V_2)$ follow the BAZILSD with pgf (2.1). Then (2.6) $\begin{array}{rl}\mathrm{E}(V_1|V_2=v)& =\frac{v{\theta }_3}{({\theta }_2+{\theta }_3)}+\frac{{\theta }_1{D}_v{R}_{v+1}({\theta }_1+\alpha )}{R_v({\theta }_1+\alpha )},\end{array}$(2.6) (2.7) $\begin{array}{rl}\mathrm{V}\mathrm{a}\mathrm{r}(V_1|V_2=v)& =\frac{v{\theta }_2{\theta }_3}{{({\theta }_2+{\theta }_3)}^2}+\frac{{\theta }_1{D}_v}{R_v^2({\theta }_1+\alpha )}[D_{v+1}{R}_v({\theta }_1+\alpha )R_{v+2}({\theta }_1+\alpha ){\theta }_1\\ & +R_v({\theta }_1+\alpha )R_{v+1}({\theta }_1+\alpha )-D_v{R}_{v+1}^2({\theta }_1+\alpha ){\theta }_1].\end{array}$(2.7)

Remark 2.2:

By a similar approach, for any non-negative integer $v$ with $P(V_1=v)>0,$ we can obtain the conditional pgf of $V_2$ given $V_1=v$ by interchanging ${\theta }_1$ and ${\theta }_2$ in (2.3). Therefore, it is evident that comments similar to those in Remark 2.1 are valid regarding conditional distribution of $V_2$ given $V_1=v$ and the explicit expression for $\mathrm{E}(V_2|V_1=v)$ and $\text{Var}(V_2|V_1=v)$ can be obtained by interchanging ${\theta }_1$ and ${\theta }_2$ in the right hand side expressions of (2.6) and (2.7) respectively.

Proposition 2.4 :

Let $V=(V_1,V_2)$ follow the BAZILSD with pgf (2.1) and let $m,n$ be any non-negative integers. The pmf $f(m,n)$ and the $(m,n)$-th factorial moment ${\mu }_{\left[m,n\right]}$ of the BAZILSD are (2.8) $\begin{array}{rl}f(m,n)& =\mathrm{\Lambda }{\theta }_1^m{\theta }_2^n\sum _{r=0}^{min(m,n)}\frac{{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }}{r!(m-r)!(n-r)!}{\left(\frac{{\theta }_3}{{\theta }_1{\theta }_2}\right)}^r,\end{array}$(2.8) (2.9) $\begin{array}{rl}{\mu }_{\left[m,n\right]}& =\mathrm{\Lambda }m!n!({\theta }_1+{\theta }_3{)}^m({\theta }_2+{\theta }_3{)}^n\sum _{r=0}^{min(m,n)}\frac{{D_r}^{\ast }{\psi }_{m+n-r}}{r!(m-r)!(n-r)!}{\xi }^r,\end{array}$(2.9) where $D_r^{{}^{\ast }}$ is defined in (1.2), for $j=1,2,\dots ,$ ${\psi }_j$, ${\text{β}}_j\left(\alpha \right)$’s are defined in (1.11) and (1.12) and $\xi =\frac{{\theta }_3}{({\theta }_1+{\theta }_3)({\theta }_2+{\theta }_3)}$.

Proof :

In order to obtain the probability mass function of the BAZILSD, we need the following derivatives of $H(t_1,t_2)$, in which $m$ is a non-negative integer. (2.10) $H^{(m,0)}(t_1,t_2)=\left(\prod _{i=0}^{m-1}{D}_i\right)({\theta }_1+{\theta }_3{t}_2{)}^m\mathrm{\Lambda }{h}_m(t_1,t_2),$(2.10) where (2.11) $h_j(t_1,t_2)=F_{2,1}(1+j,1+j;2+j;{\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha ),j=0,1,2,\dots .$(2.11) The following derivatives are needed in the sequel, in which $0\le i\le r$ and $j\ge 1$. (2.12) $\begin{array}{rl}\frac{{\mathrm{\partial }}^i{({\theta }_1+{\theta }_3{t}_2)}^m}{\mathrm{\partial }{t}_2^i}& =\frac{m!}{(m-i)!}{\theta }_3^i({\theta }_1+{\theta }_3{t}_2{)}^{m-i},\end{array}$(2.12) (2.13) $\begin{array}{rl}\frac{{\mathrm{\partial }}^j{h}_m(t_1,t_2)}{\mathrm{\partial }{t}_2^j}& =\left(\prod _{i=m}^{m+j-1}{D}_i\right)({\theta }_2+{\theta }_3{t}_1{)}^j{h}_{m+j}(t_1,t_2).\end{array}$(2.13) Differentiating both sides of (2.10) $n$ times with respect to $t_2$ and applying (2.12) and (2.13), we get the following. (2.14) $\begin{array}{rl}{H}^{(m,n)}(t_1,t_2)& =\left(\prod _{i=0}^{m-1}{D}_i\right)\mathrm{\Lambda }\sum _{r=0}^n\left(\begin{array}{c}n\\ r\end{array}\right)\frac{{\mathrm{\partial }}^r{({\theta }_1+{\theta }_3{t}_2)}^m}{\mathrm{\partial }{t}_2^r}\frac{{\mathrm{\partial }}^{n-r}{h}_m(t_1,t_2)}{\mathrm{\partial }{t}_2^{n-r}}\\ & =\left(\prod _{i=0}^{m-1}{D}_i\right)\mathrm{\Lambda }\sum _{r=0}^{min(m,n)}\left(\begin{array}{c}n\\ r\end{array}\right)\frac{m!}{(m-r)!}{\theta }_3^r({\theta }_1+{\theta }_3{t}_2{)}^{m-r}\\ & \times \left(\prod _{i=m}^{m+n-r-1}{D}_i\right)({\theta }_2+{\theta }_3{t}_1{)}^{n-r}{h}_{m+n-r}(t_1,t_2).\end{array}$(2.14) By putting $(t_1,t_2)=(0,0)$ in (2.14) and by dividing $m!n!$, we get (2.8). By putting $(t_1,t_2)=(1,1)$ in (2.14), we get (2.9).

Proposition 2.5 :

Let $V=(V_1,V_2)$ follow the BAZILSD with pgf (2.1). Then we have the following, in which ${\delta }_j=\frac{{\psi }_j}{{\psi }_0}$, (2.15) $\begin{array}{rl}\mathrm{E}\left(V_1\right)& =D_0{\delta }_1({\theta }_1+{\theta }_3),\end{array}$(2.15) (2.16) $\begin{array}{rl}\mathrm{E}\left(V_2\right)& =D_0{\delta }_1({\theta }_2+{\theta }_3),\end{array}$(2.16) and (2.17) $\mathrm{C}\mathrm{o}\mathrm{v}(V_1,V_2)=D_0(D_1{\delta }_2-D_0{\delta }_1^2)({\theta }_1+{\theta }_3)({\theta }_2+{\theta }_3)+D_0{\delta }_1{\theta }_3,$(2.17) where $D_0$ and $D_1$ are given in (2.5).

The proof follows from (2.9) in the light of the relations: $\text{E}\left(V_1\right)={\mu }_{\left[1,0\right]},\text{E}\left(V_2\right)={\mu }_{\left[0,1\right]}\text{and}\text{Cov}(V_1,V_2)={\mu }_{\left[1,1\right]}-{\mu }_{\left[1,0\right]}{\mu }_{\left[0,1\right]}.$

Proposition 2.6.

Let $V=(V_1,V_2)$ follow the BAZILSD with pgf (2.1). Then $U=V_1+V_2$ follows the modified AZILSD studied by Kumar and Riyaz (2013).

The proof follows from the fact that the pgf of $V_1+V_2$ is $H_U\left(t\right)=H(t,t)=\mathrm{\Lambda }{F}_{2,1}[1,1;2;({\theta }_1+{\theta }_2)t+{\theta }_3{t}^2+\alpha ].$

3. Recursion formulae

In this section, we develop certain recursion formulae for probabilities, raw moments and factorial moments. Let $V=(V_1,V_2)$ be a random vector with pgf (2.1). For the sake of computational simplicity, we define $\underset{\_}{u}+i=(1+i,1+i;2+i)$, for $i=0,1,2,\dots$. Now we have the following from (2.1) in which $f(m,n;\underset{\_}{u})=P(V_1=m,V_2=n)$, for $m,n\ge 0,$ (3.1) $\begin{array}{rl}H(t_1,t_2)& =\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}f(m,n;\underset{\_}{u})t_1^m{t}_2^n\\ & =\mathrm{\Lambda }{F}_{2,1}(1,1;2;{\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha ).\end{array}$(3.1) Now we obtain the following propositions.

Proposition 3.1

The probability mass function $f(m,n;\underset{\_}{u})$ of the BAZILSD satisfies the following recurrence formulae, in which ${\delta }_j$ is defined in Proposition 2.5. (3.2) $\begin{array}{rl}(m+1)f(m+1,0;\underset{\_}{u})& =D_0{\delta }_1{\theta }_{1}f(m,0;\underset{\_}{u}+1),m\ge 0,\end{array}$(3.2) (3.3) $\begin{array}{rl}(m+1)f(m+1,n;\underset{\_}{u})& =D_0{\delta }_1\left[{\theta }_{1}f(m,n;\underset{\_}{u}+1)+{\theta }_{3}f(m,n-1;\underset{\_}{u}+1)\right],m\ge 0,n\ge 1,\end{array}$(3.3) (3.4) $\begin{array}{rl}(n+1)f(0,n+1;\underset{\_}{u})& =D_0{\delta }_1{\theta }_{2}f(0,n;\underset{\_}{u}+1),n\ge 0,\end{array}$(3.4) (3.5) $\begin{array}{rl}(n+1)f(m,n+1;\underset{\_}{u})& =D_0{\delta }_1\left[{\theta }_{2}f(m,n;\underset{\_}{u}+1)+{\theta }_{3}f(m-1,n;\underset{\_}{u}+1)\right],m\ge 1,n\ge 0.\end{array}$(3.5)

Proof :

From (2.10) with $m=1$, we have the following. (3.6) $H^{(1,0)}(t_1,t_2)=\mathrm{\Lambda }{D}_0({\theta }_1+{\theta }_3{t}_2)h_1(t_1,t_2).$(3.6) On differentiating both sides of (3.1) with respect to $t_1$, we have (3.7) $\begin{array}{rl}{H}^{(1,0)}(t_1,t_2)& =\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}mf(m,n;\underset{\_}{u})t_1^{m-1}{t}_2^n\\ & =\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}(m+1)f(m+1,n;\underset{\_}{u})t_1^m{t}_2^n.\end{array}$(3.7) From (3.1), we also have the following. (3.8) $F_{2,1}(2,2;3;{\theta }_1{t}_1+{\theta }_2{t}_2+{\theta }_3{t}_1{t}_2+\alpha )={\psi }_1\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}f(m,n;\underset{\_}{u}+1)t_1^m{t}_2^n.$(3.8) Now by using (3.7) and (3.8) in (3.6) we get (3.9) $\begin{array}{rl}& \sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}(m+1)f(m+1,n;\underset{\_}{u})t_1^m{t}_2^n\\ & =D_0{\delta }_1[{\theta }_1\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}f(m,n;\underset{\_}{u}+1)t_1^m{t}_2^n+{\theta }_3\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}f(m,n;\underset{\_}{u}+1)t_1^m{t}_2^{n+1}].\end{array}$(3.9) On equating the coefficient of $t_1^m{t}_2^0$ on both sides of (3.9), we get (3.2). By equating the coefficient of $t_1^m{t}_2^n$ on both sides of (3.9), we get the relation (3.3). We omit the proof of relations (3.4) and (3.5) as it is similar to that of relations (3.2) and (3.3).

Proposition 3.2 :

Two recurrence formulae for the $(m,n)$-th raw moment ${\mu }_{m,n}\left(\underset{\_}{u}\right)$ of the BAZILSD are the following, for $m,n\ge 0$. (3.10) $\begin{array}{rl}{\mu }_{m+1,n}\left(\underset{\_}{u}\right)& =D_0{\delta }_1{\theta }_1\sum _{j=0}^m\left(\begin{array}{c}m\\ j\end{array}\right){\mu }_{m-j,n}(\underset{\_}{u}+1)+D_0{\delta }_1{\theta }_3\sum _{j=0}^m\sum _{k=0}^n\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right){\mu }_{m-j,n-k}(\underset{\_}{u}+1),\end{array}$(3.10) (3.11) $\begin{array}{rl}{\mu }_{m+1,n+1}\left(\underset{\_}{u}\right)& =D_0{\delta }_1{\theta }_2\sum _{k=0}^n\left(\begin{array}{c}n\\ j\end{array}\right){\mu }_{m,n-k}(\underset{\_}{u}+1)+D_0{\delta }_1{\theta }_3\sum _{j=0}^m\sum _{k=0}^n\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right){\mu }_{m-j,n-k}(\underset{\_}{u}+1).\end{array}$(3.11)

Proof :

The characteristic function $\phi (t_1,t_2)$ of the BAZILSD with pgf (2.1) is the following. For $(t_1,t_2)$ in ${\mathbb{R}}^2$ and $\mathrm{i}=\sqrt{-1}$, (3.12) $\begin{array}{rl}\phi (t_1,t_2)& =H({\mathrm{e}}^{\mathrm{i}{t}_1},{\mathrm{e}}^{\mathrm{i}{t}_2})\\ & =\mathrm{\Lambda }{F}_{2,1}[1;1;2;\gamma (\underset{\_}{t};\underset{\_}{\theta }\left)\right]\\ & =\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{m,n}\left(\underset{\_}{u}\right)\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!},\end{array}$(3.12) where $\gamma (\underset{\_}{t};\underset{\_}{\theta })=\gamma (t_1,t_2;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )$ $={\theta }_1{\mathrm{e}}^{\mathrm{i}{t}_1}+{\theta }_2{\mathrm{e}}^{\mathrm{i}{t}_2}+{\theta }_3{\mathrm{e}}^{\mathrm{i}(t_1+t_2)}+\alpha$.

On differentiating (3.12) with respect to $t_1$ we get, (3.13) $D_0\mathrm{\Lambda }{F}_{2,1}[2;2;3;\gamma (\underset{\_}{t};\underset{\_}{\theta }\left)\right]\left\{\mathrm{i}\right({\theta }_1+{\theta }_3{\mathrm{e}}^{\mathrm{i}{t}_2}\left){\mathrm{e}}^{\mathrm{i}{t}_1}\right\}=\sum _{m=1}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}\mathrm{i}{\mu }_{m,n}\left(\underset{\_}{u}\right)\frac{{\left(\mathrm{i}{t}_1\right)}^{m-1}{\left(\mathrm{i}{t}_2\right)}^n}{(m-1)!n!}.$(3.13) In the light of (3.12), we have the following from (3.13). $D_0{\delta }_1{\theta }_1\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!}{\mathrm{e}}^{\mathrm{i}{t}_1}+D_0{\delta }_1{\theta }_3\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!}{\mathrm{e}}^{\mathrm{i}{t}_1}{\mathrm{e}}^{\mathrm{i}{t}_2}=\sum _{m=1}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{m,n}\left(\underset{\_}{u}\right)\frac{{\left(\mathrm{i}{t}_1\right)}^{m-1}{\left(\mathrm{i}{t}_2\right)}^n}{(m-1)!n!}.$Now, on expanding exponential functions, rearranging the term and by using standard properties of double sum we obtain the following. (3.14) $\begin{array}{rl}& \sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{m+1,n}\left(\underset{\_}{u}\right)\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!}\\ & =D_0{\delta }_1\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!}\left[{\theta }_1\sum _{j=0}^m\left(\begin{array}{c}m\\ j\end{array}\right){\mu }_{m-j,n}\right(\underset{\_}{u}+1)+{\theta }_3\sum _{j=0}^m\sum _{k=0}^n\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}n\\ k\end{array}\right){\mu }_{m-j,n-k}(\underset{\_}{u}+1)].\end{array}$(3.14) On equating coefficients of $\frac{{\left(\mathrm{i}{t}_1\right)}^m{\left(\mathrm{i}{t}_2\right)}^n}{m!n!}$ on both sides of (3.14), we get the relation (3.10). A similar procedure will give (3.11).

Proposition 3.3

: The $(m,n)$-th order factorial moment ${\mu }_{\left[m,n\right]}\left(\underset{\_}{u}\right)$ of the BAZILSD satisfies the following recurrence formulae, for $m,n\ge 0$, in which ${\mu }_{\left[0,0\right]}\left(\underset{\_}{u}\right)=1$. (3.15) $\begin{array}{rl}{\mu }_{\left[m+1,n\right]}\left(\underset{\_}{u}\right)& =D_0{\delta }_1({\theta }_1+{\theta }_3){\mu }_{\left[m,n\right]}(\underset{\_}{u}+1)+D_0{\delta }_1{\theta }_{3}n{\mu }_{\left[m,n-1\right]}(\underset{\_}{u}+1),\end{array}$(3.15) (3.16) $\begin{array}{rl}{\mu }_{\left[m,n+1\right]}\left(\underset{\_}{u}\right)& =D_0{\delta }_1({\theta }_2+{\theta }_3){\mu }_{\left[m,n\right]}(\underset{\_}{u}+1)+D_0{\delta }_1{\theta }_{3}m{\mu }_{\left[m-1,n\right]}(\underset{\_}{u}+1).\end{array}$(3.16)

Proof:

Let $V=(V_1,V_2)$ be a random vector having the BAZILSD with pgf $H(t_1,t_2)$ as given in (3.1). Then the factorial moment generating function $F(t_1,t_2)$ of the BAZILSD is (3.17) $\begin{array}{rl}F(t_1,t_2)& =H(1+t_1,1+t_2)\\ & =\mathrm{\Lambda }{F}_{2,1}[1,1;2;\eta (\underset{\_}{t};\underset{\_}{\theta }\left)\right]\\ & =\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{\left[m,n\right]}\left(\underset{\_}{u}\right)\frac{t_1^m{t}_2^n}{m!n!},\end{array}$(3.17) where $\eta (\underset{\_}{t};\underset{\_}{\theta })=\eta (t_1,t_2;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )$ $={\theta }_1+{\theta }_2+{\theta }_3+({\theta }_1+{\theta }_3)t_1+({\theta }_2+{\theta }_3)t_2+{\theta }_3{t}_1{t}_3+\alpha$.

On differentiating (3.16) with respect to $t_1$, we get $\frac{\mathrm{\partial }F(t_1,t_2)}{\mathrm{\partial }{t}_1}=\left[({\theta }_1+{\theta }_3)+{\theta }_3{t}_2\right]D_0{F}_{2,1}[2,2;3;\eta (\underset{\_}{t};\underset{\_}{\theta }\left)\right].$In the light of (3.17), we can write this as (3.18) $\begin{array}{rl}& \sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{\left[m+1,n\right]}\left(\underset{\_}{u}\right)\frac{t_1^m{t}_2^n}{m!n!}\\ & =D_0{\delta }_1\left[\right({\theta }_1+{\theta }_3)\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{\left[m,n\right]}(\underset{\_}{u}+1)\frac{t_1^m{t}_2^n}{m!n!}+{\theta }_3\sum _{m=0}^{\mathrm{\infty }}\sum _{n=0}^{\mathrm{\infty }}{\mu }_{\left[m,n\right]}(\underset{\_}{u}+1)\frac{t_1^m{t}_2^{n+1}}{m!n!}].\end{array}$(3.18) Equating the coefficient of $\frac{t_1^m{t}_2^n}{m!n!}$ on both sides of (3.18), we get (3.15). Similar procedures will lead to (3.16).

4. Estimation and testing

In this section, we discuss the estimation of the parameters ${\theta }_1$, ${\theta }_2$, ${\theta }_3$ and $\alpha$ of the BAZILSD by the method of method maximum likelihood and construct certain test procedures for testing the significance of the additional parameter $\alpha$ of the BAZILSD.

4.1. Maximum likelihood estimation

Let $a(m,n)$ be the frequency of the $(m,n)$-th cell of a bivariate data. Let $y$ be the highest value of $m$ observed and $z$ be the highest value of $n$ observed. Then the likelihood function of the sample is (4.1) $L=\prod _{m=0}^y\prod _{n=0}^z[f(m,n){]}^{a(m,n)},$(4.1) where $f(m,n)$ is the pmf of the BAZILSD as given in (2.8). Taking logarithm on both sides of (4.1), we get (4.2) $\log L=\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\left[\log \mathrm{\Lambda }+\log \mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )\right],$(4.2) where $\mathrm{\Lambda }$ is given in (1.9), $\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )=\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{{}^{\ast }}\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!},$and $D_r^{\ast }$ is defined in Proposition 2.4.

Let ${\hat{\theta }}_1$, ${\hat{\theta }}_2$, ${\hat{\theta }}_3$ and $\hat{\alpha }$ denote the maximum likelihood estimators of the parameters ${\theta }_1$, ${\theta }_2$, ${\theta }_3$ and $\alpha$ of the BAZILSD. On differentiating (4.2), partially with respect to the parameters ${\theta }_1$, ${\theta }_2$, ${\theta }_3$ and $\alpha$, respectively, and equating to zero, we get the following likelihood equations, in which $\begin{array}{rl}\mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )& =\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\mathrm{\nabla }\log \mathrm{\Lambda }\\ & =\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\mathrm{\nabla }\log R_0^{-1}\left(\theta \right)\\ & =\sum _{m=0}^y\sum _{n=0}^z-a(m,n)R_0^{-2}\left(\theta \right)\mathrm{\nabla }{R}_0\left(\theta \right)\\ & =\sum _{m=0}^y\sum _{n=0}^z-a(m,n)R_0^{-2}\left(\theta \right)D_0{R}_1\left(\theta \right),\end{array}$in the light of $\mathrm{\nabla }{R}_j\left(\theta \right)=D_j{R}_{j+1}\left(\theta \right)$, where $D_j$ and $R_j\left(\theta \right)$ are defined in (2.5) and (1.8), respectively. (4.3) (4.4) $\begin{array}{rl}& \mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}=0,\\ & \mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}=0,\end{array}$(4.3) (4.4) (4.5) $\begin{array}{rl}& \mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}=0,\end{array}$(4.5) and (4.6) $\begin{array}{rl}& \mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}=0.\end{array}$(4.6) Now on solving these likelihood equations (4.3)–(4.6) by using some mathematical software such as MATHLAB, MATHCAD, MATHEMATICA, etc., one can obtain the maximum likelihood estimators of the parameters ${\theta }_1$, ${\theta }_2$, ${\theta }_3$ and $\alpha$.

4.2. Testing of the hypothesis

For testing the hypothesis $H_0:\alpha =0$ against the alternative hypothesis $H_1:\alpha \ne 0$, we construct the generalized likelihood ratio test (GLRT) and Rao’s efficient score test (REST) as follows.

In case of (GLRT), the test statistic is (4.7) $-2\log \lambda =2\left[\log L\right(\underset{\_}{\hat{\mathrm{\Omega }}};x)-\log L(\underset{\_}{{\hat{\mathrm{\Omega }}}^{\ast }};x\left)\right],$(4.7) where $\underset{\_}{\hat{\mathrm{\Omega }}}$ is the maximum likelihood estimator of $\underset{\_}{\mathrm{\Omega }}=({\theta }_1,{\theta }_2,\alpha )$ with no restrictions, and $\underset{\_}{{\hat{\mathrm{\Omega }}}^{\ast }}$ is the maximum likelihood estimator of $\underset{\_}{\mathrm{\Omega }}$ when $\alpha =0$. The test statistic $-2\log \lambda$ given in (4.7) is asymptotically distributed as Chi-square with one degree of freedom. For details, see Rao (1973).

In case of (REST), the following test statistic can be used. (4.8) $S=T^{\mathrm{\prime }}{\varphi }^{-1}T,$(4.8) where $T^{\mathrm{\prime }}=(T_1,T{}_2,T_3,T_4)$ and $\varphi ={\left(I_{rs}\right)}_{4\times 4}$ are the Fisher information matrices in which $T_i$ and $I_{rs}$ for $i=1,2,3,4$ and $r,s=1,2,3,4$ are as given in the Appendix. The test statistic given in (4.8) follows Chi-square distribution with one degree of freedom (see Rao, 1973).

5. Applications

For numerical applications, we consider two real-life data sets of which the first data set is from MitchelL and Paulson (1981), which consists of the number of aborts by 109 aircrafts in two consecutive six months of one year period and the second data set, taken from Partrat (1993), is the yearly frequency of hurricanes affecting tropical cyclones in two zones belonging to the North Atlantic coastal states in the USA. We have fitted the BZILSD, the BAZILSD and the bivariate Poisson distribution (BPD) to these data sets by the method of the maximum likelihood estimates of the parameter of the models. For the first data set, the maximum likelihood estimates (MLE_S) of the parameters in case of the BZILSD are ${\hat{\theta }}_1=0.75$, ${\hat{\theta }}_2=0.17$ and ${\hat{\theta }}_3=0.01$, those in case of the BAZILSD are ${\hat{\theta }}_1=0.65$, ${\hat{\theta }}_2=0.23$, ${\hat{\theta }}_3=0.04$ and $\hat{\alpha }=0.02$, and those in case of the BPD are ${\hat{\lambda }}_1=0.67$ ${\hat{\lambda }}_2=0.47$and ${\hat{\lambda }}_3=\mathrm{0.01.}$ For the second data set, the MLE_S of the parameters in case of the BZILSD are ${\hat{\theta }}_1=0.55$, ${\hat{\theta }}_2=0.36$ and ${\hat{\theta }}_3$ = 0.02, those in case of the BAZILSD are ${\hat{\theta }}_1=0.35$, ${\hat{\theta }}_2=0.31$, ${\hat{\theta }}_3=0.04$ and $\hat{\alpha }=0.01$, and those in case of the BPD are ${\hat{\lambda }}_1=0.62,{\hat{\lambda }}_2=0.61$ and ${\hat{\lambda }}_3=\mathrm{0.01.}$ The computed values of the expected frequencies of the BZILSD, the BAZILSD and the BPD are all presented in the Tables 1 and 2.

Table 1. Observed frequencies and computed values of expected frequencies of the BZILSD, the BAZILSD and the BPD by method of maximum likelihood for the first data set.

	0	1	2	3	4	∑
0	34	20	4	6	4	68
	34.99	12.03	7.15	4.66	3.38	65.21
	36.58	18.01	5.01	4.13	3.49	67.22
	34.51	16.22	3.81	0.59	0.07	55.20
1	17	7	0	0	0	24
	12.29	10.16	1.56	0.51	0.21	24.73
	16.53	6.01	1.23	0.12	0.03	23.92
	23.12	11.21	2.72	0.44	0.05	37.54
2	6	4	1	0	0	11
	7.56	2.17	1.01	0.25	0.01	11
	5.56	3.38	1.89	0.14	0.01	10.98
	7.75	3.87	0.97	0.16	0.02	12.77
3	0	4	2	2	0	4
	3.15	1.12	0.91	0.21	0.01	5.4
	3.11	1.70	0.56	0.09	0.01	5.53
	1.73	0.89	0.23	0.04	0.005	2.89
4	0	0	0	1	0	0
	1.40	0.71	0.11	0.09	0.01	2.32
	1.03	0.05	0.03	0.02	0.01	1.14
	0.29	0.15	0.04	0.007	0.0009	0.50
5	2	0	0	0	0	2
	0.15	0.08	0.07	0.03	0.01	0.34
	0.09	0.05	0.04	0.02	0.01	0.21
	0.04	0.02	0.005	0.001	0.04	0.11
∑	59	35	5	6	4	109
	62.54	26.27	10.81	5.75	3.63	109
	62.90	29.60	8.76	4.52	3.56	109
	67.44	32.36	7.78	1.24	0.18	109

Table 2. Observed frequencies and computed values of expected frequencies of the BZILSD, the BAZILSD and the BPD by method of maximum likelihood for the second data set.

	0	1	2	3	∑
0	27	9	3	2	41
	35.52	8.01	4.61	2.52	50.06
	28.01	7.98	3.93	1.23	41.15
	26.91	16.42	5.01	1.02	49.36
1	24	13	1	0	38
	16.01	11.21	0.91	0.11	28.24
	23.12	14.21	1.23	0.16	38.72
	16.69	10.45	3.27	0.68	31.09
2	8	2	1	0	11
	5.61	4.31	1.76	0.51	12.19
	7.62	3.02	0.76	0.28	11.68
	5.17	3.32	1.07	0.23	9.79
3	1	0	2	0	4
	1.04	1.12	0.91	0.21	5.4
	1.28	1.7	0.56	0.09	5.53
	1.07	0.7	0.23	0.76	2.76
∑	60	24	7	2	93
	58.42	24.04	7.39	3.15	93
	59.79	25.52	6.01	1.68	93
	49.84	30.89	9.58	2.69	93

(In each cell, the first row represents the observed frequency, the second row represents theoretical frequency of the BZILSD, the third row represents theoretical frequency of BAZILSD and the last row represents theoretical frequency of BPD).

(In each cell, the first row represents the observed frequency, the second row represents theoretical frequency of the BZILSD, the third row represents theoretical frequency of BAZILSD and the last row represents theoretical frequency of BPD).

The goodness of fit is applied to the first data set in case of the BAZILSD in nine categories [such as (0,0), (0,1), (0,2), (0, 3 and above); (1,0), (1, 1 and above); (2,0), (2, 1 and above) and (3,0 and above)], that in case of the BZILSD in eight categories [such as (0,0), (0,1), (0,2), (0, 3 and above); (1,0), (1, 1 and above); (2, 0 and above) and (3,0 and above)] and that in case of the BPD in seven categories [such as (0,0), (0,1 and above); (1,0), (1, 1 and above); (2, 0), (2, 1 and above); (3,0 and above)]. In the second data set, in case of the BAZILSD the goodness of fit is applied in seven categories [such as (0,0), (0,1), (0, 2 and above); (1,0), (1, 1 and above); (2, 0 and above) and (3,0 and above)], that in case of the BZILSD there are seven categories [such as (0,0), (0,1), (0, 2 and above); (1,0), (1, 1 and above) and (2, 0), (2,1 and above)] and that in case of the BPD in seven categories [such as (0,0), (0,1), (0, 2 and above); (1,0), (1, 1 and above); (2, 0), (2, 1 and above)]. The computed values of the Chi-square statistic and $P$ in case of both the models – BZILSD, BAZILSD and BPD for data set 1 and data set 2 are all presented in Table 3. Based on the values of Chi-square statistic and $P$, it can be observed that BAZILSD gives a better fit to both data sets compared to the existing models – the BZILSD and the BPD.

Table 3. The computed Chi-square value and $P$ value while fitting the models – BZILSD, BAZILSD and BPD for the Data set 1 and Data set 2.

Data set/models	Chi-square value	Degrees of freedom	P-value
Data set 1
BZILSD	12.28	4	.02
BAZILSD	1.56	4	.82
BPD	17.96	3	<.0001
Data set 2
BZILSD	8.16	3	.04
BAZILSD	0.43	2	0.81
BPD	8.55	3	0.04

Table 4 contains the computed values of $\log L(\underset{\_}{\hat{\mathrm{\Omega }}};x),\log L(\underset{\_}{{\hat{\mathrm{\Omega }}}^{\ast }};x)$ and the GLRT statistic for the BAZILSD in case of for both the data sets. We have also computed the values of $S$ based on (4.8) for the BAZILSD in the case of first data set as $S_1$ and for the BAZILSD in the case of second data set $S_2$ as given below. $\begin{array}{rl}{S}_1& =(-1.583.287.8212.57)\left[\begin{array}{cccc}0.08& -0.04& -0.05& 0.01\\ -0.04& 0.06& 0.01& -0.04\\ -0.05& 0.01& 0.06& -0.02\\ 0.01& -0.04& -0.02& 0.04\end{array}\right]\left(\begin{array}{c}-1.58\\ 3.28\\ 7.82\\ 12.57\end{array}\right)\\ & =6.26,\\ S_2& =\left(0.131.295.597.96\right)\left[\begin{array}{cccc}0.40& -0.26& -0.14& -0.03\\ -0.26& 0.18& 0.08& 0.02\\ -0.05& 0.08& 0.06& 0.008\\ 0.01& 0.02& 0.008& 0.02\end{array}\right]\left(\begin{array}{c}0.13\\ 1.29\\ 5.59\\ 7.96\end{array}\right)\\ & =\mathrm{4.98.}\end{array}$Since the critical value for the test at 5% level of significance and one degree of freedom is 3.84, the null hypothesis that $H_0:\alpha =0$ is rejected in both the above cases in respect of GLRT and REST.

Table 4. The computed the values of $\log L(\underset{\_}{\hat{\mathrm{\Omega }}};x),\log L(\underset{\_}{{\hat{\mathrm{\Omega }}}^{\ast }};x)$ and the generalized likelihood ratio test statistic under $H_0$.

	$\log L(\underset{\_}{{\hat{\mathrm{\Omega }}}^{\ast }};x)$	$\log L(\underset{\_}{\hat{\mathrm{\Omega }}};x)$	Test statistic
Data set 1	−154.03	−146.16	15.74
Data set 2	−98.42	−95.49	5.86

6. Simulation

It is quite difficult to examine the theoretical performance of the estimators of different parameters of the BAZILSD obtained by the method of maximum likelihood. So we have attempted a simulation study for assessing the performance of the estimators. We have simulated three data sets of sample size 150, 300 and 600 in both the positively correlated and negatively correlated situations of the BAZILSD by using Markov chain Monte Carlo (MCMC) procedure, and considered 200 replications in each case. We have considered the following two sets of parameters: (i) ${\theta }_1=0.4361$, ${\theta }_2=0.2679$, ${\theta }_3=0.1905$, $\alpha =0.0110$ (positively correlated) and (ii) ${\theta }_1=0.0847$, ${\theta }_2=0.0439$, ${\theta }_3=0.0216$, $\alpha =0.0112$ (negatively correlated) as initial values of the parameters while simulating the data sets. The computed values of the bias and standard errors in case of each of the estimators are given Table 5. From Table 5, it can be observed that both the bias and standard errors of the estimators of the parameters are in decreasing order as the sample size increases.

Table 5. Bias and standard errors (within parenthesis) of the estimators of the parameters ${\theta }_1$,${\theta }_2$, ${\theta }_3$ and $\alpha$ of the BAZILSD for the simulated data sets.

Parameters set	Sample size	${\hat{\theta }}_1$	${\hat{\theta }}_2$	${\hat{\theta }}_3$	$\hat{\alpha }$
Set (i) ${\theta }_1=0.4361$, ${\theta }_2=0.2679$, ${\theta }_3=0.1905$, $\alpha =0.0110$	150	0.0889 (0.1446)	0.0744 (0.1533)	0.0613 (0.1308)	0.0503 (0.1153)
	300	0.0598 (0.1281)	0.0506 (0.1273)	0.0432 (0.1005)	0.0297 (0.0927)
	600	0.0288 (0.1049)	0.0118 (0.0949)	0.0232 (0.0917)	0.0074 (0.0693)
Set (ii) ${\theta }_1=0.0847$, ${\theta }_2=0.0439$, ${\theta }_3=0.0216$, $\alpha =0.0112$	150	0.0466 (0.1456)	0.0263 (0.1122)	0.0113 (0.1204)	0.0121 (0.0917)
	300	0.0284 (0.1136)	0.0176 (0.0927)	0.0099 (0.0866)	0.0071 (0.0721)
	600	0.0084 (0.1030)	0.0097 (0.0843)	0.0044 (0.0748)	0.0026 (0.0583)

Acknowledgements

The authors are grateful to the Editor-in-Chief and the anonymous Referees for their valuable comments on an earlier version which helped to improve the quality of this article.

Appendix

The entries of $T^{\mathrm{\prime }}$ for the computations of the test statistic in case of REST are as given below. $\begin{array}{rl}{T}_1& =\frac{1}{\sqrt{n}}\frac{\mathrm{\partial }\log L}{\mathrm{\partial }{\theta }_1}\\ & =\frac{1}{\sqrt{n}}\left(\mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\right),\\ T_2& =\frac{1}{\sqrt{n}}\frac{\mathrm{\partial }\log L}{\mathrm{\partial }{\theta }_2}\\ & =\frac{1}{\sqrt{n}}\left(\mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\right),\\ T_3& =\frac{1}{\sqrt{n}}\frac{\mathrm{\partial }\log L}{\mathrm{\partial }{\theta }_2}\\ & =\frac{1}{\sqrt{n}}\left(\mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\right)\end{array}$and $\begin{array}{rl}{T}_4& =\frac{1}{\sqrt{n}}\frac{\mathrm{\partial }\log L}{\mathrm{\partial }{\theta }_4}\\ & =\frac{1}{\sqrt{n}}\left[\mathrm{\Phi }\right({\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]\\ & +\frac{1}{\sqrt{n}}\left(\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\right)\end{array}$in which $\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )$ and $\mathrm{\Phi }({\theta }_1,{\theta }_2,{\theta }_3,\alpha )$are defined in Equations (4.2) and (4.3).

The entries of $I_{rs}$ for the computations of the test statistic in case of REST are as given below. For $r,s=1,2,3$ and 4, $I_{rs}$’s are given below in which $\begin{array}{rl}\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )& =\sum _{m=0}^y\sum _{n=0}^z-a(m,n)[D_0{D}_2{R}_0^{-1}\left(\theta \right)R_2\left(\theta \right)-D_0^2\left[R_0^{-1}\right(\theta ){]}^2\left[R_1\right(\theta ){]}^2.\end{array}$ $\begin{array}{rl}{I}_{11}& =\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_1^2}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-2}}{(m-r-2)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{{\left[\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}\right]}^2}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2},\end{array}$ $\begin{array}{rl}{I}_{12}& =I_{21}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_1\mathrm{\partial }{\theta }_2}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!},\end{array}$ $\begin{array}{rl}{I}_{13}& =I_{31}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_1\mathrm{\partial }{\theta }_3}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!},\end{array}$ $\begin{array}{rl}{I}_{14}& =I_{41}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_1\mathrm{\partial }\alpha }\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r-1}}{(m-r-1)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)},\end{array}$ $\begin{array}{rl}{I}_{22}& =\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_2^2}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-2}}{(n-r-2)!}\frac{{\theta }_3^r}{r!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{{\left[\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}\right]}^2}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2},\end{array}$ $\begin{array}{rl}{I}_{23}& =I_{32}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_2\mathrm{\partial }{\theta }_3}=\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )\\ & +\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!},\end{array}$ $\begin{array}{rl}{I}_{24}& =I_{42}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_2\mathrm{\partial }\alpha }=\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )\\ & +\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r-1}}{(n-r-1)!}\frac{{\theta }_3^r}{r!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)},\end{array}$ $\begin{array}{rl}{I}_{33}& =\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_3^2}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )+\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-2}}{(r-2)!}}{\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{{\left[\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}\right]}^2}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2},\end{array}$ $\begin{array}{rl}{I}_{34}& =I_{43}=\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_3\mathrm{\partial }\alpha }=\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )\\ & +\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{(m+n-r+2)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{\text{β}}_{m+n-r}\left(\alpha \right)D^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^{r-1}}{(r-1)!}}{{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}\\ & \times \frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+1}\left(\alpha \right){D_r}^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)}\end{array}$ and $\begin{array}{rl}{I}_{44}& =\frac{\mathrm{\partial }{\log }^{2}L}{\mathrm{\partial }{\theta }_3^2}\\ & =\eta ({\theta }_1,{\theta }_2,{\theta }_3,\alpha )\\ & +\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{(m+n-r+2)}^2{\text{β}}_{m+n-r+2}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}}{(m+n-r+2)(m+n-r+3)\mathrm{\Omega }(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha )}\\ & -\sum _{m=0}^y\sum _{n=0}^{z}a(m,n)\frac{{\left(\sum _{r=0}^{min(m,n)}{(m+n-r+1)}^2{\text{β}}_{m+n-r+2}\left(\alpha \right)D_r^{\ast }\frac{{\theta }_1^{m-r}}{(m-r)!}\frac{{\theta }_2^{n-r}}{(n-r)!}\frac{{\theta }_3^r}{r!}\right)}^2}{{(m+n-r+2)}^2{\left[\mathrm{\Omega }\right(m,n;{\theta }_1,{\theta }_2,{\theta }_3,\alpha \left)\right]}^2}.\end{array}$