Bivariate exponentiated discrete Weibull distribution: statistical properties, estimation, simulation and applications

In this paper, a new bivariate discrete distribution is defined and studied in-detail, in the so-called the bivariate exponentiated discrete Weibull distribution. Several of its statistical properties including the joint cumulative distribution function, joint probability mass function, joint hazard rate function, joint moment generating function, mathematical expectation and reliability function for stress–strength model are derived. Its marginals are exponentiated discrete Weibull distributions. Hence, these marginals can be used to analyze the hazard rates in the discrete cases. The model parameters are estimated using the maximum likelihood method. Simulation study is performed to discuss the bias and mean square error of the estimators.
Finally, two real data sets are analyzed to illustrate the flexibility of the proposed model.



Introduction
Weibull (W) distribution is one of the most important and well-recognized continuous probability models in research and also in teaching. It has become important because if you have right-skewed, left-skewed or symmetric data, you can use this distribution to model it. Moreover, the hazard rate of its can be constant, increasing or decreasing. That flexibility of W distribution and its generalizations has made many researchers using it into their data analysis in different fields such as medicine, pharmacy, engineering, reliability, industry, social sciences, economics and environmental. See for example, Mudholkar and Srivastara [41], Bebbington et al. [4], Sarhan and Apaloo [49], El-Gohary et al. [12,13], El-Bassiouny et al. [9][10][11], El-Morshedy et al. [23,30], Eliwa et al. [19,21,2], El-Morshedy and Eliwa [22], among others.

The BEDsW distribution
Recently, Nekoukhou and Bidram [45] proposed a new three parameters distribution called the exponentiated discrete Weibull (EDsW) distribution. The CDF of the EDsW distribution is given by where a; b [ 0, 0\p\1 and N ¼ f0; 1; 2; . . .g. The corresponding PMF of Equation (1) can be written as For integer values of b, Eq. (2) can be represented as x 2 N : ð3Þ Table 1 presents some discrete distributions which can be obtained as special cases from the EDsW distribution. Suppose that V i ; i ¼ 1; 2; 3 are three independently distributed random variables which V i $ EDsWða; p; b i Þ. If X 1 ¼ maxfV 1 ; V 3 g and X 2 ¼ maxfV 2 ; V 3 g, then the bivariate vector X ¼ ðX 1 ; X 2 Þ has the BEDsW distribution with parameter vector X ¼ ða; p; b 1 ; b 2 ; b 3 Þ. The joint CDF of X is given by where can be written as The corresponding marginal probability mass function (PMF) to Equation (5) can be proposed as The joint PMF of the bivariate vector X can be easily obtained by using the following relation Thus, the corresponding joint PMF to Eq. (4) can be written as Mathematical Sciences with p 1 ¼ ½1 À p ðxþ1Þ a b 1 and p 2 ¼ ½1 À p x a b 1 þb 3 . Figure 1 shows the scatter plot of the joint PMF of the BEDsW distribution for various values of the parameters. As expected, the joint PMF of the BEDsW distribution can take various shapes depending on the values of its parameter vector X: Suppose ðX j1 ; X j2 Þ $ BEDsWða; p; b j1 ; b j2 ; b j3 Þ for j ¼ 1; . . .; n, and they are independently distributed. If Z 1 ¼ maxfx 11 ; . . .; x n1 g and Z 2 ¼ maxfx 12 ; . . .; x n2 g, then The joint survival function (SF) of the random vector X can be defined as Thus, the joint SF of the BEDsW distribution is given by where If X $ BEDsW X ð Þ, then the stress-strength reliability can be expressed as The joint hazard rate function (HRF) can be written as where . The scatter plot of the joint HRF of the BEDsW distribution is shown in Fig. 2.
It is observed that the joint HRF of the proposed model can take various shapes depending on the model parameters which makes this model more flexible to fit different data sets.

Statistical properties
The positive quadrant dependent (PQD) and total positivity of order two (TP2) properties Assume X $ BEDsW X ð Þ, then X 1 and X 2 are PQD for some values of x 1 and x 2 where Further, for every pair of increasing functions f X 1 ð:Þ and for all p 1 ; q 1 ; p 2 ; q 2 2 R. Assume x 11 ; x 21 ; x 12 ; x 22 2 N 0 and x 11 \x 21 \x 12 \x 22 from X $ BEDsW X ð Þ, then the joint SF of X satisfies the TP2 property for some values of x 1 and x 2 where Similarly, when

The joint probability generating function (PGF)
If Xs BEDsW(X), then the joint PGF can be expressed as where v j j\1: Hence, different moments and product moments of the BEDsW distribution can be obtained, as infinite series, using the joint PGF.
The conditional expectation (COEX) of X 1 given If Xs BEDsWðX), then the conditional PMF of X 1 j can be expressed as The conditional probability can be used in various areas, especially, in diagnostic reasoning and decision making. Table 2 lists the COEX for some specific parameter selections.

Maximum likelihood (ML) estimation
In this section, we use the method of the ML to estimate the unknown parameters a; p; b 1 ; b 2 and b 3 of the BEDsW distribution. Suppose that, we have a sample of size n, of the form ðx 11 ; x 21 Þ; ðx 12 ; x 22 Þ; . . .; ðx 1n ; x 2n Þ f g from the BEDsW distribution. We use the following notations: Based on the observations, the likelihood function is given by The log-likelihood function becomes where g 1 ðx; bÞ ¼ ½1 À p ðxþ1Þ a b À ½1 À p x a b : The ML estimation of the parameters a; p; b 1 ; b 2 and b 3 can be obtained by computing the first partial derivatives of (20) with respect to a; p; b 1 ; b 2 and b 3 , and then putting the results equal zeros. We get the likelihood equations as in the following form Table 2 The COEX for some specific parameter selections and where g 2 ðx; bÞ ¼ ½1 À p x a b lnð1 À p x a Þ; g 3 ðx; bÞ ¼ Àbx a p x a À1 ½1 À p x a bÀ1 ; g 4 ðx; bÞ ¼ Àb lnðxÞx a p x a lnðpÞ½1 À p x a bÀ1 : The ML estimation of the parameters a; p, b 1 ; b 2 and b 3 can be obtained by solving the above system of five nonlinear equations from (21) to (25). The solution of these equations is not easy to solve, so we need a numerical technique to get the ML estimators like the Newton-Raphson method.

Simulation
In this section, we estimate the bias and mean square error (MSE) for the proposed model parameters using simulations under complete sample. The population parameter is generated using software ''R'' package program. It is clear that the bias and MSE are reduced as the sample size is increased. This shows the consistency of the estimators, and therefore, the ML method is a proper for estimating the model parameters.

Data analysis
In this section, we explain the experimental importance of the BEDsW distribution using two applications to real data sets. The tested distributions are compared using some criteria, namely, the maximized log-likelihood (ÀL), Akaike information criterion (AIC) (see [29]), and Hannan-Quinn information criterion (HQIC) (see [28]).

Data set I: football data
These data are reported in Lee and Cha [39], and it represents a football match score in Italian football match (Serie A) during 1996 to 2011, between ACF Fiorentina(X 1 ) and Juventus(X 2 ). We shall compare the fits of BEDsW distribution with some competitive models like BDsE, BDsR, BDsW, bivariate Poisson with minimum operator (BPo min ), bivariate Poisson with three parameters (BPo-3P), independent bivariate Poisson (IBPo), bivariate discrete inverse exponential (BDsIE) and bivariate discrete inverse Rayleigh (BDsIR) distributions. Before trying to analyze the data by using the BEDsW distribution, we fit at first the marginals X 1 and X 2 separately and the minðX 1 ; X 2 Þ on these data. The ML estimation of the parameters a; p and b of the corresponding EDsW distribution for X 1 , X 2 and minðX 1 ;  Figure 5 shows the estimated PMF plots for the marginals X 1 , X 2 and minðX 1 ; X 2 Þ by using data set I.
From Figure 5, it is clear that EDsW distribution fits the data for the marginals. Now, we fit BEDsW distribution on these data. The ML estimators (MLEs), ÀL, AIC and HQIC values for the tested bivariate models are reported in Table 3.
From Table 3, it is clear that BEDsW distribution provides a better fit than the other tested distributions, because it has the smallest values among ÀL, AIC and HQIC.
Mathematical Sciences Figure 6 shows the profiles of the L function, which indicate that the estimators are unique. Figure 7 shows the estimated joint PMF for BEDsW, BDsW, BDsR and BDsE distributions by using data set I.
From Fig. 7, it is clear that the BEDsW model is the best among all tested models, which support the results of Table 3.

Data set II: nasal drainage severity score
These data are reported in Davis [5], and it represents the efficacy of steam inhalation in the treatment of common cold symptoms. We shall compare the fits of the BEDsW distribution with some competitive models like bivariate Poisson with four parameters (BPo-4P), IBPo, BDsE, BDsIE and BDsIR distributions. We fit at first the  Fig. 5 The estimated PMF for the marginals X 1 , X 2 and minðX 1 ; X 2 Þ by using football data set Table 3 The   Figure 8 shows the estimated PMF plots for the marginals X 1 , X 2 and minðX 1 ; X 2 Þ by using data set II. It is clear that the EDsW distribution fits the data for the marginals. Now, we fit BEDsW distribution on these data. The MLEs, ÀL, AIC and HQIC values for the tested bivariate models are reported in Table 4.  Fig. 8 The estimated PMF for the marginals X 1 , X 2 and minðX 1 ; X 2 Þ by using nasal drainage severity score

Mathematical Sciences
From Table 4, it is clear that BEDsW distribution provides a better fit than the other tested distributions. Figure 9 shows the profiles of the L function. Figure 10 shows the estimated joint PMF for BEDsW and BDsE distributions by using these data, which support the results of Table 4.

Conclusions
In this paper, we have introduced a new five parameters bivariate discrete distribution, in the so-called the bivariate exponentiated discrete Weibull distribution. Several of its statistical properties have been studied, and it is found that the marginals are positive quadrant dependent. Moreover, the joint reliability function satisfies the total positivity of order two for some values of x 1 and x 2 . The maximum likelihood method has been used to estimate the model parameters. Simulation has been performed, and it is found that the maximum likelihood method works quite well in estimating the model parameters. Finally, two real data sets have been analyzed, and it is found that the proposed model provides better fit than other well-known discrete distributions like bivariate discrete Weibull, bivariate discrete Rayleigh, bivariate discrete exponential, bivariate Poisson with minimum operator, bivariate Poisson with three parameters, independent bivariate Poisson, bivariate discrete inverse exponential and bivariate discrete inverse Rayleigh distributions.