# Multivariate Inverted Kumaraswamy Distribution: Derivation and Estimation

Authors: Hanan M. Aly; Ola A. Abuelamayem
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. Introduction

Global competition in combination with using higher manufacturing technologies results in producing two or multicomponent products. This 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. This leads to the use of multivariate distribution. For more details, see Li and Li [10]. There 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–15, 18] applied Bayesian approach for estimating certain bivariate lifetime distributions.The 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.The 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. This work could be applied to all Marshall–Olkin multivariate distributions.The paper is organized as follows. In Section2, 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. The 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.


## 2. 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.


### 2.1. Derivation of the Bivariate Inverted Kumaraswamy Distribution

The 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]):(1)fIKx,β,α=αβ1+x−α+11−1+x−αβ−1,x>0,α,β>0,(2)FIKx,β,α=1−1+x−αβ,x>0,α,β>0.Assume thatUi,i=1,2,3, are three random variables, such that Ui follows IK (βi,α). Define X1=maxU1,U2 and X2=maxU2,U3. Hence, the bivariate vector (X1, X2) follows BIK with shape parameters βi,α,i=1,2,3. The joint cdf of (X1,X2) has the following form:(3)FX1,X2x1,x2=PX1≤x1,X2≤x2=PmaxU1,U2≤x1,maxU2,U3≤x2,=PU1≤x1,U2≤x2,U3≤minx1,x2,=PU1≤x1,U2≤x2,U3≤z,=FIKx1,β1,αFIKx2,β2,αFIKz,β3,α,where z=minx1,x2. The joint cdf can also be written as follows:(4)FX1,X2x1,x2=FIKx1,β1+β3,αFIKx2,β2,α,ifx1<x2,FIKx1,β1,αFIKx2,β2+β3,α,ifx2<x1,FIKx,β1+β2+β3,α,ifx1=x2=x,where FIK is as illustrated in equation (2).Proposition 1.
The joint pdf ofX1,X2 has the following form:(5)fX1,X2x1,x2=f1x1,x2,ifx1<x2,f2x1,x2,ifx2<x1,f3x,ifx1=x2=x,where(6)f1x1,x2=fIKx1,β1+β3,αfIKx2,β2,α,f2x1,x2=fIKx1,β1,αfIKx2,β2+β3,α,f3x=β3β1+β2+β3fIKx,β1+β2+β3,α,in which fIK is as illustrated in equation (1).Proof.
See AppendixA.1.
The joint pdf ofX1,X2 can be expressed as a mixture of absolutely continuous part and singular part as follows.Proposition 2.
IfX1,X2 follows BIK (β1,β2,β3,α), then(7)FX1,X2x1,x2=β1+β2β1+β2+β3Fax1,x2+β3β1+β2+β3Fsz,where(8)Fax1,x2=β1+β2+β3β1+β21−1+x1−αβ11−1+x2−αβ21−1+z−αβ3−β3β1+β2+β31−1+z−αβ1+β2+β3,Fsz=1−1+z−αβ1+β2+β3,where z=minx1,x2,Fs. and Fa.,. are the singular and absolute parts, respectively.Proof.
See AppendixA.1.Corollary 1.
The joint pdf ofX1,X2 can be written as follows:(9)fX1,X2x1,x2=β1+β2β1+β2+β3fax1,x2+β3β1+β2+β3fsz,where(10)fax1,x2=β1+β2+β3β1+β2f1x1,x2ifx1<x2f2x1,x2ifx2<x1,fsz=fIKz,β1+β2+β3.The absolutely continuous part of BIK can be unimodal depending on the values ofβ1,β2,β3, and α. fax1,x2 is unimodal if β1+β3<β2 (under the case x1<x2) or β2+β3<β1 (under the case x1>x2). The respective modes are(11)αβ1+β3+1α+11/α−1,αβ2+1α+11/α−1,αβ1+1α+11/α−1,αβ2+β3+1α+11/α−1.


### 2.2. 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.Proposition 3 (Marginal and conditional distributions).
IfX1,X2 follows BIK (β1,β2,β3,α), then(a)
X1 follows IK (β1+β3,α) and X2 follows IK (β2+β3,α).(b)
max {X1,X2} follows IK (β1+β2+β3,α).(c)
FX1X2≤x2x1=1−1+x1−αβ1+β31−1+x2−α−β3,ifx1≤x2,1−1+x1−αβ1,ifx2<x1.Proof.
See AppendixA.1.Proposition 4 (Moment generating function).
IfX1,X2 follows BIK (β1,β2,β3,α), then the joint moment generating function of X1 and X2 is given by Mt1,t2=l1+l2+l3, where(12)l1=β2β1+β3∑i=0∞t2ii!∑k2=0iik2−1i−k2∑j=0∞t1ij!∑k1=0jjk1−1i−k1XB1−k2α,β2B1−k1α,β1+β3−αα−k1B2−k2+k1α,β2X3F22−k2+k1α,1−k1α,1−β1−β3;2−k1α;2+β2−k2+k1α;1,l2=β1β2+β3∑j=0∞t1jj!∑k1=0jjk1−1j−k1∑i=0∞t2ii!∑k2=0iik2−1i−k2XB1−k1α,β1B1−k2α,β2+β3−αα−k2B2−k1+k2α,β1X3F22−k1+k2α,1−k2α,1−β2−β3;2−k2α;2+β1−k1+k2α;1,l3=β1∑i=0∞t1+t2ii!∑k=0iik−1i−kB1−kα,β1+β2+β3,where Fpqb1,…,bq;c1,…,cq;u=∑k=0∞b1k,…,bpk/c1k,…,cqkuk/k!, bi=Γb+i/Γb,i=1,…,p,cj=Γc+j/Γc,j=1,…,q, p and q are nonnegative integers, and B (.,.) is beta function.Proof.
See AppendixA.1.Proposition 5 (Product moments).
IfX1,X2 follows BIK (β1,β2,β3,α), then the product moments of X1 and X2 are given by EX1sX2r=J1+J2+J3, where(13)J1=∑i=0r−1r−iri∑j=0s−1s−jsjβ2β1+β3B1−jα,β1+β3B1−iα,β2−αα−jB2−i+jα,β23F22−i+jα,1−jα,1−β1−β3;2−jα,2−i+jα+β2;1,J2=∑i=0r−1r−iri∑j=0s−1s−jsjβ1β2+β3B1−iα,β2+β3B1−jα,β1−αα−iB2−i+jα,β1F232−j+iα,1−iα,1−β2−β3;2−iα,2−j+iα+β1;1,J3=∑k=0r+s−1r+s−kr+skβ3B1−kα,β1+β2+β3B1−kα,β2+β3.Proof.
See AppendixA.1.


## 2.1. Derivation of the Bivariate Inverted Kumaraswamy Distribution

The 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]):(1)fIKx,β,α=αβ1+x−α+11−1+x−αβ−1,x>0,α,β>0,(2)FIKx,β,α=1−1+x−αβ,x>0,α,β>0.Assume thatUi,i=1,2,3, are three random variables, such that Ui follows IK (βi,α). Define X1=maxU1,U2 and X2=maxU2,U3. Hence, the bivariate vector (X1, X2) follows BIK with shape parameters βi,α,i=1,2,3. The joint cdf of (X1,X2) has the following form:(3)FX1,X2x1,x2=PX1≤x1,X2≤x2=PmaxU1,U2≤x1,maxU2,U3≤x2,=PU1≤x1,U2≤x2,U3≤minx1,x2,=PU1≤x1,U2≤x2,U3≤z,=FIKx1,β1,αFIKx2,β2,αFIKz,β3,α,where z=minx1,x2. The joint cdf can also be written as follows:(4)FX1,X2x1,x2=FIKx1,β1+β3,αFIKx2,β2,α,ifx1<x2,FIKx1,β1,αFIKx2,β2+β3,α,ifx2<x1,FIKx,β1+β2+β3,α,ifx1=x2=x,where FIK is as illustrated in equation (2).Proposition 1.
The joint pdf ofX1,X2 has the following form:(5)fX1,X2x1,x2=f1x1,x2,ifx1<x2,f2x1,x2,ifx2<x1,f3x,ifx1=x2=x,where(6)f1x1,x2=fIKx1,β1+β3,αfIKx2,β2,α,f2x1,x2=fIKx1,β1,αfIKx2,β2+β3,α,f3x=β3β1+β2+β3fIKx,β1+β2+β3,α,in which fIK is as illustrated in equation (1).Proof.
See AppendixA.1.
The joint pdf ofX1,X2 can be expressed as a mixture of absolutely continuous part and singular part as follows.Proposition 2.
IfX1,X2 follows BIK (β1,β2,β3,α), then(7)FX1,X2x1,x2=β1+β2β1+β2+β3Fax1,x2+β3β1+β2+β3Fsz,where(8)Fax1,x2=β1+β2+β3β1+β21−1+x1−αβ11−1+x2−αβ21−1+z−αβ3−β3β1+β2+β31−1+z−αβ1+β2+β3,Fsz=1−1+z−αβ1+β2+β3,where z=minx1,x2,Fs. and Fa.,. are the singular and absolute parts, respectively.Proof.
See AppendixA.1.Corollary 1.
The joint pdf ofX1,X2 can be written as follows:(9)fX1,X2x1,x2=β1+β2β1+β2+β3fax1,x2+β3β1+β2+β3fsz,where(10)fax1,x2=β1+β2+β3β1+β2f1x1,x2ifx1<x2f2x1,x2ifx2<x1,fsz=fIKz,β1+β2+β3.The absolutely continuous part of BIK can be unimodal depending on the values ofβ1,β2,β3, and α. fax1,x2 is unimodal if β1+β3<β2 (under the case x1<x2) or β2+β3<β1 (under the case x1>x2). The respective modes are(11)αβ1+β3+1α+11/α−1,αβ2+1α+11/α−1,αβ1+1α+11/α−1,αβ2+β3+1α+11/α−1.


## 2.2. 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.Proposition 3 (Marginal and conditional distributions).
IfX1,X2 follows BIK (β1,β2,β3,α), then(a)
X1 follows IK (β1+β3,α) and X2 follows IK (β2+β3,α).(b)
max {X1,X2} follows IK (β1+β2+β3,α).(c)
FX1X2≤x2x1=1−1+x1−αβ1+β31−1+x2−α−β3,ifx1≤x2,1−1+x1−αβ1,ifx2<x1.Proof.
See AppendixA.1.Proposition 4 (Moment generating function).
IfX1,X2 follows BIK (β1,β2,β3,α), then the joint moment generating function of X1 and X2 is given by Mt1,t2=l1+l2+l3, where(12)l1=β2β1+β3∑i=0∞t2ii!∑k2=0iik2−1i−k2∑j=0∞t1ij!∑k1=0jjk1−1i−k1XB1−k2α,β2B1−k1α,β1+β3−αα−k1B2−k2+k1α,β2X3F22−k2+k1α,1−k1α,1−β1−β3;2−k1α;2+β2−k2+k1α;1,l2=β1β2+β3∑j=0∞t1jj!∑k1=0jjk1−1j−k1∑i=0∞t2ii!∑k2=0iik2−1i−k2XB1−k1α,β1B1−k2α,β2+β3−αα−k2B2−k1+k2α,β1X3F22−k1+k2α,1−k2α,1−β2−β3;2−k2α;2+β1−k1+k2α;1,l3=β1∑i=0∞t1+t2ii!∑k=0iik−1i−kB1−kα,β1+β2+β3,where Fpqb1,…,bq;c1,…,cq;u=∑k=0∞b1k,…,bpk/c1k,…,cqkuk/k!, bi=Γb+i/Γb,i=1,…,p,cj=Γc+j/Γc,j=1,…,q, p and q are nonnegative integers, and B (.,.) is beta function.Proof.
See AppendixA.1.Proposition 5 (Product moments).
IfX1,X2 follows BIK (β1,β2,β3,α), then the product moments of X1 and X2 are given by EX1sX2r=J1+J2+J3, where(13)J1=∑i=0r−1r−iri∑j=0s−1s−jsjβ2β1+β3B1−jα,β1+β3B1−iα,β2−αα−jB2−i+jα,β23F22−i+jα,1−jα,1−β1−β3;2−jα,2−i+jα+β2;1,J2=∑i=0r−1r−iri∑j=0s−1s−jsjβ1β2+β3B1−iα,β2+β3B1−jα,β1−αα−iB2−i+jα,β1F232−j+iα,1−iα,1−β2−β3;2−iα,2−j+iα+β1;1,J3=∑k=0r+s−1r+s−kr+skβ3B1−kα,β1+β2+β3B1−kα,β2+β3.Proof.
See AppendixA.1.


## 3. Estimation of Bivariate Inverted Kumaraswamy Distribution

In this section, we estimate the unknown parameters using both maximum likelihood and Bayesian approaches.


### 3.1. Using Maximum Likelihood Approach

In this section, we derive the maximum likelihood (ML) estimators of the unknown parameters of BIK distribution. Suppose {(x11,x21), …, (x1n,x2n)} is a random sample of size n from BIK (β1,β2,β3,α); then, the ML estimators of the unknown parameters are obtained as follows.The log likelihood function of the sample of sizen is given by logLθ¯=∑iϵI1logf1x1i,x2i+∑iϵI2logf2x1i,x2i+∑iϵI3logf3xi, where(14)I1=i;x1i<x2i,I2=i;x1i>x2i,I3=i;x1i=x2i=xi,I=I1∪I2∪I3,I1=n1,I2=n2,I3=n3,n1+n2+n3=n,where Ij denotes the cardinality of the set Ij, for j=1,2,3.Thus,(15)logLθ¯=2n1+2n2+n3logα+n1logβ1+β3+n1logβ2+n2logβ1+n2logβ2+β3+n3logβ3−α+1∑iϵI1∪I2log1+x1i+∑iϵI1∪I2log1+x2i+∑iϵI3log1+xi+β1+β3−1∑i∈I1log1−1+x1i−α+β2−1∑i∈I1log1−1+x2i−α+β1−1∑i∈I2log1−1+x1i−α+β2+β3−1∑i∈I2log1−1+x2i−α+β1+β2+β3−1∑i∈I3log1−1+xi−α,where θ¯=α,β1,β2,β3.The first derivatives of the log likelihood with respect to the unknown parameters and also the observed Fisher information matrix are obtained in AppendixA.2. The ML estimates θ¯^=α^,β^1,β^2,β^3 of θ¯=α,β1,β2,β3 are numerically obtained in Section 6.


### 3.2. Using Bayesian Approach

Let (X1,X2) follow BIK (θ¯), where θ¯=α,β1,β2,β3 is the vector of unknown parameters. The posterior pdf of parameters can be obtained as follows:(16)Pθ¯X1,X2∝LX1,X2;θ¯Pθ¯,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. They 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:(17)Pβi,ai,bi∝βiai−1exp−biβi,Pα,a4,b4∝αa4−1exp−b4α,where i=1,2,3 and ai,bi,a4,b4 are the hyperparameters. The posterior pdf has the following form:(18)Pθ¯X1,X2∝∏i=1n1α2β1+β3β21+x1i−α+11+x2i−α+1X1−1+x1i−αβ1+β3−11−1+x2i−αβ2−1∏i=1n2α2β2+β3β11+x1i−α+11+x2i−α+1X1−1+x1i−αβ1−11−1+x2i−αβ2+β3−1∏i=1n3αβ31+xi−α+11−1+xi−αβ1+β2+β3−1X∏i=13βiai−1exp−biβiαa4−1exp−b4α.It is observed that Bayesian estimators under square error loss function cannot be obtained in explicit forms. Therefore, we obtain the posterior mean using MCMC technique which is illustrated in Section6.


## 3.1. Using Maximum Likelihood Approach

In this section, we derive the maximum likelihood (ML) estimators of the unknown parameters of BIK distribution. Suppose {(x11,x21), …, (x1n,x2n)} is a random sample of size n from BIK (β1,β2,β3,α); then, the ML estimators of the unknown parameters are obtained as follows.The log likelihood function of the sample of sizen is given by logLθ¯=∑iϵI1logf1x1i,x2i+∑iϵI2logf2x1i,x2i+∑iϵI3logf3xi, where(14)I1=i;x1i<x2i,I2=i;x1i>x2i,I3=i;x1i=x2i=xi,I=I1∪I2∪I3,I1=n1,I2=n2,I3=n3,n1+n2+n3=n,where Ij denotes the cardinality of the set Ij, for j=1,2,3.Thus,(15)logLθ¯=2n1+2n2+n3logα+n1logβ1+β3+n1logβ2+n2logβ1+n2logβ2+β3+n3logβ3−α+1∑iϵI1∪I2log1+x1i+∑iϵI1∪I2log1+x2i+∑iϵI3log1+xi+β1+β3−1∑i∈I1log1−1+x1i−α+β2−1∑i∈I1log1−1+x2i−α+β1−1∑i∈I2log1−1+x1i−α+β2+β3−1∑i∈I2log1−1+x2i−α+β1+β2+β3−1∑i∈I3log1−1+xi−α,where θ¯=α,β1,β2,β3.The first derivatives of the log likelihood with respect to the unknown parameters and also the observed Fisher information matrix are obtained in AppendixA.2. The ML estimates θ¯^=α^,β^1,β^2,β^3 of θ¯=α,β1,β2,β3 are numerically obtained in Section 6.


## 3.2. Using Bayesian Approach

Let (X1,X2) follow BIK (θ¯), where θ¯=α,β1,β2,β3 is the vector of unknown parameters. The posterior pdf of parameters can be obtained as follows:(16)Pθ¯X1,X2∝LX1,X2;θ¯Pθ¯,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. They 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:(17)Pβi,ai,bi∝βiai−1exp−biβi,Pα,a4,b4∝αa4−1exp−b4α,where i=1,2,3 and ai,bi,a4,b4 are the hyperparameters. The posterior pdf has the following form:(18)Pθ¯X1,X2∝∏i=1n1α2β1+β3β21+x1i−α+11+x2i−α+1X1−1+x1i−αβ1+β3−11−1+x2i−αβ2−1∏i=1n2α2β2+β3β11+x1i−α+11+x2i−α+1X1−1+x1i−αβ1−11−1+x2i−αβ2+β3−1∏i=1n3αβ31+xi−α+11−1+xi−αβ1+β2+β3−1X∏i=13βiai−1exp−biβiαa4−1exp−b4α.It is observed that Bayesian estimators under square error loss function cannot be obtained in explicit forms. Therefore, we obtain the posterior mean using MCMC technique which is illustrated in Section6.


## 4. 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.


### 4.1. Derivation of MIK

In this section, we will derive the cdf and pdf of MIK.LetU1,…,Um+1 be m+1 independent random variables such that Ui∼IKβi,α, for i=1,…,m+1. Define Xj=maxUj,Um+1, j=1,…,m. Then, X=X1,…,XmT is an m-variate IK with parameters (β1,…,βm+1,α), and it will be denoted by MIK m,β1,…,βm+1,α. We have the following results regarding MIK distribution.Proposition 6.
IfX¯=X1,…,Xm∼MIKm,β1,…,βm+1,α, then the joint cdf of X¯ for x1>0,…,xm>0 is(19)FX¯x¯=∏i=1m+1FIKxi,βi,α,where x¯=x1,…,xm and xm+1=minx1,…,xm.Proof.
We prove by generalizing the same method illustrated in Section2.
Similar to the bivariate case, the MIK distribution can be written as, form>1,(20)FX¯x¯=kFax¯+1−kFsx¯,where 0<k<1 and Fa and Fs denote the absolute and continuous parts, respectively. The corresponding pdf can also be written as(21)fX¯x¯=kfax¯+1−kfsx¯.
The absolutely continuous part offax¯ can be obtained from(22)∂mFX¯x1,…,xm∂x1,…,∂xm,where x¯=x1,…,xmT belongs to the set where FX¯x¯ is absolutely continuous, if and only if xi′s are different. For each x¯,where xi are different, there exists a permutation ρ=i1,…,im such that xi1<,⋯,<xim.
Definefρx¯=fIKxi1,βi1+βm+1,αfIKxi2,βi2,α,…,fIKxim,βim,α. Then, for xi1<,⋯,<xim, ∂mFX¯x1,…,xm/∂x1,…,∂xm=kfax1,…,xm=fρxi1,…,xim, where k can be obtained as follows:(23)k=∑ρJρ,where(24)∫Rmfax1,…,xmdx1,…,dxm=∑ρ∫xim=0∞∫xim−1=0xim,…,∫xi1=0xi2fρx1,…,xmdxi1,…,dxim=∑ρJρ.
Since(25)∫xi1=0xi2fρx1,…,xmdxi1=FIKxi2;βi1+βm+1,α∏j=2mfIKxij,βij,α,∫xi2=0xi3∫xi1=0xi2fρx1,…,xmdxi1dxi2=βi2βi1+βi2+βm+1FIKxi3;βi1+βi2+βm+1,α∏j=3mfIKxij,βij,α,⋮Jρ=βi2βi1+βi2+βm+1Xβi3βi1+βi2+βi3+βm+1X,…,Xβimβi1+…+βim+βm+1,then(26)k=∑ρβi2βi1+βi2+βm+1Xβi3βi1+βi2+βi3+βm+1X,…,Xβimβi1+…+βim+βm+1,and for all xi1<,⋯,<xim, fax¯=1/kfρx¯.
Now, letIl=i1,…,il⊂I=1,…,m such that i1<,⋯,<il. fX¯x¯ can be written as(27)fX¯x¯=kfax¯+∑l=2m∑Il⊂IklfIlx¯,where fIl is a pdf with respect to m−l+1 dimensional Lebesgue measure on the hyperplane AIl=x¯∈ℜm:xi1=,⋯,=xil. The exact meaning of fX¯x¯ is as follows.
For any Borel measurable setB∈ℜm,(28)px∈B=k∫Bfax¯+∑l=2m∑Il⊂IkIl∫BIlfIlx¯,where BIl=B∩AIl is the projection of set B onto m−l+1 dimensional hyperplane AIl.
Now, we providekIl and fIlx¯. Note that if x¯ϵAIl, then x¯ has the following form:(29)x¯=x1…,xi1−1,x∗,xi1+1,…,xi2−1,x∗,xi2+1,…,xil−1,x∗,xil+1,…,xm.
For a givenx¯∈ℜm, we define gIl from the m−l+1 dimensional hyperplane AIl to ℜ as follows:(30)gIlx¯=fIKx∗,βm+1,αFIKx∗,∑i∈Ilβi,α∏iϵI−IlfIkxi,βi,α,if xi>x∗ for i∈I−Il and zero otherwise. Similar to k, we can obtain kIl as follows:(31)kIl=∫AIlgIlx¯dx¯∑ρI−Il∫xjm−l=0∞∫xjm−l−1=0xjm−l∫xj1=0xj2gIlx¯dx∗dxj1,…,dxjm−l=∑ρI−Ilβm+1∑iϵIlβi+βm+1∗βj1∑iϵIlβi+βm+1+βj1∗βjm−l∑iϵIβi+βm+1,fIlx¯=1kIlgIlx¯.


### 4.2. Properties of MIK

In this section, we will get the marginal and conditional distributions of MIK.Proposition 7.
IfX¯=X1,…,Xm∼MIKm,β1,…,βm,βm+1,α, then(a)
X1∼IKα,β1+βm+1,…,Xm∼IKα,βm+βm+1.(b)
The conditional distribution ofX1,…,Xs given Xs+1≤xs+1,…,Xm≤xm is(32)PX1≤x1,…,Xs≤xs|Xs+1≤xs+1,…,Xm≤xm=∏j=1sFIKxj,βj,αifz=υ∏j=1sFIKxj,βj,αFIKz,βm+1,αFIKυ,βm+1,αifz<υ,wherez=minx1,…,xs, υ=minxs+1,…,xm.(c)
For2≤s≤m,X1,…,Xs∼MIKs,β1,…,βs,βm+1,α.(d)
IfT=maxX1,…,Xm, then T∼IKβ1+,⋯,+βm+1,α.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.


## 4.1. Derivation of MIK

In this section, we will derive the cdf and pdf of MIK.LetU1,…,Um+1 be m+1 independent random variables such that Ui∼IKβi,α, for i=1,…,m+1. Define Xj=maxUj,Um+1, j=1,…,m. Then, X=X1,…,XmT is an m-variate IK with parameters (β1,…,βm+1,α), and it will be denoted by MIK m,β1,…,βm+1,α. We have the following results regarding MIK distribution.Proposition 6.
IfX¯=X1,…,Xm∼MIKm,β1,…,βm+1,α, then the joint cdf of X¯ for x1>0,…,xm>0 is(19)FX¯x¯=∏i=1m+1FIKxi,βi,α,where x¯=x1,…,xm and xm+1=minx1,…,xm.Proof.
We prove by generalizing the same method illustrated in Section2.
Similar to the bivariate case, the MIK distribution can be written as, form>1,(20)FX¯x¯=kFax¯+1−kFsx¯,where 0<k<1 and Fa and Fs denote the absolute and continuous parts, respectively. The corresponding pdf can also be written as(21)fX¯x¯=kfax¯+1−kfsx¯.
The absolutely continuous part offax¯ can be obtained from(22)∂mFX¯x1,…,xm∂x1,…,∂xm,where x¯=x1,…,xmT belongs to the set where FX¯x¯ is absolutely continuous, if and only if xi′s are different. For each x¯,where xi are different, there exists a permutation ρ=i1,…,im such that xi1<,⋯,<xim.
Definefρx¯=fIKxi1,βi1+βm+1,αfIKxi2,βi2,α,…,fIKxim,βim,α. Then, for xi1<,⋯,<xim, ∂mFX¯x1,…,xm/∂x1,…,∂xm=kfax1,…,xm=fρxi1,…,xim, where k can be obtained as follows:(23)k=∑ρJρ,where(24)∫Rmfax1,…,xmdx1,…,dxm=∑ρ∫xim=0∞∫xim−1=0xim,…,∫xi1=0xi2fρx1,…,xmdxi1,…,dxim=∑ρJρ.
Since(25)∫xi1=0xi2fρx1,…,xmdxi1=FIKxi2;βi1+βm+1,α∏j=2mfIKxij,βij,α,∫xi2=0xi3∫xi1=0xi2fρx1,…,xmdxi1dxi2=βi2βi1+βi2+βm+1FIKxi3;βi1+βi2+βm+1,α∏j=3mfIKxij,βij,α,⋮Jρ=βi2βi1+βi2+βm+1Xβi3βi1+βi2+βi3+βm+1X,…,Xβimβi1+…+βim+βm+1,then(26)k=∑ρβi2βi1+βi2+βm+1Xβi3βi1+βi2+βi3+βm+1X,…,Xβimβi1+…+βim+βm+1,and for all xi1<,⋯,<xim, fax¯=1/kfρx¯.
Now, letIl=i1,…,il⊂I=1,…,m such that i1<,⋯,<il. fX¯x¯ can be written as(27)fX¯x¯=kfax¯+∑l=2m∑Il⊂IklfIlx¯,where fIl is a pdf with respect to m−l+1 dimensional Lebesgue measure on the hyperplane AIl=x¯∈ℜm:xi1=,⋯,=xil. The exact meaning of fX¯x¯ is as follows.
For any Borel measurable setB∈ℜm,(28)px∈B=k∫Bfax¯+∑l=2m∑Il⊂IkIl∫BIlfIlx¯,where BIl=B∩AIl is the projection of set B onto m−l+1 dimensional hyperplane AIl.
Now, we providekIl and fIlx¯. Note that if x¯ϵAIl, then x¯ has the following form:(29)x¯=x1…,xi1−1,x∗,xi1+1,…,xi2−1,x∗,xi2+1,…,xil−1,x∗,xil+1,…,xm.
For a givenx¯∈ℜm, we define gIl from the m−l+1 dimensional hyperplane AIl to ℜ as follows:(30)gIlx¯=fIKx∗,βm+1,αFIKx∗,∑i∈Ilβi,α∏iϵI−IlfIkxi,βi,α,if xi>x∗ for i∈I−Il and zero otherwise. Similar to k, we can obtain kIl as follows:(31)kIl=∫AIlgIlx¯dx¯∑ρI−Il∫xjm−l=0∞∫xjm−l−1=0xjm−l∫xj1=0xj2gIlx¯dx∗dxj1,…,dxjm−l=∑ρI−Ilβm+1∑iϵIlβi+βm+1∗βj1∑iϵIlβi+βm+1+βj1∗βjm−l∑iϵIβi+βm+1,fIlx¯=1kIlgIlx¯.


## 4.2. Properties of MIK

In this section, we will get the marginal and conditional distributions of MIK.Proposition 7.
IfX¯=X1,…,Xm∼MIKm,β1,…,βm,βm+1,α, then(a)
X1∼IKα,β1+βm+1,…,Xm∼IKα,βm+βm+1.(b)
The conditional distribution ofX1,…,Xs given Xs+1≤xs+1,…,Xm≤xm is(32)PX1≤x1,…,Xs≤xs|Xs+1≤xs+1,…,Xm≤xm=∏j=1sFIKxj,βj,αifz=υ∏j=1sFIKxj,βj,αFIKz,βm+1,αFIKυ,βm+1,αifz<υ,wherez=minx1,…,xs, υ=minxs+1,…,xm.(c)
For2≤s≤m,X1,…,Xs∼MIKs,β1,…,βs,βm+1,α.(d)
IfT=maxX1,…,Xm, then T∼IKβ1+,⋯,+βm+1,α.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.


## 5. 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. Therefore, in this section, we will consider the process of estimation for MIK parameters. The proposed techniques could be applied for any multivariate MO distribution. Here, we will apply both maximum likelihood and Bayesian approaches.


### 5.1. Using Maximum Likelihood Approach

In this section, for simplification, we consider the case when we have three random variablesX1,X2, and X3. Applying Proposition 6, we have the following cdf:(33)FX¯x¯=∏i=14FIKxi,βi,α.The cdf can be rewritten as(34)FX1,X2,X3x1,x2,x3=F1x1,x2,x3,ifx1<x2<x3orx1<x3<x2,F2x1,x2,x3,ifx2<x1<x3orx2<x3<x1,F3x1,x2,x3,ifx3<x1<x2orx3<x2<x1F4x,x3,ifx1=x2=x<x3,F5x,x2,ifx1=x3=x<x2,F6x,x2,ifx2=x3=x<x1,F7x,ifx1=x2=x3=x,F8x,x3,ifx1=x2=x>x3,F9x,x2,ifx1=x3=x>x2,F10x,x1,ifx2=x3=x>x1,where(35)F1x1,x2,x3=FIKx1,β1+β4,αFIKx2,β2,αFIKx3,β3,α,F2x1,x2,x3=FIKx1,β1,αFIKx2,β2+β4,αFIKx3,β3,α,F3x1,x2,x3=FIKx1,β1,αFIKx2,β2,αFIKx3,β3+β4,α,F4x,x3=FIKx,β1+β2+β4,αFIKx3,β3,α,F5x,x2=FIKx,β1+β3+β4,αFIKx2,β2,α,F6x,x2=FIKx,β2+β3+β4,αFIKx1,β1,α,F7x=FIKx,β1+β2+β3+β4,α,F8x,x3=FIKx,β1+β2,αFIKx3,β3+β4,α,F9x,x2=FIKx,β1+β3,αFIKx2,β2+β4,α,F10x,x1=FIKx,β2+β3,αFIKx1,β1+β4,α.The pdf can be obtained by taking the derivative, except forf10 where we should take into consideration that the sum of all probabilities equals one.Now, supposex11,x21,x31,…,x1n,x2n,x3n is a random sample of size n from MIKα,β1,β2,β3,β4; the problem is to find the ML estimators of the unknown parameters. Consider the following notation:(36)I1=i;x1i<x2i<x3i,I2=i;x1i<x3i<x2i,I3=i;x2i<x1i<x3i,I4=i;x2i<x3i<x1i,I5=i;x3i<x1i<x2i,I6=i;x3i<x2i<x1i,I7=i;x1i=x2i=xi<x3i,I8=i;x1i=x3i=xi<x2i,I9=i;x2i=x3i=xi<x1i,I10=i;x1i=x2i=x3i=xi,I11=i;x1i=x2i=xi>x3i,I12=i;x1i=x3i=xi>x2i,I13=i;x2i=x3i=xi>x1i,Ij=nj,j=1,2,…,13,∑j=113nj=n,where Ij denotes the cardinality of the set Ij, for j=1,2,…,13.The log likelihood function of the sample of sizen is given by(37)logLθ¯=∑i∈I1∪I2logα3β1+β4β2β31+x1i−α+11−1+x1i−αβ1+β4−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3−1+∑i∈I3∪I4logα3β2+β4β1β31+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2+β4−11+x3i−α+1X1−1+x3i−αβ3−1+∑i∈I5∪I6logα3β3+β4β1β21+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3+β4−1+∑i∈I7logα2β1+β2+β4β31+xi−α+11−1+xi−αβ1+β2+β4−1X1+x3i−α+11−1+x3i−αβ3−1+∑i∈I8logα2β1+β3+β4β21+xi−α+11−1+xi−αβ1+β3+β4−1X1+x2i−α+11−1+x2i−αβ2−1+∑i∈I9logα2β2+β3+β4β11+xi−α+11−1+xi−αβ2+β3+β4−1X1+x1i−α+11−1+x1i−αβ1−1+∑i∈I10logαS1+xi−α+11−1+xi−αβ1+β2+β3+β4−1+∑i∈I11logα2β1+β2β3+β41+xi−α+11−1+xi−αβ1+β2−1X1+x3i−α+11−1+x3i−αβ3+β4−1+∑i∈I12logα2β1+β3β2+β41+xi−α+11−1+xi−αβ1+β3−1X1+x2i−α+11−1+x2i−αβ2+β4−1+∑i∈I13logα2β2+β3β1+β41+xi−α+11−1+xi−αβ2+β3−1X1+x1i−α+11−1+x1i−αβ1+β4−1,where θ¯=α,β1,β2,β3,β4.(38)S=β1+β2+β3+β4−β3β1+β2β1+β2+β4−β2β1+β3β1+β3+β4−β1β2+β3β3+β2+β4−3β1+β2+β3.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.


### 5.2. Using Bayesian Approach

LetX1,X2,X3 be three random variables from MIK (θ¯), where θ¯=α,β1,β2,β3,β4 is the vector of unknown parameters. The posterior pdf can be obtained as follows:(39)Pθ¯X1,X2,X3∝LX1,X2,X3;θ¯Pθ¯,where Pθ¯ is the prior distribution.We considered a gamma prior distribution with the following pdf:(40)Pβi,ai,bi∝βiai−1exp−biβi,Pα,a4,b4∝αa4−1exp−b4α,where i=1,2,3,4 and ai,bi,a4,b4 are the hyperparameters. The posterior pdf has the following form:(41)Pθ¯X¯∝∏i∈I1∪I2α3β1+β4β2β31+x1i−α+11−1+x1i−αβ1+β4−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3−1X∏i∈I3∪I4α3β2+β4β1β31+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2+β4−11+x3i−α+1X1−1+x3i−αβ3−1X∏i∈I5∪I6α3β3+β4β1β21+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3+β4−1X∏i∈I7α2β1+β2+β4β31+xi−α+11−1+xi−αβ1+β2+β4−1X1+x3i−α+11−1+x3i−αβ3−1X∏i∈I8α2β1+β3+β4β21+xi−α+11−1+xi−αβ1+β3+β4−1X1+x2i−α+11−1+x2i−αβ2−1X∏i∈I9α2β2+β3+β4β11+xi−α+11−1+xi−αβ2+β3+β4−1X1+x1i−α+11−1+x1i−αβ1−1X∏i∈I10αS1+xi−α+11−1+xi−αβ1+β2+β3+β4−1X∏i∈I11α2β1+β2β3+β41+xi−α+11−1+xi−αβ1+β2−1X1+x3i−α+11−1+x3i−αβ3+β4−1X∏i∈I12α2β1+β3β2+β41+xi−α+11−1+xi−αβ1+β3−1X1+x2i−α+11−1+x2i−αβ2+β4−1X∏i∈I13α2β2+β3β1+β41+xi−α+11−1+xi−αβ2+β3−1X1+x1i−α+11−1+x1i−αβ1+β4−1X∏i=14βiai−1exp−biβiαa4−1exp−b4α.As in the bivariate case, Bayesian estimation will be obtained numerically using MCMC which is illustrated in the next section.


## 5.1. Using Maximum Likelihood Approach

In this section, for simplification, we consider the case when we have three random variablesX1,X2, and X3. Applying Proposition 6, we have the following cdf:(33)FX¯x¯=∏i=14FIKxi,βi,α.The cdf can be rewritten as(34)FX1,X2,X3x1,x2,x3=F1x1,x2,x3,ifx1<x2<x3orx1<x3<x2,F2x1,x2,x3,ifx2<x1<x3orx2<x3<x1,F3x1,x2,x3,ifx3<x1<x2orx3<x2<x1F4x,x3,ifx1=x2=x<x3,F5x,x2,ifx1=x3=x<x2,F6x,x2,ifx2=x3=x<x1,F7x,ifx1=x2=x3=x,F8x,x3,ifx1=x2=x>x3,F9x,x2,ifx1=x3=x>x2,F10x,x1,ifx2=x3=x>x1,where(35)F1x1,x2,x3=FIKx1,β1+β4,αFIKx2,β2,αFIKx3,β3,α,F2x1,x2,x3=FIKx1,β1,αFIKx2,β2+β4,αFIKx3,β3,α,F3x1,x2,x3=FIKx1,β1,αFIKx2,β2,αFIKx3,β3+β4,α,F4x,x3=FIKx,β1+β2+β4,αFIKx3,β3,α,F5x,x2=FIKx,β1+β3+β4,αFIKx2,β2,α,F6x,x2=FIKx,β2+β3+β4,αFIKx1,β1,α,F7x=FIKx,β1+β2+β3+β4,α,F8x,x3=FIKx,β1+β2,αFIKx3,β3+β4,α,F9x,x2=FIKx,β1+β3,αFIKx2,β2+β4,α,F10x,x1=FIKx,β2+β3,αFIKx1,β1+β4,α.The pdf can be obtained by taking the derivative, except forf10 where we should take into consideration that the sum of all probabilities equals one.Now, supposex11,x21,x31,…,x1n,x2n,x3n is a random sample of size n from MIKα,β1,β2,β3,β4; the problem is to find the ML estimators of the unknown parameters. Consider the following notation:(36)I1=i;x1i<x2i<x3i,I2=i;x1i<x3i<x2i,I3=i;x2i<x1i<x3i,I4=i;x2i<x3i<x1i,I5=i;x3i<x1i<x2i,I6=i;x3i<x2i<x1i,I7=i;x1i=x2i=xi<x3i,I8=i;x1i=x3i=xi<x2i,I9=i;x2i=x3i=xi<x1i,I10=i;x1i=x2i=x3i=xi,I11=i;x1i=x2i=xi>x3i,I12=i;x1i=x3i=xi>x2i,I13=i;x2i=x3i=xi>x1i,Ij=nj,j=1,2,…,13,∑j=113nj=n,where Ij denotes the cardinality of the set Ij, for j=1,2,…,13.The log likelihood function of the sample of sizen is given by(37)logLθ¯=∑i∈I1∪I2logα3β1+β4β2β31+x1i−α+11−1+x1i−αβ1+β4−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3−1+∑i∈I3∪I4logα3β2+β4β1β31+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2+β4−11+x3i−α+1X1−1+x3i−αβ3−1+∑i∈I5∪I6logα3β3+β4β1β21+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3+β4−1+∑i∈I7logα2β1+β2+β4β31+xi−α+11−1+xi−αβ1+β2+β4−1X1+x3i−α+11−1+x3i−αβ3−1+∑i∈I8logα2β1+β3+β4β21+xi−α+11−1+xi−αβ1+β3+β4−1X1+x2i−α+11−1+x2i−αβ2−1+∑i∈I9logα2β2+β3+β4β11+xi−α+11−1+xi−αβ2+β3+β4−1X1+x1i−α+11−1+x1i−αβ1−1+∑i∈I10logαS1+xi−α+11−1+xi−αβ1+β2+β3+β4−1+∑i∈I11logα2β1+β2β3+β41+xi−α+11−1+xi−αβ1+β2−1X1+x3i−α+11−1+x3i−αβ3+β4−1+∑i∈I12logα2β1+β3β2+β41+xi−α+11−1+xi−αβ1+β3−1X1+x2i−α+11−1+x2i−αβ2+β4−1+∑i∈I13logα2β2+β3β1+β41+xi−α+11−1+xi−αβ2+β3−1X1+x1i−α+11−1+x1i−αβ1+β4−1,where θ¯=α,β1,β2,β3,β4.(38)S=β1+β2+β3+β4−β3β1+β2β1+β2+β4−β2β1+β3β1+β3+β4−β1β2+β3β3+β2+β4−3β1+β2+β3.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.


## 5.2. Using Bayesian Approach

LetX1,X2,X3 be three random variables from MIK (θ¯), where θ¯=α,β1,β2,β3,β4 is the vector of unknown parameters. The posterior pdf can be obtained as follows:(39)Pθ¯X1,X2,X3∝LX1,X2,X3;θ¯Pθ¯,where Pθ¯ is the prior distribution.We considered a gamma prior distribution with the following pdf:(40)Pβi,ai,bi∝βiai−1exp−biβi,Pα,a4,b4∝αa4−1exp−b4α,where i=1,2,3,4 and ai,bi,a4,b4 are the hyperparameters. The posterior pdf has the following form:(41)Pθ¯X¯∝∏i∈I1∪I2α3β1+β4β2β31+x1i−α+11−1+x1i−αβ1+β4−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3−1X∏i∈I3∪I4α3β2+β4β1β31+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2+β4−11+x3i−α+1X1−1+x3i−αβ3−1X∏i∈I5∪I6α3β3+β4β1β21+x1i−α+11−1+x1i−αβ1−1X1+x2i−α+11−1+x2i−αβ2−11+x3i−α+1X1−1+x3i−αβ3+β4−1X∏i∈I7α2β1+β2+β4β31+xi−α+11−1+xi−αβ1+β2+β4−1X1+x3i−α+11−1+x3i−αβ3−1X∏i∈I8α2β1+β3+β4β21+xi−α+11−1+xi−αβ1+β3+β4−1X1+x2i−α+11−1+x2i−αβ2−1X∏i∈I9α2β2+β3+β4β11+xi−α+11−1+xi−αβ2+β3+β4−1X1+x1i−α+11−1+x1i−αβ1−1X∏i∈I10αS1+xi−α+11−1+xi−αβ1+β2+β3+β4−1X∏i∈I11α2β1+β2β3+β41+xi−α+11−1+xi−αβ1+β2−1X1+x3i−α+11−1+x3i−αβ3+β4−1X∏i∈I12α2β1+β3β2+β41+xi−α+11−1+xi−αβ1+β3−1X1+x2i−α+11−1+x2i−αβ2+β4−1X∏i∈I13α2β2+β3β1+β41+xi−α+11−1+xi−αβ2+β3−1X1+x1i−α+11−1+x1i−αβ1+β4−1X∏i=14βiai−1exp−biβiαa4−1exp−b4α.As in the bivariate case, Bayesian estimation will be obtained numerically using MCMC which is illustrated in the next section.


## 6. Numerical Analysis

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


### 6.1. 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.


#### 6.1.1. A Simulation Study

Here, we present an algorithm to generate BIK distribution (Algorithm1). To perform a simulation study, we first need to select initial values for the parameters. Here, the following eight different populations are considered:(i)
Case 1:α=0.2,β1=0.2,β2=0.2,β3=0.2.(ii)
Case 2:α=0.5,β1=0.5,β2=0.5,β3=0.5.(iii)
Case 3:α=0.8,β1=0.8,β2=0.8,β3=0.8.(iv)
Case 4:α=1.5,β1=1.2,β2=1.3,β3=2.5.(v)
Case 5:α=2,β1=2,β2=2,β3=2.(vi)
Case 6:α=2,β1=2.5,β2=3,β3=3.5.(vii)
Case 7:α=2.5,β1=2.5,β2=3,β3=2.5.(viii)
Case 8:α=3.5,β1=4,β2=3.5,β3=3.5.Algorithm 1: An Algorithm to Generate BIK Distribution.
The following algorithm is to generateX1,X2 from BIK distribution.Step 1: generateU1,U2, and U3 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, and T3=1−U31/β31/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, and Z3=1/T3−1.Step 4: obtainX1=maxZ1,Z3,X2=maxZ2,Z3.The parameters are selected to cover different shapes of the distribution. It can be seen from Figure1 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.Figure 1
Different shapes of the absolutely continuous part of BIK pdf. (a)α = 0.2, β1 = 0.2, β2 = 0.2, β3 = 0.2; (b) α = 0.5, β1 = 0.5, β2 = 0.5, β3 = 0.5; (c) α = 0.8, β1 = 0.8, β2 = 0.8, β3 = 0.8; (d) α = 2, β1 = 2, β2 = 2, β3 = 2; (e) α = 1.5, β1 = 1.2, β2 = 1.3, β3 = 2.5; (f) α = 2, β1 = 2.5, β2 = 3, β3 = 3.5; (g) α = 2.5, β1 = 2.5, β2 = 3, β3 = 2.5; (h) α = 3.5, β1 = 4, β2 = 3.5, β3 = 3.5.
(a)(b)(c)(d)(e)(f)(g)(h)(1) Maximum Likelihood Approach. The 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 and presented in Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.Table 1
The results of MLEs and Bayesian estimates.
MLEBayesianαβ1β2β3αβ1β2β3α=0.2,β1=0.2, β2=0.2,β3=0.2n = 30ABias0.0200.0150.0150.0100.0160.0120.0080.012MSE0.0030.0040.0040.0040.0030.0040.0040.003CW (CL)0.2140.2480.2480.2480.2130.2430.2390.230CP0.9500.9380.9570.9480.9430.9320.9490.951n = 50ABias0.0110.0090.0080.0070.0100.0090.0040.009MSE0.0020.0020.0020.0020.0020.0020.0020.002CW (CL)0.1760.1760.1760.1760.1620.1880.1850.179CP0.9540.9400.9380.9510.9560.9430.9490.942n = 70ABias0.0100.0060.0060.0060.0080.0060.0040.007MSE0.0010.0020.0020.0010.0010.0020.0020.001CW (CL)0.1240.1760.1760.1240.1360.1570.1560.143CP0.9520.9570.9460.9510.9580.9450.9490.948α=0.5,β1=0.5, β2=0.5,β3=0.5n = 30ABias0.0320.0480.0460.0320.0280.0410.0270.038MSE0.0110.0340.0340.0260.0110.0310.0300.025CW (CL)0.3920.7020.7020.6200.3940.6640.6590.587CP0.9470.9450.9560.9510.9420.9300.9450.951n = 50ABias0.0180.0280.0260.0220.0200.0270.0260.025MSE0.0060.0180.0180.0140.0060.0170.0180.014CW (CL)0.3040.5120.5120.4640.2990.5010.5040.448CP0.9480.9440.9450.9570.9540.9490.9460.937n = 70ABias0.0170.0190.0210.0190.0130.0210.0170.019MSE0.0040.0120.0120.0100.0040.0120.0120.010CW (CL)0.2480.4300.4300.3920.2500.4170.4180.375CP0.9450.9580.9480.9500.9510.9290.9510.946α=0.8,β1=0.8, β2=0.8,β3=0.8n = 30ABias0.0450.0880.0860.0660.0370.0720.0730.060MSE0.0220.1010.0990.0760.0220.0900.0910.068CW (CL)0.5541.1961.1881.0520.5571.1291.1330.979CP0.9410.9510.9540.9460.9400.9390.9500.951n = 50ABias0.0260.0510.0490.0400.0260.0590.0490.040MSE0.0130.0520.0510.0400.0130.0520.0500.038CW (CL)0.4300.8680.8680.7640.4250.8540.8500.744CP0.9500.9440.9410.9560.9560.9480.9480.940n = 70ABias0.0240.0360.0390.0350.0190.0370.0300.030MSE0.0090.0340.0360.0280.0090.0340.0330.026CW (CL)0.3500.7120.7220.6440.3580.7030.7030.619CP0.9440.9560.9490.9490.9580.9400.9520.946α=1.5,β1=1.2, β2=1.3,β3=2.5n = 30ABias0.0680.1420.1580.3070.0230.1200.1620.137MSE0.0620.3910.4420.9130.0590.3550.4470.664CW (CL)0.9362.3882.5323.5480.9242.2072.4173.063CP0.9500.9400.9500.9500.9400.9500.9300.960n = 50ABias0.0460.0990.1110.1910.0180.0760.1030.091MSE0.0350.2030.2310.4440.0340.1890.2190.369CW (CL)0.7121.7221.8342.5040.7031.6331.7382.043CP0.9600.9500.9500.9600.9600.9200.9500.970n = 70ABias0.0390.0680.0820.1570.0130.0680.0510.090MSE0.0260.1340.1530.3010.0230.1310.1390.258CW (CL)0.6081.4081.4981.6270.5861.3641.4161.941CP0.9300.9400.9400.9600.9100.9400.9500.960α=2,β1=2, β2=2,β3=2n = 30ABias0.0940.2960.2950.2370.0600.2640.2470.147MSE0.0990.8900.8920.6800.0940.8110.7900.519CW (CL)1.1763.5103.5183.0961.1663.3123.2372.723CP0.9430.9490.9500.9520.9490.9310.9480.952n = 50ABias0.0550.1690.1660.1360.0400.1850.1690.101MSE0.0550.4270.4280.3260.0550.4320.4180.290CW (CL)1.1763.5103.5183.0960.8982.4432.3882.050CP0.9530.9460.9410.9610.9520.9460.9450.952n = 70ABias0.0510.1230.1190.0510.0310.1290.1010.063MSE0.0400.2810.2850.2260.0380.2770.2630.195CW (CL)0.7542.0222.0301.8040.7531.9791.9481.695CP0.9420.9530.9450.9420.9610.9420.9490.942α=2,β1=2.5, β2=3,β3=3.5n = 30ABias0.0860.3830.4930.5070.0370.2680.4260.178MSE0.0891.8322.5832.4760.2191.4972.4110.908CW (CL)1.1225.0885.9965.8400.8434.5435.5633.387CP0.9500.9500.9500.9500.9420.9500.8800.938n = 50ABias0.0530.2260.2790.2970.0190.1970.3040.095MSE0.0510.8701.1941.1320.0490.8151.1970.865CW (CL)0.8583.5484.1414.0060.8423.3763.9963.541CP0.9500.9400.9400.9700.9300.9100.9300.950n = 70ABias0.0490.1730.2360.2440.0140.1660.1490.091MSE0.0360.5740.8000.7660.0330.5430.6880.607CW (CL)0.7122.8943.3823.2940.6972.7703.1162.997CP0.9400.9400.9500.9600.9500.9300.9300.960α=2.5,β1=2.5, β2=3,β3=2.5n = 30ABias0.1110.3900.4890.3260.0670.3740.4180.168MSE0.1401.5352.1911.2340.1341.4351.9380.891CW (CL)1.4024.6105.4764.1641.3954.3755.0163.545CP0.9430.9500.9540.9470.9460.9030.9150.955n = 50ABias0.0680.2320.2790.1830.0410.2360.2620.108MSE0.0790.7371.0170.5770.0780.7250.9770.496CW (CL)1.0663.2403.7982.8921.0743.1673.6572.697CP0.9500.9400.9490.9610.9490.9170.8860.952n = 70ABias0.0600.1640.2140.1600.0280.1660.1660.069MSE0.0570.4790.6690.3990.0370.4720.6220.333CW (CL)0.9022.6363.0942.3940.9032.5812.9722.217CP0.9480.9500.9500.9430.9660.8870.8830.950α=3.5,β1=4, β2=3.5,β3=3.5n = 30ABias0.1520.7170.6070.5330.0680.6630.5360.115MSE0.2574.5743.5522.8880.2374.1683.1271.682CW (CL)1.8967.8986.9946.3221.8737.4036.3564.974CP0.9380.9490.9490.9540.8360.9280.9550.956n = 50ABias0.0910.4140.3520.2890.0530.4380.3570.133MSE0.1432.0751.6371.2770.1412.0801.6061.038CW (CL)1.4405.4104.8224.2821.4485.3114.6413.905CP0.9520.9470.9340.9660.7610.9190.9320.963n = 70ABias0.0820.3090.2680.2590.0310.2830.2320.060MSE0.1031.3521.0660.8830.0991.2971.0120.701CW (CL)1.2144.3943.9083.5421.2224.2753.7733.240CP0.9400.9510.9450.9530.6960.9280.9570.953(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.The 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θ¯X1,X2. The 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. The 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 statisticzn converges to normal distribution for large sample sizes. Hence, large absolute values of zn 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 Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.From Table1, it can be seen that under different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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. This is clear from Figure2.Figure 2
MSE of (a)α=2, (b) β1=2, (c) β2=2, and (d) β3=2.
(a)(b)(c)(d)Also, it can be seen that as the sample sizen 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.Figure 3
(a) MSE and (b) ABias forα.
(a)(b)


#### 6.1.2. 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. The dataset has been obtained from Meintanis [23]; he used the bivariate MO exponential distribution (BE) to analyze the data. The 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 (X1) and the time of the first goal of any type scored by the home team (X2).The bivariate dataset has the following structure:X1 < X2, X1 = X2, and X1 > X2. Since, X1=X2 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. This is not going to make any difference in the statistical inference.First, before analyzing the data, we fit inverted Kumaraswamy distribution toX1, X2 and max (X1,X2). To guess the initial values for the parameters of BIK model, the MLEs of the shape parameters α,β of the respective inverted Kumaraswamy distribution for X1, X2 and max (X1,X2) 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 variableX in a population. The 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 F0x, then the null hypothesis is(42)H0:Fx=F0x.One approach for applying GOF tests is based on the empirical distribution function (EDF)Fnx which is defined as follows:(43)Fnx=number of elements in the sample≤xtotal sample sizen.This approach is based on defining a statistic to measure the discrepancy betweenFnx and F0x. The most used statistics are modified Cramér–von Mises statistic W∗ and Anderson–Darling statistic A∗ which have the following formulas:(44)W∗=112n+∑i=1nF0xi−2i−12n2,A∗=−n−2∑i=1n2i−12nlnF0xi+1−2i−12nln1−F0xi,where xi is the value of the ith order statistics in the sample.Large values of test statistics (or small correspondingp 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, and max (X1,X2) are illustrated in Table 2.Table 2
Cramér–von Mises statisticW∗ and Anderson–Darling statistic A∗ for IK.
X1X2maxX1,X2A∗0.966 (0.375)0.549 (0.696)1.395 (0.204)W∗0.139 (0.426)0.083 (0.679)0.233 (0.212)Based on the values of these statistics and the correspondingp 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.Table 3
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE5.2581.7630.5601.411Confidence interval width2.8242.0260.8861.374BayesianPosterior mean4.7271.6730.6791.442Credible interval length2.1681.2860.4911.061Now, 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:(45)IC=−2logL+kp,where k is the penalized term and p refers to number of parameters in the model.The 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 Table4.Table 4
Information criteria.
Log (L)AICBICCAICHQICk—2lnnlnn+1lnlnnBIK−23.53∗55.0661.561.6157.33BE−44.56∗∗95.1299.95100.0396.82BIW−30.25∗∗∗68.574.9475.0570.77∗Based on the estimates in Table 3. ∗∗Based on the estimates obtained by Meintanis [23]. ∗∗∗Based on the estimates presented by Muhammed [9].By analyzing equation (45), we can see that the first term −2logL 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. The model with the lowest IC is the best (for more details about IC, see Vrieze [25]).The log likelihood value, AIC, BIC, CAIC, and HQIC are represented in Table4. All of the criteria suggest that BIK provides the best fit compared with BE and BIW models. This may show the importance of BIK.(2) Motor Data. The dataset has been obtained from [26]. The data are about failure time in days for a parallel system containing two motors. The variables are time to failure for the first motor X1 and time to failure for the second motor (X2). 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.Table 5
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗1.133 (0.294)0.294 (0.942)0.488 (0.756)W∗0.204 (0.261)0.036 (0.958)0.081 (0.692)IWA∗1.491 (0.179)0.874 (0.429)7.165 (0.0003)W∗0.267 (0.168)0.119 (0.505)1.509 (9.18 ×10−5)EA∗3.902 (0.009)3.363 (0.014)94.4 (3.3 ×10−5)W∗0.778 (0.007)0.672 (0.014)5.964 (2.2 ×10−16)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 Table6.Table 6
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE16.97712.60313.0670.01Confidence interval width10.50819.5920.04618.730BayesianPosterior mean16.95012.90013.2100.034Credible interval length9.02019.15220.7200.594(3) Diabetic Retinopathy Data. The dataset has been obtained from [27]. The 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. The variables are time to failure for treated eye (X1) and time to failure for untreated eye (X2). 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.Table 7
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗0.662 (0.591)0.570 (0.676)0.707 (0.551)W∗0.092 (0.627)0.099 (0.594)0.090 (0.635)IWA∗2.538 (0.048)73.841 (8.4 ×10−6)166.45 (8.4 ×10−6)W∗0.369 (0.087)11.41 (2.2 ×10−16)18.25 (2.2 ×10−16)EA∗3.886 (0.009)214.7 (8.4 ×10−6)271.7 (8.4 ×10−6)W∗0.671 (0.015)22.4 (2.2 ×10−16)23.1 (2.2 ×10−16)In Table8, 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.Table 8
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE3.8302.8892.5700.900Confidence interval width1.2882.1141.9521.044BayesianPosterior mean3.5862.7392.37960.907Credible interval length0.9301.3261.1330.700From these three datasets, we can conclude that the derived BIK distribution will be of great importance.


### 6.2. For MIK Distribution

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


#### 6.2.1. A Simulation Study

Here, we present an algorithm to generate MIK distribution. Also, we illustrate the simulation results for both ML and Bayesian approaches (Algorithm2).Algorithm 2: An algorithm to generate MIK distribution.
The following algorithm is to generateX1,X2,X3 from MIK distribution.Step 1: generateU1,U2,U3, and U4 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, T3=1−U31/β31/α, and T4=1−U41/β41/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, Z3=1/T3−1, and Z4=1/T4−1.Step 4: obtainX1=maxZ1,Z4,X2=maxZ2,Z4, and X3=maxZ3,Z4.(1) Maximum Likelihood Approach. To obtain the maximum likelihood estimates, a Monte Carlo simulation is performed using R package with 1000 replications and three different sample sizes n=30,50, and 70. Different initial values of the parameters are arbitrarily chosen, varying between small and large values to cover different cases of the distribution. The parameter β4 is assumed known for simplicity. The following twelve populations are considered.(i)
α=0.9,β1=0.9,β2=0.9,β3=0.9.(ii)
α=1.8,β1=0.8,β2=1.8,β3=0.8.(iii)
α=1.8,β1=1.8,β2=1.8,β3=1.8.(iv)
α=2,β1=1.5,β2=2,β3=2.(v)
α=2.4,β1=2.4,β2=2.4,β3=2.4.(vi)
α=2.5,β1=1.5,β2=1.5,β3=1.5.(vii)
α=3,β1=3,β2=3,β3=3.(viii)
α=3,β1=2.8,β2=2.8,β3=3.5.(ix)
α=3,β1=3,β2=4,β3=3.(x)
α=3.5,β1=3.5,β2=3.5,β3=4.(xi)
α=3.5,β1=3.5,β2=3.5,β3=3.5.(xii)
α=4,β1=1.5,β2=2,β3=2.5.To check the behavior of the estimates, ABias, MSE, CW, and CP are computed in Table9. The algorithm is explained in detail in Appendix B.2.Table 9
The results of MLEs and Bayesian estimates for MIK.
MLEBayesianαβ1β2β3αβ1β2β3α=0.9,β1=0.9, β2=0.9,β3=0.9n = 30ABias0.0460.0630.0500.0610.0450.0370.0200.032MSE0.0200.0450.0420.0450.0190.0400.0370.039CW (CL)0.5260.7940.7740.7940.5050.7610.7350.749CP0.9530.9730.9540.9490.9390.9250.9480.933n = 50ABias0.0350.0330.0460.0400.0390.0050.0200.019MSE0.0110.0230.0250.0250.0120.0210.0220.023CW (CL)0.3920.5820.5940.5940.3910.5670.5730.578CP0.9600.9680.9500.9600.9450.9460.9290.943n = 70ABias0.0200.0260.0210.0320.0290.0010.0050.006MSE0.0070.0170.0150.0170.0080.0150.0150.015CW (CL)0.3280.4960.4800.4960.3280.4770.4720.479CP0.9450.9510.9580.9640.9300.9310.9330.950α=1.8,β1=0.8, β2=1.8,β3=0.8n = 30ABias0.1350.0830.2240.0890.1190.0430.1990.046MSE0.0840.0440.2820.0460.0780.0330.2480.034CW (CL)1.0080.7541.8880.7640.9770.6851.7560.687CP0.9180.9610.9730.9520.9270.9400.9480.935n = 50ABias0.1110.0670.1550.0620.0850.0150.1670.022MSE0.0510.0250.1480.0250.0440.0180.1480.018CW (CL)0.7740.5681.3800.5680.7490.5111.3350.516CP0.9390.9620.9500.9580.9350.9560.9250.938n = 70ABias0.0860.0570.1550.0490.0720.0110.1200.015MSE0.0340.0170.1130.0160.0310.0120.0940.013CW (CL)0.6440.4641.1700.4640.6280.4321.0190.433CP0.9140.9690.9300.9770.9190.9510.9280.948α=1.8,β1=1.8, β2=1.8,β3=1.8n = 30ABias0.1120.2100.1730.2190.0900.1950.1690.167MSE0.0640.2620.2360.2690.0600.2450.2290.224CW (CL)0.8861.8301.7801.8420.8841.7631.7251.718CP0.9210.9620.9720.9620.9170.9420.9530.934n = 50ABias0.0810.1370.1500.1470.0740.1300.1620.115MSE0.0370.1340.1410.1390.0360.1280.1430.123CW (CL)0.6781.3301.3461.3400.6771.2991.3241.288CP0.9510.9720.9440.9650.9440.9640.9220.952n = 70ABias0.0590.1130.1160.1220.0540.0970.1280.091MSE0.0240.0920.0930.0950.0240.0860.0950.084CW (CL)0.5681.1021.1081.3680.5621.0751.0901.074CP0.9320.9520.9490.9590.9410.9410.9280.946α=2,β1=1.5, β2=2,β3=2n = 30ABias0.1190.1590.2080.2730.0840.1300.1890.196MSE0.0760.1790.3100.3670.0660.1430.2800.286CW (CL)0.9761.4782.0262.1180.9401.3711.9001.924CP0.9230.9640.9790.9640.9230.9450.9520.945n = 50ABias0.0920.1030.1730.1760.0750.0920.1570.124MSE0.0440.0860.1810.0670.0410.0800.1650.154CW (CL)0.7441.0741.5241.5240.7261.0361.4441.435CP0.9480.9690.9440.9650.9450.9630.9270.954n = 70ABias0.0650.0860.1220.1350.0450.0620.1040.097MSE0.0290.0550.1150.1190.0270.0530.1050.104CW (CL)0.6200.8941.2401.2560.6060.8561.1881.189CP0.9270.9500.9570.9670.9330.9450.9500.950α=2.4,β1=2.4, β2=2.4,β3=2.4n = 30ABias0.2840.2450.2980.1330.1060.2640.2400.209MSE0.4940.4560.5080.0980.0930.4630.4410.410CW (CL)2.5202.4662.5381.1081.1102.4242.3662.340CP0.9670.9730.9670.9340.9430.9480.9540.950n = 50ABias0.1840.2040.1980.0950.0700.1560.2000.132MSE0.2510.2540.2600.0560.0520.2320.2570.219CW (CL)1.8261.8561.8420.8500.8461.7791.8041.751CP0.9750.9520.9640.9520.9550.9660.9250.955n = 70ABias0.1460.1360.1570.0700.0590.1490.1390.117MSE0.1700.1660.1760.0380.0370.1680.1650.156CW (CL)1.5141.5081.5240.7120.7131.4911.4811.469CP0.9440.9520.9620.9390.9380.9280.9500.956α=2.5,β1=1.5, β2=1.5,β3=1.5n = 30ABias0.1680.1680.1510.1750.1460.1510.1410.138MSE0.1350.1730.1630.1770.1240.1550.1490.149CW (CL)1.2821.4921.4661.4981.2401.4081.3881.387CP0.9310.9610.9790.9630.9190.9430.9550.950n = 50ABias0.1250.1200.1190.1250.0960.1020.1020.094MSE0.0780.0920.1150.0950.0690.0830.1520.082CW (CL)0.9761.0941.1021.1020.9471.0471.0431.044CP0.9460.9690.9520.9660.9490.9560.9230.953n = 70ABias0.0920.1000.0940.0940.0860.0920.1050.062MSE0.0510.0640.0620.0620.0500.0600.0630.054CW (CL)0.8120.9100.9020.9020.7990.8760.8780.862CP0.9130.9420.9500.9710.9270.9490.9210.966α=3,β1=3, β2=3,β3=3n = 30ABias0.1570.3770.3180.3740.1170.3610.2930.249MSE0.1410.8370.7950.8330.1310.7960.7160.658CW (CL)1.3363.2683.1803.2641.3283.1503.0252.975CP0.9360.9660.9740.9690.9440.9320.9610.956n = 50ABias0.1210.2420.2830.2580.0960.2140.2810.174MSE0.0830.4200.4560.4340.0790.3970.4510.367CW (CL)1.0222.3562.4042.3741.0242.3062.3452.248CP0.9410.9750.9510.9630.9390.9620.9290.958n = 70ABias0.0740.1650.1590.1880.0580.1500.1580.137MSE0.0520.2690.2660.2810.0510.2590.2640.249CW (CL)0.8501.9281.9241.9440.8551.8931.8921.876CP0.9370.9510.9610.9630.9430.9380.9480.956α=3,β1=2.8,β2=2.8,β3=3.5n = 30ABias0.1540.3390.2830.4550.1030.2840.2430.257MSE0.1390.7000.6331.2210.1260.6210.5820.896CW (CL)1.3302.9982.9163.9481.3142.8342.7553.513CP0.9370.9640.9770.9700.8640.9380.8470.914n = 50ABias0.1210.2180.2520.3240.0870.1690.2330.210MSE0.0830.3540.3810.6420.0760.3200.3620.529CW (CL)1.0222.1682.2082.8721.0182.0982.1362.697CP0.9470.9730.9550.9650.8050.9740.7970.857n = 70ABias0.0810.1560.1480.2380.0530.1210.1330.137MSE0.0540.2300.2270.4140.0500.2110.2170.344CW (CL)0.8501.7801.7742.3420.8491.7261.7332.220CP0.9330.9510.9650.9630.8000.9470.7060.839α=3,β1=3,β2=4,β3=3n = 30ABias0.1540.3740.4770.3700.1090.3270.4550.230MSE0.1360.8271.5580.8220.1200.7261.4410.617CW (CL)1.3123.2504.5203.2441.2843.0484.3072.914CP0.9370.9630.9730.9710.9410.9510.9700.963n = 50ABias0.1160.2320.4060.2560.0910.2060.3880.171MSE0.0790.4090.9150.4300.0730.3800.8510.354CW (CL)1.0082.3363.3942.3660.9892.2543.2542.213CP0.9420.9730.9600.9650.9420.9530.9650.952n = 70ABias0.0740.1630.2340.1820.0580.1560.2280.123MSE0.0510.2660.5300.2760.0490.2570.5130.241CW (CL)0.8401.9162.7021.9320.8341.8692.6361.847CP0.9430.9560.9600.9600.9550.9650.9670.958α=3.5,β1=3.5,β2=3.5,β3=4n = 30ABias0.1710.4500.3780.5360.1000.4020.3890.238MSE0.1761.1901.0741.6660.1571.0951.0711.149CW (CL)1.5023.8943.7824.6041.4863.7023.6604.027CP0.9390.9630.9750.9690.9490.9520.9310.937n = 50ABias0.1320.2810.3300.3710.0800.2530.2680.166MSE0.1040.5870.6390.8590.0930.5620.5730.653CW (CL)1.1562.7942.8543.3281.1462.7392.7153.059CP0.9440.9720.9590.9660.9360.9610.9000.903n = 70ABias0.0850.1970.1850.2610.0550.1450.1940.159MSE0.0680.3810.3720.5450.0650.3510.3780.465CW (CL)0.9682.2922.2802.7080.9672.2302.2582.583CP0.9400.9520.9630.9700.9500.9590.8710.864α=3.5,β1=3.5, β2=3.5,β3=3.5n = 30ABias0.1720.4500.3780.4450.1070.3690.3230.250MSE0.1791.1941.0791.1880.1621.0560.9920.885CW (CL)1.5143.9023.7923.9001.5033.6893.5933.488CP0.9410.9660.9750.9710.9550.9450.9590.955n = 50ABias0.1330.2840.3310.3080.1000.2630.3080.199MSE0.1060.5920.6430.6170.0990.5730.6100.516CW (CL)1.1622.3832.8622.8321.1602.7552.7552.665CP0.9460.9710.9590.9650.9460.9640.9290.949n = 70ABias0.0870.2040.1860.2240.0590.1630.1820.144MSE0.0690.3870.3750.3980.0660.3620.3740.346CW (CL)0.9682.3022.2862.3120.9782.2522.2592.220CP0.9370.9540.9630.9630.9410.9510.9560.956α=4,β1=1.5, β2=2,β1=2.5n = 30ABias0.2340.1520.1990.3300.1770.1430.1850.236MSE0.3510.1620.3020.5900.2720.1550.2840.474CW (CL)1.9281.4622.0062.7221.9061.4151.9122.486CP0.9320.9630.9740.9710.9500.9410.9450.954n = 50ABias0.1780.1000.1670.2310.1410.0920.1550.167MSE0.1730.0850.1770.3060.1620.0830.1690.260CW (CL)1.4721.0741.5141.9721.4681.0591.4681.871CP0.9440.9720.9500.9630.9600.9580.9530.929n = 70ABias0.1360.0820.1200.1930.1210.0730.1370.126MSE0.1170.0590.1130.2090.1160.0570.1200.177CW (CL)1.2340.8941.2341.6261.2410.8831.2361.559CP0.9290.9500.9620.9710.9430.9390.9460.935(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. The algorithm is explained in detail in Appendix B.2.From Table9, it can be seen that for different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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 Figure4. Also, for the majority of cases, Bayesian estimates have smaller ABias than ML ones.Figure 4
MSE of (a)α=3.5, (b) β1=3.5, (c) β2=3.5, and (d) β3=3.5.
(a)(b)(c)(d)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.Figure 5
(a) MSE and (b) ABias forα.
(a)(b)


#### 6.2.2. A Real Dataset

Here, we analyze a real dataset to show the applicability of MIK. The 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. The variables are as follows:X1: first systolic blood pressure measurement.X2: second systolic blood pressure measurement.X3: third systolic blood pressure measurement.All data points have been divided by 1000. This 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, X3, and maxX1,X2,X3 are illustrated in Table 10.Table 10
W∗ and A∗ for IK, IW, and E.
ModelX1X2X3maxX1,X2IKA∗0.31 (0.93)0.43 (0.82)0.49 (0.76)0.50 (0.75)W∗0.04 (0.94)0.07 (0.76)0.08 (0.68)0.08 (0.72)IWA∗13.60 (7 ×10−6)14.48 (7 ×10−6)15.80 (7 ×10−6)14.20 (7 ×10−6)W∗2.79 (2 ×10−7)2.99 (5 ×10−8)3.30 (9 ×10−9)2.91 (7 ×10−8)EA∗23.95 (7 ×10−6)24.18 (7 ×10−6)24.75 (7 ×10−6)24.29 (7 ×10−6)W∗5.06 (2 ×10−16)5.12 (2 ×10−16)5.26 (2 ×10−16)5.14 (2 ×10−16)Based on the values of these statistics and the correspondingp 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.Table 11
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE46.949264.199245.852245.519Confidence interval width8.126269.684256.424263.866BayesianPosterior mean48.697270.528226.877257.231Credible interval length4.27069.21561.11064.305


## 6.1. 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.


### 6.1.1. A Simulation Study

Here, we present an algorithm to generate BIK distribution (Algorithm1). To perform a simulation study, we first need to select initial values for the parameters. Here, the following eight different populations are considered:(i)
Case 1:α=0.2,β1=0.2,β2=0.2,β3=0.2.(ii)
Case 2:α=0.5,β1=0.5,β2=0.5,β3=0.5.(iii)
Case 3:α=0.8,β1=0.8,β2=0.8,β3=0.8.(iv)
Case 4:α=1.5,β1=1.2,β2=1.3,β3=2.5.(v)
Case 5:α=2,β1=2,β2=2,β3=2.(vi)
Case 6:α=2,β1=2.5,β2=3,β3=3.5.(vii)
Case 7:α=2.5,β1=2.5,β2=3,β3=2.5.(viii)
Case 8:α=3.5,β1=4,β2=3.5,β3=3.5.Algorithm 1: An Algorithm to Generate BIK Distribution.
The following algorithm is to generateX1,X2 from BIK distribution.Step 1: generateU1,U2, and U3 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, and T3=1−U31/β31/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, and Z3=1/T3−1.Step 4: obtainX1=maxZ1,Z3,X2=maxZ2,Z3.The parameters are selected to cover different shapes of the distribution. It can be seen from Figure1 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.Figure 1
Different shapes of the absolutely continuous part of BIK pdf. (a)α = 0.2, β1 = 0.2, β2 = 0.2, β3 = 0.2; (b) α = 0.5, β1 = 0.5, β2 = 0.5, β3 = 0.5; (c) α = 0.8, β1 = 0.8, β2 = 0.8, β3 = 0.8; (d) α = 2, β1 = 2, β2 = 2, β3 = 2; (e) α = 1.5, β1 = 1.2, β2 = 1.3, β3 = 2.5; (f) α = 2, β1 = 2.5, β2 = 3, β3 = 3.5; (g) α = 2.5, β1 = 2.5, β2 = 3, β3 = 2.5; (h) α = 3.5, β1 = 4, β2 = 3.5, β3 = 3.5.
(a)(b)(c)(d)(e)(f)(g)(h)(1) Maximum Likelihood Approach. The 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 and presented in Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.Table 1
The results of MLEs and Bayesian estimates.
MLEBayesianαβ1β2β3αβ1β2β3α=0.2,β1=0.2, β2=0.2,β3=0.2n = 30ABias0.0200.0150.0150.0100.0160.0120.0080.012MSE0.0030.0040.0040.0040.0030.0040.0040.003CW (CL)0.2140.2480.2480.2480.2130.2430.2390.230CP0.9500.9380.9570.9480.9430.9320.9490.951n = 50ABias0.0110.0090.0080.0070.0100.0090.0040.009MSE0.0020.0020.0020.0020.0020.0020.0020.002CW (CL)0.1760.1760.1760.1760.1620.1880.1850.179CP0.9540.9400.9380.9510.9560.9430.9490.942n = 70ABias0.0100.0060.0060.0060.0080.0060.0040.007MSE0.0010.0020.0020.0010.0010.0020.0020.001CW (CL)0.1240.1760.1760.1240.1360.1570.1560.143CP0.9520.9570.9460.9510.9580.9450.9490.948α=0.5,β1=0.5, β2=0.5,β3=0.5n = 30ABias0.0320.0480.0460.0320.0280.0410.0270.038MSE0.0110.0340.0340.0260.0110.0310.0300.025CW (CL)0.3920.7020.7020.6200.3940.6640.6590.587CP0.9470.9450.9560.9510.9420.9300.9450.951n = 50ABias0.0180.0280.0260.0220.0200.0270.0260.025MSE0.0060.0180.0180.0140.0060.0170.0180.014CW (CL)0.3040.5120.5120.4640.2990.5010.5040.448CP0.9480.9440.9450.9570.9540.9490.9460.937n = 70ABias0.0170.0190.0210.0190.0130.0210.0170.019MSE0.0040.0120.0120.0100.0040.0120.0120.010CW (CL)0.2480.4300.4300.3920.2500.4170.4180.375CP0.9450.9580.9480.9500.9510.9290.9510.946α=0.8,β1=0.8, β2=0.8,β3=0.8n = 30ABias0.0450.0880.0860.0660.0370.0720.0730.060MSE0.0220.1010.0990.0760.0220.0900.0910.068CW (CL)0.5541.1961.1881.0520.5571.1291.1330.979CP0.9410.9510.9540.9460.9400.9390.9500.951n = 50ABias0.0260.0510.0490.0400.0260.0590.0490.040MSE0.0130.0520.0510.0400.0130.0520.0500.038CW (CL)0.4300.8680.8680.7640.4250.8540.8500.744CP0.9500.9440.9410.9560.9560.9480.9480.940n = 70ABias0.0240.0360.0390.0350.0190.0370.0300.030MSE0.0090.0340.0360.0280.0090.0340.0330.026CW (CL)0.3500.7120.7220.6440.3580.7030.7030.619CP0.9440.9560.9490.9490.9580.9400.9520.946α=1.5,β1=1.2, β2=1.3,β3=2.5n = 30ABias0.0680.1420.1580.3070.0230.1200.1620.137MSE0.0620.3910.4420.9130.0590.3550.4470.664CW (CL)0.9362.3882.5323.5480.9242.2072.4173.063CP0.9500.9400.9500.9500.9400.9500.9300.960n = 50ABias0.0460.0990.1110.1910.0180.0760.1030.091MSE0.0350.2030.2310.4440.0340.1890.2190.369CW (CL)0.7121.7221.8342.5040.7031.6331.7382.043CP0.9600.9500.9500.9600.9600.9200.9500.970n = 70ABias0.0390.0680.0820.1570.0130.0680.0510.090MSE0.0260.1340.1530.3010.0230.1310.1390.258CW (CL)0.6081.4081.4981.6270.5861.3641.4161.941CP0.9300.9400.9400.9600.9100.9400.9500.960α=2,β1=2, β2=2,β3=2n = 30ABias0.0940.2960.2950.2370.0600.2640.2470.147MSE0.0990.8900.8920.6800.0940.8110.7900.519CW (CL)1.1763.5103.5183.0961.1663.3123.2372.723CP0.9430.9490.9500.9520.9490.9310.9480.952n = 50ABias0.0550.1690.1660.1360.0400.1850.1690.101MSE0.0550.4270.4280.3260.0550.4320.4180.290CW (CL)1.1763.5103.5183.0960.8982.4432.3882.050CP0.9530.9460.9410.9610.9520.9460.9450.952n = 70ABias0.0510.1230.1190.0510.0310.1290.1010.063MSE0.0400.2810.2850.2260.0380.2770.2630.195CW (CL)0.7542.0222.0301.8040.7531.9791.9481.695CP0.9420.9530.9450.9420.9610.9420.9490.942α=2,β1=2.5, β2=3,β3=3.5n = 30ABias0.0860.3830.4930.5070.0370.2680.4260.178MSE0.0891.8322.5832.4760.2191.4972.4110.908CW (CL)1.1225.0885.9965.8400.8434.5435.5633.387CP0.9500.9500.9500.9500.9420.9500.8800.938n = 50ABias0.0530.2260.2790.2970.0190.1970.3040.095MSE0.0510.8701.1941.1320.0490.8151.1970.865CW (CL)0.8583.5484.1414.0060.8423.3763.9963.541CP0.9500.9400.9400.9700.9300.9100.9300.950n = 70ABias0.0490.1730.2360.2440.0140.1660.1490.091MSE0.0360.5740.8000.7660.0330.5430.6880.607CW (CL)0.7122.8943.3823.2940.6972.7703.1162.997CP0.9400.9400.9500.9600.9500.9300.9300.960α=2.5,β1=2.5, β2=3,β3=2.5n = 30ABias0.1110.3900.4890.3260.0670.3740.4180.168MSE0.1401.5352.1911.2340.1341.4351.9380.891CW (CL)1.4024.6105.4764.1641.3954.3755.0163.545CP0.9430.9500.9540.9470.9460.9030.9150.955n = 50ABias0.0680.2320.2790.1830.0410.2360.2620.108MSE0.0790.7371.0170.5770.0780.7250.9770.496CW (CL)1.0663.2403.7982.8921.0743.1673.6572.697CP0.9500.9400.9490.9610.9490.9170.8860.952n = 70ABias0.0600.1640.2140.1600.0280.1660.1660.069MSE0.0570.4790.6690.3990.0370.4720.6220.333CW (CL)0.9022.6363.0942.3940.9032.5812.9722.217CP0.9480.9500.9500.9430.9660.8870.8830.950α=3.5,β1=4, β2=3.5,β3=3.5n = 30ABias0.1520.7170.6070.5330.0680.6630.5360.115MSE0.2574.5743.5522.8880.2374.1683.1271.682CW (CL)1.8967.8986.9946.3221.8737.4036.3564.974CP0.9380.9490.9490.9540.8360.9280.9550.956n = 50ABias0.0910.4140.3520.2890.0530.4380.3570.133MSE0.1432.0751.6371.2770.1412.0801.6061.038CW (CL)1.4405.4104.8224.2821.4485.3114.6413.905CP0.9520.9470.9340.9660.7610.9190.9320.963n = 70ABias0.0820.3090.2680.2590.0310.2830.2320.060MSE0.1031.3521.0660.8830.0991.2971.0120.701CW (CL)1.2144.3943.9083.5421.2224.2753.7733.240CP0.9400.9510.9450.9530.6960.9280.9570.953(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.The 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θ¯X1,X2. The 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. The 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 statisticzn converges to normal distribution for large sample sizes. Hence, large absolute values of zn 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 Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.From Table1, it can be seen that under different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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. This is clear from Figure2.Figure 2
MSE of (a)α=2, (b) β1=2, (c) β2=2, and (d) β3=2.
(a)(b)(c)(d)Also, it can be seen that as the sample sizen 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.Figure 3
(a) MSE and (b) ABias forα.
(a)(b)


### 6.1.2. 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. The dataset has been obtained from Meintanis [23]; he used the bivariate MO exponential distribution (BE) to analyze the data. The 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 (X1) and the time of the first goal of any type scored by the home team (X2).The bivariate dataset has the following structure:X1 < X2, X1 = X2, and X1 > X2. Since, X1=X2 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. This is not going to make any difference in the statistical inference.First, before analyzing the data, we fit inverted Kumaraswamy distribution toX1, X2 and max (X1,X2). To guess the initial values for the parameters of BIK model, the MLEs of the shape parameters α,β of the respective inverted Kumaraswamy distribution for X1, X2 and max (X1,X2) 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 variableX in a population. The 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 F0x, then the null hypothesis is(42)H0:Fx=F0x.One approach for applying GOF tests is based on the empirical distribution function (EDF)Fnx which is defined as follows:(43)Fnx=number of elements in the sample≤xtotal sample sizen.This approach is based on defining a statistic to measure the discrepancy betweenFnx and F0x. The most used statistics are modified Cramér–von Mises statistic W∗ and Anderson–Darling statistic A∗ which have the following formulas:(44)W∗=112n+∑i=1nF0xi−2i−12n2,A∗=−n−2∑i=1n2i−12nlnF0xi+1−2i−12nln1−F0xi,where xi is the value of the ith order statistics in the sample.Large values of test statistics (or small correspondingp 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, and max (X1,X2) are illustrated in Table 2.Table 2
Cramér–von Mises statisticW∗ and Anderson–Darling statistic A∗ for IK.
X1X2maxX1,X2A∗0.966 (0.375)0.549 (0.696)1.395 (0.204)W∗0.139 (0.426)0.083 (0.679)0.233 (0.212)Based on the values of these statistics and the correspondingp 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.Table 3
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE5.2581.7630.5601.411Confidence interval width2.8242.0260.8861.374BayesianPosterior mean4.7271.6730.6791.442Credible interval length2.1681.2860.4911.061Now, 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:(45)IC=−2logL+kp,where k is the penalized term and p refers to number of parameters in the model.The 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 Table4.Table 4
Information criteria.
Log (L)AICBICCAICHQICk—2lnnlnn+1lnlnnBIK−23.53∗55.0661.561.6157.33BE−44.56∗∗95.1299.95100.0396.82BIW−30.25∗∗∗68.574.9475.0570.77∗Based on the estimates in Table 3. ∗∗Based on the estimates obtained by Meintanis [23]. ∗∗∗Based on the estimates presented by Muhammed [9].By analyzing equation (45), we can see that the first term −2logL 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. The model with the lowest IC is the best (for more details about IC, see Vrieze [25]).The log likelihood value, AIC, BIC, CAIC, and HQIC are represented in Table4. All of the criteria suggest that BIK provides the best fit compared with BE and BIW models. This may show the importance of BIK.(2) Motor Data. The dataset has been obtained from [26]. The data are about failure time in days for a parallel system containing two motors. The variables are time to failure for the first motor X1 and time to failure for the second motor (X2). 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.Table 5
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗1.133 (0.294)0.294 (0.942)0.488 (0.756)W∗0.204 (0.261)0.036 (0.958)0.081 (0.692)IWA∗1.491 (0.179)0.874 (0.429)7.165 (0.0003)W∗0.267 (0.168)0.119 (0.505)1.509 (9.18 ×10−5)EA∗3.902 (0.009)3.363 (0.014)94.4 (3.3 ×10−5)W∗0.778 (0.007)0.672 (0.014)5.964 (2.2 ×10−16)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 Table6.Table 6
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE16.97712.60313.0670.01Confidence interval width10.50819.5920.04618.730BayesianPosterior mean16.95012.90013.2100.034Credible interval length9.02019.15220.7200.594(3) Diabetic Retinopathy Data. The dataset has been obtained from [27]. The 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. The variables are time to failure for treated eye (X1) and time to failure for untreated eye (X2). 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.Table 7
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗0.662 (0.591)0.570 (0.676)0.707 (0.551)W∗0.092 (0.627)0.099 (0.594)0.090 (0.635)IWA∗2.538 (0.048)73.841 (8.4 ×10−6)166.45 (8.4 ×10−6)W∗0.369 (0.087)11.41 (2.2 ×10−16)18.25 (2.2 ×10−16)EA∗3.886 (0.009)214.7 (8.4 ×10−6)271.7 (8.4 ×10−6)W∗0.671 (0.015)22.4 (2.2 ×10−16)23.1 (2.2 ×10−16)In Table8, 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.Table 8
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE3.8302.8892.5700.900Confidence interval width1.2882.1141.9521.044BayesianPosterior mean3.5862.7392.37960.907Credible interval length0.9301.3261.1330.700From these three datasets, we can conclude that the derived BIK distribution will be of great importance.


## 6.1.1. A Simulation Study

Here, we present an algorithm to generate BIK distribution (Algorithm1). To perform a simulation study, we first need to select initial values for the parameters. Here, the following eight different populations are considered:(i)
Case 1:α=0.2,β1=0.2,β2=0.2,β3=0.2.(ii)
Case 2:α=0.5,β1=0.5,β2=0.5,β3=0.5.(iii)
Case 3:α=0.8,β1=0.8,β2=0.8,β3=0.8.(iv)
Case 4:α=1.5,β1=1.2,β2=1.3,β3=2.5.(v)
Case 5:α=2,β1=2,β2=2,β3=2.(vi)
Case 6:α=2,β1=2.5,β2=3,β3=3.5.(vii)
Case 7:α=2.5,β1=2.5,β2=3,β3=2.5.(viii)
Case 8:α=3.5,β1=4,β2=3.5,β3=3.5.Algorithm 1: An Algorithm to Generate BIK Distribution.
The following algorithm is to generateX1,X2 from BIK distribution.Step 1: generateU1,U2, and U3 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, and T3=1−U31/β31/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, and Z3=1/T3−1.Step 4: obtainX1=maxZ1,Z3,X2=maxZ2,Z3.The parameters are selected to cover different shapes of the distribution. It can be seen from Figure1 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.Figure 1
Different shapes of the absolutely continuous part of BIK pdf. (a)α = 0.2, β1 = 0.2, β2 = 0.2, β3 = 0.2; (b) α = 0.5, β1 = 0.5, β2 = 0.5, β3 = 0.5; (c) α = 0.8, β1 = 0.8, β2 = 0.8, β3 = 0.8; (d) α = 2, β1 = 2, β2 = 2, β3 = 2; (e) α = 1.5, β1 = 1.2, β2 = 1.3, β3 = 2.5; (f) α = 2, β1 = 2.5, β2 = 3, β3 = 3.5; (g) α = 2.5, β1 = 2.5, β2 = 3, β3 = 2.5; (h) α = 3.5, β1 = 4, β2 = 3.5, β3 = 3.5.
(a)(b)(c)(d)(e)(f)(g)(h)(1) Maximum Likelihood Approach. The 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 and presented in Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.Table 1
The results of MLEs and Bayesian estimates.
MLEBayesianαβ1β2β3αβ1β2β3α=0.2,β1=0.2, β2=0.2,β3=0.2n = 30ABias0.0200.0150.0150.0100.0160.0120.0080.012MSE0.0030.0040.0040.0040.0030.0040.0040.003CW (CL)0.2140.2480.2480.2480.2130.2430.2390.230CP0.9500.9380.9570.9480.9430.9320.9490.951n = 50ABias0.0110.0090.0080.0070.0100.0090.0040.009MSE0.0020.0020.0020.0020.0020.0020.0020.002CW (CL)0.1760.1760.1760.1760.1620.1880.1850.179CP0.9540.9400.9380.9510.9560.9430.9490.942n = 70ABias0.0100.0060.0060.0060.0080.0060.0040.007MSE0.0010.0020.0020.0010.0010.0020.0020.001CW (CL)0.1240.1760.1760.1240.1360.1570.1560.143CP0.9520.9570.9460.9510.9580.9450.9490.948α=0.5,β1=0.5, β2=0.5,β3=0.5n = 30ABias0.0320.0480.0460.0320.0280.0410.0270.038MSE0.0110.0340.0340.0260.0110.0310.0300.025CW (CL)0.3920.7020.7020.6200.3940.6640.6590.587CP0.9470.9450.9560.9510.9420.9300.9450.951n = 50ABias0.0180.0280.0260.0220.0200.0270.0260.025MSE0.0060.0180.0180.0140.0060.0170.0180.014CW (CL)0.3040.5120.5120.4640.2990.5010.5040.448CP0.9480.9440.9450.9570.9540.9490.9460.937n = 70ABias0.0170.0190.0210.0190.0130.0210.0170.019MSE0.0040.0120.0120.0100.0040.0120.0120.010CW (CL)0.2480.4300.4300.3920.2500.4170.4180.375CP0.9450.9580.9480.9500.9510.9290.9510.946α=0.8,β1=0.8, β2=0.8,β3=0.8n = 30ABias0.0450.0880.0860.0660.0370.0720.0730.060MSE0.0220.1010.0990.0760.0220.0900.0910.068CW (CL)0.5541.1961.1881.0520.5571.1291.1330.979CP0.9410.9510.9540.9460.9400.9390.9500.951n = 50ABias0.0260.0510.0490.0400.0260.0590.0490.040MSE0.0130.0520.0510.0400.0130.0520.0500.038CW (CL)0.4300.8680.8680.7640.4250.8540.8500.744CP0.9500.9440.9410.9560.9560.9480.9480.940n = 70ABias0.0240.0360.0390.0350.0190.0370.0300.030MSE0.0090.0340.0360.0280.0090.0340.0330.026CW (CL)0.3500.7120.7220.6440.3580.7030.7030.619CP0.9440.9560.9490.9490.9580.9400.9520.946α=1.5,β1=1.2, β2=1.3,β3=2.5n = 30ABias0.0680.1420.1580.3070.0230.1200.1620.137MSE0.0620.3910.4420.9130.0590.3550.4470.664CW (CL)0.9362.3882.5323.5480.9242.2072.4173.063CP0.9500.9400.9500.9500.9400.9500.9300.960n = 50ABias0.0460.0990.1110.1910.0180.0760.1030.091MSE0.0350.2030.2310.4440.0340.1890.2190.369CW (CL)0.7121.7221.8342.5040.7031.6331.7382.043CP0.9600.9500.9500.9600.9600.9200.9500.970n = 70ABias0.0390.0680.0820.1570.0130.0680.0510.090MSE0.0260.1340.1530.3010.0230.1310.1390.258CW (CL)0.6081.4081.4981.6270.5861.3641.4161.941CP0.9300.9400.9400.9600.9100.9400.9500.960α=2,β1=2, β2=2,β3=2n = 30ABias0.0940.2960.2950.2370.0600.2640.2470.147MSE0.0990.8900.8920.6800.0940.8110.7900.519CW (CL)1.1763.5103.5183.0961.1663.3123.2372.723CP0.9430.9490.9500.9520.9490.9310.9480.952n = 50ABias0.0550.1690.1660.1360.0400.1850.1690.101MSE0.0550.4270.4280.3260.0550.4320.4180.290CW (CL)1.1763.5103.5183.0960.8982.4432.3882.050CP0.9530.9460.9410.9610.9520.9460.9450.952n = 70ABias0.0510.1230.1190.0510.0310.1290.1010.063MSE0.0400.2810.2850.2260.0380.2770.2630.195CW (CL)0.7542.0222.0301.8040.7531.9791.9481.695CP0.9420.9530.9450.9420.9610.9420.9490.942α=2,β1=2.5, β2=3,β3=3.5n = 30ABias0.0860.3830.4930.5070.0370.2680.4260.178MSE0.0891.8322.5832.4760.2191.4972.4110.908CW (CL)1.1225.0885.9965.8400.8434.5435.5633.387CP0.9500.9500.9500.9500.9420.9500.8800.938n = 50ABias0.0530.2260.2790.2970.0190.1970.3040.095MSE0.0510.8701.1941.1320.0490.8151.1970.865CW (CL)0.8583.5484.1414.0060.8423.3763.9963.541CP0.9500.9400.9400.9700.9300.9100.9300.950n = 70ABias0.0490.1730.2360.2440.0140.1660.1490.091MSE0.0360.5740.8000.7660.0330.5430.6880.607CW (CL)0.7122.8943.3823.2940.6972.7703.1162.997CP0.9400.9400.9500.9600.9500.9300.9300.960α=2.5,β1=2.5, β2=3,β3=2.5n = 30ABias0.1110.3900.4890.3260.0670.3740.4180.168MSE0.1401.5352.1911.2340.1341.4351.9380.891CW (CL)1.4024.6105.4764.1641.3954.3755.0163.545CP0.9430.9500.9540.9470.9460.9030.9150.955n = 50ABias0.0680.2320.2790.1830.0410.2360.2620.108MSE0.0790.7371.0170.5770.0780.7250.9770.496CW (CL)1.0663.2403.7982.8921.0743.1673.6572.697CP0.9500.9400.9490.9610.9490.9170.8860.952n = 70ABias0.0600.1640.2140.1600.0280.1660.1660.069MSE0.0570.4790.6690.3990.0370.4720.6220.333CW (CL)0.9022.6363.0942.3940.9032.5812.9722.217CP0.9480.9500.9500.9430.9660.8870.8830.950α=3.5,β1=4, β2=3.5,β3=3.5n = 30ABias0.1520.7170.6070.5330.0680.6630.5360.115MSE0.2574.5743.5522.8880.2374.1683.1271.682CW (CL)1.8967.8986.9946.3221.8737.4036.3564.974CP0.9380.9490.9490.9540.8360.9280.9550.956n = 50ABias0.0910.4140.3520.2890.0530.4380.3570.133MSE0.1432.0751.6371.2770.1412.0801.6061.038CW (CL)1.4405.4104.8224.2821.4485.3114.6413.905CP0.9520.9470.9340.9660.7610.9190.9320.963n = 70ABias0.0820.3090.2680.2590.0310.2830.2320.060MSE0.1031.3521.0660.8830.0991.2971.0120.701CW (CL)1.2144.3943.9083.5421.2224.2753.7733.240CP0.9400.9510.9450.9530.6960.9280.9570.953(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.The 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θ¯X1,X2. The 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. The 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 statisticzn converges to normal distribution for large sample sizes. Hence, large absolute values of zn 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 Table1. The numerical steps and the corresponding equations are explained in detail in Appendix B.1.From Table1, it can be seen that under different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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. This is clear from Figure2.Figure 2
MSE of (a)α=2, (b) β1=2, (c) β2=2, and (d) β3=2.
(a)(b)(c)(d)Also, it can be seen that as the sample sizen 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.Figure 3
(a) MSE and (b) ABias forα.
(a)(b)


## 6.1.2. 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. The dataset has been obtained from Meintanis [23]; he used the bivariate MO exponential distribution (BE) to analyze the data. The 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 (X1) and the time of the first goal of any type scored by the home team (X2).The bivariate dataset has the following structure:X1 < X2, X1 = X2, and X1 > X2. Since, X1=X2 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. This is not going to make any difference in the statistical inference.First, before analyzing the data, we fit inverted Kumaraswamy distribution toX1, X2 and max (X1,X2). To guess the initial values for the parameters of BIK model, the MLEs of the shape parameters α,β of the respective inverted Kumaraswamy distribution for X1, X2 and max (X1,X2) 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 variableX in a population. The 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 F0x, then the null hypothesis is(42)H0:Fx=F0x.One approach for applying GOF tests is based on the empirical distribution function (EDF)Fnx which is defined as follows:(43)Fnx=number of elements in the sample≤xtotal sample sizen.This approach is based on defining a statistic to measure the discrepancy betweenFnx and F0x. The most used statistics are modified Cramér–von Mises statistic W∗ and Anderson–Darling statistic A∗ which have the following formulas:(44)W∗=112n+∑i=1nF0xi−2i−12n2,A∗=−n−2∑i=1n2i−12nlnF0xi+1−2i−12nln1−F0xi,where xi is the value of the ith order statistics in the sample.Large values of test statistics (or small correspondingp 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, and max (X1,X2) are illustrated in Table 2.Table 2
Cramér–von Mises statisticW∗ and Anderson–Darling statistic A∗ for IK.
X1X2maxX1,X2A∗0.966 (0.375)0.549 (0.696)1.395 (0.204)W∗0.139 (0.426)0.083 (0.679)0.233 (0.212)Based on the values of these statistics and the correspondingp 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.Table 3
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE5.2581.7630.5601.411Confidence interval width2.8242.0260.8861.374BayesianPosterior mean4.7271.6730.6791.442Credible interval length2.1681.2860.4911.061Now, 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:(45)IC=−2logL+kp,where k is the penalized term and p refers to number of parameters in the model.The 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 Table4.Table 4
Information criteria.
Log (L)AICBICCAICHQICk—2lnnlnn+1lnlnnBIK−23.53∗55.0661.561.6157.33BE−44.56∗∗95.1299.95100.0396.82BIW−30.25∗∗∗68.574.9475.0570.77∗Based on the estimates in Table 3. ∗∗Based on the estimates obtained by Meintanis [23]. ∗∗∗Based on the estimates presented by Muhammed [9].By analyzing equation (45), we can see that the first term −2logL 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. The model with the lowest IC is the best (for more details about IC, see Vrieze [25]).The log likelihood value, AIC, BIC, CAIC, and HQIC are represented in Table4. All of the criteria suggest that BIK provides the best fit compared with BE and BIW models. This may show the importance of BIK.(2) Motor Data. The dataset has been obtained from [26]. The data are about failure time in days for a parallel system containing two motors. The variables are time to failure for the first motor X1 and time to failure for the second motor (X2). 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.Table 5
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗1.133 (0.294)0.294 (0.942)0.488 (0.756)W∗0.204 (0.261)0.036 (0.958)0.081 (0.692)IWA∗1.491 (0.179)0.874 (0.429)7.165 (0.0003)W∗0.267 (0.168)0.119 (0.505)1.509 (9.18 ×10−5)EA∗3.902 (0.009)3.363 (0.014)94.4 (3.3 ×10−5)W∗0.778 (0.007)0.672 (0.014)5.964 (2.2 ×10−16)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 Table6.Table 6
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE16.97712.60313.0670.01Confidence interval width10.50819.5920.04618.730BayesianPosterior mean16.95012.90013.2100.034Credible interval length9.02019.15220.7200.594(3) Diabetic Retinopathy Data. The dataset has been obtained from [27]. The 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. The variables are time to failure for treated eye (X1) and time to failure for untreated eye (X2). 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.Table 7
W∗ and A∗ for IK, IW, and E.
ModelX1X2maxX1,X2IKA∗0.662 (0.591)0.570 (0.676)0.707 (0.551)W∗0.092 (0.627)0.099 (0.594)0.090 (0.635)IWA∗2.538 (0.048)73.841 (8.4 ×10−6)166.45 (8.4 ×10−6)W∗0.369 (0.087)11.41 (2.2 ×10−16)18.25 (2.2 ×10−16)EA∗3.886 (0.009)214.7 (8.4 ×10−6)271.7 (8.4 ×10−6)W∗0.671 (0.015)22.4 (2.2 ×10−16)23.1 (2.2 ×10−16)In Table8, 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.Table 8
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE3.8302.8892.5700.900Confidence interval width1.2882.1141.9521.044BayesianPosterior mean3.5862.7392.37960.907Credible interval length0.9301.3261.1330.700From these three datasets, we can conclude that the derived BIK distribution will be of great importance.


## 6.2. For MIK Distribution

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


### 6.2.1. A Simulation Study

Here, we present an algorithm to generate MIK distribution. Also, we illustrate the simulation results for both ML and Bayesian approaches (Algorithm2).Algorithm 2: An algorithm to generate MIK distribution.
The following algorithm is to generateX1,X2,X3 from MIK distribution.Step 1: generateU1,U2,U3, and U4 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, T3=1−U31/β31/α, and T4=1−U41/β41/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, Z3=1/T3−1, and Z4=1/T4−1.Step 4: obtainX1=maxZ1,Z4,X2=maxZ2,Z4, and X3=maxZ3,Z4.(1) Maximum Likelihood Approach. To obtain the maximum likelihood estimates, a Monte Carlo simulation is performed using R package with 1000 replications and three different sample sizes n=30,50, and 70. Different initial values of the parameters are arbitrarily chosen, varying between small and large values to cover different cases of the distribution. The parameter β4 is assumed known for simplicity. The following twelve populations are considered.(i)
α=0.9,β1=0.9,β2=0.9,β3=0.9.(ii)
α=1.8,β1=0.8,β2=1.8,β3=0.8.(iii)
α=1.8,β1=1.8,β2=1.8,β3=1.8.(iv)
α=2,β1=1.5,β2=2,β3=2.(v)
α=2.4,β1=2.4,β2=2.4,β3=2.4.(vi)
α=2.5,β1=1.5,β2=1.5,β3=1.5.(vii)
α=3,β1=3,β2=3,β3=3.(viii)
α=3,β1=2.8,β2=2.8,β3=3.5.(ix)
α=3,β1=3,β2=4,β3=3.(x)
α=3.5,β1=3.5,β2=3.5,β3=4.(xi)
α=3.5,β1=3.5,β2=3.5,β3=3.5.(xii)
α=4,β1=1.5,β2=2,β3=2.5.To check the behavior of the estimates, ABias, MSE, CW, and CP are computed in Table9. The algorithm is explained in detail in Appendix B.2.Table 9
The results of MLEs and Bayesian estimates for MIK.
MLEBayesianαβ1β2β3αβ1β2β3α=0.9,β1=0.9, β2=0.9,β3=0.9n = 30ABias0.0460.0630.0500.0610.0450.0370.0200.032MSE0.0200.0450.0420.0450.0190.0400.0370.039CW (CL)0.5260.7940.7740.7940.5050.7610.7350.749CP0.9530.9730.9540.9490.9390.9250.9480.933n = 50ABias0.0350.0330.0460.0400.0390.0050.0200.019MSE0.0110.0230.0250.0250.0120.0210.0220.023CW (CL)0.3920.5820.5940.5940.3910.5670.5730.578CP0.9600.9680.9500.9600.9450.9460.9290.943n = 70ABias0.0200.0260.0210.0320.0290.0010.0050.006MSE0.0070.0170.0150.0170.0080.0150.0150.015CW (CL)0.3280.4960.4800.4960.3280.4770.4720.479CP0.9450.9510.9580.9640.9300.9310.9330.950α=1.8,β1=0.8, β2=1.8,β3=0.8n = 30ABias0.1350.0830.2240.0890.1190.0430.1990.046MSE0.0840.0440.2820.0460.0780.0330.2480.034CW (CL)1.0080.7541.8880.7640.9770.6851.7560.687CP0.9180.9610.9730.9520.9270.9400.9480.935n = 50ABias0.1110.0670.1550.0620.0850.0150.1670.022MSE0.0510.0250.1480.0250.0440.0180.1480.018CW (CL)0.7740.5681.3800.5680.7490.5111.3350.516CP0.9390.9620.9500.9580.9350.9560.9250.938n = 70ABias0.0860.0570.1550.0490.0720.0110.1200.015MSE0.0340.0170.1130.0160.0310.0120.0940.013CW (CL)0.6440.4641.1700.4640.6280.4321.0190.433CP0.9140.9690.9300.9770.9190.9510.9280.948α=1.8,β1=1.8, β2=1.8,β3=1.8n = 30ABias0.1120.2100.1730.2190.0900.1950.1690.167MSE0.0640.2620.2360.2690.0600.2450.2290.224CW (CL)0.8861.8301.7801.8420.8841.7631.7251.718CP0.9210.9620.9720.9620.9170.9420.9530.934n = 50ABias0.0810.1370.1500.1470.0740.1300.1620.115MSE0.0370.1340.1410.1390.0360.1280.1430.123CW (CL)0.6781.3301.3461.3400.6771.2991.3241.288CP0.9510.9720.9440.9650.9440.9640.9220.952n = 70ABias0.0590.1130.1160.1220.0540.0970.1280.091MSE0.0240.0920.0930.0950.0240.0860.0950.084CW (CL)0.5681.1021.1081.3680.5621.0751.0901.074CP0.9320.9520.9490.9590.9410.9410.9280.946α=2,β1=1.5, β2=2,β3=2n = 30ABias0.1190.1590.2080.2730.0840.1300.1890.196MSE0.0760.1790.3100.3670.0660.1430.2800.286CW (CL)0.9761.4782.0262.1180.9401.3711.9001.924CP0.9230.9640.9790.9640.9230.9450.9520.945n = 50ABias0.0920.1030.1730.1760.0750.0920.1570.124MSE0.0440.0860.1810.0670.0410.0800.1650.154CW (CL)0.7441.0741.5241.5240.7261.0361.4441.435CP0.9480.9690.9440.9650.9450.9630.9270.954n = 70ABias0.0650.0860.1220.1350.0450.0620.1040.097MSE0.0290.0550.1150.1190.0270.0530.1050.104CW (CL)0.6200.8941.2401.2560.6060.8561.1881.189CP0.9270.9500.9570.9670.9330.9450.9500.950α=2.4,β1=2.4, β2=2.4,β3=2.4n = 30ABias0.2840.2450.2980.1330.1060.2640.2400.209MSE0.4940.4560.5080.0980.0930.4630.4410.410CW (CL)2.5202.4662.5381.1081.1102.4242.3662.340CP0.9670.9730.9670.9340.9430.9480.9540.950n = 50ABias0.1840.2040.1980.0950.0700.1560.2000.132MSE0.2510.2540.2600.0560.0520.2320.2570.219CW (CL)1.8261.8561.8420.8500.8461.7791.8041.751CP0.9750.9520.9640.9520.9550.9660.9250.955n = 70ABias0.1460.1360.1570.0700.0590.1490.1390.117MSE0.1700.1660.1760.0380.0370.1680.1650.156CW (CL)1.5141.5081.5240.7120.7131.4911.4811.469CP0.9440.9520.9620.9390.9380.9280.9500.956α=2.5,β1=1.5, β2=1.5,β3=1.5n = 30ABias0.1680.1680.1510.1750.1460.1510.1410.138MSE0.1350.1730.1630.1770.1240.1550.1490.149CW (CL)1.2821.4921.4661.4981.2401.4081.3881.387CP0.9310.9610.9790.9630.9190.9430.9550.950n = 50ABias0.1250.1200.1190.1250.0960.1020.1020.094MSE0.0780.0920.1150.0950.0690.0830.1520.082CW (CL)0.9761.0941.1021.1020.9471.0471.0431.044CP0.9460.9690.9520.9660.9490.9560.9230.953n = 70ABias0.0920.1000.0940.0940.0860.0920.1050.062MSE0.0510.0640.0620.0620.0500.0600.0630.054CW (CL)0.8120.9100.9020.9020.7990.8760.8780.862CP0.9130.9420.9500.9710.9270.9490.9210.966α=3,β1=3, β2=3,β3=3n = 30ABias0.1570.3770.3180.3740.1170.3610.2930.249MSE0.1410.8370.7950.8330.1310.7960.7160.658CW (CL)1.3363.2683.1803.2641.3283.1503.0252.975CP0.9360.9660.9740.9690.9440.9320.9610.956n = 50ABias0.1210.2420.2830.2580.0960.2140.2810.174MSE0.0830.4200.4560.4340.0790.3970.4510.367CW (CL)1.0222.3562.4042.3741.0242.3062.3452.248CP0.9410.9750.9510.9630.9390.9620.9290.958n = 70ABias0.0740.1650.1590.1880.0580.1500.1580.137MSE0.0520.2690.2660.2810.0510.2590.2640.249CW (CL)0.8501.9281.9241.9440.8551.8931.8921.876CP0.9370.9510.9610.9630.9430.9380.9480.956α=3,β1=2.8,β2=2.8,β3=3.5n = 30ABias0.1540.3390.2830.4550.1030.2840.2430.257MSE0.1390.7000.6331.2210.1260.6210.5820.896CW (CL)1.3302.9982.9163.9481.3142.8342.7553.513CP0.9370.9640.9770.9700.8640.9380.8470.914n = 50ABias0.1210.2180.2520.3240.0870.1690.2330.210MSE0.0830.3540.3810.6420.0760.3200.3620.529CW (CL)1.0222.1682.2082.8721.0182.0982.1362.697CP0.9470.9730.9550.9650.8050.9740.7970.857n = 70ABias0.0810.1560.1480.2380.0530.1210.1330.137MSE0.0540.2300.2270.4140.0500.2110.2170.344CW (CL)0.8501.7801.7742.3420.8491.7261.7332.220CP0.9330.9510.9650.9630.8000.9470.7060.839α=3,β1=3,β2=4,β3=3n = 30ABias0.1540.3740.4770.3700.1090.3270.4550.230MSE0.1360.8271.5580.8220.1200.7261.4410.617CW (CL)1.3123.2504.5203.2441.2843.0484.3072.914CP0.9370.9630.9730.9710.9410.9510.9700.963n = 50ABias0.1160.2320.4060.2560.0910.2060.3880.171MSE0.0790.4090.9150.4300.0730.3800.8510.354CW (CL)1.0082.3363.3942.3660.9892.2543.2542.213CP0.9420.9730.9600.9650.9420.9530.9650.952n = 70ABias0.0740.1630.2340.1820.0580.1560.2280.123MSE0.0510.2660.5300.2760.0490.2570.5130.241CW (CL)0.8401.9162.7021.9320.8341.8692.6361.847CP0.9430.9560.9600.9600.9550.9650.9670.958α=3.5,β1=3.5,β2=3.5,β3=4n = 30ABias0.1710.4500.3780.5360.1000.4020.3890.238MSE0.1761.1901.0741.6660.1571.0951.0711.149CW (CL)1.5023.8943.7824.6041.4863.7023.6604.027CP0.9390.9630.9750.9690.9490.9520.9310.937n = 50ABias0.1320.2810.3300.3710.0800.2530.2680.166MSE0.1040.5870.6390.8590.0930.5620.5730.653CW (CL)1.1562.7942.8543.3281.1462.7392.7153.059CP0.9440.9720.9590.9660.9360.9610.9000.903n = 70ABias0.0850.1970.1850.2610.0550.1450.1940.159MSE0.0680.3810.3720.5450.0650.3510.3780.465CW (CL)0.9682.2922.2802.7080.9672.2302.2582.583CP0.9400.9520.9630.9700.9500.9590.8710.864α=3.5,β1=3.5, β2=3.5,β3=3.5n = 30ABias0.1720.4500.3780.4450.1070.3690.3230.250MSE0.1791.1941.0791.1880.1621.0560.9920.885CW (CL)1.5143.9023.7923.9001.5033.6893.5933.488CP0.9410.9660.9750.9710.9550.9450.9590.955n = 50ABias0.1330.2840.3310.3080.1000.2630.3080.199MSE0.1060.5920.6430.6170.0990.5730.6100.516CW (CL)1.1622.3832.8622.8321.1602.7552.7552.665CP0.9460.9710.9590.9650.9460.9640.9290.949n = 70ABias0.0870.2040.1860.2240.0590.1630.1820.144MSE0.0690.3870.3750.3980.0660.3620.3740.346CW (CL)0.9682.3022.2862.3120.9782.2522.2592.220CP0.9370.9540.9630.9630.9410.9510.9560.956α=4,β1=1.5, β2=2,β1=2.5n = 30ABias0.2340.1520.1990.3300.1770.1430.1850.236MSE0.3510.1620.3020.5900.2720.1550.2840.474CW (CL)1.9281.4622.0062.7221.9061.4151.9122.486CP0.9320.9630.9740.9710.9500.9410.9450.954n = 50ABias0.1780.1000.1670.2310.1410.0920.1550.167MSE0.1730.0850.1770.3060.1620.0830.1690.260CW (CL)1.4721.0741.5141.9721.4681.0591.4681.871CP0.9440.9720.9500.9630.9600.9580.9530.929n = 70ABias0.1360.0820.1200.1930.1210.0730.1370.126MSE0.1170.0590.1130.2090.1160.0570.1200.177CW (CL)1.2340.8941.2341.6261.2410.8831.2361.559CP0.9290.9500.9620.9710.9430.9390.9460.935(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. The algorithm is explained in detail in Appendix B.2.From Table9, it can be seen that for different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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 Figure4. Also, for the majority of cases, Bayesian estimates have smaller ABias than ML ones.Figure 4
MSE of (a)α=3.5, (b) β1=3.5, (c) β2=3.5, and (d) β3=3.5.
(a)(b)(c)(d)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.Figure 5
(a) MSE and (b) ABias forα.
(a)(b)


### 6.2.2. A Real Dataset

Here, we analyze a real dataset to show the applicability of MIK. The 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. The variables are as follows:X1: first systolic blood pressure measurement.X2: second systolic blood pressure measurement.X3: third systolic blood pressure measurement.All data points have been divided by 1000. This 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, X3, and maxX1,X2,X3 are illustrated in Table 10.Table 10
W∗ and A∗ for IK, IW, and E.
ModelX1X2X3maxX1,X2IKA∗0.31 (0.93)0.43 (0.82)0.49 (0.76)0.50 (0.75)W∗0.04 (0.94)0.07 (0.76)0.08 (0.68)0.08 (0.72)IWA∗13.60 (7 ×10−6)14.48 (7 ×10−6)15.80 (7 ×10−6)14.20 (7 ×10−6)W∗2.79 (2 ×10−7)2.99 (5 ×10−8)3.30 (9 ×10−9)2.91 (7 ×10−8)EA∗23.95 (7 ×10−6)24.18 (7 ×10−6)24.75 (7 ×10−6)24.29 (7 ×10−6)W∗5.06 (2 ×10−16)5.12 (2 ×10−16)5.26 (2 ×10−16)5.14 (2 ×10−16)Based on the values of these statistics and the correspondingp 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.Table 11
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE46.949264.199245.852245.519Confidence interval width8.126269.684256.424263.866BayesianPosterior mean48.697270.528226.877257.231Credible interval length4.27069.21561.11064.305


## 6.2.1. A Simulation Study

Here, we present an algorithm to generate MIK distribution. Also, we illustrate the simulation results for both ML and Bayesian approaches (Algorithm2).Algorithm 2: An algorithm to generate MIK distribution.
The following algorithm is to generateX1,X2,X3 from MIK distribution.Step 1: generateU1,U2,U3, and U4 from uniform (0, 1).Step 2: computeT1=1−U11/β11/α,T2=1−U21/β21/α, T3=1−U31/β31/α, and T4=1−U41/β41/α.Step 3: defineZ1=1/T1−1,Z2=1/T2−1, Z3=1/T3−1, and Z4=1/T4−1.Step 4: obtainX1=maxZ1,Z4,X2=maxZ2,Z4, and X3=maxZ3,Z4.(1) Maximum Likelihood Approach. To obtain the maximum likelihood estimates, a Monte Carlo simulation is performed using R package with 1000 replications and three different sample sizes n=30,50, and 70. Different initial values of the parameters are arbitrarily chosen, varying between small and large values to cover different cases of the distribution. The parameter β4 is assumed known for simplicity. The following twelve populations are considered.(i)
α=0.9,β1=0.9,β2=0.9,β3=0.9.(ii)
α=1.8,β1=0.8,β2=1.8,β3=0.8.(iii)
α=1.8,β1=1.8,β2=1.8,β3=1.8.(iv)
α=2,β1=1.5,β2=2,β3=2.(v)
α=2.4,β1=2.4,β2=2.4,β3=2.4.(vi)
α=2.5,β1=1.5,β2=1.5,β3=1.5.(vii)
α=3,β1=3,β2=3,β3=3.(viii)
α=3,β1=2.8,β2=2.8,β3=3.5.(ix)
α=3,β1=3,β2=4,β3=3.(x)
α=3.5,β1=3.5,β2=3.5,β3=4.(xi)
α=3.5,β1=3.5,β2=3.5,β3=3.5.(xii)
α=4,β1=1.5,β2=2,β3=2.5.To check the behavior of the estimates, ABias, MSE, CW, and CP are computed in Table9. The algorithm is explained in detail in Appendix B.2.Table 9
The results of MLEs and Bayesian estimates for MIK.
MLEBayesianαβ1β2β3αβ1β2β3α=0.9,β1=0.9, β2=0.9,β3=0.9n = 30ABias0.0460.0630.0500.0610.0450.0370.0200.032MSE0.0200.0450.0420.0450.0190.0400.0370.039CW (CL)0.5260.7940.7740.7940.5050.7610.7350.749CP0.9530.9730.9540.9490.9390.9250.9480.933n = 50ABias0.0350.0330.0460.0400.0390.0050.0200.019MSE0.0110.0230.0250.0250.0120.0210.0220.023CW (CL)0.3920.5820.5940.5940.3910.5670.5730.578CP0.9600.9680.9500.9600.9450.9460.9290.943n = 70ABias0.0200.0260.0210.0320.0290.0010.0050.006MSE0.0070.0170.0150.0170.0080.0150.0150.015CW (CL)0.3280.4960.4800.4960.3280.4770.4720.479CP0.9450.9510.9580.9640.9300.9310.9330.950α=1.8,β1=0.8, β2=1.8,β3=0.8n = 30ABias0.1350.0830.2240.0890.1190.0430.1990.046MSE0.0840.0440.2820.0460.0780.0330.2480.034CW (CL)1.0080.7541.8880.7640.9770.6851.7560.687CP0.9180.9610.9730.9520.9270.9400.9480.935n = 50ABias0.1110.0670.1550.0620.0850.0150.1670.022MSE0.0510.0250.1480.0250.0440.0180.1480.018CW (CL)0.7740.5681.3800.5680.7490.5111.3350.516CP0.9390.9620.9500.9580.9350.9560.9250.938n = 70ABias0.0860.0570.1550.0490.0720.0110.1200.015MSE0.0340.0170.1130.0160.0310.0120.0940.013CW (CL)0.6440.4641.1700.4640.6280.4321.0190.433CP0.9140.9690.9300.9770.9190.9510.9280.948α=1.8,β1=1.8, β2=1.8,β3=1.8n = 30ABias0.1120.2100.1730.2190.0900.1950.1690.167MSE0.0640.2620.2360.2690.0600.2450.2290.224CW (CL)0.8861.8301.7801.8420.8841.7631.7251.718CP0.9210.9620.9720.9620.9170.9420.9530.934n = 50ABias0.0810.1370.1500.1470.0740.1300.1620.115MSE0.0370.1340.1410.1390.0360.1280.1430.123CW (CL)0.6781.3301.3461.3400.6771.2991.3241.288CP0.9510.9720.9440.9650.9440.9640.9220.952n = 70ABias0.0590.1130.1160.1220.0540.0970.1280.091MSE0.0240.0920.0930.0950.0240.0860.0950.084CW (CL)0.5681.1021.1081.3680.5621.0751.0901.074CP0.9320.9520.9490.9590.9410.9410.9280.946α=2,β1=1.5, β2=2,β3=2n = 30ABias0.1190.1590.2080.2730.0840.1300.1890.196MSE0.0760.1790.3100.3670.0660.1430.2800.286CW (CL)0.9761.4782.0262.1180.9401.3711.9001.924CP0.9230.9640.9790.9640.9230.9450.9520.945n = 50ABias0.0920.1030.1730.1760.0750.0920.1570.124MSE0.0440.0860.1810.0670.0410.0800.1650.154CW (CL)0.7441.0741.5241.5240.7261.0361.4441.435CP0.9480.9690.9440.9650.9450.9630.9270.954n = 70ABias0.0650.0860.1220.1350.0450.0620.1040.097MSE0.0290.0550.1150.1190.0270.0530.1050.104CW (CL)0.6200.8941.2401.2560.6060.8561.1881.189CP0.9270.9500.9570.9670.9330.9450.9500.950α=2.4,β1=2.4, β2=2.4,β3=2.4n = 30ABias0.2840.2450.2980.1330.1060.2640.2400.209MSE0.4940.4560.5080.0980.0930.4630.4410.410CW (CL)2.5202.4662.5381.1081.1102.4242.3662.340CP0.9670.9730.9670.9340.9430.9480.9540.950n = 50ABias0.1840.2040.1980.0950.0700.1560.2000.132MSE0.2510.2540.2600.0560.0520.2320.2570.219CW (CL)1.8261.8561.8420.8500.8461.7791.8041.751CP0.9750.9520.9640.9520.9550.9660.9250.955n = 70ABias0.1460.1360.1570.0700.0590.1490.1390.117MSE0.1700.1660.1760.0380.0370.1680.1650.156CW (CL)1.5141.5081.5240.7120.7131.4911.4811.469CP0.9440.9520.9620.9390.9380.9280.9500.956α=2.5,β1=1.5, β2=1.5,β3=1.5n = 30ABias0.1680.1680.1510.1750.1460.1510.1410.138MSE0.1350.1730.1630.1770.1240.1550.1490.149CW (CL)1.2821.4921.4661.4981.2401.4081.3881.387CP0.9310.9610.9790.9630.9190.9430.9550.950n = 50ABias0.1250.1200.1190.1250.0960.1020.1020.094MSE0.0780.0920.1150.0950.0690.0830.1520.082CW (CL)0.9761.0941.1021.1020.9471.0471.0431.044CP0.9460.9690.9520.9660.9490.9560.9230.953n = 70ABias0.0920.1000.0940.0940.0860.0920.1050.062MSE0.0510.0640.0620.0620.0500.0600.0630.054CW (CL)0.8120.9100.9020.9020.7990.8760.8780.862CP0.9130.9420.9500.9710.9270.9490.9210.966α=3,β1=3, β2=3,β3=3n = 30ABias0.1570.3770.3180.3740.1170.3610.2930.249MSE0.1410.8370.7950.8330.1310.7960.7160.658CW (CL)1.3363.2683.1803.2641.3283.1503.0252.975CP0.9360.9660.9740.9690.9440.9320.9610.956n = 50ABias0.1210.2420.2830.2580.0960.2140.2810.174MSE0.0830.4200.4560.4340.0790.3970.4510.367CW (CL)1.0222.3562.4042.3741.0242.3062.3452.248CP0.9410.9750.9510.9630.9390.9620.9290.958n = 70ABias0.0740.1650.1590.1880.0580.1500.1580.137MSE0.0520.2690.2660.2810.0510.2590.2640.249CW (CL)0.8501.9281.9241.9440.8551.8931.8921.876CP0.9370.9510.9610.9630.9430.9380.9480.956α=3,β1=2.8,β2=2.8,β3=3.5n = 30ABias0.1540.3390.2830.4550.1030.2840.2430.257MSE0.1390.7000.6331.2210.1260.6210.5820.896CW (CL)1.3302.9982.9163.9481.3142.8342.7553.513CP0.9370.9640.9770.9700.8640.9380.8470.914n = 50ABias0.1210.2180.2520.3240.0870.1690.2330.210MSE0.0830.3540.3810.6420.0760.3200.3620.529CW (CL)1.0222.1682.2082.8721.0182.0982.1362.697CP0.9470.9730.9550.9650.8050.9740.7970.857n = 70ABias0.0810.1560.1480.2380.0530.1210.1330.137MSE0.0540.2300.2270.4140.0500.2110.2170.344CW (CL)0.8501.7801.7742.3420.8491.7261.7332.220CP0.9330.9510.9650.9630.8000.9470.7060.839α=3,β1=3,β2=4,β3=3n = 30ABias0.1540.3740.4770.3700.1090.3270.4550.230MSE0.1360.8271.5580.8220.1200.7261.4410.617CW (CL)1.3123.2504.5203.2441.2843.0484.3072.914CP0.9370.9630.9730.9710.9410.9510.9700.963n = 50ABias0.1160.2320.4060.2560.0910.2060.3880.171MSE0.0790.4090.9150.4300.0730.3800.8510.354CW (CL)1.0082.3363.3942.3660.9892.2543.2542.213CP0.9420.9730.9600.9650.9420.9530.9650.952n = 70ABias0.0740.1630.2340.1820.0580.1560.2280.123MSE0.0510.2660.5300.2760.0490.2570.5130.241CW (CL)0.8401.9162.7021.9320.8341.8692.6361.847CP0.9430.9560.9600.9600.9550.9650.9670.958α=3.5,β1=3.5,β2=3.5,β3=4n = 30ABias0.1710.4500.3780.5360.1000.4020.3890.238MSE0.1761.1901.0741.6660.1571.0951.0711.149CW (CL)1.5023.8943.7824.6041.4863.7023.6604.027CP0.9390.9630.9750.9690.9490.9520.9310.937n = 50ABias0.1320.2810.3300.3710.0800.2530.2680.166MSE0.1040.5870.6390.8590.0930.5620.5730.653CW (CL)1.1562.7942.8543.3281.1462.7392.7153.059CP0.9440.9720.9590.9660.9360.9610.9000.903n = 70ABias0.0850.1970.1850.2610.0550.1450.1940.159MSE0.0680.3810.3720.5450.0650.3510.3780.465CW (CL)0.9682.2922.2802.7080.9672.2302.2582.583CP0.9400.9520.9630.9700.9500.9590.8710.864α=3.5,β1=3.5, β2=3.5,β3=3.5n = 30ABias0.1720.4500.3780.4450.1070.3690.3230.250MSE0.1791.1941.0791.1880.1621.0560.9920.885CW (CL)1.5143.9023.7923.9001.5033.6893.5933.488CP0.9410.9660.9750.9710.9550.9450.9590.955n = 50ABias0.1330.2840.3310.3080.1000.2630.3080.199MSE0.1060.5920.6430.6170.0990.5730.6100.516CW (CL)1.1622.3832.8622.8321.1602.7552.7552.665CP0.9460.9710.9590.9650.9460.9640.9290.949n = 70ABias0.0870.2040.1860.2240.0590.1630.1820.144MSE0.0690.3870.3750.3980.0660.3620.3740.346CW (CL)0.9682.3022.2862.3120.9782.2522.2592.220CP0.9370.9540.9630.9630.9410.9510.9560.956α=4,β1=1.5, β2=2,β1=2.5n = 30ABias0.2340.1520.1990.3300.1770.1430.1850.236MSE0.3510.1620.3020.5900.2720.1550.2840.474CW (CL)1.9281.4622.0062.7221.9061.4151.9122.486CP0.9320.9630.9740.9710.9500.9410.9450.954n = 50ABias0.1780.1000.1670.2310.1410.0920.1550.167MSE0.1730.0850.1770.3060.1620.0830.1690.260CW (CL)1.4721.0741.5141.9721.4681.0591.4681.871CP0.9440.9720.9500.9630.9600.9580.9530.929n = 70ABias0.1360.0820.1200.1930.1210.0730.1370.126MSE0.1170.0590.1130.2090.1160.0570.1200.177CW (CL)1.2340.8941.2341.6261.2410.8831.2361.559CP0.9290.9500.9620.9710.9430.9390.9460.935(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. The algorithm is explained in detail in Appendix B.2.From Table9, it can be seen that for different combinations of the parameters and for different sample sizes, ABias and MSE are relatively small. This 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 Figure4. Also, for the majority of cases, Bayesian estimates have smaller ABias than ML ones.Figure 4
MSE of (a)α=3.5, (b) β1=3.5, (c) β2=3.5, and (d) β3=3.5.
(a)(b)(c)(d)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.Figure 5
(a) MSE and (b) ABias forα.
(a)(b)


## 6.2.2. A Real Dataset

Here, we analyze a real dataset to show the applicability of MIK. The 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. The variables are as follows:X1: first systolic blood pressure measurement.X2: second systolic blood pressure measurement.X3: third systolic blood pressure measurement.All data points have been divided by 1000. This 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 statisticW∗ and Anderson–Darling statistic A∗. The values of these statistics and the corresponding p values (in brackets) for X1, X2, X3, and maxX1,X2,X3 are illustrated in Table 10.Table 10
W∗ and A∗ for IK, IW, and E.
ModelX1X2X3maxX1,X2IKA∗0.31 (0.93)0.43 (0.82)0.49 (0.76)0.50 (0.75)W∗0.04 (0.94)0.07 (0.76)0.08 (0.68)0.08 (0.72)IWA∗13.60 (7 ×10−6)14.48 (7 ×10−6)15.80 (7 ×10−6)14.20 (7 ×10−6)W∗2.79 (2 ×10−7)2.99 (5 ×10−8)3.30 (9 ×10−9)2.91 (7 ×10−8)EA∗23.95 (7 ×10−6)24.18 (7 ×10−6)24.75 (7 ×10−6)24.29 (7 ×10−6)W∗5.06 (2 ×10−16)5.12 (2 ×10−16)5.26 (2 ×10−16)5.14 (2 ×10−16)Based on the values of these statistics and the correspondingp 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.Table 11
The estimates based on ML and Bayesian approaches.
Approachαβ1β2β3MLMLE46.949264.199245.852245.519Confidence interval width8.126269.684256.424263.866BayesianPosterior mean48.697270.528226.877257.231Credible interval length4.27069.21561.11064.305


## 7. Conclusion

In this paper, bivariate inverted Kumaraswamy (BIK) distribution is derived as a new member of bivariate Marshall–Olkin family. Its properties are also studied. Estimation is performed using both maximum likelihood (ML) and Bayesian approaches. To see the applicability of BIK distribution, three real datasets in different fields like sports, engineering, and biology have been analyzed. It is observed that the BIK model provides the best fit. Due to the wide applicability and great performance of the BIK model, a generalization to multivariate inverted Kumaraswamy (MIK) distribution is performed. To the best of our knowledge, estimation of multivariate Marshall–Olkin family was not studied before. Here, estimation of MIK is performed using both ML and Bayesian approaches. A very convenient algorithm is proposed for both approaches. This algorithm could be applied for any multivariate Marshall–Olkin distributions. Finally, a real dataset has been analyzed to illustrate the applicability of MIK distribution, and it is observed that the MIK model provides good fit, while multivariate inverse Weibull and multivariate exponential distributions failed to fit this dataset.



---

*Source: 6349523-2020-10-23.xml*