Bivariate Burr X Generator of Distributions: Properties and Estimation Methods with Applications to Complete and Type-II Censored Samples

: Burr proposed twelve different forms of cumulative distribution functions for modeling data. Among those twelve distribution functions is the Burr X distribution. In statistical literature, a ﬂexible family called the Burr X-G (BX-G) family is introduced. In this paper, we propose a bivariate extension of the BX-G family, in the so-called bivariate Burr X-G (BBX-G) family of distributions based on the Marshall–Olkin shock model. Important statistical properties of the BBX-G family are obtained, and a special sub-model of this bivariate family is presented. The maximum likelihood and Bayesian methods are used for estimating the bivariate family parameters based on complete and Type II censored data. A simulation study was carried out to assess the performance of the family parameters. Finally, two real data sets are analyzed to illustrate the importance and the ﬂexibility of the proposed bivariate distribution, and it is found that the proposed model provides better ﬁt than the competitive bivariate distributions.


Introduction
The Burr X (BX) model, as one of twelve models, was explored by utilizing the method of differential equation (see [1]). The random variable X is said to have the BX if its cumulative distribution function (CDF) is given by where γ > 0 is the shape parameter. This model has found many applications in many areas such as reliability study, agricultural, biological, health, the lifetime of random phenomenon and engineering, see for example, [2][3][4][5][6][7][8].

1.
Bivariate compounded distributions and families: Reference [20] obtained four bivariate extended exponential-geometric distributions from the extended exponential-geometric model introduced by [21]. Reference [22] compounded two discrete distributions and proposed bivariate geometric-Poisson distribution. Reference [23] proposed the bivariate Weibull-geometric distribution and discussed some of its properties and estimation methods. Reference [24] proposed the bivariate exponentiated generalized Weibull-Gompertz distribution. Reference [25] proposed the bivariate exponentiated modified Weibull extension distribution. Reference [26] introduced and studied complementary bivariate generalized linear failure rate-power series family of distributions.
The aim of our paper was to introduce a new bivariate family, the bivariate Burr X-G (BBX-G) family based on the Marshall-Olkin shock model (see [35]), whose marginal distributions are BX-G families. The structure of the proposed paper follows similarly to that of [32,33]. A random vector X = (X 1 , X 2 ) follows the bivariate Marshall-Olkin model if and only if there exist three independent random variables U 1 , U 2 and U 3 such that (X 1 = max(U 1 , U 3 ) and X 2 = max(U 2 , U 3 )) or (X 1 = min(U 1 , U 3 ) . Equation (9) can be obtained by substituting from Equations (5) and (8) in the relation f X i |X j (x i | f X j (x j ) , (i = j = 1, 2). The PDF and CDF marginals can be represented as a linear representation as follows and Mathematics 2020, 8,264 5 of 32 respectively, where Υ (i) 2(m+1)+l (x i ; ω) = (2(m + 1) + l)g(x i ; ω)G(x i ; ω) 2m+l+1 , and Λ (i) 2(m+1)+l is the CDF of the exponential-G (exp-G) family with power parameter 2(m + 1) + l. For more detail around exp-G family of distributions (see [36]).

The Moments, Product Moment, Covariance, Skewness and Kurtosis
The rth moment of X i , say M (r) i , can be defined as M Hence, by using Equation (10), we get where Z i,2(m+1)+l ; i = 1, 2 be a random variable having the exp-G CDF with power parameter 2(m + 1) + l. The moments of the exp-G distributions are given by [41]. Setting r = 1 in Equation (18), we get the mean of X i ; i = 1, 2. Thus, the nth central moment of X i , say L (n) i , is given by The sth incomplete moment of X i , say ϕ (s) Then, the sth incomplete moment can be expressed as follows where ϕ * (s) i 2(m+1)+l (x i ; ω) dx i . Therefore, the mean deviations of X i about the mean and median are given by ρ i = 2M (1) respectively. The sth incomplete moment has more applications in various fields, for more details, see [42]. The product moment can be expressed as follows Mathematics 2020, 8, 264 Using Equations (18) and (21) when r = 1, we get the covariance of the bivariate distribution as follows where Cov(X 1 , . Moreover, the correlation of X 1 and X 2 is the number , where 0 ≤ ρ ≤ 1 and Var(X i ) = M (2) i − M (1) i 2 ; i = 1, 2. By using [43] measures of multivariate and bivariate skewness and kurtosis, we get where

The Joint RF, Joint Reversed (Hazard) Rate Functions and Stress-Strength Reliability
Assume (X 1 , X 2 ) be two dimensional random variable with joint CDF F X 1 ,X 2 (x 1 , x 2 ) and the marginal functions are F X 1 (x 1 ) and F X 2 (x 2 ), then the joint RF can be defined as R X 1 ,X 2 (x 1 , . So, the joint RF of the BBX-G family can be expressed as follows where Reference [44] defined the bivariate hazard rate function (BHRF) as follows h X 1 ,X 2 (x 1 , . So, the BHRF of the BBX-G family can be written as follows where The marginal hazard rate functions h i (x i ); i = 1, 2 can be expressed as follows Reference [45] defined the bivariate reversed hazard rate function (BRHRF) as a scalar, given by F X 1 ,X 2 (x 1 ,x 2 ) . So, the BRHRF for the random vector (X 1 , X 2 ) can be expressed as follows where The marginal reversed hazard rate functions r i (x i ); i = 1, 2 can be expressed as follows On the other hand, the proposed bivariate family has a nice interpretation, namely, the stress-strength model does not depend on the baseline function G(x; ω), by another way P[X 1 < X 2 ] = γ 2 +γ 3 γ 1 +γ 2 +2γ 3 and P[X 2 < X 1 ] = γ 1 +γ 3 γ 1 +γ 2 +2γ 3 .

Special Case of BBX-G Family: Bivariate Burr X-Exponential Distribution with Properties
The random variable X is said to have the exponential (Ex) distribution if its CDF is given by The joint CDF of the bivariate Burr X-exponential (BBXEx) distribution can be expressed as follows where z = min(x 1 , x 2 ). By substituting from Equation (30) in Equations (5), (25) and (26), we get the joint PDF, joint RF and BHRF of the BBXEx distribution, respectively. Figures 1-3 show the surface plots of those functions for γ 1 = γ 2 = γ 3 = 0.3 and a = 0.1, 0.3 and 0.5, respectively.

Special Case of BBX-G Family: Bivariate Burr X-Exponential Distribution with Properties
The random variable X is said to have the exponential (Ex) distribution if its CDF is given by The joint CDF of the bivariate Burr X-exponential (BBXEx) distribution can be expressed as follows where z = min(x 1 , x 2 ). By substituting from Equation (30) in Equations (5), (25) and (26), we get the joint PDF, joint RF and BHRF of the BBXEx distribution, respectively. Figures 1-3 show the surface plots of those functions for γ 1 = γ 2 = γ 3 = 0.3 and a = 0.1, 0.3 and 0.5, respectively.

103
Although the baseline distribution G(x; ω) can be presented by several distributions, we choose the exponential (Ex) distribution for an example. The random variable X is said to have the Ex distribution if its CDF is given by The joint CDF of the bivariate Burr X-exponential (BBXEx) distribution can be expressed as follows where z = min(x 1 , x 2 ). By substituting from Equation (30) in Equations (5), (25) and (26), we get the joint PDF, joint RF and BHRF of the BBXEx distribution, respectively.   The marginal reversed hazard rate functions r i (x i ); i = 1, 2 can be expressed as follows On the other hand, the proposed bivariate family has a nice interpretation, namely, the stress-strength 100 model does not depend on the baseline function G(x; ω), by another way P[X 1 < X 2 ] = γ 2 +γ 3 γ 1 +γ 2 +2γ 3 101 and P[X 2 < X 1 ] = γ 1 +γ 3 γ 1 +γ 2 +2γ 3 .

103
Although the baseline distribution G(x; ω) can be presented by several distributions, we choose the exponential (Ex) distribution for an example. The random variable X is said to have the Ex distribution if its CDF is given by The joint CDF of the bivariate Burr X-exponential (BBXEx) distribution can be expressed as follows where z = min(x 1 , x 2 ). By substituting from Equation (30) in Equations (5), (25) and (26), we get the joint PDF, joint RF and BHRF of the BBXEx distribution, respectively.         It is clear that the joint density has a long left tail as compared to its right tail. Moreover, the BBXEx distribution presents different shapes for the BHRF. Furthermore, the joint RF decreases for fixed values of γ 1 , γ 2 and γ 3 with a → ∞. Thus, this model can be used to discuss several phenomena in different fields. [45] defined the local dependence function, say η(x 1 , x 2 ), in order to study the dependence between X 1 and X 2 , where For the BBXEx distribution, it can be verified that 106 η(x 1 , x 2 ) > 0, and then X 1 and X 2 are PT2. As a consequence, 107 1. The linear correlation coefficient between X 1 and X 2 is always positive.
Recall, Equations (23) and (24), the correlation, skewness and kurtosis measures of the BBXEx distribution are listed in Table 1 for (a, γ 1 , γ 2 , γ 3 ) = (1.5, 0.6, γ 2 , 1.5). In this section, we compute the maximum likelihood estimation (MLE) for the unknown parameters Θ =(ω, γ 1 , γ 2 , γ 3 ) based on complete and Type-II censored data. Suppose that (x 11 , x 21 ), (x 12 , x 22 ),..., (x 1n , x 2n ) be the observed values from the BBX-G family. We use the following notation It is clear that the joint density has a long left tail as compared to its right tail. Moreover, the BBXEx distribution presents different shapes for the BHRF. Furthermore, the joint RF decreases for fixed values of γ 1 , γ 2 and γ 3 with a → ∞. Thus, this model can be used to discuss several phenomena in different fields. [45] defined the local dependence function, say η(x 1 , x 2 ), in order to study the dependence between X 1 and X 2 , where For the BBXEx distribution, it can be verified that η(x 1 , x 2 ) > 0, and then X 1 and X 2 are PT2. As a consequence,

1.
The linear correlation coefficient between X 1 and X 2 is always positive.

2.
The conditional hazard rate of The conditional hazard rate of Recall, Equations (23) and (24), the correlation, skewness and kurtosis measures of the BBXEx distribution are listed in Table 1 for (a, γ 1 , γ 2 , γ 3 ) = (1.5, 0.6, γ 2 , 1.5). From Table 1, it is observed that the value of correlation increases with γ 2 −→ ∞ for fixed values of a, γ 1 and γ 3 . Moreover, this distribution can be used to model skewed as well as symmetric data sets. It is clear that the joint density has a long left tail as compared to its right tail. Moreover, the BBXEx distribution presents different shapes for the BHRF. Furthermore, the joint RF decreases for fixed values of γ 1 , γ 2 and γ 3 with a → ∞. Thus, this model can be used to discuss several phenomena in different fields. Reference [46] defined the local dependence function, say η(x 1 , x 2 ), in order to study the dependence between X 1 and X 2 , where

Estimation Based on Complete and Type-II Censored Samples
For the BBXEx distribution, it can be verified that η(x 1 , x 2 ) > 0, and then X 1 and X 2 are PT2. As a consequence,

1.
The linear correlation coefficient between X 1 and X 2 is always positive.

2.
The conditional hazard rate of The conditional hazard rate of X 2 |X 1 = x 1 is decreasing in x 1 .
Recall, Equations (23) and (24), the correlation, skewness and kurtosis measures of the BBXEx distribution are listed in Table 1 for (a, γ 1 , γ 2 , γ 3 ) = (1.5, 0.6, γ 2 , 1.5). From Table 1, it is observed that the value of correlation increases with γ 2 −→ ∞ for fixed values of a, γ 1 and γ 3 . Moreover, this distribution can be used to model skewed as well as symmetric data sets.

Maximum Likelihood Estimation
In this section, we compute the maximum likelihood estimation (MLE) for the unknown parameters Θ =(ω, γ 1 , γ 2 , γ 3 ) based on complete and Type-II censored data. Suppose that (x 11 , x 21 ), (x 12 , x 22 ),..., (x 1n , x 2n ) are the observed values from the BBX-G family. We use the following notation The total likelihood function for Θ based on complete data can be defined as follows Substituting Equation (5) into Equation (33), the log-likelihood function L(Θ) is given by The first partial derivatives of Equation (34) with respect to γ 1 , γ 2 , γ 3 and ω k (k = 1, 2, 3, ...) are where [U(.)] (ω) means the derivative of the function U (.) with respect to ω. By equating the Equations (35)-(38) by zeros, we get the non-linear normal equations. The likelihood function for the bivariate distribution based on Type-II censored data can be written as follows (see [47]). The log-likelihood function L * (Θ) can be expressed as follows Substituting from Equations (5) and (7) into Equation (40), and then differentiation the result equation with respect to γ 1 , γ 2 , γ 3 and ω k (k = 1, 2, 3, ...). The MLEs of the parameters can be obtained by solving the normal equations simultaneously.
Thus, the Bayesian estimators of the parameters ω and γ l under square error loss function can be calculated through the following equations as follows and respectively, where j = 1, ..., k and l = 1, 2, 3. Generally, the ratio of k + l integrals given by Equations (44) and (45) [48,49]. Similarly to acceptance-rejection sampling, the M-H algorithm consider that to each iteration of the algorithm, a candidate value can be generated from a proposal distribution. So, the candidate value is accepted according to an adequate acceptance probability. This procedure guarantees the convergence of the Markov chain for the target density. For more details regarding the implementation of M-H algorithm, the readers may refer to [50][51][52].
Regarding to the Type-II censored data, Equation (39) can be used instead of Equation (33) to get the Bayes estimates of the unknown parameters ω, γ 1 , γ 2 and γ 3 . At the end of this section, we can conclude that the advantage of using the MCMC method over the MLE method is that we can always obtain a reasonable interval estimate of the parameters by constructing the probability intervals based on empirical posterior distribution. This is often unavailable in MLE.

Percentile bootstrap confidence interval
The following algorithm shows how to calculate the percentile bootstrap confidence interval (P-BCI) for the model parameters: 1.
Compute the MLE of Θ k where k =length(Θ) for BBXEx model.

2.
Generate the bootstrap samples using Θ k to obtain the bootstrap estimate of Θ k , sayΘ b k , using the bootstrap sample.

5.
A two side 100(1 − α)% P-BCI for the unknown parameters Θ k is given by

Percentile Bootstrap-t Confidence Interval
The following algorithm shows how to calculate the percentile bootstrap-t confidence interval (B-TCI) for the model parameters: 1.
Same as steps 1 and 2 in P-BCI. 2.
k and it can be obtained using the Fisher information matrix.

5.
A two side 100(1 − α)% B-TCI for the unknown parameters Θ k is given by

Simulation Results Based on Complete Data
In this section, the MLE, Bayesian estimation (BSE) and bootstrap confidence interval (BCI) methods are used to estimate the parameters a, γ 1 , γ 2 and γ 3 of the BBXEx distribution by using different sample sizes n = [50,100,150,200,300] from N =1000 replications. The population parameters are generated using the software R package. For more details around the R package, see [50] and [51]. This study presents an assessment of the properties for both MLE and BSE in terms of bias and mean square error (MSE) as well as the BCI for the parameters. The following algorithm shows how to generate data from the BBXEx distribution.
The MLEs and BSEs as well as the BCI values are listed in Table 2 for the BBXEx distribution when (a, γ 1 , γ 2 , γ 3 ) = (5, 0.7, 0.8, 0.9) based on complete data. The magnitude of bias in general always close to zero when n grows.

3.
Based on the MSE, the performance of the BSE method is better than the MLE method.

4.
The confidence in the results increases as the sample size increases where the BCI decreases when n grows.

Simulation Results Based on Type-II Censored Samples
The following algorithm shows how to generate Type-II censored bivariate samples from the BBXEx distribution: 1.

5.
Type-II censored data are obtained by keeping the first r pairs of ordered observations (X 1 j:n , X 2 [j:n] ); i = 1, 2, ..., n and dropping the remaining n − r observations.
The MLEs and BSEs as well as the BCI values are reported in Tables 3 and 4 for the BBXEx distribution when (a, γ 1 , γ 2 , γ 3 ) = (5, 0.7, 0.8, 0.9) based on Type-II censored data for different sample sizes n = 100 and 200, respectively.  Based on the simulation results, it is clear that: 1. The biases and MSEs of both MLEs and BSEs decrease when the sampling r increases for a fixed sample size n.

2.
The MLE and BSE methods provide a fit for estimating the model parameters.

3.
The ACI, BT and BP decrease when the sampling r increases for a fixed sample size n. So, confidence in the results increases as the sample size increases where the results approaching the real average.

Real Data
In this section, we illustrate the empirical importance of the BBXEx distribution using two applications to real data. The fitted distributions are compared using some criteria, namely, the maximized log-likelihood (L), Akaike information criterion (AIC), corrected AIC (CAIC), Bayesian IC (BIC) and Hannan-Quinn IC (HQIC); in addition to the Kolmogorov-Smirnov (KS) statistic and its p-value for the marginals. For more details regarding these criteria, see [53][54][55][56].

Data Set I: Football Data
Here, consider the data obtained by [57], which represent football (soccer) data. This data describes the games where at least one kick goal scored by any team has been considered, and the home team must have scored at least one goal. This data was analyzed by several authors, see for example, [24,25,58,59]. We consider the BBXEx model to analyze this data, comparing with other famous bivariate models, such as bivariate generalized exponential (BGEx), bivariate exponential (BEx), bivariate Gumbel exponential (BGuEx), bivariate generalized linear failure rate (BGLFR), bivariate Weibull (BW), bivariate exponentiated Weibull (BEW), bivariate generalized power Weibull (BGPW) and bivariate Gompertz (BGz) distributions. Figure 4 shows that the scatter plot for data set I.

Data Set I: Football Data
Here, consider the data obtained by [56], which represent football (soccer) data. This data describes the games which at least one kick goal scored by any team have been considered, and the home team must be scored at least one goal. This data was analyzed by several authors, see for example, [57][58][59][60]. We consider the BBXEx model to analyze this data comparing with other famous bivariate models, such as bivariate generalized exponential (BGEx), bivariate exponential (BEx), bivariate Gumbel exponential (BGuEx), bivariate generalized linear failure rate (BGLFR), bivariate Weibull (BW), bivariate exponentiated Weibull (BEW), bivarite generalized power Weibull (BGPW) and bivarite Gompertz (BGz) distributions. Figure 4 shows that the scatter plot for data set I.
Here, consider the data obtained by [56], which represent football (soccer) data. This data describes the games which at least one kick goal scored by any team have been considered, and the home team must be scored at least one goal. This data was analyzed by several authors, see for example [57][58][59][60]. We consider the BBXEx model to analyze this data comparing with other famous bivariate models, such as bivariate generalized exponential (BGEx), bivariate exponential (BEx), bivariate Gumbel exponential (BGuEx), bivariate generalized linear failure rate (BGLFR), bivariate Weibull (BW) bivariate exponentiated Weibull (BEW), bivarite generalized power Weibull (BGPW) and bivarite Gompertz (BGz) distributions. Figure 4 shows that the scatter plot for data set I. Before trying to analyze the data using the BBXEx model, we fit at first the marginals X 1 , X 2 and min(X 1 , X 2 ) separately on the UEFA Champion's League data. The MLEs of the parameters (γ, a) of the corresponding Burr X-exponential (BXEx) model for X 1 , X 2 and min(X 1 , X 2 ) are (0.724, 0.013) We fit at first the marginals X 1 , X 2 and min(X 1 , X 2 ) separately on the UEFA Champion's League data. The MLEs of the parameters (γ, a) of the corresponding Burr X-exponential (BXEx) model for X 1 , X 2 and min(X 1 , X 2 ) are (0.724, 0.013), (0.445, 0.012) and (0.459, 0.014), respectively with standard error (STER) (0.137, 0.001), (0.080, 0.001) and (0.083, 0.001). The −L, KS distance and its p-value for the marginals are listed in Table 5. Table 5. The-L, KS and p-values for the marginals using data set I. It is clear that the BXEx model fits the data for the marginals. The fitted PDF, estimated CDF and PP plots displayed in Figures 5-7 which support our results in Table 5. The fitted PDF, estimated CDF and Probability-Probability (PP) plots displayed in Figures 5-7 which support our results in Table 5. We fit at first the marginals X 1 , X 2 and min(X 1 , X 2 ) separately on the UEFA Champion's League data. The MLEs of the parameters (γ, a) of the corresponding Burr X-exponential (BXEx) model for X 1 , X 2 and min(X 1 , X 2 ) are (0.724, 0.013), (0.445, 0.012) and (0.459, 0.014), respectively with standard error (STER) (0.137, 0.001), (0.080, 0.001) and (0.083, 0.001). The −L, KS distance and its p-value for the marginals are listed in Table 5. Table 5. The log-likelihood (L), Kolmogorov-Smirnov (KS) and p-values for the marginals using data set I. It is clear that the BXEx model fits the data for the marginals. The fitted PDF, estimated CDF and PP plots displayed in Figures 5-7 which support our results in Table 5. The fitted PDF, estimated CDF and probability-probability (PP) plots displayed in Figures 5-7 which support our results in Table 5. Mathematics 2020, xx, 5 18 of 32 fitted PDF, estimated CDF and Probability-Probability (PP) plots displayed in Figures 5, 6 and 7 which 214 support our results in Table 5.  Figure 5. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 6. The estimated CDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I. Figure 5. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 5. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 6. The estimated CDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 7. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set I.
From Figures 5, 6 and 7, it is quite apparent that the marginals can be used to discuss this data. Therefore, the BBXEx model may be used for this purpose. Now, we fit the BBXEx model on this data. In the enclosed Table 6, we provide the MLEs with its corresponding standard error (STER), −L, AIC, CAIC, BIC and HQIC values for tested distributions. Figure 7. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set I. Figure 5. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 5. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 6. The estimated CDF for X 1 , X 2 and min(X 1 , X 2 ) for data set I.  Figure 7. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set I.
From Figures 5, 6 and 7, it is quite apparent that the marginals can be used to discuss this data. Therefore, the BBXEx model may be used for this purpose. Now, we fit the BBXEx model on this data. In the enclosed Table 6, we provide the MLEs with its corresponding standard error (STER), −L, AIC, CAIC, BIC and HQIC values for tested distributions. Figure 7. The probability-probability (PP) plots for X 1 , X 2 and min(X 1 , X 2 ) for data set I.
From Figures 5-7, it is quite apparent that the marginals can be used to discuss this data. Therefore, the BBXEx model may be used for this purpose. Now, we fit the BBXEx model on this data. In the enclosed Table 6, we provide the MLEs with its corresponding standard error (STER), −L, AIC, CAIC, BIC and HQIC values for tested distributions. From Table 6, it is observed that, the BBXEx model provides a better fit than the other competitive models, because it has the smallest value among −L, AIC, CAIC, BIC and HQIC. The  Table 7. The results presented in Table 7 are very similar to the MLE results. Regarding to the hyper-parameter elicitation, the elicitation of the hyper-parameters will rely on the informative priors. These informative priors will be obtained from the maximum likelihood estimates for (a, γ 1 , γ 2 , γ 3 ) by equating the mean and variance with the mean and variance of the considered priors (Gamma priors). . Now, in regards to solving the above two equations, the estimated hyper-parameters are a 1 = 118.249, a 2 = 18.288, a 3 = 7.206 and a 4 = 7.206 whereas b 1 = 11367.16, b 2 = 120.0362, b 3 = 77.288 and b 4 = 110.049. For more details around credible interval algorithm, see [60,61]. The MCMC plots for data set I are displayed in Figure 8.  Figure 8. The MCMC plots for data set I using the BBXEx model. Table 8 shows some descriptive statistics for data set I utilizing the BBXEx distribution and its marginals. According to Table 8, it is clear that the bivariate data has positively skewed with platykurtic. Moreover, the correlation between the two random variables is positive and strong. Positive correlation  Table 8 shows some descriptive statistics for data set I utilizing the BBXEx distribution and its marginals. According to Table 8, it is clear that the bivariate data has positively skewed with platykurtic. Moreover, the correlation between the two random variables is positive and strong. Positive correlation is a relationship between two variables in which both variables move in tandem that is, in the same direction.
Tables 9-11 list estimation summaries for the BBXEx model and the competitive models based on Type-II censored data using data set I. Regarding Tables 9-11, it is clear that both BW and BGPW models are better than the BBXEx model in case of small values of r as seen in Table 9, whereas the BBXEx model provides better fit than other competitive models when the value of r grows as seen in Tables 10 and 11.

Data Set II: Motor Data
This data is reported in [62], and it represents the failure times of a parallel system constituted by two identical motors in days. We consider the BBXEx model to analyze the censored samples. We fit at first the marginals X 1 , X 2 and max(X 1 , X 2 ) separately on the motor data. The MLEs of the parameters (γ, a) of the BXEx model for X 1 , X 2 and min(X 1 , X 2 ) are (1.548, 0.004), (1.233, 0.003) and (1.343, 0.004), respectively with STER (0.465, 0.0003), (0.359, 0.0003) and (0.394, 0.0003). The −L, KS distance and its p-value for the marginals are reported in Table 12.  Figure 9. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 10. The estimated CDF for X 1 , X 2 and min(X 1 , X 2 ) for data set II. Figure 9. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 9. The fitted PDF for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 10. The estimated CDF for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 11. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 11. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  Figure 11. The PP plots for X 1 , X 2 and min(X 1 , X 2 ) for data set II.  From Figure 12, it is clear that the marginals have increasing HRF. Now, we fit the BBXEx model based on a complete sample. In the enclosed Table 13, we provide the MLEs with its corresponding STER, −L, AIC, CAIC, BIC and HQIC values for tested distributions.  Table 14. The results presented in Table 14 are very similar to the MLE results. For the BSE of the BBXEx parameters, the estimated hyper-parameters are a 1 = 354.948, a 2 = 7.899, a 3 = 9.601 and a 4 = 21.055 whereas b 1 = 106405.318, b 2 = 21.849, b 3 = 22.623 and b 4 = 23.216. The MCMC plots for data set II based on complete sample are displayed in Figure 13.  Figure 13. The MCMC plots for data set II using the BBXEx model based on complete sample. Here, we fit the BBXEx model on data set II based on censored samples. In the enclosed Tables 15-17, we provide the MLEs, BSEs, AIC, CAIC, BIC and HQIC values for all tested models. From Tables 15-17 it is observed that, the BBXEx model provides a better fit than the other competitive models. Table 18 shows some descriptive statistics for data set II utilizing the BBXEx distribution and its marginals. According to Table 18, it is clear that the bivariate data has positively skewed with platykurtic. Moreover, the correlation between the two random variables is positive and strong. Positive correlation is a relationship between two variables in which both variables move in tandem that is, in the same direction.

Conclusions
In this paper, we have proposed a bivariate BBX-G family of distributions, whose marginal distributions are BX-G families. It was found that the BBX-G family is suitable of modeling positive skewness and symmetric data sets with leptokurtic phenomena. Moreover, the stress-strength reliability does not depend on the baseline function, but only on the family parameters. The family parameters have been estimated using Bayesian and maximum likelihood methods based on complete and Type-II censored samples, and it was found that the two methods performed quite well in estimating the family parameters. The usefulness of the proposed family is illustrated by two real data sets and it was found that the new family provides a better fit than others sub models and non-nested models. Finally, we can say that the new family will serve as an alternative model to other models available in the literature for modeling positive real data in many areas.