Bivariate Discrete Generalized Exponential Distribution

TL;DR

This paper introduces a flexible bivariate discrete generalized exponential distribution, extending existing models, with properties analysis, an EM algorithm for parameter estimation, and application to real data.

Contribution

It develops a new bivariate discrete distribution with a novel EM algorithm for efficient parameter estimation, expanding the modeling toolkit for discrete data.

Findings

Distribution includes bivariate geometric as a special case

Efficient EM algorithm for maximum likelihood estimation

Successful application to real data set

Abstract

In this paper we develop a bivariate discrete generalized exponential distribution, whose marginals are discrete generalized exponential distribution as proposed by Nekoukhou, Alamatsaz and Bidram ("Discrete generalized exponential distribution of a second type", Statistics, 47, 876 - 887, 2013). It is observed that the proposed bivariate distribution is a very flexible distribution and the bivariate geometric distribution can be obtained as a special case of this distribution. The proposed distribution can be seen as a natural discrete analogue of the bivariate generalized exponential distribution proposed by Kundu and Gupta ("Bivariate generalized exponential distribution", Journal of Multivariate Analysis, 100, 581 - 593, 2009). We study different properties of this distribution and explore its dependence structures. We propose a new EM algorithm to compute the maximum likelihood estimators of the unknown parameters which can be implemented very efficiently, and discuss some inferential issues also. The analysis of one data set has been performed to show the effectiveness of the proposed model. Finally we propose some open problems and conclude the paper.

Tables

Table 1. Table 1: UEFA Champion’s League data

Obs.	ACF	Juventus	Obs.	ACF	Juventus
	Firontina			Firontina
	($X_1$)	($X_2$)		($X_1$)	($X_2$)
1	1	2	14	1	2
2	0	0	15	1	1
3	1	1	16	1	3
4	2	2	17	3	3
5	1	1	18	0	1
6	0	1	19	1	1
7	1	1	20	1	2
8	3	2	21	1	0
9	1	1	22	3	0
10	2	1	23	1	2
11	1	2	24	1	1
12	3	3	25	0	1
13	0	1	26	0	1

Table 2. Table 2: UEFA Champion’s League data

$X_1\downarrow$ $X_2\to$	0	1	2	3	Total
0	1	5	0	0	6
1	1	7	5	1	14
2	0	1	1	0	2
3	1	0	1	2	4
Total	3	13	7	3	26

Table 3. Table 3: MLEs, chi-square and associated p 𝑝 p -values for X 1 subscript 𝑋 1 X_{1} , X 2 subscript 𝑋 2 X_{2} and max { X 1 , X 2 } subscript 𝑋 1 subscript 𝑋 2 \{X_{1},X_{2}\} .

Data	$\hat{\alpha }$	$\hat{p}$	${\chi }^2$	$p$-value
$X_1$	4.6681	0.2617	3.9322	0.2689
$X_2$	8.4382	0.2311	0.0993	0.9619
max$\{X_1,X_2\}$	12.2939	0.2283	1.064	0.7857

Table 4. Table 4: Expected cell frequencies based on fitted BDGE

$X_1\downarrow$ $X_2\to$	0	1	2	3
0	1.28	2.60	1.43	0.57
1	1.32	5.46	2.50	0.99
2	0.58	1.79	2.47	0.52
3	0.21	0.66	0.43	0.78

Table 5. Table 5: Expected cell frequencies based on fitted BDGE when α 1 = α 2 subscript 𝛼 1 subscript 𝛼 2 \alpha_{1}=\alpha_{2} .

$X_1\downarrow$ $X_2\to$	0	1	2	3
0	1.33	1.99	0.92	0.33
1	1.99	6.06	2.23	0.79
2	0.92	2.22	2.52	0.44
3	0.33	0.79	0.44	0.73

Table 6. Table 6: Expected cell frequencies based on fitted bivariate Poisson with joint PMF ( 47 ).

$X_1\downarrow$ $X_2\to$	0	1	2	3
0	2.48	2.42	1.18	0.38
1	2.01	3.36	2.32	0.98
2	0.81	1.92	1.89	1.05
3	0.22	0.67	0.87	0.64

Equations

$$f_{GE}(x;\alpha ,\lambda )$$

$$F_{GE}(x;\alpha ,\lambda )$$

$$f_{DGE}(x;\alpha ,p)=P(X=x)=(1-p^{x+1}{)}^{\alpha }-(1-p^x{)}^{\alpha };forx\in {\mathcal{N}}_0=\{0,1,2,\dots ,\}.$$

$$F_{DGE}(x;\alpha,p)=P(X\leq x)=\left\{\matrix{0&\hbox{if}&x<0\cr(1-p^{[x]+1})^{\alpha}&\hbox{if}&x\geq 0.}\right.$$

$$G_X(z)$$

$$P(X=i)=(1-p^{i+1}{)}^{\alpha }-(1-p^i{)}^{\alpha }=\sum _{j=1}^{\infty }(-1{)}^{j+1}(\genfrac{}{}{0ex}{}{\alpha }{j})p^{ji}(1-p^j).$$

$$Y\sim GE(\alpha ,\lambda )⟹X=[Y]\sim DGE(\alpha ,p).$$

$$P(X_1<X_2)=\sum _{j=0}^{\infty }\left\{(1-p^{j+2}{)}^{{\alpha }_2}-(1-p^{j+1}{)}^{{\alpha }_2}\right\}(1-p^{j+1}{)}^{{\alpha }_1}\le \frac{{\alpha }_2}{{\alpha }_1+{\alpha }_2}.$$

$$P(X_1<X_2)$$

$$P(X_1<X_2)\le P(Y_1<Y_2)=\frac{{\alpha }_2}{{\alpha }_1+{\alpha }_2}.$$

$$f_{U|X=j}(u)=\frac{\alpha \lambda (1-p^{j+u}{)}^{\alpha -1}{p}^{j+u}}{(1-p^{j+1}{)}^{\alpha }-(1-p^j{)}^{\alpha }};0<u<1,$$

$$f_U(u)=\sum _{j=1}^{\infty }\alpha \lambda (1-p^{j+u}{)}^{\alpha -1}{p}^{j+u};0<u<1.$$

$$F_{BGE}(y_{1},y_{2})=\left\{\matrix{F_{GE}(y_{1}:\alpha_{1}+\alpha_{3},\lambda)F_{GE}(y_{2};\alpha_{2},\lambda)&\hbox{if}&y_{1}<y_{2}\cr F_{GE}(y_{1};\alpha_{1},\lambda)F_{GE}(y_{2}:\alpha_{2}+\alpha_{3},\lambda)&\hbox{if}&y_{1}>y_{2}\cr F_{GE}(y_{1}:\alpha_{1}+\alpha_{2}+\alpha_{3},\lambda)&\hbox{if}&y_{1}=y_{2}.\cr}\right.$$

$$f_{BGE}(y_{1},y_{2})=\left\{\matrix{f_{1}(y_{1},y_{2})&\hbox{if}&y_{1}<y_{2}\cr f_{2}(y_{1},y_{2})&\hbox{if}&y_{1}>y_{2}\cr f_{0}(y)&\hbox{if}&y_{1}=y_{2}=y,\cr}\right.$$

$$f_1(y_1,y_2)$$

$$f_2(y_1,y_2)$$

$$f_0(y)$$

$$F_{X_1,X_2}(x_1,x_2)$$

$$\displaystyle f_{X_{1},X_{2}}(x_{1},x_{2})=\left\{\begin{array}[]{ll}f_{1}(x_{1},x_{2})&\hbox{if~{}~{}~{}~{} $0\leq x_{1}<x_{2}$}\\ f_{2}(x_{1},x_{2})&\hbox{if~{}~{}~{}~{} $0\leq x_{2}<x_{1}$}\\ f_{0}(x)&\hbox{if~{}~{}~{}~{}$0\leq x_{1}=x_{2}=x$},\\ \end{array}\right.$$

$$f_1(x_1,x_2)=f_{DGE}(x_1;{\alpha }_1+{\alpha }_3,p)f_{DGE}(x_2;{\alpha }_2,p),$$

$$f_2(x_1,x_2)=f_{DGE}(x_1;{\alpha }_1,p)f_{DGE}(x_2;{\alpha }_2+{\alpha }_3,p),$$

$$f_0(x)=p_1{f}_{DGE}(x;{\alpha }_1+{\alpha }_3,p)-p_2{f}_{DGE}(x;{\alpha }_1,p),$$

$$f_{X_1,X_2}(x_1,x_2)=F_{X_1,X_2}(x_1,x_2)-F_{X_1,X_2}(x_1-1,x_2)-F_{X_1,X_2}(x_1,x_2-1)+F_{X_1,X_2}(x_1-1,x_2-1).$$

$$S_{X_1,X_2}(x_1,x_2)=1-F_{X_1}(x_1)-F_{X_2}(x_2)+F_{X_1,X_2}(x_1,x_2).$$

$$G_{X_1,X_2}(z_1,z_2)$$

$$\displaystyle f_{X_{1}|X_{2}=x_{2}}(x_{1}|x_{2})=\left\{\begin{array}[]{ll}f_{1}(x_{1}|x_{2})&\hbox{if~{}~{}~{}~{} $0\leq x_{1}<x_{2}$}\\ f_{2}(x_{1}|x_{2})&\hbox{if~{}~{}~{}~{} $0\leq x_{2}<x_{1}$}\\ f_{0}(x_{1}|x_{2})&\hbox{if~{}~{}~{}~{}$0\leq x_{1}=x_{2}=x$},\\ \end{array}\right.$$

$$f_i(x_1|x_2)=\frac{f_i(x_1,x_2)}{f_{DGE}(x_2;{\alpha }_2+{\alpha }_3)},i=1,2$$

$$f_0(x_1|x_2)=\frac{f_0(x)}{f_{DGE}(x_2;{\alpha }_2+{\alpha }_3)}.$$

$$F_{X_1|X_2\le x_2}(x_1)$$

$$F_{X_1|X_2=x_2}(x_1)$$

$$F_{X_1,X_2}(x_1,x_2)=(1-p^{x_1+1}{)}^{1-\alpha }(1-p^{x_2+1}{)}^{1-\alpha }(1-p^{z+1}{)}^{\alpha }.$$

$$F_{X_1,X_2}(x_1,x_2)\ge F_{X_1}(x_1)F_{X_2}(x_2).$$

Full text

Key Words and Phrases: Discrete bivariate model; generalized exponential distribution; maximum likelihood estimators; positive dependence; joint probability mass function; EM algorithm.

AMS 2000 Subject Classification: Primary 62F10; Secondary: 62H10

1Department of Statistics, Khansar Faculty of Mathematics and Computer Science, Khansar, Iran.

2Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur, Pin 208016, India.   e-mail: kundu@iitk.ac.in.

1 Introduction

The Generalized exponential (GE) distribution originally introduced by Gupta and Kundu [6] has received considerable attention in recent years. It is an absolutely continuous univariate distribution with several interesting properties. It has been used quite successfully as an alternative to a gamma or a Weibull distribution. Recently, Nekoukhou et al. [12] introduced a discrete generalized exponential (DGE) distribution, which can be considered as the discrete analogue of the absolutely continuous GE distribution of Gupta and Kundu [6]. The DGE distribution proposed by Nekoukhou et al. [12] is a very flexible two-parameter distribution. The probability mass function of the DGE distribution can be a decreasing or a unimodal function. Similarly, the hazard function of the DGE distribution can be an increasing, decreasing or a constant function depending on the shape parameter. Hence, the geometric distribution can be obtained as a special case of the DGE distribution. It has been used to analyze various discrete data sets, and the performances are quite satisfactory.

Discrete bivariate data also occur quite naturally in practice. For example, the number of goals scored by two competing teams or the number of insurance claims for two different causes is purely discrete in nature. Recently, Lee and Cha [10] proposed two very general methods namely (i) minimization and (ii) maximization methods, to generate a class of discrete bivariate distributions. They have discussed in details some special cases namely bivariate Poisson, bivariate geometric, bivariate negative binomial and bivariate binomial distributions. Although, the method proposed by Lee and Cha [10] can produce a very flexible class of discrete bivariate distributions, the joint probability mass function (PMF) may not be in a very convenient form in many cases. Due to this reason, developing the inference procedure of the unknown parameters becomes difficult in many cases. Another point may be mentioned that the bivariate discrete distributions proposed by Lee and Cha [10] may not have the same corresponding univariate marginals. For example, the bivariate discrete Poisson distribution proposed by them does not have Poisson marginals, which may not be very desirable.

The main aim of this paper is to consider the bivariate discrete generalized exponential (BDGE) distribution which can be obtained from three independent DGE distributions by using the maximization method as suggested by Lee and Cha [10]. It is observed that the BDGE distribution is a natural discrete analogue of the bivariate generalized exponential distribution (BGE) proposed by Kundu and Gupta [9]. The BDGE distribution is a very flexible bivariate distribution, and its joint PMF can take various shapes depending on the parameter values. The generation from a BDGE distribution is straight forward, hence the simulation experiments can be performed quite conveniently. It has some interesting physical interpretations also. In addition, its marginals are DGE distributions and the bivariate geometric distribution can be obtained as a special case of this model.

The estimation and the construction of confidence intervals of the unknown parameters play an important role in any statistical problem. The BDGE model has four parameters. The maximum likelihood estimators (MLEs) of the unknown parameters cannot be obtained in explicit forms. They can be obtained by solving a four dimensional optimization problem. Algorithms like Newton-Raphson or Gauss-Newton method may be used to solve this problem. But they have the standard problem of convergence to a local optimum rather than the global optimum and choice of the initial guesses. To avoid those problems, we propose to use an EM algorithm to compute the MLEs of the unknown parameters. We treat this problem as a missing value problem. In each E-Step we use the maximum likelihood predictor method to estimate the missing values and it avoids computation of the explicit expectation. At each M-step of the EM algorithm, the maximization of the ’pseudo’ log-likelihood function can be performed by solving one non-linear equation only. Hence, the implementation of the EM algorithm is quite simple in practice. Moreover, at the last step of the EM algorithm using the idea of Louis [11] the observed Fisher information matrix also can be obtained, and it will be used for construction of the confidence intervals of the unknown parameters. One real data set has been analyzed to see the effectiveness of the proposed model and the EM algorithm. The performances are quite satisfactory. Finally we propose some open problems.

Rest of the paper is organized as follows. In Section 2 we provide the preliminaries. The BDGE distribution is proposed and its properties are discussed in Section 3. In Section 4, we provide statistical inference procedures of the unknown parameters of a BDGE model. In Section 5 we provide the analysis of a real data set. Finally, in Section 6 we propose some open problems and conclude the paper.

2 Preliminaries

2.1 The DGE Distribution

The absolutely continuous GE distribution was proposed by Gupta and Kundu [6] as an alternative to the well known gamma and Weibull distributions. The two-parameter GE distribution has the following probability density function (PDF) and cumulative distribution function (CDF), respectively;

[TABLE]

Here $\alpha>0$ and $\lambda>0$ are the shape and the scale parameters, respectively. From now on a GE distribution with the shape parameter $\alpha$ and the scale parameter $\lambda$ will be denoted by GE$(\alpha,\lambda)$.

Recently, the DGE distribution was proposed by Nekoukhou et al. [12]. It has been defined as follows. A discrete random variable $X$ is said to have a DGE distribution with parameters $\alpha$ and $\displaystyle p$ $\displaystyle(=e^{-\lambda})$, if the probability mass function (PMF) of $X$ can be written as follows:

[TABLE]

The corresponding CDF becomes

[TABLE]

Here $[x]$ denotes the largest integer less than or equal to $x$. From now on a DGE distribution with parameters $\alpha$ and $p$ will be denoted by DGE$(\alpha,p)$. The PMF and the hazard function (HF) of a DGE distribution can take various shapes. The PMF can be a decreasing or a unimodal function, and the HF can be an increasing or a decreasing function. A DGE model is appropriate for modelling both over and under-dispersed data, since, in this model, the variance can be larger or smaller than the mean which is not the case for some of the standard classical discrete distributions. The authors did not provide the probability generating function (PGF) of DGE, which is quite useful for a discrete distribution. We provide the PGF of DGE$(\alpha,p)$ for completeness purposes. Let us consider the PGF of $\displaystyle G_{X}(z)=E(z^{X})$, for $|z|<1$. Therefore,

[TABLE]

here $\displaystyle{\alpha\choose{j}}=\alpha(\alpha-1)\ldots(\alpha-j+1)/j!$. Note that to compute $G_{X}(z)$, we have used the identity for $i=0,1,\ldots,$

[TABLE]

The following representation of a DGE random variable becomes very useful. Suppose $X\sim DGE(\alpha,p)$, then for $\displaystyle\lambda=-\ln p$,

[TABLE]

Using (6), the generation of a random sample from a DGE$(\alpha,p)$ becomes very simple. For example, first we can generate a random sample $Y$ from a GE$(\alpha,\lambda)$, and then considering $X=[Y]$, we can obtain a generated sample from DGE$(\alpha,p)$. Suppose $Y_{1}\sim$ GE$(\alpha_{1},\lambda)$, $Y_{2}\sim$ GE$(\alpha_{2},\lambda)$, $X_{1}=[Y_{1}]$, $X_{2}=[Y_{2}]$, then using (6) the following results can be easily obtained.

RESULT 1: If $Y_{1}$ and $Y_{2}$ are independently distributed then

[TABLE]

Proof: The first part (equality) of the above result follows from the definition as given below:

[TABLE]

The second part follows from the following observation: Since $X_{1}<X_{2}$ implies $Y_{1}<Y_{2}$, therefore we have

[TABLE]

RESULT 2: If $X_{1}\sim$ DGE$(\alpha_{1},p)$, $X_{2}\sim$ DGE$(\alpha_{2},p)$ and $\alpha_{1}<\alpha_{2}$, then $X_{2}>_{st}X_{1}$, i.e. $X_{2}$ is stochastically larger than $X_{1}$.

Proof: Obvious.

RESULT 3: If $X_{1}\sim$ DGE$(\alpha_{1},p)$, $\ldots$, $X_{n}\sim$ DGE$(\alpha_{n},p)$, and they are independently distributed then $\max\{X_{1},\ldots,X_{n}\}\sim$ DGE$\displaystyle\left(\sum_{i=1}^{n}\alpha_{i},p\right)$.

Proof: It follows from the CDF of the DGE distribution.

RESULT 4: Let $Y\sim$ GE$(\alpha,\lambda)$ with $\lambda=-\ln p$, $X=[Y]$ and $U=Y-X$. Then the conditional PDF of $U$ given $X=j$, for $\lambda=-\ln p$ and $j=1,2,\ldots$, is

[TABLE]

and the PDF of $U$ is

[TABLE]

Proof: It is simple, hence, the details are avoided.

2.2 The BGE Distribution

Kundu and Gupta [9] introduced the bivariate generalized exponential (BGE) distribution, whose marginals are GE distributions. The model has also some interesting physical interpretations. The joint CDF of the BGE model is as follows:

[TABLE]

The corresponding joint PDF becomes:

[TABLE]

where

[TABLE]

Kundu and Gupta [9] provided several properties and discussed inferential issues of the above mentioned model in details. For some recent work on the BGE distribution one is referred to Ashour et al. [1], Dey and Kundu [4], Dewan and Nandi [3], Genc [5] and the references cited there in.

3 The BDGE Distribution and its Properties

3.1 Definition and Interpretations

Definition: Suppose $U_{1}\sim$ DGE$(\alpha_{1},p)$, $U_{2}\sim$ DGE$(\alpha_{2},p)$ and $U_{3}\sim$ DGE$(\alpha_{3},p)$ and they are independently distributed. If $X_{1}=\max\{U_{1},U_{3}\}$ and $X_{2}=\max\{U_{2},U_{3}\}$, then we say that the bivariate vector $(X_{1},X_{2})$ has a BDGE distribution with the parameter vector $\boldsymbol{\theta}=(\alpha_{1},\alpha_{2},\alpha_{3},p)^{T}$. From now on we will denote this discrete bivariate distribution by BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$.

If $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then the joint CDF of $(X_{1},X_{2})$ for $x_{1}\in{\cal N}_{0}$, $x_{2}\in{\cal N}_{0}$ and for $z=\min\{x_{1},x_{2}\}$ is

[TABLE]

The corresponding joint PMF of $(X_{1},X_{2})$ for $x_{1},x_{2}\in{\cal N}_{0}$ is given by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $\displaystyle p_{1}=(1-p^{x+1})^{\alpha_{2}}$, $\displaystyle p_{2}=(1-p^{x})^{\alpha_{2}+\alpha_{3}}$. Note that the expressions $f_{1}(x_{1},x_{2})$, $f_{2}(x_{1},x_{2})$ and $f_{3}(x_{1},x_{2})$ for $x_{1},x_{2}\in{\cal N}_{0}$ can be easily obtained by using the relation

[TABLE]

If $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then the joint survival function (SF) of the vector $(X_{1},X_{2})$ can also be expressed in a compact form due to the relation

[TABLE]

Now we provide the joint PGF of $X_{1}$ and $X_{2}$. The joint PGF of $X_{1}$ and $X_{2}$ for $|z_{1}|<1$ and $|z_{2}|<1$, can be written as

[TABLE]

Using the joint PGF, different moments and product moments can be obtained as infinite series.

The following shock model and maintenance model interpretations can be provided for the BDGE distribution.

Shock Model: Suppose a system has two components, and it is assumed that the amount of shocks is measured in a digital (discrete) unit. Each component is subjected to individual shocks say $U_{1}$ and $U_{2}$, respectively. The system faces an overall shock $U_{3}$, which is transmitted to both the component equally, independent of their individual shocks. Therefore, the observed shocks at the two components are $X_{1}=\max\{U_{1},U_{3}\}$ and $X_{2}=\max\{U_{2},U_{3}\}$.

Maintenance Model: Suppose a system has two components and it is assumed that each component has been maintained independently and also there is an overall maintenance. Due to component maintenance, suppose the lifetime of the individual component is increased by $U_{i}$ amount and because of the overall maintenance, the lifetime of each component is increased by $U_{3}$ amount. Here, $U_{1}$, $U_{2}$ and $U_{3}$ are all measured in a discrete unit. Therefore, the increased lifetimes of the two components are $X_{1}=\max\{U_{1},U_{3}\}$ and $X_{2}=\max\{U_{2},U_{3}\}$, respectively.

3.2 Properties

RESULT 5: If $(Y_{1},Y_{2})\sim$ BGE$(\alpha_{1},\alpha_{2},\alpha_{3},\lambda)$, then $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, where $X_{1}=[Y_{1}]$, $X_{2}=[Y_{2}]$ and $p=e^{-\lambda}$.

Proof: It can be easily obtained from the joint CDF of $X_{1}$ and $X_{2}$.

The Result 5 indicates that the proposed BDGE distribution is a natural discrete version of BGE distribution. In addition, the marginals are DGE distributions. More precisely, we see that $X_{1}\sim$ DGE$(\alpha_{1}+\alpha_{3},p)$ and $X_{2}\sim$ DGE$(\alpha_{2}+\alpha_{3},p)$. The following algorithm can be used to generate a random sample from a BDGE distribution using Result 4.

Algorithm:

• •

Generate $U_{1}\sim$ GE$(\alpha_{1},\lambda)$, $U_{2}\sim$ GE$(\alpha_{2},\lambda)$, $U_{3}\sim$ GE($\alpha_{3},\lambda$) by using inverse transformation method.

• •

Obtain $Y_{1}=\max\{U_{1},U_{3}\}$ and $Y_{2}=\max\{U_{2},U_{3}\}$.

• •

$(X_{1},X_{2})$, where $X_{1}=[Y_{1}]$ and $X_{2}=[Y_{2}]$, is the desired random sample.

RESULT 6: We have the following results regarding the conditional distribution of $X_{1}$ given $X_{2}$ when $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$. The proofs are quite standard and the details are avoided.

(a) The conditional PMF of $X_{1}$ given $X_{2}=x_{2}$, say $f_{X_{1}|X_{2}=x_{2}}(x_{1}|x_{2})$, is given by

[TABLE]

where

[TABLE]

and

[TABLE]

(b) The conditional CDF of $X_{1}$ given $X_{2}\leq x_{2}$, say $F_{X_{1}|X_{2}\leq x_{2}}(x_{1})$, is given by

[TABLE]

(c) The conditional CDF of $X_{1}$ given $X_{2}=x_{2}$, say $F_{X_{1}|X_{2}=x_{2}}(x_{1})$, is given by

[TABLE]

Note that if $\alpha_{3}=\alpha$, $0<\alpha<1$, and also $\alpha_{1}=\alpha_{2}=1-\alpha$, then $(X_{1},X_{2})$ has geometric marginals. Therefore, we have a new bivariate geometric distribution with parameters $0<p<1$ and $0<\alpha<1$, and whose joint CDF is

[TABLE]

Here $x_{1},x_{2}\in{\cal N}_{0}$ and $z=\min\{x_{1},x_{2}\}$ as before. Moreover in this case $X_{1}$ and $X_{2}$ both have geometric distributions with parameter $p$.

From (12) it follows that for all values of $x_{1}\geq 0$ and $x_{2}\geq 0$,

[TABLE]

Therefore, $X_{1}$ and $X_{2}$ are positive quadrant dependent, i.e., for every pair of increasing functions $m_{1}(.)$ and $m_{2}(.)$, it follows that $Cov(m_{1}(X_{1}),m_{2}(X_{2}))\geq 0$, see for example Nelsen [13]. Moreover, in view of the fact that

[TABLE]

we see that

[TABLE]

Moreover, $X_{1}$ and $X_{2}$ are independent when $\alpha_{3}$ = 0. Therefore, $Corr\{X_{1},X_{2}\}=0$. For fixed $\alpha_{1}$ and $\alpha_{2}$

[TABLE]

where $U_{1},U_{2}$ and $U_{3}$ are same as defined in Section 3.1. Now, for any two arbitrary random variables $Y_{1}$ and $Y_{2}$ with finite second moments, let us define a new random variable

[TABLE]

It is clear that $E|g(Y_{1},Y_{2})|\leq 1$, for any two arbitrary random variables $Y_{1}$ and $Y_{2}$ with finite second moments. Therefore,

[TABLE]

Therefore, $Corr\{X_{1},X_{2}\}\rightarrow 1$. Hence, for a BDGE distribution the correlation coefficient has the range $[0,1)$.

Let us recall the following two definitions. Let $(X,Y)$ be a pair of random variables, then (a) $Y$ said to be left-tail decreasing in $X$, if and only if $P(Y\leq y|X\leq x)$ is a non-increasing function of $x$ for every $y$, and (b) $Y$ said to be stochastically increasing in $X$ if and only if $P(Y\leq y|X=x)$ is a non-increasing function of $x$ for every $y$, see Nelsen [13]. We have the following result for a BDGE distribution.

RESULT 7: Suppose $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then $X_{2}$ is left-tail decreasing in $X_{1}$ and $X_{2}$ is stochastically increasing in $X_{1}$.

Proof: Both the assertions follow from Result 6.

Suppose $(X,Y)$ is a pair of discrete random variables having support on ${\cal N}_{0}\times{\cal N}_{0}$, then it is said to have a total positivity of order two (TP2) property if the joint probability mass function $f(x,y)$ satisfies

[TABLE]

RESULT 8: Suppose $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then it has the TP2 property.

Proof: Suppose $x_{11},x_{21},x_{12},x_{22}\in{\cal N}_{0}$ and $x_{11}<x_{21}<x_{12}<x_{22}$, then observe that

[TABLE]

Similarly considering all other cases such as $x_{11}=x_{21}<x_{12}<x_{22}$, $x_{21}<x_{11}<x_{12}<x_{22}$ etc. it can be shown that it satisfies (30). Hence, the result is proved.

The following result is an extension of Result 1 to the bivariate case.

RESULT 9: If $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then

[TABLE]

Proof: It follows exactly the same way as the proof of Result 1, hence the details are avoided.

The following simple result is useful for goodness of fitting purposes. Suppose $(X_{1},X_{2})\sim$ BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$, then $\max\{X_{1},X_{2}\}\sim$ DGE$(\alpha_{1}+\alpha_{2}+\alpha_{3},p)$. The proof can be easily obtained as follows.

[TABLE]

In addition, the following result is an extension of Result 3 to the bivariate case.

RESULT 10: Suppose $(X_{i1},X_{i2})\sim$ BDGE$(\alpha_{i1},\alpha_{i2},\alpha_{i3},p)$, for $i=1,\ldots,n$, and they are independently distributed. If $\displaystyle Y_{1}=\max\{X_{11},\ldots,X_{n1}\}$ and $\displaystyle Y_{2}=\max\{X_{12},\ldots,X_{n2}\}$, then $(Y_{1},Y_{2})\sim$ BDGE$\displaystyle\left(\sum_{i=1}^{n}\alpha_{i1},\sum_{i=1}^{n}\alpha_{i2},\sum_{i=1}^{n}\alpha_{i3},p\right)$.

Proof: The proof can be easily obtained from the joint CDF and hence, the details are avoided.

4 Statistical Inference

4.1 Maximum Likelihood Estimation

In this section we consider the maximum likelihood estimation of the unknown parameters of a BDGE$(\alpha_{1},\alpha_{2},\alpha_{3},p)$ model based on a sample of size $n$. It is assumed that we have the following bivariate sample ${\cal D}=\{(x_{11},x_{21}),\ldots,(x_{1n},x_{2n})\}$. We use the following notations

[TABLE]

and $\displaystyle n_{1}=|I_{1}|$, $\displaystyle n_{2}=|I_{2}|$, $\displaystyle n_{0}=|I_{0}|$. Here $|I_{j}|$ denotes the number of elements in set $I_{j}$, for $j=0,1,2$. Based on the above data (31), the log-likelihood function can be written as

[TABLE]

Therefore, the MLEs of the unknown parameters can be obtained by maximizing (32) with respect to the unknown parameters. It can be obtained by solving four non-linear equations simultaneously. Newton-Raphson method may be used to solve these four non-linear equations.

We propose to use an EM algorithm technique to compute the MLEs of the unknown parameters, which is very specific to a discrete distribution, mainly to avoid solving four dimensional optimization problem. We treat this problem as a missing value problem as it is usually being done for any implementation of an EM algorithm. We estimate the missing values by maximum likelihood predictor method similar to Karlis [8]. The basic idea comes from the fact that if we know $\{(u_{1i},u_{2i},u_{3i});i=1,\ldots,n\}$ then the log-likelihood function of the complete observations becomes

[TABLE]

where

[TABLE]

The following result will be useful for further development.

RESULT 11: For any $0<p<1$, $g_{1}(\alpha_{1},p)$, $g_{2}(\alpha_{2},p)$ and $g_{3}(\alpha_{3},p)$ as defined in (34), (35) and (36), respectively are unimodal functions of $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$, respectively.

Proof: See in the Appendix.

Therefore, the maximization of (33) can be easily performed using profile likelihood method. For a fixed $p$ maximize $g_{1}(\alpha_{1},p)$, $g_{2}(\alpha_{2},p)$ and $g_{3}(\alpha_{3},p)$, with respect to $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$, respectively, and they are unique due to Result 11. If they are $\widehat{\alpha}_{1}(p)$, $\widehat{\alpha}_{2}(p)$ and $\widehat{\alpha}_{3}(p)$, respectively, then obtain the MLE of $p$ as $\widehat{p}$ by maximizing

[TABLE]

and the MLEs of $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ as

[TABLE]

Hence, it is clear that if we have the complete observations, then the MLEs of the unknown parameters can be obtained by solving three one dimensional optimization problems.

Therefore, to implement the EM algorithm at each E-Step, we obtain the missing $\{u_{ji}\}$ for $j=1,2,3$, say $\{(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i});i=1,\ldots,n\}$, by using the maximum likelihood predictor method, and then at the corresponding M-Step, we maximize the complete log-likelihood function to compute the estimates for the next iterate. It is assumed that at the $k$-th step the values of the unknown parameters are $\Theta^{(k)}=(\alpha_{1}^{(k)},\alpha_{2}^{(k)},\alpha_{3}^{(k)},p^{(k)})$ and the available data are ${\cal D}^{(k)}=\{(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i});i=1,\ldots,n\}$. We will provide the methodology how to compute $\Theta^{(k+1)}$ from $\Theta^{(k)}$ based on ${\cal D}^{(k)}$ by using EM algorithm.

E-Step: We will mention how to obtain $(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i})$ from $(x_{1i},x_{2i})$ by maximum likelihood prediction method as follows:

Case I:  $x_{1i}<x_{2i}$

Since $x_{1i}=\max\{u_{1i},u_{3i}\}$ and $x_{2i}=\max\{u_{2i},u_{3i}\}$, it is clear $u_{1i}<u_{2i}$ and $u_{3i}<u_{2i}$. Therefore, $\widetilde{u}_{2i}=x_{2i}$, and $(\widetilde{u}_{1i},\widetilde{u}_{3i})$ are chosen so that it maximizes the corresponding probability, i.e.

[TABLE]

Case II:  $x_{2i}<x_{1i}$

In this case similarly as before, $\widetilde{u}_{1i}=x_{1i}$, and $(\widetilde{u}_{2i},\widetilde{u}_{3i})$ are chosen so that

[TABLE]

Case III:  $x_{1i}=x_{2i}=x_{i}$

In this case $(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i})$ are chosen so that

[TABLE]

where

[TABLE]

M-Step: In this step we maximize $l_{complete}(\alpha_{1},\alpha_{2},\alpha_{3},p|{\cal D}^{(k)})$ as defined in (33) with respect to $\alpha_{1},\alpha_{2},\alpha_{3}$ and $p$ to obtain $\Theta^{(k+1)}$ by profile likelihood method as described before.

The choice of initial estimates of $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ and $p$ are important. Based on the marginal $\{x_{1i};i=1,\ldots,n\}$ we can obtain estimates of $\alpha_{1}+\alpha_{3}$ and $p$. Similarly, from $\{x_{2i};i=1,\ldots,n\}$, we can obtain estimates of $\alpha_{2}+\alpha_{3}$ and $p$, and from $\{z_{i}=\max\{x_{1i},x_{2i}\};i=1,\ldots,n\}$, we can obtain initial estimates of $\alpha_{1}+\alpha_{2}+\alpha_{3}$ and $p$. From these estimates we can obtain initial estimates of $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ and $p$.

Therefore, the EM algorithm can be implemented as follows:

Algorithm:

Initial Estimate: Get initial estimates of $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ and $p$ from the marginals and from the maximum of the two marginals. Let it be denoted by $\Theta^{(0)}=(\alpha_{1}^{(0)},\alpha_{2}^{(0)},\alpha_{3}^{(0)},p^{(0)})$. 2. 2.

E-Step: Obtain $\{(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i});i=1,\ldots,n\}$, when $\Theta=\Theta^{(k)}$. 3. 3.

M-Step: From ${\cal D}^{(k)}=\{(\widetilde{u}_{1i},\widetilde{u}_{2i},\widetilde{u}_{3i});i=1,\ldots,n\}$, using profile likelihood method obtain $\Theta^{(k+1)}$ which maximizes the complete log-likelihood function based on ${\cal D}^{(k)}$. 4. 4.

Check the convergence, if not satisfied then go back to E-Step and continue the process.

4.2 Testing of Hypotheses

In this section we discuss different testing of hypotheses problems which have some practical relevance. We mainly use the likelihood ratio test (LRT) for the different problems. In each case under the null hypothesis the MLE of any arbitrary parameter $\delta$ will be denoted by $\widetilde{\delta}$.

Test 1: We want to test the following hypothesis for unknown $\alpha>0$

[TABLE]

The problem is of interest as it tests whether all the $U_{i}$’s have the same distribution or not. Under $H_{0}$, the MLEs of $\alpha$ and $p$ can be obtained using the same EM algorithm as it has been described before. Here at the M-Step we need to maximize the function

[TABLE]

with respect to $\alpha>0$ and $0<p<1$, where $g_{1}$, $g_{2}$ and $g_{3}$ are the same as defined by (34), (35) and (36), respectively. Under $H_{0}$,

[TABLE]

Test 2: For $0<\alpha<1$, we want to test

[TABLE]

This is an important problem as it tests whether the marginals have geometric distribution or not. We can use the EM algorithm to compute the MLEs of the unknown parameters under $H_{0}$. In this case at the M-Step we need to maximize

[TABLE]

with respect to $0<\alpha<0$ and $0<p<1$, where $g_{1}$, $g_{2}$ and $g_{3}$ are same as defined (34), (35) and (36), respectively. Under $H_{0}$

[TABLE]

Test 3: Finally we want to test the following hypothesis

[TABLE]

This testing problem is of interest as it tests whether the marginals are independent or not. In this case we do not use the EM algorithm to compute the MLEs of $\alpha_{1}$, $\alpha_{2}$ and $p$ under $H_{0}$. However it is an iterative process. We obtain $\widetilde{\alpha}_{1}$, $\widetilde{\alpha}_{2}$ and $\widetilde{p}$ as follows.

[TABLE]

here $g_{1}(\cdot)$ and $g_{2}(\cdot)$ are the same defined before, and

[TABLE]

Since $\alpha_{3}$ is in the boundary, the standard result does not work. Using Theorem 3 of Self and Liang (1987), it follows that

[TABLE]

5 Data Analysis

In this section we present the analysis of a data set for illustrative purposes. The data set represents Italian Series A football match score data between ‘ACF Firontina’ ($X_{1}$) and ‘Juventus’ ($X_{2}$) during 1996 to 2011. The data set is presented in Table 1. It is presented in the contingency table form in Table 2. Before progressing further first we have fitted the DGE to $X_{1}$, $X_{2}$ and to max$\{X_{1},X_{2}\}$. The MLEs of the unknown parameters and the fitted chi-square values and the associated $p$-values are reported in Table 3. From Table 3 it is clear that the DGE fits quite well to $X_{1}$, $X_{2}$ and max$\{X_{1},X_{2}\}$.

Now we would like to fit the BDGE model to the above data set. We use the EM algorithm to compute the MLEs of the unknown parameters. We start the EM algorithm with the initial estimates as suggested in the previous section. If we denote the initial estimates of $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ and $p$ as $\alpha_{1}^{(0)}$, $\alpha_{2}^{(0)}$, $\alpha_{3}^{(0)}$ and $p^{(0)}$, respectively, then $\alpha_{1}^{(0)}$, $\alpha_{2}^{(0)}$ and $\alpha_{3}^{(0)}$ satisfy the following equations.

[TABLE]

Hence,

[TABLE]

We have taken $p^{(0)}$ = (0.2617+0.2311+0.2283)/3.0 = 0.2404. With these initial estimates we have started the EM algorithm and we stop the EM algorithm when the absolute difference between the two consecutive log-likelihood values is less than $10^{-4}$. The EM algorithm stops after 24 steps and the final estimates of the unknown parameters and the associated 95% confidence intervals are reported within brackets as follows: $\widehat{\alpha}_{1}$ = 1.2836 (0.5317,2.0355), $\widehat{\alpha}_{2}$ = 3.7705 (2.0187,5.5223), $\widehat{\alpha}_{3}$ = 1.0358 (0.4887,1.5829) and $\widehat{p}$ = 0.3410 (0.2863,0.3957). The associated log-likelihood value becomes -51.0549. The expected frequencies for each cell based on the fitted BDGE distribution are provided in Table 4.

The observed chi-square value is 15.8524 with the $p$-value greater than 0.30, for the $\chi^{2}$ distribution with 14 degrees of freedom. Hence, it implies that the BDGE provides a very good fit to the bivariate data set.

It is observed that the final estimates are quite far away from the initial estimates. So to check, whether they are actually MLEs or not we have performed a grid search on four dimensions with the range of $\alpha$ values between 0.1 to 10.0 with a grid size 0.0001, and the $p$-values between 0.01 to 0.99 with the same grid size. We obtain the global optimum values of $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ and $p$ as 1.2827, 3.7783, 1.0401 and 0.3428, respectively with the log-likelihood value as -51.0538. Therefore, it is clear that the obtained estimates using EM algorithm are very close to the true MLEs. It may be mentioned here that the grid search took more than 6 hours, where as in the same machine the EM algorithm took less than 50 seconds.

Now we would like to test the following hypothesis:

[TABLE]

It mainly indicates whether the performances of the two teams against each other are the same or not. Under the null hypothesis we obtain the MLEs as follows: $\widehat{\alpha}_{10}=\widehat{\alpha}_{20}$ = 3.3025, $\widehat{\alpha}_{30}$ = 1.1423, $\widehat{p}_{0}$ = 0.3175 and the associated log-likelihood value is -51.9978. We would like to use the likelihood ratio test, and the value of the test statistic is 2(-51.0538 + 51.9978) = 1.888. The associated $p$-value for $\chi^{2}_{1}$ degrees of freedom is less than 0.17. Hence, we cannot reject the null hypothesis. The expected frequencies for each cell based on the fitted BDGE distribution with $\alpha_{1}=\alpha_{2}$, are provided in Table 5.

The observed chi-square value is 18.1072 with the $p$-value 0.20. Hence, it implies that the BDGE with $\alpha_{1}=\alpha_{2}$ also provides a good fit to the bivariate data set. Therefore, based on the available data we make the conclusion that the performances of the two teams against each other are not significantly different.

Now we would like to test whether bivariate geometric distribution can be used or not to analyze this data set, as it has been considered in Test 2 in the previous section. We obtain the MLEs of $\alpha$ and $p$ under the null hypothesis as $\widehat{\alpha}$ = 0.5334 and $\widehat{p}$ = 0.2879. The associated log-likelihood value becomes -93.3893. It is observed that the $p$-value of the test statistic, is less than 0.001. Hence, bivariate geometric distribution cannot be used in this case.

For comparison purposes we would like to examine whether bivariate Poisson distribution provides a better fit or not to this data set. We have used the following joint PMF of a bivariate Poisson distribution with parameters $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$.

[TABLE]

The MLEs of $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are as follows: $\widehat{\lambda}_{1}$ = 0.8089, $\widehat{\lambda}_{2}$ = 0.9737, $\widehat{\lambda_{3}}$ = 0.5643, respectively, and the associated log-likelihood value is -53.3251. The expected frequencies for each cell based on the fitted bivariate Poisson distribution are provided in Table 6.

The observed chi-square value and the associated $p$-value are 21.8381 and 0.08, respectively. Therefore, it is clear that the bivariate Poisson distribution (47) does not provide a good fit to this data set.

6 Conclusions

In this paper we have introduced a bivariate discrete generalized exponential (BDGE) distribution which is a natural discrete version of the continuous bivariate generalized exponential distribution of Kundu and Gupta [9]. The proposed BDGE distribution has the marginals which are discrete generalized exponential distribution. The BDGE distribution has four parameters and it is a flexible bivariate model. We have derived several properties of the distribution. It is observed that the MLEs of the BDGE distribution cannot be obtained in explicit forms and they have to be obtained by solving four non-linear equations simultaneously. We have proposed a new EM algorithm which is very specific to this particular problem. The proposed EM algorithm avoids solving four non-linear equations. It only needs to solve four one-dimensional optimization problems. Hence, it is quite easy to implement in practice. One data set has been analyzed and it is observed that the proposed model and the EM algorithm work quite well in practice.

In this paper we have developed mainly the classical inference. It will be interesting to develop the Bayesian inference also and compare their performances. Moreover, here we have considered the bivariate model, it is important to see how it can be generalized to the multivariate case. More work is needed along these directions.

Acknowledgements:

The authors would like to thank two unknown referees and the associate editor for their constructive comments which had helped to improve the manuscript in a significant manner.

Appendix

PROOF OF RESULT 11:

First we will show that the function $h(\alpha|p)$ for $0<p<1$, and for any $j=0,1,\ldots$,

[TABLE]

is a log-concave function. By straight forward calculation it can be seen that

[TABLE]

where

[TABLE]