# Classical and Bayesian Inference of a Mixture of Bivariate Exponentiated Exponential Model

Authors: Refah Alotaibi; Mervat Khalifa; Ehab M. Almetwally; Indranil Ghosh; Rezk. H.
Journal: Journal of Mathematics
 (2021)

---

## Abstract

Exponentiated exponential (EE) model has been used effectively in reliability, engineering, biomedical, social sciences, and other applications. In this study, we introduce a new bivariate mixture EE model with two parameters assuming two cases, independent and dependent random variables. We develop a bivariate mixture starting from two EE models assuming two cases, two independent and two dependent EE models. We study some useful statistical properties of this distribution, such as marginals and conditional distributions and product moments and conditional moments. In addition, we study a dependent case, a new mixture of the bivariate model based on EE distribution marginal with two parameters and with a bivariate Gaussian copula. Different methods of estimation for the model parameters are used both under the classical and under the Bayesian paradigm. Some simulation studies are presented to verify the performance of the estimation methods of the proposed model. To illustrate the flexibility of the proposed model, a real dataset is reanalyzed.


---

## Body

## 1. Introduction

In the last two decades or so, a major point of interest for statisticians and practitioners was to study populations that exhibit similar behaviors with respect to some predetermined criteria. The earliest evidence regarding the study of heterogeneous populations was mostly due to Newcomb [1] and Pearson [2] who utilized/developed an approach, commonly known as finite mixture distributions. With modern days’ stellar advancement on long data-related computation facilities, studies focusing on heterogeneous populations became more popular in the modern era; for some useful references, see the work of Titterington et al. [3], Everitt and Hand [4], McLachlan and Basford [5], and AL-Hussaini and Sultan [6] and the references cited therein. In recent times, there is a growing trend to study and explore the application of finite mixture models; for more details, see the work of Al-Hussaini and Sultan [6].In several studies concerning heterogeneous population, the EE probability model appears to be really useful. Two-parameter EE distribution is a right skewed unimodal distribution. The behaviors of the probability density function and the hazard function of the EE distribution is quite close to the behavior of the pdf and the hazard function of the gamma or Weibull model. The two-parameter EE distribution has received an increasing amount of interest in recent times. The efficacy of the EE distribution in modeling lifetime data can be found in the works of Gupta and Kundu [7–12]. Several studies have demonstrated that, in specific real-life scenarios, the EE distribution provides a better fit (based on several well-known goodness-of-fit measures) as compared to the gamma or the Weibull model. Kundu and Gupta [13] introduced a bivariate generalized exponential (BGE) distribution and constructed BGE distribution which has three parameters. Some studies obtained a mixture of bivariate inverse Weibull and gamma models; for more details, see the work of Jones et al [14], Sarhan and Balakrishnan [15], Chen and Tan [16], Khosravi [17], and AL-Moisheer et al. [18] and the references cited therein.The main objective of this paper is to develop and study the mixture of a new bivariate absolutely continuous distribution via a mixture of two independent two-parameter EE distributions. We call this new bivariate distribution as the bivariate mixture of exponentiated exponential distribution (henceforth, in short, the BMEE). The proposed model is constructed under two mechanisms. In the first case, letX and Y be two random variables where each variable is independent and distributed as EE distribution with parameters λ1,θ1 and λ2,θ2, respectively. In the second case, we construct the BMEE distribution via copula approach using the well-known bivariate Gaussian copula (see, for details, the work of Nelson [19]). Several useful mathematical properties of the proposed model are derived. Classical and Bayesian estimation methods are discussed. In addition, the performance of the suggested BMEE model is examined using simulation in estimating the model parameters and a real dataset. The rest of the paper is organized as follows. In Section 2, we introduce the BMEE (type I) distribution, discuss its construction via two mechanisms, and provide some contour plots. In Section 3, we provide some useful mathematical properties and obtain expressions for the bivariate survival function, hazard rate function, bivariate moment generating functions, conditional moments, joint moments, stochastic ordering, etc., for the BMEE distribution constructed starting from two independent EE distributions. In Section 4, we discuss the estimation strategy of the model parameters via EM algorithm. In Section 5, we discuss the estimation strategy for the BMEE (type II) distribution constructed via the bivariate Gaussian copula. In Section 6, we study and explore the estimation of the model parameters under the Bayesian paradigm. Simulation results are presented in Section 7. In Section 8, a well-known motor dataset has been reanalyzed to exhibit the efficacy of the proposed BMEE-type models. Finally, we conclude the paper by providing some final remarks in Section 9.


## 2. Mixture Bivariate Independent EE Model

Here, we begin our discussion with two independent univariate EE distributions with parametersλ1,θ1 and λ2,θ2, respectively. The central idea of compounding is to consider that θ1andθ2 are indeed random variables not constant, and the observed (marginal) distribution of Xin and Yin can be obtained from the joint distribution of θ1andθ2 which is as follows:(1)hxin,yinxin,yin=∬hxin,yin,θ1,θ2dθ1dθ2.Next, we construct a bivariate of EE mixture distribution by assuming two cases. In the first case,Xin and Yin are independent EE distributions with the scale parameters having a generalized bivariate Bernoulli distribution. In the second case, Xd and Yd are dependent. A random variable with an EE distribution has a cumulative distribution function (cdf) and a probability density function (pdf) for Xin>0, given by(2)Fxin,θ,λ=1−e−λxinθ,xin>0,(3)fxin,θ,λ=θλe−λxin1−e−λxinθ−1,xin>0,where θ>0 and λ>0 are the shape and scale parameters.In the bivariate case, letXin and Yin be two random variables with parameters θ1, and θ2, respectively. For given fixed values of θ1 and θ2, Xin and Yin are independent. The pdf of BMEE distribution is defined as(4)fxin,θ,=∑i=12pifixini,θi,where pi are the mixing proportions which must satisfy ∑i=12pi=1 and pi≥0, and all parameters are unknowns. The pdf of the first component of EE is given by (2), with fixed shape parameter θ>0, and a random scale parameter λ>0 that takes two distinct values λ1 and λ2. Likewise, for fixed shape parameter θ2, let Y have an EE mixture density, and the pdf of second component (EE) is given by(5)gyin,ϕ,β=ϕβie−βiyin1−e−βiyinφ−1,yin>0,where β is a random scale parameter β>0 that takes two distinct values β1 and β2. For given values of λ,β, we assume that Xin and Yin are independent, but λ and β are correlated through their generalized bivariate distribution with the following probability matrix:(6)β1β2,P=λ1λ2Pλ1β1Pλ1β2Pλ2β1Pλ2β2,where P is the mixture components and Pλ1β1+Pλ1β2+Pλ2β1+Pλ2β2=1. Let hxin,yinxin,yin be the joint pdf of (Xin, Yin); then,(7)hxin,yinxin,yin=fxin|θ,λ1gyiny|ϕ,β1Pλ1β1+fxin|θ,λ1gyin|ϕ,β2Pλ1β2+fxin|θ,λ2gyin|ϕ,β1Pλ2β1+fxin|θ,λ2gyin|ϕ,β2Pλ2β2.For simplicity in the independent case, we usexin=xandyin=y as follows:(8)hx,y=θλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1Pλ1β1+θλ1e−λ1x1−e−λ1xθ−1ϕβ2e−β2y1−e−β2yϕ−1Pλ1β2+θλ2e−λ2x1−e−λ2xθ−1ϕβ1e−β1y1−e−β1yφ−1Pλ2β1+θλ2e−λ2x1−e−λ2xθ−1ϕβ2e−β2y1−e−β2yφ−1Pλ2β2.For simplification, leta=pλ1β1, b=pλ1β2, c=pλ2β1, and d=pλ2β2.From Figures1 and 2, it is evident that the joint pdf in equation (4) can produce various shapes corresponding to several parameter choices. The joint pdf mixture of four univariate EE mixture distributions involves a total of 9 parameters for its specification. In application to real-life datasets, as one can imagine, not all four components might be necessary. Consequently, one may put some restrictions, such as b = c = 0, a = d = 0, or a = b = c = 0. These restrictions result in correlation values (among scale parameters) of +1, −1, and 0, respectively.Figure 1
The plot of the BMEE pdf model for varying parameter choices in equation (4).Figure 2
The contour plots of BMEE distribution for varying parameter choices in equation (4).The marginal densities ofX and Y, respectively, are given as follows:(9)hxx=π1f1x+1−π1f2x,where π1=a+b.(10)hyY=π2g1Y+1−π2g2Y,where π2=a+c.The joint cdf will be(11)Fx,y=∫0y∫0xht,sdtds=∫0y∫0xaθλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1+bθλ1e−λ1x1−e−λ1xθ−1ϕβ2e−β2y1−e−β2yϕ−1+cθλ2e−λ2x1−e−λ2xθ−1ϕβ1e−β1y1−e−β1yφ−1+dθλ2e−λ2x1−e−λ2xθ−1ϕβ2e−β2y1−e−β2yφ−1dtds=a1−e−λ1xθ1−e−β1yφ+b1−e−λ1xθ1−e−β2yφ+c1−e−λ2xθ1−e−β1yφ+d1−e−λ2xθ1−e−β2yφ.The associated survival function of BMEE distribution will be(12)Rx,y=PX>x,Y>y=1−PX<x−PY<y+PX<x,Y<y=1−FXx−FYy+FX,Yx,y,where FXx and FYy are the marginal density functions of X and Y, respectively. The hazard rate function (hrf) is given as follows:(13)hrfx,y=AB,where A is given in (8) and B is given in (12). The conditional pdf of X, given Y, for each fixed Y = y will be(14)hx|y=hx,yhyy.The conditional pdf ofY, given X, will be(15)hy|x=hx,yhxx.


## 3. Structural Properties

Here, we derive some properties of the BMEE distribution which are as follows.Proposition 1.
LetX,Y∼BMEEa,b,c,d,λ1,λ2,β1,β2,θ,ϕ. Then, X,Y has a total positivity of order 2 positive association (TP2 property).Proof.
Observe that an absolute continuous bivariate random vector, sayU1,U2, has TP2 property if and only if, for any u11, u12,u21,u22, whenever u11<u12 and u21>u22, we have(16)fU1,U2u11,u21fU1,U2u11,u12−fU1,U2u12,u21fU1,U2u11,u22≥0,where fU1,U2, is the joint PDF of U1,U2. Observe that, by taking different ordered u11, u12,u21,u22, such that u11<u12, and u21>u22, our result immediately follows.
As a consequence, the positive quadrant dependence property (alternatively, theTP2 property) will indicate several other nonincreasing properties related to conditional survival, conditional cdf of X, given Y, and Y, given X, including this result: g1. and g2., Covg1X,g2Y≥0.Proposition 2 (moments).
The product moments,μx,y′, about zero is(17)μX,Y′=∫0∞∫0∞xyhx,ydxdy=∫0∞∫0∞xyaθλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1+bθλ1e−λ1x1−e−λ1xθ−1ϕβ2e−β2y1−e−β2yϕ−1+cθλ2e−λ2x1−e−λ2xθ−1ϕβ1e−β1y1−e−β1yφ−1+dθλ2e−λ2x1−e−λ2xθ−1ϕβ2e−β2y1−e−β2yφ−1dxdy=Harmonic  numberθHarmonic numberφaλ1β1+bλ1β2+cλ2β1+dλ2β2.where Harmonic  number. is the sum of the reciprocals of the first n natural numbers.
The expected value ofX and Y are, respectively, given by(18)EX=∫0∞xhxxdx=Harmonic  numberφθaλ1β1e−λ1x1−e−λ1xθ−1+bλ2β1e−λ2x1−e−λ2xθ−1+λ1β2e−λ1x1−e−λ1xθ−1+dλ2β2e−λ2x1−e−λ2xθ−1,(19)EY=∫0∞yhxydy=Harmonic  numberθφaβ1λ1e−β1y1−e−β1yϕ−1+bβ2λ1e−β2y1−e−β2yϕ−1+cβ1λ2e−β1y1−e−β1yϕ−1+dβ2λ2e−β2y1−e−β2yϕ−1.
A Mathematica 11.2 is used to obtain the integral in (17)–(19).Proposition 3.
The joint moment generating function of BMEE will be(20)Mx,yt,s=∫0∞∫0∞etx+syhx,ydxdy=∫0∞∫0∞etx+syaθλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1+bθλ1e−λ1x1−e−λ1xθ−1ϕβ2e−β2y1−e−β2yϕ−1+cθλ2e−λ2x1−e−λ2xθ−1ϕβ1e−β1y1−e−β1yφ−1+dθλ2e−λ2x1−e−λ2xθ−1ϕβ2e−β2y1−e−β2yφ−1dxdy=θφΓθΓφaΓ1−s/β1Γ1−t/λ1Γ1+θ−t/λ1Γ1+φ−s/β1+bΓ1−s/β2Γ1−t/λ2Γ1+θ−t/λ2Γ1+φ−s/β2+cΓ1−s/β1Γ1−t/λ1Γ1+θ−t/λ1Γ1+φ−s/β1+dΓ1−s/β2Γ1−t/λ2Γ1+θ−t/λ2Γ1+φ−s/β2,provided 1+θ−t/λ2≤1 and 1+φ−s/β2≤1.
Using the joint pdf ofX and Y and/or using the MGF expression given above, the different product moments of the order XmXn,m,n≥1, can be obtained from the above.Proposition 4 (shape of the distribution).
A critical point of a function with two variables is a point where the partial derivatives of first order are equal to zero. The most compelling two reasons to study the critical points for a bivariate distribution are as follows:(1)
A real-life dataset(s) can have several different shapes. The flexibility of any proposed model can be well determined from such a study.(2)
In dealing with bivariate distributions, quite often, it is imperative to study the tails of the joint pdf as well as the point of inflection. Knowledge on critical point(s) will help to better understand these properties.
Let us now consider the shape of the BMEE distribution:(21)∂fx,y∂xx=θϕ−aβ1λ12eλ1x−θ1−e−λ1xθ1−e−β1yφe−λ1x−12eβ1y−1−bβ2λ12eλ1x−θ1−e−λ1xθ1−e−β2yφ/e−λ1x−12+λ22eλ2x−θ1−e−λ2xθcβ1eβ2y−11−e−β1yφ+dβ2eβ1y−11−e−β2yφ/e−λ1x−12eβ1y−1eβ2y−1.∂fx,y∂y=θϕ−aβ1λ1eβ1y−ϕ1−e−λ1xθ1−e−β1yφeβ1y−12e−λ1x−1−bλ1β22eβ2y−ϕ1−e−λ1xθ1−e−β2yφ/e−λ1x−1+λ21−e−λ2xθcβ12eβ2y−12eβ1y−φ1−e−β1yφ+dβ22eβ1y−121−e−β2yφeβ2y−φ/e−λ1x−12eβ1y−1eβ2y−1−2.
Consequently, for BMEE distribution, there may be several critical points. For specific choices of the model parameters, a numerical study can be made here.


## 4. Estimation (Independence Case Scenario)

In this section, we discuss the estimation of the model parameters of the BMEE distribution assuming the independence ofX and Y under the EM algorithm.


### 4.1. EM Algorithm

Mclachlan and Krishnan [20] introduced expectation-maximization (EM) algorithm which is an iterative method to find the maximum likelihood estimator of a parameter Φ of a parametric probability distribution. To invoke the idea of EM algorithm, we augment the data (xk,yk,k=1,…,n, with the group membership variables Φk=ak,bk,ck, k = 1,…., n, where ak is an indicator variable, and one if the kth observation is in fx,λ1,β1, and zero otherwise. Similarly, for bk,ck, we have four groups Gij,i,j=1,2, for which the densities are(22)fijx,y=fixfjy=θλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1.The corresponding mixing proportions arePG11=a, PG12=b, PG21=c, and PG22=1−a−b−c.The corresponding log-likelihood function to the complete sample,ℓijx,y=logfijx,y, is represented by(23)ℓ=∑k=1nakℓ11xk,yk+=∑k=1nbkℓ12xk,yk+∑k=1nckℓ21xk,yk+∑k=1n1−ak−bk−ckℓ22xk,yk.Each iteration of the EM algorithm involves two steps: stepE (expectation) and step M (maximization), defined by (16). This is linear in the group membership variables Φk, so in the E-step, we replace in (9) the associated expected values, given the current estimates (θ^,φ,^λ^1,λ^2,β^1,β^2,a^,b,^c^ of the parameters are calculated as(24)a^k=a^f11xk,yka^f11xk,yk+b^f12xk,yk+c^f21xk,yk+1−a^−b^−c^f22xk,yk,and follow the same procedure similarly for bkandck.Next, in the M-step, we need to maximize (23) over θ,φ,β1,β2,λ1,λ2 for fixed values of Φk. It is achieved by the conditional independence of X and Y, given the group membership. Differentiating (23), we obtain(25)∂ι∂λ1=nλ1∑k=1nak+bk−∑k=1nak+bkxk+θ−1∑k=1nxkak+bke−λ1xk1−e−λ1xk,(26)∂ι∂λ2=nλ2∑k=1nck+dk−∑k=1nck+dkxk+θ−1∑k=1nxkck+dke−λ2xk1−e−λ2xk,(27)∂ι∂β1=nβ1∑k=1nak+ck−∑k=1nak+ckyk+φ−1∑k=1nykak+cke−β1yk1−e−β1yk,(28)∂ι∂β2=nβ2∑k=1nbk+dk−∑k=1nbk+dkyk+φ−1∑k=1nykbk+dke−β2yk1−e−β2yk,(29)∂ι∂θ=nθ∑k=1nak+bk+ck+dk+∑k=1nak+ck1−e−λ1xkθ−1log1−e−λ1xk+∑k=1nbk+dk1−e−λ2xkθ−1log1−e−λ2xk,(30)∂ι∂φ=nφ∑k=1nak+bk+ck+dk+∑k=1nak+ck1−e−β1ykφ−1log1−e−β1yk+∑k=1nbk+dk1−e−β2ykφ−1log1−e−β2yk.Next, theM-step is completed by setting(31)a^=1n∑k=1nak,etc.This algorithm requires an initial value of the model parameters, designated byΦ0. A judicious choice of these initial values requires special attention on the fact that the rate of convergence of the assumed EM algorithm may not become quite slow. Another point of concern is that the maximum likelihood equation may have multiple solutions corresponding to local maxima; therefore, the selection of the starting values is indeed very important. A comparative study of various strategies in the choice of initial values can be found in Karlis and Xekalaki [21]. We use the copula R package to solve these equations numerically. After the maximum likelihood estimators for θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22 are obtained, we substitute these estimates in (ak,bk,ck). We complete the M-step by setting a^=1/n∑k=1nak, etc.Initial values for the mixing proportions are obtained by the moment's method of the marginal univariate EE parameter separately. Next, we take the resulting estimates of the BEE parameters as starting values for the EM algorithm. After that, we merge the moment estimators of the marginal mixing parameters to obtain initial values for the bivariate mixing parameters, assuming that the independence between two variablesX and Y. We apply this method in application as mentioned in Section 8; later on, specifically, in Tables 1–4, we equate (25)–(30) to zero, to obtain estimates of the parameter of the distribution.Table 1
MLE and Bayesian estimation with different sample sizes in the independent case.
θ=1.8,λ=3.7,β=3.8,δ=2,ρ=0.5nθλβδρ50MLEBias−0.7241−0.4413−0.3785−0.52160.4645MSE0.80390.49550.20150.61520.2158L.CI2.07372.15070.94612.29740.0305BayesBias−0.4340−0.2691−0.2255−0.31900.4434MSE0.34490.23710.09880.30740.1967L.CI1.55171.59130.85871.77830.0469100MLEBias−0.7098−0.4190−0.3405−0.51760.4565MSE0.79000.42150.19770.55520.2159L.CI1.53871.66930.71951.85400.0231BayesBias−0.3439−0.2276−0.2151−0.30210.4251MSE0.26930.13210.08870.20930.1920L.CI1.08751.11110.63021.27740.0265200MLEBias−0.6857−0.4056−0.3143−0.53640.4264MSE0.72830.41370.18200.50350.2106L.CI1.19001.27440.51331.37150.0169BayesBias−0.3431−0.2142−0.2055−0.29310.4153MSE0.23660.09490.07890.15640.1912L.CI0.88180.86860.46120.97030.0181Table 2
MLE and Bayesian estimation with different sample sizes in the independent case.
θ=2.8,λ=3.7,β=3.7,δ=1.5,ρ=0.5nθλβδρ50MLEBias−0.4552−0.4981−0.5148−0.41160.4647MSE0.30380.59740.38280.52880.2160L.CI1.21872.31791.34602.35100.0295BayesBias−0.0552−0.0646−0.0605−0.04860.3255MSE0.02620.03100.02750.03110.1142L.CI0.59720.64230.60500.66510.3554100MLEBias−0.3504−0.4569−0.4568−0.37500.4064MSE0.28300.49720.34850.36410.1922L.CI0.83951.63510.97771.70620.0224BayesBias−0.0547−0.0617−0.0607−0.03560.3188MSE0.02180.02620.02310.02330.1514L.CI0.51720.57110.51900.58230.1170200MLEBias−0.3053−0.4262−0.4160−0.35330.3846MSE0.23060.35030.30940.39270.1522L.CI0.64461.32350.72251.29130.0162BayesBias−0.0490−0.0611−0.0600−0.02660.3041MSE0.02040.02180.02410.01720.1723L.CI0.43650.51900.46590.50400.0492Table 3
MLE and Bayesian estimation with different sample sizes in dependent case (Case I).
θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.2,ρ12=0.15,ρ21=0.05, and ρ22=0.1nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.5957−0.0916−0.2343−0.0876−0.2470−0.05430.70930.77310.55380.7953MSE0.57300.17290.13360.05520.28410.08360.50850.60640.57130.6626L.CI1.83201.59091.10010.85471.85261.11370.28670.36592.01730.6805BayesBias0.03220.01080.0748−0.0658−0.27860.01420.11080.10940.04330.0975MSE0.02980.03500.03340.02860.23240.02700.04980.04790.03340.0468L.CI0.66560.73260.65370.61111.32760.64250.76010.74390.69660.7571100MLEBias−0.5762−0.0823−0.2263−0.0867−0.2029−0.04990.71200.77320.56030.7975MSE0.53420.13150.10060.03970.20270.05990.50810.60940.57190.6610L.CI1.76371.38540.87170.70381.57650.93990.13280.32041.99210.6190BayesBias0.02150.02070.0349−0.0399−0.76710.02630.04490.04660.02220.0397MSE0.02790.02900.02350.01680.56800.02080.03260.03250.03200.0333L.CI0.64950.66310.58480.48421.18560.55600.68570.68340.69630.6984200MLEBias−0.5672−0.0468−0.2164−0.0944−0.1768−0.04820.71080.77500.52950.7978MSE0.47200.10200.08140.03740.16430.03950.50610.60480.55340.6601L.CI1.52041.23920.72940.66261.43030.75660.10790.25232.04940.6131BayesBias0.01600.0190−0.0367−0.0215−0.68830.02280.25140.23270.09730.1968MSE0.02610.02550.01860.00880.54630.01670.04230.04620.03420.0427L.CI0.65400.62220.51450.35871.05620.49840.76730.82750.77020.8576Table 4
MLE and Bayesian estimation with different sample sizes in the dependent case (Case II).
θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.5,ρ12=0.4,ρ21=0.2, and ρ22=0.3nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.4624−0.0246−0.2569−0.0415−0.2156−0.09370.40200.55330.43520.6387MSE0.41180.13160.12870.04070.21980.07200.16540.30660.38930.4132L.CI1.74521.41930.98230.77421.63260.98600.24070.07831.75380.2861BayesBias0.0165−0.00310.0702−0.0744−0.78250.02610.15960.17940.07000.1406MSE0.03260.03720.03270.02750.71740.02570.04840.06610.04300.0566L.CI0.70530.75670.65330.58181.27120.62050.59410.72270.76550.7529100MLEBias−0.4742−0.0217−0.2463−0.0580−0.2125−0.08400.39710.55140.44870.6365MSE0.37380.09990.11260.03140.18600.05760.16050.30430.38360.4105L.CI1.51371.23670.89410.65711.47210.88210.20970.06121.67440.2871BayesBias0.00050.01470.0232−0.0402−0.27290.03130.08570.09350.02980.0569MSE0.02950.02940.02190.01520.16190.02260.03250.04200.03300.0349L.CI0.67360.67000.57290.45721.16040.57630.62180.71520.70260.6974200MLEBias−0.4674−0.0309−0.2227−0.0681−0.1844−0.06430.39630.55120.40060.6404MSE0.34120.07880.08600.02660.14460.04300.15900.30400.36390.4081L.CI1.37401.09430.74880.58151.30440.77310.17230.04791.76880.1133BayesBias−0.01750.0344−0.0357−0.0153−0.59800.03240.26800.31210.13970.2743MSE0.02360.02290.01700.00880.14220.01820.04820.06120.04160.0511L.CI0.59820.57760.49240.36350.99840.51300.38960.58910.76770.7211


## 4.1. EM Algorithm

Mclachlan and Krishnan [20] introduced expectation-maximization (EM) algorithm which is an iterative method to find the maximum likelihood estimator of a parameter Φ of a parametric probability distribution. To invoke the idea of EM algorithm, we augment the data (xk,yk,k=1,…,n, with the group membership variables Φk=ak,bk,ck, k = 1,…., n, where ak is an indicator variable, and one if the kth observation is in fx,λ1,β1, and zero otherwise. Similarly, for bk,ck, we have four groups Gij,i,j=1,2, for which the densities are(22)fijx,y=fixfjy=θλ1e−λ1x1−e−λ1xθ−1ϕβ1e−β1y1−e−β1yφ−1.The corresponding mixing proportions arePG11=a, PG12=b, PG21=c, and PG22=1−a−b−c.The corresponding log-likelihood function to the complete sample,ℓijx,y=logfijx,y, is represented by(23)ℓ=∑k=1nakℓ11xk,yk+=∑k=1nbkℓ12xk,yk+∑k=1nckℓ21xk,yk+∑k=1n1−ak−bk−ckℓ22xk,yk.Each iteration of the EM algorithm involves two steps: stepE (expectation) and step M (maximization), defined by (16). This is linear in the group membership variables Φk, so in the E-step, we replace in (9) the associated expected values, given the current estimates (θ^,φ,^λ^1,λ^2,β^1,β^2,a^,b,^c^ of the parameters are calculated as(24)a^k=a^f11xk,yka^f11xk,yk+b^f12xk,yk+c^f21xk,yk+1−a^−b^−c^f22xk,yk,and follow the same procedure similarly for bkandck.Next, in the M-step, we need to maximize (23) over θ,φ,β1,β2,λ1,λ2 for fixed values of Φk. It is achieved by the conditional independence of X and Y, given the group membership. Differentiating (23), we obtain(25)∂ι∂λ1=nλ1∑k=1nak+bk−∑k=1nak+bkxk+θ−1∑k=1nxkak+bke−λ1xk1−e−λ1xk,(26)∂ι∂λ2=nλ2∑k=1nck+dk−∑k=1nck+dkxk+θ−1∑k=1nxkck+dke−λ2xk1−e−λ2xk,(27)∂ι∂β1=nβ1∑k=1nak+ck−∑k=1nak+ckyk+φ−1∑k=1nykak+cke−β1yk1−e−β1yk,(28)∂ι∂β2=nβ2∑k=1nbk+dk−∑k=1nbk+dkyk+φ−1∑k=1nykbk+dke−β2yk1−e−β2yk,(29)∂ι∂θ=nθ∑k=1nak+bk+ck+dk+∑k=1nak+ck1−e−λ1xkθ−1log1−e−λ1xk+∑k=1nbk+dk1−e−λ2xkθ−1log1−e−λ2xk,(30)∂ι∂φ=nφ∑k=1nak+bk+ck+dk+∑k=1nak+ck1−e−β1ykφ−1log1−e−β1yk+∑k=1nbk+dk1−e−β2ykφ−1log1−e−β2yk.Next, theM-step is completed by setting(31)a^=1n∑k=1nak,etc.This algorithm requires an initial value of the model parameters, designated byΦ0. A judicious choice of these initial values requires special attention on the fact that the rate of convergence of the assumed EM algorithm may not become quite slow. Another point of concern is that the maximum likelihood equation may have multiple solutions corresponding to local maxima; therefore, the selection of the starting values is indeed very important. A comparative study of various strategies in the choice of initial values can be found in Karlis and Xekalaki [21]. We use the copula R package to solve these equations numerically. After the maximum likelihood estimators for θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22 are obtained, we substitute these estimates in (ak,bk,ck). We complete the M-step by setting a^=1/n∑k=1nak, etc.Initial values for the mixing proportions are obtained by the moment's method of the marginal univariate EE parameter separately. Next, we take the resulting estimates of the BEE parameters as starting values for the EM algorithm. After that, we merge the moment estimators of the marginal mixing parameters to obtain initial values for the bivariate mixing parameters, assuming that the independence between two variablesX and Y. We apply this method in application as mentioned in Section 8; later on, specifically, in Tables 1–4, we equate (25)–(30) to zero, to obtain estimates of the parameter of the distribution.Table 1
MLE and Bayesian estimation with different sample sizes in the independent case.
θ=1.8,λ=3.7,β=3.8,δ=2,ρ=0.5nθλβδρ50MLEBias−0.7241−0.4413−0.3785−0.52160.4645MSE0.80390.49550.20150.61520.2158L.CI2.07372.15070.94612.29740.0305BayesBias−0.4340−0.2691−0.2255−0.31900.4434MSE0.34490.23710.09880.30740.1967L.CI1.55171.59130.85871.77830.0469100MLEBias−0.7098−0.4190−0.3405−0.51760.4565MSE0.79000.42150.19770.55520.2159L.CI1.53871.66930.71951.85400.0231BayesBias−0.3439−0.2276−0.2151−0.30210.4251MSE0.26930.13210.08870.20930.1920L.CI1.08751.11110.63021.27740.0265200MLEBias−0.6857−0.4056−0.3143−0.53640.4264MSE0.72830.41370.18200.50350.2106L.CI1.19001.27440.51331.37150.0169BayesBias−0.3431−0.2142−0.2055−0.29310.4153MSE0.23660.09490.07890.15640.1912L.CI0.88180.86860.46120.97030.0181Table 2
MLE and Bayesian estimation with different sample sizes in the independent case.
θ=2.8,λ=3.7,β=3.7,δ=1.5,ρ=0.5nθλβδρ50MLEBias−0.4552−0.4981−0.5148−0.41160.4647MSE0.30380.59740.38280.52880.2160L.CI1.21872.31791.34602.35100.0295BayesBias−0.0552−0.0646−0.0605−0.04860.3255MSE0.02620.03100.02750.03110.1142L.CI0.59720.64230.60500.66510.3554100MLEBias−0.3504−0.4569−0.4568−0.37500.4064MSE0.28300.49720.34850.36410.1922L.CI0.83951.63510.97771.70620.0224BayesBias−0.0547−0.0617−0.0607−0.03560.3188MSE0.02180.02620.02310.02330.1514L.CI0.51720.57110.51900.58230.1170200MLEBias−0.3053−0.4262−0.4160−0.35330.3846MSE0.23060.35030.30940.39270.1522L.CI0.64461.32350.72251.29130.0162BayesBias−0.0490−0.0611−0.0600−0.02660.3041MSE0.02040.02180.02410.01720.1723L.CI0.43650.51900.46590.50400.0492Table 3
MLE and Bayesian estimation with different sample sizes in dependent case (Case I).
θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.2,ρ12=0.15,ρ21=0.05, and ρ22=0.1nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.5957−0.0916−0.2343−0.0876−0.2470−0.05430.70930.77310.55380.7953MSE0.57300.17290.13360.05520.28410.08360.50850.60640.57130.6626L.CI1.83201.59091.10010.85471.85261.11370.28670.36592.01730.6805BayesBias0.03220.01080.0748−0.0658−0.27860.01420.11080.10940.04330.0975MSE0.02980.03500.03340.02860.23240.02700.04980.04790.03340.0468L.CI0.66560.73260.65370.61111.32760.64250.76010.74390.69660.7571100MLEBias−0.5762−0.0823−0.2263−0.0867−0.2029−0.04990.71200.77320.56030.7975MSE0.53420.13150.10060.03970.20270.05990.50810.60940.57190.6610L.CI1.76371.38540.87170.70381.57650.93990.13280.32041.99210.6190BayesBias0.02150.02070.0349−0.0399−0.76710.02630.04490.04660.02220.0397MSE0.02790.02900.02350.01680.56800.02080.03260.03250.03200.0333L.CI0.64950.66310.58480.48421.18560.55600.68570.68340.69630.6984200MLEBias−0.5672−0.0468−0.2164−0.0944−0.1768−0.04820.71080.77500.52950.7978MSE0.47200.10200.08140.03740.16430.03950.50610.60480.55340.6601L.CI1.52041.23920.72940.66261.43030.75660.10790.25232.04940.6131BayesBias0.01600.0190−0.0367−0.0215−0.68830.02280.25140.23270.09730.1968MSE0.02610.02550.01860.00880.54630.01670.04230.04620.03420.0427L.CI0.65400.62220.51450.35871.05620.49840.76730.82750.77020.8576Table 4
MLE and Bayesian estimation with different sample sizes in the dependent case (Case II).
θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.5,ρ12=0.4,ρ21=0.2, and ρ22=0.3nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.4624−0.0246−0.2569−0.0415−0.2156−0.09370.40200.55330.43520.6387MSE0.41180.13160.12870.04070.21980.07200.16540.30660.38930.4132L.CI1.74521.41930.98230.77421.63260.98600.24070.07831.75380.2861BayesBias0.0165−0.00310.0702−0.0744−0.78250.02610.15960.17940.07000.1406MSE0.03260.03720.03270.02750.71740.02570.04840.06610.04300.0566L.CI0.70530.75670.65330.58181.27120.62050.59410.72270.76550.7529100MLEBias−0.4742−0.0217−0.2463−0.0580−0.2125−0.08400.39710.55140.44870.6365MSE0.37380.09990.11260.03140.18600.05760.16050.30430.38360.4105L.CI1.51371.23670.89410.65711.47210.88210.20970.06121.67440.2871BayesBias0.00050.01470.0232−0.0402−0.27290.03130.08570.09350.02980.0569MSE0.02950.02940.02190.01520.16190.02260.03250.04200.03300.0349L.CI0.67360.67000.57290.45721.16040.57630.62180.71520.70260.6974200MLEBias−0.4674−0.0309−0.2227−0.0681−0.1844−0.06430.39630.55120.40060.6404MSE0.34120.07880.08600.02660.14460.04300.15900.30400.36390.4081L.CI1.37401.09430.74880.58151.30440.77310.17230.04791.76880.1133BayesBias−0.01750.0344−0.0357−0.0153−0.59800.03240.26800.31210.13970.2743MSE0.02360.02290.01700.00880.14220.01820.04820.06120.04160.0511L.CI0.59820.57760.49240.36350.99840.51300.38960.58910.76770.7211


## 5. BMEE Model (Type II) Distribution Using Gaussian Copula

The BMEE model proposed in this paper involves EE marginals, with greatest flexibility with its marginals as well as in the correlation structure. On the contrary, this proposed distribution has several fields of applicability. Usually, in the dependence study, copulas play a vital role. Copulas are a general tool to construct bivariate and multivariate distributions and to study dependence structure between random variables. Several bivariate and multivariate lifetime distributions are suggested using several methods of constructing bivariate and multivariate distributions, and copula functions have been proposed by Nelsen [22], Trivedi and Zimmer [23], Adham and Walker [24], Kundu et al. [25], Kundu and Gupta [26], Kundu [27], El-Morshedy et al. [28], and Alotaibi et al. [29].


### 5.1. BMEE Distribution Based on Gaussian Copula

The concept of copula, suggested and derived by Sklar [30], states that any multivariate distribution can be disintegrated to a copula and its continuous marginal. In a bivariate setup, copulas are used to link two marginal distributions with a joint distribution such that, for every bivariate distribution function Fx,y with continuous marginal Fx,Fy, there exists a unique copula function C given by(32)Fx,y=CFx,Fy,x,y∈−∞,∞X−∞,∞.The associated density function of bivariate distribution will be(33)fx,y=fxfycFx,Fy,where cFx,Fy is the density function of a copula; for further details, see the work of Nelsen [19, 22]. A plethora of choices are available to construct BMEE distributions via copula using EE marginals as given in (1). Here, the Gaussian copula is utilized to construct BMEE distribution.The Gaussian copula has the following form:(34)Cu,v=φΣφ−1u,φ−1v,where φΣ denotes the distribution function of a bivariate standard normal random variable and φ−1 represents its inverse. The joint pdf of Xd and Yd based on Gaussian copula becomes(35)fxd,yd=fxdfyd11−ρ2exp−ρ21−ρ2ρz12+z22−2z1z2,where ρ∈−1,1 is a dependence parameter and fxdandfyd are the density function of EE distributions given in (1). Suppose that the marginals are EE distribution; then, the bivariate exponentiated exponential (BEE) distribution pdf is(36)fxd,yd,θ,λ,β,δ=θλe−λxd1−exp−λxdθ−1δβe−βyd1−exp−βydδ−111−ρ2exp−ρ21−ρ2ρz12+z22−2z1z2,xd,yd,θ,λ,β,δ>0.The pdf of BMEE distribution is defined as(37)fxdx,θ,=∑i=12pifixdi,θi,where pi, are the mixing proportions, and it must satisfy ∑i=12pi=1 and pi≥0, and all of them are unknown. The pdf of first component of EE is given by (1):(38)fxd,θ,λi=θλie−λixd1−e−λixdθ−1,xd>0,with fixed shape parameter θ>0 and random scale parameter λ>0, which take two distinct values λ1 and λ2, respectively. Similarly, for fixed shape parameter θ2, let Yd have an EE mixture density, and the pdf of second component EE is given by(39)gyd,ϕ,β=δβie−βiyd1−e−βiydδ−1,yd>0,with β being a random scale parameter taking values β1 and β2.For a given valuesλ,β, we assume that Xd and Yd are dependent, and λ and β are correlated through their generalized bivariate distribution with the probability matrix given by(40)β1β2p=λ1λ2pλ1β1pλ1β2pλ2β1pλ2β2.Letf (xd, yd) be the joint pdf of (Xd, Yd); then,(41)fxd,yd=fxd|θ,λ1gyd|δ,β1cFxd|θ,λ1,Gyd|δ,β1pλ1β1+fxd|θ,λ1gyd|δ,β2cFxd|θ,λ1,Gyd|δ,β2pλ1β2+fxd|θ,λ2gyd|δ,β1cFxd|θ,λ2,Gyd|δ,β1pλ2β1+fxd|θ,λ2gyd|δ,β2cFxd|θ,λ2,Gyd|δ,β2pλ2β2.Like the joint pdf in (8), the joint pdf in (41) can assume several different shapes as well. Let xd=xandyd=y. Consequently, the associated BMEE distribution pdf will be(42)fx,y,θ,λ,β,δ=θλ1e−λ1x1−exp−λ1xθ−1δβ1e−β1y1−exp−β1yδ−111−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2pλ1β1+θλ1e−λ1x1−exp−λ1xθ−1δβ2e−β2y1−exp−β2yδ−111−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2pλ1β2+λ2e−λ2x1−exp−λ2xθ−1δβ1e−β1y1−exp−β1yδ−111−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2pλ2β1+λ2e−λ2x1−exp−λ2xθ−1δβ2e−β2y1−exp−β2yδ−111−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2pλ2β2,x,y,θ,λ,β,δ>0,where ρij∈−1,1 is the dependence parameter.


### 5.2. EM Algorithm under Gaussian Copula

The EM algorithm is introduced as a method of estimation. To apply the EM algorithm, as before, we augment the data (xk,yk), k = 1,…, n, with the group membership variables (ak,bk,ck), k = 1,…, n, where ak with the group membership variables (ak,bk,ck), k = 1,…, n, where ak is one if the kth observation is in fijx,y,θ,λ1,β1,δ,ρ11, and zero otherwise. Similarly, for bk and ck, we have four groups Gij, i, j = 1, 2, for which the densities are(43)fijx,y,θ,λi,βj,δ,ρij=θλie−λix1−exp−λixθ−1δβje−βjy1−exp−βjyδ−111−ρij2exp−ρij21−ρij2ρijz12+z22−2z1z2,x,y,θ,λi,βj,δ>0,where ρij∈−1,1 is a dependence parameter.The mixing proportions are given as follows:P (G11) = a, P (G12) = b, P (G21) = c, and P (G22) = 1 − a − b − c. We define ℓijx,y=logfijx,y,θ,λi,βj,δ,ρij; then, the EM algorithm as the method of estimation is given by finding the complete log likelihood, ℓ, as follows:(44)ℓ=∑k=1nakℓ11xk,yk+∑k=1nbkℓ12xk,yk+∑k=1nckℓ12xk,yk+∑k=1n1−ak−bk−ckℓ22xk,yk.This is linear in the group membership variables (ak,bk,ck); consequently, in the E step, we enter into (26); their expected values, given the current estimates θ^,λ1^,λ2^,β1^,β2^,δ^,ρ11^,ρ12^,ρ21^,ρ22^,a,b,c of the parameter, are calculated as(45)ak^=a^f11xk,yka^f11xk,yk+b^f12xk,yk+c^f21xk,yk+1−a^−b^−c^f22xk,yk.Similarly, forbk and ck, we follow the same strategies. Note that the algebraic simplification of the above might be necessary to avoid numerical problems. For the M-step, we need to maximize (27) over θ,λ1,β1,δ,ρ11,ρ12,ρ21,ρ22, for fixed values of (ak,bk,ck). This is achieved by the conditional dependence of X and Y, given the group membership. We can essentially deal with the univariates and the Gaussian copula parameter separately. Differentiating (25) gives(46)∂l∂λ1=nλ1∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+∑k=1nxθexp−λ1x−11−exp−λ1x∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2=0,(47)nλ2∑k=1nck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2+∑k=1nxθexp−λ2xk−11−exp−λ2xk∑k=1nck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(48)∂l∂β1=nβ1∑k=1nak11−ρ112exp−ρ2121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+∑k=1nxθexp−β1yk−11−exp−β1yk∑k=1nak11−ρ112exp−ρ2121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2=0,(49)∂l∂β2=nβ2∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2+∑k=1nykθexp−β2yk−11−exp−β2yk∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(50)∂l∂θ=nθ+∑k=1nlog1−exp−λ1xk∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+nθ+∑k=1nlog1−exp−λ2xk∑k=1nck11−ρ112exp−ρ2121−ρ212ρ11z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(51)∂l∂δ=∂l∂δnδ+∑k=1nlog1−exp−β1yk∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+nδ+∑k=1nlog1−exp−β2yk∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0.The approach involves a two-step procedure in estimating the marginal ofX and Y and the copula function independently that gives the maximum likelihood estimation of θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22. The solution of the nonlinear (46)–(51) gives the MLE of θ,λ1,β1,λ2,β2,andδ and θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22.Then, copula density is estimated as(52)logLijρ=∑k=1nlogcFi^xk,F^jyk,where Fi^xk and ,F^jyk denote the maximum likelihood estimates of the pdf from the first step. The solution of the nonlinear equation (34) gives the MLE of ρ11,ρ12,ρ21, and ρ22.The M-step is completed by setting(53)a^=1n∑k=1nak.We use the copulaR package to solve these equations numerically. After the maximum likelihood estimators for θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22 are obtained, next, we substitute these estimates in (ak,bk,ck). We complete the M-step by setting a^=1/n∑k=1nak, etc.Initial values of the parameters for the mixing proportions are obtained by the method of matching moments that are obtained from the marginal univariate EEM and the Gaussian copula parameter separately, Then, we take the resulting estimates of the BEE parameters as starting values for the EM algorithm. Next, we merge the moment estimators of the marginal mixing parameters to obtain initial values for the bivariate mixing parameters, assuming the dependence between two variablesX and Y. We apply this method in application as mentioned in Section 9, specifically, in Tables 1–4. For more details, see the work of Kosmidis and Karlis [31].Next, we provide the estimation procedure of the unknown parameters for the density in (44). In the copula-based estimation, we adopt two approaches that are termed as parametric and semiparametric.


### 5.3. Maximum Likelihood Estimation (MLE)

Here, we discuss the estimation of the unknown parameters of BEE distributions by the approach of the maximum likelihood, by using the two-step estimation. It involves a two-step procedure by which we estimate the marginal and the copula function separately.The log-likelihood function is expressed as(54)logL=∑i=1nlogf1xi+logf2yi+logcF1xi,F2yi.The log-likelihood function in (35) can be re-expressed as(55)logL=∑i=1nlogf1xi+∑i=1nlogf2x2i+∑i=1nlogcF1x1i,F2x2i.The first step involves estimating the parameters of marginals distributionF1 and F2 by MLE, separately given as follows:(56)logL1=∑i=1nloglogf1xi,logL2=∑i=1nlogf2yi.Then, estimating copula parameters by maximizing the copula density, we will obtain(57)logL=∑i=1nlogcF1xi,F2yi.By considering the first step with EE distributions, the parameters of each marginal distribution will be estimated by the MLE. Ifx1,…,xn is a random sample from EEθ,λ and y1,…,yn is a random sample from EEδ,β, then the log-likelihood functions are, respectively, given by(58)logL1x,θ,λ=nlogθ+nlogλ−∑i=1nxi+θ−1∑i=1nlog1−exp−λxi,(59)logL2y,δ,β=nlogδ+nlogβ−∑i=1nyi+δ−1∑i=1nlog1−exp−βyi.So, the maximum likelihood equations are(60)∂l∂λ=nλ+∑i=1nxiθexp−λxi−11−exp−λxi=0,(61)∂l∂β=nβ+∑i=1nyiδexp−βyi−11−exp−βyi=0,(62)∂l∂θ=nθ+∑i=1nlog1−exp−λxi=0,(63)∂l∂δ=nδ+∑i=1nlog1−exp−βyi=0,The solution of the system of nonlinear equations (60)–(63) gives the MLE of θ,λ,β, and δ. Then, copula density is estimated as follows:(64)logLγ=∑i=1nlogcF1^xi,F^2yi,where F1^x and F^2y denote the ML estimates of the parameters from the first step.The solution of the nonlinear equation (64) gives the MLE of γ.


### 5.4. Semiparametric Methods of Estimation

Two semiparametric methods are used to estimate the copula parameter in the copula models and are compared with the two methods of moments approaches which are the inversion Kendall'sτ and inversion of Spearman's ρ, respectively.


### 5.5. Methods of Moments

From the moment’s method of inversion of Kendall'sτ and the inversion of Spearman's ρ mentioned in Kojadinovic and Yan [32], we provide a brief details which are given as follows. Let c be a bivariate random sample from a cdf CγF1x,F2y, where F1 and F2 are continuous cdf’s and Cγ is an absolutely continuous copula such that γ∈O, where O is an open subset of R2. Furthermore, let R1,…,Rn are the vectors of ranks associated with x1,…,xn unless otherwise stated. In what follows, all vectors are row vectors. Moment’s approaches are based on the inversion of a consistent estimator of a moment of the copula Cγ. The two best-known moments, Spearman’s rho and Kendall’s tau, are, respectively, given by(65)ργ=12∫∫0,12uvdCγu,v−3,(66)τγ=4∫∫0,12Cγu,vdCγu,v−1.Consistent estimators of these two moments can be expressed as(67)ρn=12nn+1n−1∑i=1nRi,1Ri,2−3n+1n−1,(68)τn=4nn−1∑i=1n1xi,1≤xj,11xi,2≤xj,2−1.Ifρ and τ are one-to-one, consistent estimators of γ will be γn,ρ=ρ−1ρn and γn,τ=τ−1τn, respectively.It can be called inversion of Kendall'sτ and inversion of Spearman's ρ, respectively. For more information, see the work of Kojadinovic and Yan [32] and the references cited therein.As explained above, the moment’s method ofτ and ρ estimation for copula may be considered under the umbrella of semiparametric approach estimation.


### 5.6. Goodness of Fit Tests for Copula

We want to compare the empirical copula with the parametric estimator derived under the null hypothesis; for details, see the work of Fermanian [33]. Theory suggests a test if C is well-represented by a specific copula Cγ:(69)H0:C=CγVs.H1:C≠Cγ.Several well-known approaches are available in the literature; for example, see the work of Genest and Rémillard [34], or the fast multiplier approach, Genest et al. [35], and Kojadinovic et al. [36]. The goodness of fit tests based on the empirical process is given as(70)ℂnu,v=nCnu,v−Cγnu,v,where Cnu,v is the empirical copula of the data of X and YCnu,v=1/n∑i=1n1Ui,n≤u,Vi,n≤v,u,v∈0,1, and Ui,nandVi,n are pseudoobservations from C calculated from data as follows.Ui,n=R1i/n+1, Vi,n=R2i/n+1, and R1i and R2i are, respectively, the ranks of XiandYi.Here,Cnu,v is a consistent estimator, and θn is an estimator of γ obtained using the pseudoobservations. According to Genest et al. [35], the appropriate test statistics is the Cramer-von Miss and is defined as(71)Sn=∑i=1nCnUi,n,Vi,n−CγnUi,n,Vi,n2.


## 5.1. BMEE Distribution Based on Gaussian Copula

The concept of copula, suggested and derived by Sklar [30], states that any multivariate distribution can be disintegrated to a copula and its continuous marginal. In a bivariate setup, copulas are used to link two marginal distributions with a joint distribution such that, for every bivariate distribution function Fx,y with continuous marginal Fx,Fy, there exists a unique copula function C given by(32)Fx,y=CFx,Fy,x,y∈−∞,∞X−∞,∞.The associated density function of bivariate distribution will be(33)fx,y=fxfycFx,Fy,where cFx,Fy is the density function of a copula; for further details, see the work of Nelsen [19, 22]. A plethora of choices are available to construct BMEE distributions via copula using EE marginals as given in (1). Here, the Gaussian copula is utilized to construct BMEE distribution.The Gaussian copula has the following form:(34)Cu,v=φΣφ−1u,φ−1v,where φΣ denotes the distribution function of a bivariate standard normal random variable and φ−1 represents its inverse. The joint pdf of Xd and Yd based on Gaussian copula becomes(35)fxd,yd=fxdfyd11−ρ2exp−ρ21−ρ2ρz12+z22−2z1z2,where ρ∈−1,1 is a dependence parameter and fxdandfyd are the density function of EE distributions given in (1). Suppose that the marginals are EE distribution; then, the bivariate exponentiated exponential (BEE) distribution pdf is(36)fxd,yd,θ,λ,β,δ=θλe−λxd1−exp−λxdθ−1δβe−βyd1−exp−βydδ−111−ρ2exp−ρ21−ρ2ρz12+z22−2z1z2,xd,yd,θ,λ,β,δ>0.The pdf of BMEE distribution is defined as(37)fxdx,θ,=∑i=12pifixdi,θi,where pi, are the mixing proportions, and it must satisfy ∑i=12pi=1 and pi≥0, and all of them are unknown. The pdf of first component of EE is given by (1):(38)fxd,θ,λi=θλie−λixd1−e−λixdθ−1,xd>0,with fixed shape parameter θ>0 and random scale parameter λ>0, which take two distinct values λ1 and λ2, respectively. Similarly, for fixed shape parameter θ2, let Yd have an EE mixture density, and the pdf of second component EE is given by(39)gyd,ϕ,β=δβie−βiyd1−e−βiydδ−1,yd>0,with β being a random scale parameter taking values β1 and β2.For a given valuesλ,β, we assume that Xd and Yd are dependent, and λ and β are correlated through their generalized bivariate distribution with the probability matrix given by(40)β1β2p=λ1λ2pλ1β1pλ1β2pλ2β1pλ2β2.Letf (xd, yd) be the joint pdf of (Xd, Yd); then,(41)fxd,yd=fxd|θ,λ1gyd|δ,β1cFxd|θ,λ1,Gyd|δ,β1pλ1β1+fxd|θ,λ1gyd|δ,β2cFxd|θ,λ1,Gyd|δ,β2pλ1β2+fxd|θ,λ2gyd|δ,β1cFxd|θ,λ2,Gyd|δ,β1pλ2β1+fxd|θ,λ2gyd|δ,β2cFxd|θ,λ2,Gyd|δ,β2pλ2β2.Like the joint pdf in (8), the joint pdf in (41) can assume several different shapes as well. Let xd=xandyd=y. Consequently, the associated BMEE distribution pdf will be(42)fx,y,θ,λ,β,δ=θλ1e−λ1x1−exp−λ1xθ−1δβ1e−β1y1−exp−β1yδ−111−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2pλ1β1+θλ1e−λ1x1−exp−λ1xθ−1δβ2e−β2y1−exp−β2yδ−111−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2pλ1β2+λ2e−λ2x1−exp−λ2xθ−1δβ1e−β1y1−exp−β1yδ−111−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2pλ2β1+λ2e−λ2x1−exp−λ2xθ−1δβ2e−β2y1−exp−β2yδ−111−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2pλ2β2,x,y,θ,λ,β,δ>0,where ρij∈−1,1 is the dependence parameter.


## 5.2. EM Algorithm under Gaussian Copula

The EM algorithm is introduced as a method of estimation. To apply the EM algorithm, as before, we augment the data (xk,yk), k = 1,…, n, with the group membership variables (ak,bk,ck), k = 1,…, n, where ak with the group membership variables (ak,bk,ck), k = 1,…, n, where ak is one if the kth observation is in fijx,y,θ,λ1,β1,δ,ρ11, and zero otherwise. Similarly, for bk and ck, we have four groups Gij, i, j = 1, 2, for which the densities are(43)fijx,y,θ,λi,βj,δ,ρij=θλie−λix1−exp−λixθ−1δβje−βjy1−exp−βjyδ−111−ρij2exp−ρij21−ρij2ρijz12+z22−2z1z2,x,y,θ,λi,βj,δ>0,where ρij∈−1,1 is a dependence parameter.The mixing proportions are given as follows:P (G11) = a, P (G12) = b, P (G21) = c, and P (G22) = 1 − a − b − c. We define ℓijx,y=logfijx,y,θ,λi,βj,δ,ρij; then, the EM algorithm as the method of estimation is given by finding the complete log likelihood, ℓ, as follows:(44)ℓ=∑k=1nakℓ11xk,yk+∑k=1nbkℓ12xk,yk+∑k=1nckℓ12xk,yk+∑k=1n1−ak−bk−ckℓ22xk,yk.This is linear in the group membership variables (ak,bk,ck); consequently, in the E step, we enter into (26); their expected values, given the current estimates θ^,λ1^,λ2^,β1^,β2^,δ^,ρ11^,ρ12^,ρ21^,ρ22^,a,b,c of the parameter, are calculated as(45)ak^=a^f11xk,yka^f11xk,yk+b^f12xk,yk+c^f21xk,yk+1−a^−b^−c^f22xk,yk.Similarly, forbk and ck, we follow the same strategies. Note that the algebraic simplification of the above might be necessary to avoid numerical problems. For the M-step, we need to maximize (27) over θ,λ1,β1,δ,ρ11,ρ12,ρ21,ρ22, for fixed values of (ak,bk,ck). This is achieved by the conditional dependence of X and Y, given the group membership. We can essentially deal with the univariates and the Gaussian copula parameter separately. Differentiating (25) gives(46)∂l∂λ1=nλ1∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+∑k=1nxθexp−λ1x−11−exp−λ1x∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2=0,(47)nλ2∑k=1nck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2+∑k=1nxθexp−λ2xk−11−exp−λ2xk∑k=1nck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(48)∂l∂β1=nβ1∑k=1nak11−ρ112exp−ρ2121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+∑k=1nxθexp−β1yk−11−exp−β1yk∑k=1nak11−ρ112exp−ρ2121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2=0,(49)∂l∂β2=nβ2∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2+∑k=1nykθexp−β2yk−11−exp−β2yk∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(50)∂l∂θ=nθ+∑k=1nlog1−exp−λ1xk∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+bk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+nθ+∑k=1nlog1−exp−λ2xk∑k=1nck11−ρ112exp−ρ2121−ρ212ρ11z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0,(51)∂l∂δ=∂l∂δnδ+∑k=1nlog1−exp−β1yk∑k=1nak11−ρ112exp−ρ1121−ρ112ρ11z12+z22−2z1z2+ck11−ρ212exp−ρ2121−ρ212ρ21z12+z22−2z1z2+nδ+∑k=1nlog1−exp−β2yk∑k=1nbk11−ρ122exp−ρ1221−ρ122ρ12z12+z22−2z1z2+1−ak−bk−ck11−ρ222exp−ρ2221−ρ222ρ22z12+z22−2z1z2=0.The approach involves a two-step procedure in estimating the marginal ofX and Y and the copula function independently that gives the maximum likelihood estimation of θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22. The solution of the nonlinear (46)–(51) gives the MLE of θ,λ1,β1,λ2,β2,andδ and θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22.Then, copula density is estimated as(52)logLijρ=∑k=1nlogcFi^xk,F^jyk,where Fi^xk and ,F^jyk denote the maximum likelihood estimates of the pdf from the first step. The solution of the nonlinear equation (34) gives the MLE of ρ11,ρ12,ρ21, and ρ22.The M-step is completed by setting(53)a^=1n∑k=1nak.We use the copulaR package to solve these equations numerically. After the maximum likelihood estimators for θ,λ1,β1,λ2,β2,δ,ρ11,ρ12,ρ21, and ρ22 are obtained, next, we substitute these estimates in (ak,bk,ck). We complete the M-step by setting a^=1/n∑k=1nak, etc.Initial values of the parameters for the mixing proportions are obtained by the method of matching moments that are obtained from the marginal univariate EEM and the Gaussian copula parameter separately, Then, we take the resulting estimates of the BEE parameters as starting values for the EM algorithm. Next, we merge the moment estimators of the marginal mixing parameters to obtain initial values for the bivariate mixing parameters, assuming the dependence between two variablesX and Y. We apply this method in application as mentioned in Section 9, specifically, in Tables 1–4. For more details, see the work of Kosmidis and Karlis [31].Next, we provide the estimation procedure of the unknown parameters for the density in (44). In the copula-based estimation, we adopt two approaches that are termed as parametric and semiparametric.


## 5.3. Maximum Likelihood Estimation (MLE)

Here, we discuss the estimation of the unknown parameters of BEE distributions by the approach of the maximum likelihood, by using the two-step estimation. It involves a two-step procedure by which we estimate the marginal and the copula function separately.The log-likelihood function is expressed as(54)logL=∑i=1nlogf1xi+logf2yi+logcF1xi,F2yi.The log-likelihood function in (35) can be re-expressed as(55)logL=∑i=1nlogf1xi+∑i=1nlogf2x2i+∑i=1nlogcF1x1i,F2x2i.The first step involves estimating the parameters of marginals distributionF1 and F2 by MLE, separately given as follows:(56)logL1=∑i=1nloglogf1xi,logL2=∑i=1nlogf2yi.Then, estimating copula parameters by maximizing the copula density, we will obtain(57)logL=∑i=1nlogcF1xi,F2yi.By considering the first step with EE distributions, the parameters of each marginal distribution will be estimated by the MLE. Ifx1,…,xn is a random sample from EEθ,λ and y1,…,yn is a random sample from EEδ,β, then the log-likelihood functions are, respectively, given by(58)logL1x,θ,λ=nlogθ+nlogλ−∑i=1nxi+θ−1∑i=1nlog1−exp−λxi,(59)logL2y,δ,β=nlogδ+nlogβ−∑i=1nyi+δ−1∑i=1nlog1−exp−βyi.So, the maximum likelihood equations are(60)∂l∂λ=nλ+∑i=1nxiθexp−λxi−11−exp−λxi=0,(61)∂l∂β=nβ+∑i=1nyiδexp−βyi−11−exp−βyi=0,(62)∂l∂θ=nθ+∑i=1nlog1−exp−λxi=0,(63)∂l∂δ=nδ+∑i=1nlog1−exp−βyi=0,The solution of the system of nonlinear equations (60)–(63) gives the MLE of θ,λ,β, and δ. Then, copula density is estimated as follows:(64)logLγ=∑i=1nlogcF1^xi,F^2yi,where F1^x and F^2y denote the ML estimates of the parameters from the first step.The solution of the nonlinear equation (64) gives the MLE of γ.


## 5.4. Semiparametric Methods of Estimation

Two semiparametric methods are used to estimate the copula parameter in the copula models and are compared with the two methods of moments approaches which are the inversion Kendall'sτ and inversion of Spearman's ρ, respectively.


## 5.5. Methods of Moments

From the moment’s method of inversion of Kendall'sτ and the inversion of Spearman's ρ mentioned in Kojadinovic and Yan [32], we provide a brief details which are given as follows. Let c be a bivariate random sample from a cdf CγF1x,F2y, where F1 and F2 are continuous cdf’s and Cγ is an absolutely continuous copula such that γ∈O, where O is an open subset of R2. Furthermore, let R1,…,Rn are the vectors of ranks associated with x1,…,xn unless otherwise stated. In what follows, all vectors are row vectors. Moment’s approaches are based on the inversion of a consistent estimator of a moment of the copula Cγ. The two best-known moments, Spearman’s rho and Kendall’s tau, are, respectively, given by(65)ργ=12∫∫0,12uvdCγu,v−3,(66)τγ=4∫∫0,12Cγu,vdCγu,v−1.Consistent estimators of these two moments can be expressed as(67)ρn=12nn+1n−1∑i=1nRi,1Ri,2−3n+1n−1,(68)τn=4nn−1∑i=1n1xi,1≤xj,11xi,2≤xj,2−1.Ifρ and τ are one-to-one, consistent estimators of γ will be γn,ρ=ρ−1ρn and γn,τ=τ−1τn, respectively.It can be called inversion of Kendall'sτ and inversion of Spearman's ρ, respectively. For more information, see the work of Kojadinovic and Yan [32] and the references cited therein.As explained above, the moment’s method ofτ and ρ estimation for copula may be considered under the umbrella of semiparametric approach estimation.


## 5.6. Goodness of Fit Tests for Copula

We want to compare the empirical copula with the parametric estimator derived under the null hypothesis; for details, see the work of Fermanian [33]. Theory suggests a test if C is well-represented by a specific copula Cγ:(69)H0:C=CγVs.H1:C≠Cγ.Several well-known approaches are available in the literature; for example, see the work of Genest and Rémillard [34], or the fast multiplier approach, Genest et al. [35], and Kojadinovic et al. [36]. The goodness of fit tests based on the empirical process is given as(70)ℂnu,v=nCnu,v−Cγnu,v,where Cnu,v is the empirical copula of the data of X and YCnu,v=1/n∑i=1n1Ui,n≤u,Vi,n≤v,u,v∈0,1, and Ui,nandVi,n are pseudoobservations from C calculated from data as follows.Ui,n=R1i/n+1, Vi,n=R2i/n+1, and R1i and R2i are, respectively, the ranks of XiandYi.Here,Cnu,v is a consistent estimator, and θn is an estimator of γ obtained using the pseudoobservations. According to Genest et al. [35], the appropriate test statistics is the Cramer-von Miss and is defined as(71)Sn=∑i=1nCnUi,n,Vi,n−CγnUi,n,Vi,n2.


## 6. Bayesian Estimation

In this section, the Bayes estimates of the model parameters of the joint pdf (51) are obtained under the assumption that the random variables, Φ=θ,λ1,β1,λ2,β2,δ, have an independent gamma prior distributions with hyperparameter wkandmk,k=1,2,3,4,5,6, given by(72)fΦ;w,m=mkwkΓwkΦwk−1e−mkΦ,Φ>0,and ρ11,ρ12,ρ21, and ρ22 have a noninformative prior.By multiplying (23) or (44) with (72), the joint posterior density for the vector Φ, given the data, becomes(73)πΦ|x¯∝Lx|ΦfΦ;wk,sk.Marginal distributions ofΦ can be obtained by integrating out the (nuisance) hyperparameters. Thus, the Bayesian estimators of the parameters Φ under square error loss function can be calculated as follows:(74)Φ^∝∫Φ^∞ΦπΦ|x¯dΦ.The integrals in (74) cannot be obtained in a closed form, so the Markov chain Monte Carlo (MCMC) technique is used. In MCMC methods, the posterior distribution and the intractable integrals using simulated samples from the posterior distribution are obtained. Also, Gibbs sampling and the Metropolis–Hastings (MH) algorithm as a MCMC technique are used. For more details, see Metropolis et al. [37], Hastings [38], and Mohsin et al. [39]. The M-H algorithm considers that, to each iteration of the algorithm, an applicant value can be generated from a proposed distribution. Thus, the applicant value is allowed according to a sufficient approval probability. This technique assurances the convergence of the Markov chain for the target density. Finally, we can investigate that the advantage of 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.


### 6.1. Credible Intervals

In this section, a symmetric1001−ε% two-sided Bayes probability interval estimate of Φ, denoted by LΦ,UΦ, can be obtained as satisfying the following expression:(75)pLt¯<Φ<Ut¯=∫Lt¯Ut¯πθ,β,λ|t¯dΦ=1−ε,Since it is difficult to find the intervalLΦ and UΦ analytically, therefore, we apply suitable numerical techniques to solve this nonlinear equation.


## 6.1. Credible Intervals

In this section, a symmetric1001−ε% two-sided Bayes probability interval estimate of Φ, denoted by LΦ,UΦ, can be obtained as satisfying the following expression:(75)pLt¯<Φ<Ut¯=∫Lt¯Ut¯πθ,β,λ|t¯dΦ=1−ε,Since it is difficult to find the intervalLΦ and UΦ analytically, therefore, we apply suitable numerical techniques to solve this nonlinear equation.


## 7. Simulation Study

Here, a simulation study is conducted to see the efficacy of the proposed model in two cases, independent case and the dependent case. Monte Carlo simulation is done for comparison between maximum likelihood and Bayesian estimation methods, for estimating parameters of BEEM distribution usingR language such as (bbmle). The MLE estimation methods are done based on Newton–Raphson algorithm by using “maxLik” package. The Bayesian estimation by using Markov Chain Monte Carlo (MCMC) approach and Metropolis–Hastings (MH) algorithm is carried out by using R-program. Monte Carlo simulations are considered based on a data that is generated from a Gaussian copula by using copula package in R. We generate 10000 random samples of sizes n=50,100,and200, and different cases of actual values of the parameters are listed below:Case I:θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.2,ρ12=0.15,ρ21=0.05, and ρ22=0.1Case II:θ=2.8,λ1=2.5,β1=2.3,λ2=1.2,β2=1.4,δ=1.5,ρ11=0.5,ρ12=0.4,ρ21=0.2, and ρ22=0.3Case III:θ=1.8,λ1=1.5,β1=1.3,λ2=2.2,β2=2.5,δ=2,ρ11=0.2,ρ12=0.15,ρ21=0.05, and ρ22=0.1Case IV:θ=1.8,λ1=1.5,β1=1.3,λ2=2.2,β2=2.5,δ=2,ρ11=0.5,ρ12=0.4,ρ21=0.2, and ρ22=0.3We claim the best performance as the method which minimizes the mean squared error (MSE), bias of estimation, and length of confidence interval (L.CI) of the estimator. The two-sided confidence limit with confidence level γ=0.95 of the parameters is constructed as well.From the reported bias, MSE, and L.CI, in Tables3–6, it appears that the efficiency of the estimation under the Bayesian paradigm is quite evident. In particular, the MSE values in all parametric combinations tried in this article support in favor of this statement. We have the following observation on the simulation study as follows:(i)
As sample size (n) increases and for the same case (fixed actual parameters), the bias, MSE, and L.CI associated with the parameter estimates decrease for both methods of estimation.(ii)
In majority of the situations, for fixedn and for the correlation matrix parameters ρ11,ρ12,ρ21, and ρ22, the bias, MSE, and L.CI associated with them decrease under both the classical approach (MLE method) as well as under the Bayesian paradigm. As expected, when the sample size increases, a greater efficiency is observed in the overall estimation procedure (i.e., low bias and smaller values of the MSE).(iii)
From the reported bias, MSE, and L.CI, it appears that the efficiency of the estimation under the Bayesian paradigm is quite evident for both independent and dependent cases. In particular, the MSE values in all parametric combinations tried in this article support in favor of this statement. For reference, see Tables1–6 in the revised manuscript.(iv)
Credible interval(s) constructed appear to be the best indicator on the estimation of the model parameters as expected.(v)
Bias measures whether, over many replications, the estimator yields result that is correct on an average. A statistic is positively biased if it tends to overestimate the parameter; a statistic is negatively biased if it tends to underestimate the parameter. Negative bias means that the estimator is too small on an average compared to the true value.Table 5
MLE and Bayesian estimation with different sample sizes in dependent case (Case III).
θ=1.8,λ1=1.5,β1=1.3,λ2=2.2,β2=2.5,δ=2,ρ11=0.2,ρ12=0.15,ρ21=0.05, and ρ22=0.1nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.2789−0.2044−0.4097−0.1416−0.0470−0.28710.68150.78380.66380.8068MSE0.19680.13730.31260.20710.08980.29260.50990.61650.61760.6765L.CI1.35291.21231.49221.69641.16101.79790.83710.18371.65020.6260BayesBias0.00480.02820.2110−0.2575−0.46320.03920.18020.26370.12780.2249MSE0.03040.02750.08290.13010.28320.03230.08310.11510.05760.0965L.CI0.68310.64100.76810.99071.02790.68820.88250.83690.79720.8401100MLEBias−0.2655−0.1979−0.4074−0.1566−0.0580−0.28960.69330.78160.71660.8126MSE0.15890.11010.26510.17870.06150.23220.50800.61140.60430.6766L.CI1.16621.04431.23451.53970.94541.51060.64930.09011.41280.5019BayesBias0.00840.03870.1821−0.2361−0.33870.05040.08820.12050.05210.1102MSE0.02510.02630.06800.10580.15400.03190.04820.05280.04030.0508L.CI0.62040.61760.73210.87750.77720.67210.78820.76680.75980.7716200MLEBias−0.2401−0.1699−0.4124−0.1019−0.0545−0.27060.69200.78190.71750.8152MSE0.12060.07080.23750.13950.04050.18490.49450.60120.60140.6727L.CI0.98440.80351.01891.40950.76021.31050.49010.08421.39410.3565BayesBias−0.02050.01730.1272−0.2147−0.21060.05450.18020.26370.12780.2249MSE0.02030.01660.04490.08440.06040.02740.03740.05670.02260.0359L.CI0.55360.50100.66410.76730.49630.61270.03660.03710.03120.0327Table 6
MLE and Bayesian estimation with different sample sizes in dependent case (Case IV).
θ=1.8,λ1=1.5,β1=1.3,λ2=2.2,β2=2.5,δ=2,ρ11=0.5,ρ12=0.4,ρ21=0.2, and ρ22=0.3nθλ1β1λ2β2δρ11ρ12ρ21ρ2250MLEBias−0.2722−0.2106−0.3576−0.1175−0.0275−0.22510.40880.55470.52580.6412MSE0.17220.11530.23110.15480.06020.17930.17050.30810.41370.4138L.CI1.22871.04451.26001.47280.95651.40680.22890.08051.45250.2035BayesBias0.00910.02730.2060−0.2488−0.26230.03370.23940.32420.14400.2928MSE0.02660.03160.08000.11520.05930.03100.07140.12440.06330.1180L.CI0.63920.68850.76080.90560.95790.67770.46590.54480.80910.7048100MLEBias−0.2827−0.2017−0.3537−0.1082−0.0122−0.21400.40680.55260.55480.6393MSE0.14220.08800.20320.11930.04210.14310.16780.30560.41970.4115L.CI0.97890.85371.09581.28630.80351.22350.19030.05841.31200.2063BayesBias−0.00050.02520.1543−0.2224−0.34110.04870.13350.19670.07100.1565MSE0.02210.02430.05670.09150.05820.02980.04430.07420.04520.0664L.CI0.58270.60370.71110.80410.74560.64990.63790.73960.78560.8025200MLEBias−0.2756−0.1894−0.3556−0.0970−0.0193−0.19490.40450.55240.57700.6418MSE0.12810.07420.18400.09280.03170.12300.16590.30530.43600.4124L.CI0.89560.76810.94071.13260.69461.14330.18770.04881.25920.0907BayesBias−0.01400.00090.0934−0.1745−0.19330.04610.05720.09560.02620.0594MSE0.01830.01690.03500.05900.05350.02520.02840.04470.03240.0390L.CI0.52830.51040.63590.66280.49780.59520.62170.73990.69800.7383


## 8. Real Dataset (Motor Data)

The data represent the failure times of a parallel system constituted by two identical motors in days. These data are reported in Relia Soft (2003), whereX = (102, 84, 88, 156, 148, 139, 245, 235, 220, 207, 250, 212, 213, 220, 243, 300, 257, 263). Y = (65, 148, 202, 121, 123, 150, 156, 172, 192, 214, 212, 220, 265, 275, 300, 248, 330, 350). We fit at first the marginals of X and Y separately on the motor data. The MLE and Bayesian estimation for the BMEE distribution using dataset are shown in Table 7. The MLE of the parameters Kolmogorov–Smirnov distance (KS-d) and its p value for the marginals are listed in Table 8.Table 7
MLE and Bayesian estimation for the BMEE distribution using the dataset.
θλβδρMLECoefficient6.51950.01288.45800.01290.9791S.E0.07550.00420.11700.01100.1975BayesianCoefficient6.58820.01308.57620.01330.9576S.E0.06500.00160.07030.00150.0205Table 8
KS distance and itsp value and MLE.
θ^λ1^KS-dp valueLX11.47210.01510.26640.1553102.16545.58960.0029Y9.85070.01380.10010.9853103.87034.71740.0028The empirical cdf, the histogram of the pdf, PP plots, and QQ plots are displayed for the first variable in Figure3. The empirical cdf, the histogram of the pdf, PP plots, and QQ plots are displayed for the second variable in Figure 4. The EE distribution is fitted to real data using Kolmogorov–Smirnov goodness-of-fit test. The estimated cdf with empirical cdf, the histogram of the pdf, PP plots, and QQ plots are displayed for the first variable in Figure 3. The estimated cdf with empirical cdf, the histogram of the pdf, PP plots, and QQ plots are displayed for the second variable in Figure 4. These figures show that the tow variables are fitted for the marginal EE distribution. Also, Table 8 confirmed this conclusion by using Kolmogorov–Smirnov goodness-of-fit test, where p values are more than 0.05.Figure 3
The fitted pdf, cdf, PP plot, and QQ plot for EE distribution for the first variable.Figure 4
The fitted pdf, cdf, PP plot, and QQ plot for EE distribution for the second variable.In majority of the cases, one may observe that the estimates under the Bayesian method are preferable than the measures of MLE estimates as the standard error of the estimates are smaller in all estimates of the parameters. History plots, approximate marginal posterior density, and MCMC convergence ofθ,λ,δ,β, and ρ are represented in Figure 5.Figure 5
The Markov Chain Monte Carlo (MCMC) plots for motor data using the BMEE model.


## 9. Conclusions

In this paper, we have proposed and studied a new class of BEE whose marginals are EE distributions. The proposed class of distribution is constructed via two different types of mixture: (a) type I: starting with two independent EE distributions and (b) type II: using a bivariate Gaussian copula. Estimation of the model parameters for both types of MBEE distribution are conducted using classical (the method of moments and the method of maximum likelihood) and under the Bayesian paradigm using independent gamma priors. Since the joint distribution function and the joint density function are in closed forms, consequently, this distribution can be used in practice for nonnegative and positively correlated random variables. Since the maximum likelihood estimators of the unknown parameters cannot be obtained in the closed form, we consider the EM algorithm that works quite well, and it can be effectively used to compute the MLEs. Since the choice of hyperparameters for a prior in a Bayesian paradigm is of paramount importance, as a continuation of this work in future, we will be focusing on (including but not limited to) the following:(i)
Exploring various strategies (for example, matching conditional moments, or conditional percentiles information to be provided by our expert with the corresponding theoretical moments and percentiles and subsequently assuming, say, Euclidean distance) to estimate/guess best choice(s) of the hyperparameters for the priors.(ii)
In the current work, our prior choices are mostly conjugate in nature. However, in a real-life scenario, we might not have such an information on the prior always. Also, in the presence of more concrete information, one might consider a more precise prior for the model parameter(s), possibly a partially informative improper prior. It would be interesting to see the effect on the overall efficiency of the MCMC and Gibbs sampling in such a setting. We are currently working on it, and it will be reported somewhere else.



---

*Source: 5200979-2021-10-16.xml*