Multivariate Inverted Kumaraswamy Distribution: Derivation and Estimation

Industrial revolution leads to the manufacturing of multicomponent products; to guarantee the sufficiency of the product and consumer satisfaction, the producer has to study the lifetime of the products. This leads to the use of bivariate and multivariate lifetime distributions in reliability engineering. The most popular and applicable is Marshall–Olkin family of distributions. In this paper, a new bivariate lifetime distribution which is the bivariate inverted Kumaraswamy (BIK) distribution is found and its properties are illustrated. Estimation using both maximum likelihood and Bayesian approaches is accomplished. Using different selection criteria, it is found that BIK provides the best performance compared with other bivariate distributions like bivariate exponential and bivariate inverse Weibull distributions. As a generalization, the multivariate inverted Kumaraswamy (MIK) distribution is derived. Few studies have been conducted on the multivariate Marshall–Olkin lifetime distributions. To the best of our knowledge, none of them handle estimation process. In this paper, we developed an algorithm to show how to estimate the unknown parameters of MIK using both maximum likelihood and Bayesian approaches. This algorithm could be applied in estimating other Marshall–Olkin multivariate lifetime distributions.


Introduction
Global competition in combination with using higher manufacturing technologies results in producing two or multicomponent products. is led to the use of bivariate and multivariate distributions in reliability engineering. Different families of distributions were constructed. One of the most commonly used is the Marshall-Olkin (MO) family. It is widely used due to its flexibility in considering different situations of failures (i.e., the first component has lifetime smaller, greater, or equal to the lifetime of other components).
In the literature, several lifetime distributions were derived as members of the bivariate Marshall-Olkin family. Marshall and Olkin [1] presented a bivariate exponential distribution with exponential marginals and loss of memory property. Using the same strategy, Kundu and Dey [2], Bareto-Souza and Lemonte [3], Muhammed [4], and Alqallaf and Kundu [5] introduced the bivariate Weibull, bivariate Kumaraswamy, bivariate generalized Burr, and bivariate inverse generalized exponential distributions, respectively. Using the maximum instead of the minimum in the MO scheme, Kundu and Gupta [6,7] introduced the bivariate generalized exponential and bivariate proportional reversed hazard distributions, respectively. Moreover, Sarhan et al. [8] presented bivariate generalized linear failure rate distribution. Recently, Muhammed [9] introduced bivariate inverse Weibull (BIW) distribution.
Sometimes, the use of bivariate distributions may not be sufficient and there exists a need for multivariate distributions. For example, in air fighter jets, the natural lifetime since being manufactured and the flying time since being put into service are recorded. Studying the reliability of the air fighter jets using only two variables may not be good enough. One should take into consideration the lifetime of the engine, the wing, and the fuselage. is leads to the use of multivariate distribution. For more details, see Li and Li [10].
ere is no much work performed in the multivariate case. Sarhan et al. [8] and Kundu and Gupta [11] derived the multivariate generalized linear failure rate and multivariate inverse Weibull distributions, respectively. To the best of our knowledge, there is no work dealing with estimating the unknown parameters for multivariate Marshall-Olkin family.
Several authors tackled the estimation problem for bivariate MO distributions. For example, Kundu and Gupta [6], Muhammed [9], Aly et al. [12], Eliwa and El-Morshedy [13], El-Morshedy et al. [14], and Sarhan [4,5,15] estimated the unknown parameters using maximum likelihood approach for different bivariate lifetime distributions. On the other hand, Hanagal and Ahmadi [16], Kundu and Gupta [17], and Lin et al. [11,[13][14][15]18] applied Bayesian approach for estimating certain bivariate lifetime distributions. e univariate inverted Kumaraswamy (IK) distribution has several applications in different fields (see Abd Al-Fattah et al. [19] and Abdul Hammed et al. [20]), for example, in medical research, life testing problems, and stress-strength analysis. Also, in reliability and biological studies, IK distribution may be used to model failure rates. Due to its expected wide applicability, we are interested in deriving bivariate inverted Kumaraswamy (BIK) distribution. BIK could be applied in different fields like sports, engineering, and medicine as will be explained using three different real datasets. We expect better performance of BIK than other bivariate distributions. No one has derived the distribution before or found its mathematical properties. e main purpose of this paper is to introduce BIK as a new Marshall-Olkin bivariate distribution in order to be applied efficiently in several fields. As a generalization, the multivariate inverted Kumaraswamy (MIK) distribution is derived. To the best of our knowledge, there is no work dealing with estimating the unknown parameters for multivariate Marshall-Olkin family. Here, estimation of MIK parameters is found using both maximum likelihood and Bayesian approaches.
is work could be applied to all Marshall-Olkin multivariate distributions. e paper is organized as follows. In Section 2, the bivariate inverted Kumaraswamy distribution is derived, and the cumulative distribution function and probability density function are presented. Also, the marginals and conditional distributions of the proposed model are obtained. Moreover, the product moments and the moment generating function are derived. In Section 3, the maximum likelihood estimators of the model parameters, asymptotic Fisher information matrix, and Bayesian estimators are obtained. Multivariate inverted Kumaraswamy distribution and its properties are illustrated in Section 4. e maximum likelihood and Bayesian estimators of the parameters under multivariate case are obtained in Section 5. Numerical analysis using both simulation, and real datasets are presented in Section 6. Finally, the paper is concluded in Section 7.

Bivariate Inverted Kumaraswamy Distribution
In this section, we will derive the bivariate inverted Kumaraswamy distribution as a new member in the MO family. Its properties such as the marginal and conditional distributions, joint moment generating function, and product moments are studied.

Derivation of the Bivariate Inverted Kumaraswamy
Distribution. e probability density function (pdf ) and the cumulative distribution function (cdf ) of the univariate inverted Kumaraswamy distribution (IK), respectively, are as follows (for more details, see [19]): Assume that U i , i � 1, 2, 3, are three random variables, such that U i follows IK (β i , α). Define X 1 � max(U 1 , U 2 ) and X 2 � max(U 2 , U 3 ). Hence, the bivariate vector (X 1 , X 2 ) follows BIK with shape parameters β i , α, i � 1, 2, 3. e joint cdf of (X 1 , X 2 ) has the following form: where z � min(x 1 , x 2 ). e joint cdf can also be written as follows: where where z � min x 1 , x 2 , F s (.) and F a (., .) are the singular and absolute parts, respectively.
□ Corollary 1. e joint pdf of (X 1 , X 2 ) can be written as follows: where f a x 1 , e absolutely continuous part of BIK can be unimodal depending on the values of β 1 , β 2 , β 3 , and α.
e respective modes are

Properties of Bivariate Inverted Kumaraswamy
Distribution. In this section, we illustrate different properties of BIK distribution. We provide marginal, conditional distributions, joint moment generating function, and product moments.
, then the product moments of X 1 and X 2 are Proof. See Appendix A.1.

Estimation of Bivariate Inverted Kumaraswamy Distribution
In this section, we estimate the unknown parameters using both maximum likelihood and Bayesian approaches.
e posterior pdf of parameters can be obtained as follows: where P(θ) is the prior distribution.
Pena and Gupta [21] considered Bayesian estimation of the parameters for Marshall-Olkin bivariate exponential distribution (BVE), in series and parallel systems. ey obtained posterior mode using gamma Dirichlet distribution as prior distribution. Angali et al. [22] considered Bayesian estimation for BVE using gamma prior.
Similar to [22], we considered a gamma prior distribution with the following pdf: where i � 1, 2, 3 and a i , b i , a 4 , b 4 are the hyperparameters. e posterior pdf has the following form: Mathematical Problems in Engineering It is observed that Bayesian estimators under square error loss function cannot be obtained in explicit forms. erefore, we obtain the posterior mean using MCMC technique which is illustrated in Section 6.

Multivariate Inverted Kumaraswamy Distribution
In literature, there is no much work that dealt with multivariate MO distributions. Sarhan et al. [8] derived the multivariate generalized linear failure rate distribution. Kundu and Gupta [11] derived the multivariate inverse Weibull distribution. Here, we will introduce a new multivariate MO distribution, which is the multivariate inverted Kumaraswamy (MIK) distribution. It is a generalization of the BIK considered in Section 2. We expect that MIK will be of great importance for several applied fields. In the following two sections, the derivation of MIK is explained and some of its properties are studied.
Proof. We prove by generalizing the same method illustrated in Section 2. Similar to the bivariate case, the MIK distribution can be written as, for m > 1, where 0 < k < 1 and F a and F s denote the absolute and continuous parts, respectively. e corresponding pdf can also be written as en, for where Since 6 Mathematical Problems in Engineering and for all can be written as where f I l is a pdf with respect to (m − l + 1) dimensional Lebesgue measure on the hyperplane e exact meaning of f X (x) is as follows.
For any Borel measurable set B ∈ R m , where B I l � B ∩ A I l is the projection of set B onto (m − l + 1) dimensional hyperplane A I l . Now, we provide k I l and f I l (x). Note that if x ϵA I l , then x has the following form: (29) For a given x ∈ R m , we define g I l from the (m − l + 1) dimensional hyperplane A I l to R as follows: if x i > x * for i ∈ I − I l and zero otherwise. Similar to k, we can obtain k I l as follows: x j2

Properties of MIK.
In this section, we will get the marginal and conditional distributions of MIK.
Proof. (a, c) B taking the limit for the joint cdf. (b, d) Could be directly obtained from the definition of the multivariate inverted Kumaraswamy distribution.

Estimation of the Multivariate Inverted Kumaraswamy Distribution
Although estimating the unknown parameters of a certain multivariate MO distribution is very important, no one in the literature was interested in it. erefore, in this section, we will consider the process of estimation for MIK parameters. e proposed techniques could be applied for any multivariate MO distribution. Here, we will apply both maximum likelihood and Bayesian approaches.

Using Maximum Likelihood Approach.
In this section, for simplification, we consider the case when we have three random variables X 1 , X 2 , and X 3 . Applying Proposition 6, we have the following cdf: e cdf can be rewritten as e pdf can be obtained by taking the derivative, except for f 10 where we should take into consideration that the sum of all probabilities equals one.
; the problem is to find the ML estimators of the unknown parameters. Consider the following notation: where |I j | denotes the cardinality of the set I j , for j � 1, 2, . . . , 13. e log likelihood function of the sample of size n is given by 8 Mathematical Problems in Engineering where θ � (α, β 1 , β 2 , β 3 , β 4 ).
It is seen that the ML estimates could not be obtained in explicit forms, and hence we need to use numerical analysis to obtain them.

Using Bayesian
is the vector of unknown parameters. e posterior pdf can be obtained as follows: where P(θ) is the prior distribution. We considered a gamma prior distribution with the following pdf: where i � 1, 2, 3, 4 and a i , b i , a 4 , b 4 are the hyperparameters. e posterior pdf has the following form: Mathematical Problems in Engineering As in the bivariate case, Bayesian estimation will be obtained numerically using MCMC which is illustrated in the next section.

Numerical Analysis
In this section, both simulation and MCMC techniques are carried out to investigate the performance of the derived BIK and MIK distributions. e estimation is performed using both ML and Bayesian approaches. ree real datasets are analyzed in case of BIK and another one for the case of MIK.

For BIK Distribution.
In this section, we perform a simulation study to get the estimates of the unknown parameters of BIK distribution. Also, three real datasets are analyzed.

A Simulation Study.
Here, we present an algorithm to generate BIK distribution (Algorithm 1). To perform a simulation study, we first need to select initial values for the parameters. Here, the following eight different populations are considered: e parameters are selected to cover different shapes of the distribution. It can be seen from Figure 1 that (i) For cases 1 and 2, the surface plot of the absolutely continuous part of the joint probability density function is decreasing. (ii) For cases 3 to 8, the surface plot of the absolutely continuous part of the joint probability density function is increasing till it reaches the mode; then, it is decreasing.
e maximum likelihood estimates of the model parameters are obtained by maximizing the log likelihood function given by equation (15). Monte Carlo simulation is performed using R package with 1000 replications and three different sample sizes n � 30, 50 , and 70 and eight different populations.
Absolute bias (ABias), mean square error (MSE), confidence width (CW), and coverage probability (CP) are obtained e following algorithm is to generate (X 1 , X 2 ) from BIK distribution.
Step 2: compute Step 4: obtain ALGORITHM 1: An Algorithm to Generate BIK Distribution.   (2) Bayesian Approach. Using Bayesian approach under square error loss function, the Bayesian estimator is the posterior mean. However, it is hard to obtain the posterior mean theoretically as we have four parameters to estimate. One can use Markov Chain Monte Carlo (MCMC) simulation method to obtain it numerically.
e MCMC method uses simulation techniques to obtain a Markov sequence such that it has a limiting distribution. In the Bayesian approach, the limiting distribution is the posterior pdf as it includes all needed information about the parameters θ.
Here, the MCMC method can be used to set up a Markov chain of parameters θ with distribution P(θ | (X 1 , X 2 )). e mean of the sequence can be considered as the posterior mean. To perform MCMC, we used both R and WinBugs packages. Gamma prior is used with the same three sample sizes and eight populations used in ML approach. e R package with 1000 replications is used, and for each replication, WinBugs is used with 1000 replications to generate the sequence of Markov chain.
We used the Geweke test to examine the convergence of the generated Markov chain sequence. Geweke statistic (z n ) converges to normal distribution for large sample sizes. Hence, large absolute values of z n are considered as a rejection for convergence. Only those converged sequences are used in the analysis. For more details about the Geweke test, see [16].
ABias, MSE, confidence length (CL), and CP are obtained and presented in Table 1. e numerical steps and the corresponding equations are explained in detail in Appendix B.1.
From Table 1, it can be seen that under different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. is indicates that both Bayesian and ML approaches work efficiently in estimating the parameters of BIK.
Comparing ML and Bayesian estimates, it is found that Bayesian estimates have less than or equal mean square error (MSE) than ML ones.
is is clear from Figure 2.
Also, it can be seen that as the sample size (n) increases, the ABias, MSE, CW, and CL decrease for both ML and Bayesian as seen from Figure 3. Moreover, it can be seen that for most cases, the CP is around 0.95.

Real Datasets.
Here, we analyze three real datasets to show the applicability of BIK in several fields like sports, engineering, and medicine.
(1) Football Data. e dataset has been obtained from Meintanis [23]; he used the bivariate MO exponential distribution (BE) to analyze the data. e data are about football (soccer) where at least one goal scored by the home team and at least one goal scored directly from a penalty kick, a foul kick, or any other direct kick (all of them together will be called as kick goal) by any team have been considered. Here, the variables are the time in minutes of the first kick goal scored by any team (X 1 ) and the time of the first goal of any type scored by the home team (X 2 ). e bivariate dataset has the following structure: X 1 < X 2 , X 1 � X 2 , and X 1 > X 2 . Since, X 1 � X 2 has a positive probability, we need a singular distribution to analyze this dataset. Here, we analyze the data using BIK distribution defined by equation (4). All the data points have been divided by 100. is is not going to make any difference in the statistical inference.
First, before analyzing the data, we fit inverted Kumaraswamy distribution to X 1 , X 2 and max (X 1 , X 2 ). To guess the initial values for the parameters of BIK model, the MLEs of the shape parameters (α, β) of the respective inverted Kumaraswamy distribution for X 1 , X 2 and max (X 1 , X 2 ) are obtained. To check the modelʼs fitness, we first need to illustrate goodness-of-fit tests.
Goodness-of-fit (GOF) tests are hypothesis tests regarding the distribution of some random variable (X) in a population. e objective of applying GOF tests is to measure how well the data agree with a given distribution as its population. For example, if we want to examine if the random variable (X) follows distribution F 0 (x), then the null hypothesis is One approach for applying GOF tests is based on the empirical distribution function (EDF) (F n (x)) which is defined as follows: is approach is based on defining a statistic to measure the discrepancy between F n (x) and F 0 (x). e most used statistics are modified Cramér-von Mises statistic (W * ) and Anderson-Darling statistic (A * ) which have the following formulas: where x (i) is the value of the i th order statistics in the sample. Large values of test statistics (or small corresponding p value) lead to the rejection of the null hypothesis. For more details about GOF tests, see DʼAgostnio and Stephens [24].
Here, we apply GOF tests in order to see whether the fit based on univariate inverted Kumaraswamy distributions is reasonable for this dataset. We computed the modified Cramér-von Mises statistic (W * ) and Anderson-Darling statistic (A * ). e values of these statistics and the corresponding p values (in brackets) for X 1 , X 2 , and max (X 1 , X 2 ) are illustrated in Table 2.
Based on the values of these statistics and the corresponding p values, the inverted Kumaraswamy distribution cannot be rejected for modeling the marginals and the maximum. In Table 3, the ML estimates and the posterior mean using gamma priors are obtained for the parameters of BIK. Also, credible interval length and confidence width are illustrated. Now, we try to compare the performance of BIK, bivariate exponential (BE), and bivariate inverted Weibull (BIW) to fit this dataset. To select between models, several information criteria (IC) were presented; the main idea behind IC is to afford a balance between good fit and complexity of the model as follows:  where k is the penalized term and p refers to number of parameters in the model. e most commonly used IC in model selections are Akaike information criteria (AIC), Bayesian information criteria (BIC), the consistent Akaike information criteria (CAIC), and Hannan-Quinn information criteria (HQIC). Each IC has a different penalty term illustrated in the first row of Table 4.
By analyzing equation (45), we can see that the first term (− 2log(L)) tends to decrease when the model provides good fit. But, the second term tends to increase as the number of parameters in the model increases.
e model with the lowest IC is the best (for more details about IC, see Vrieze [25]). e log likelihood value, AIC, BIC, CAIC, and HQIC are represented in Table 4. All of the criteria suggest that BIK provides the best fit compared with BE and BIW models.
is may show the importance of BIK.
(2) Motor Data. e dataset has been obtained from [26]. e data are about failure time in days for a parallel system containing two motors. e variables are time to failure for the first motor (X 1 ) and time to failure for the second motor (X 2 ). All data points have been divided by 1000. We applied GOF tests on IK, IW, and E distributions. From Table 5, based on the values of (W * ), (A * ), and the corresponding p values, only IK distribution can be used for modeling the marginals and the maximum. Hence, only BIK can be used for modeling this dataset.
ML estimates and the posterior mean using gamma priors are obtained for the parameters of BIK. Also, credible interval length and confidence width are illustrated in Table 6.
(3) Diabetic Retinopathy Data. e dataset has been obtained from [27]. e data are used by the National Eye Institute to study the effect of laser treatment on the blindness in patients with diabetic retinopathy. At the beginning of clinical trial, for each patient, one eye is randomly selected for laser treatment. e variables are time to failure for treated eye (X 1 ) and time to failure for untreated eye (X 2 ). All data points have been divided by 1000. We applied GOF on IK, IW, and E distributions. From Table 7, it can be seen that only IK distribution can be used for modeling the marginals and the maximum. Hence, only BIK can be used for modeling this dataset.
In Table 8, ML estimates and the posterior mean using gamma priors are obtained for the parameters of BIK. Also, credible interval length and confidence width are illustrated.
From these three datasets, we can conclude that the derived BIK distribution will be of great importance.

For MIK Distribution.
In this section, we present numerical results of estimation using a simulation study and a real dataset.

A Simulation Study.
Here, we present an algorithm to generate MIK distribution. Also, we illustrate the simulation results for both ML and Bayesian approaches (Algorithm 2).
(1) Maximum Likelihood Approach. To obtain the maximum likelihood estimates, a Monte Carlo simulation is performed    Table 3. * * Based on the estimates obtained by Meintanis [23]. * * * Based on the estimates presented by Muhammed [9].    e following algorithm is to generate (X 1 , X 2 , X 3 ) from MIK distribution.
ALGORITHM 2: An algorithm to generate MIK distribution. Table 9: e results of MLEs and Bayesian estimates for MIK.

MLE
Bayesian To check the behavior of the estimates, ABias, MSE, CW, and CP are computed in Table 9. e algorithm is explained in detail in Appendix B.2.
(2) Bayesian Approach. Similar to the bivariate case, MCMC simulation is used to obtain the posterior mean numerically. Absolute bias (ABias), mean square error (MSE), confidence length (CL), and coverage probability (CP) are obtained and presented in Table 9. e algorithm is explained in detail in Appendix B.2.
From Table 9, it can be seen that for different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. is indicates that both Bayesian and ML approaches work efficiently in estimating the parameters of MIK.
Comparing ML and Bayesian estimates, it is found that Bayesian estimates have less than or equal mean square error (MSE) than ML ones as seen from Figure 4. Also, for the majority of cases, Bayesian estimates have smaller ABias than ML ones.
Also, it can be seen that as the sample size (n) increases, the ABias, MSE, CW, and CL decrease for both ML and Bayesian as seen from Figure 5. Moreover, it can be seen that for most cases, the CP is around 0.95.

A Real Dataset.
Here, we analyze a real dataset to show the applicability of MIK. e dataset has been obtained from Bland and Altman [28]. It represents a set of systolic blood pressure measurement for 85 patients made by a semiautomatic blood pressure monitor; three readings were made for each patient. e variables are as follows: X 1 : first systolic blood pressure measurement. X 2 : second systolic blood pressure measurement. X 3 : third systolic blood pressure measurement.
All data points have been divided by 1000. is is not going to make any difference in the statistical inference. We applied GOF tests in order to check if the fit based on IK, IW, and E distributions is reasonable in this case. We computed the modified Cramér-von Mises statistic (W * ) and Anderson-Darling statistic (A * ). e values of these statistics and the corresponding p values (in brackets) for X 1 , X 2 , X 3 , and max(X 1 , X 2 , X 3 ) are illustrated in Table 10.
Based on the values of these statistics and the corresponding p values, only IK distribution can be used for modeling the marginals and the maximum. Hence, only   MIK can be used for modeling these data. In Table 11, the ML estimates and the posterior mean using gamma priors are obtained for the parameters of MIK (β 4 is considered known for simplicity). Also, credible interval length and confidence width are illustrated.
It can be seen that [1 − (1 + z) − α ] β 1 +β 2 +β 3 is the singular part as its second partial derivative is zero when x 1 ≠ x 2 . us, P( ) is the absolutely continuous part as its mixed partial derivative is a density function.
(c) Conditional distribution of X 1 given X 2 ≤ x 2 is given by that is, substituting for f(x 1 , x 2 ) by the corresponding formula, and then using change of variable technique and the following facts, the formula is derived.
Step 2: the maximum likelihood estimates (MLEs) and the corresponding variance (Var) are obtained and stored.
Step 4: the three previous steps are repeated 1000 times.
Step 5: for converged datasets (this is obtained from R using optim function), the average of MLE, variance, and CI is obtained as follows: CI i number of converged datasets . (B.1) Step 6: absolute bias (ABias), mean square error (MSE), and confidence width (CW) are obtained as follows: ABias � |AMLE − intial value|, MSE � AVar + Bias 2 , CW � upper limit of ACI − lower limit of ACI.

(B.2)
Step 7: the coverage probability (CP) is obtained as follows: CP � number of times the intial value falls inside the converged samples' CI number of converged samples . (B.3)

B.1.2. MCMC Simulation Algorithm for BIK.
e following algorithm is used to perform Markov Chain Monte Carlo (MCMC) simulation study using R and WinBugs packages and get Bayesian estimates for BIK distribution.
Step 1: using R package: n independent samples (X 1 , X 2 ) from BIK distribution are generated using the same procedure in Appendix B.1.1.
Step 2: the generated dataset is sent to WinBugs package and gamma priors are defined.
Step 3: WinBugs is used with 1000 replications to generate the sequence of Markov chain.
Step 5: results are sent back to R package and stored.
Step 6: the five previous steps are repeated 1000 times.
Step 7: for converged datasets (i.e., Geweke test results <1.96, it is obtained using coda function in R package), the average of posterior mean (PM), variance (Var), and credible interval (CI) is obtained as follows:
Step 9: the coverage probability (CP) is the same as step 7 in Appendix B.1.1.

B.2. Algorithms for MIK
B.2.1. Monte Carlo Simulation Algorithm for MIK. e following algorithm is used to perform Monte Carlo simulation study using R package and get ML estimates for MIK distribution.

B.2.2. MCMC Simulation Algorithm for MIK.
e following algorithm is used to perform Markov Chain Monte Carlo (MCMC) simulation study using R and WinBugs packages and get Bayesian estimates for MIK distribution.
Step 1: using R package: n independent samples (X 1 , X 2 , X 3 ) from MIK distribution are generated using the same procedure in Appendix B.2.1. Steps 2 till 9 are the same as in Appendix B.1.2.

Conflicts of Interest
e authors declare that they have no conflicts of interest.