In this article, we estimate the parameters of the bivariate Inverse Lomax distribution of Marshall-Olkin based on right censored sample. Utilizing EM algorithm is a priority because the vector of the observed data is not complete but viewed as an observable function of complete data. After that, the EM algorithm makes use of the simplicity of maximum likelihood estimation for complete data. In addition, normal deviations of the estimates of bivariate Lomax distribution are derived.

A comparison is conducted via a simulation study between estimates obtained by using the EM algorithm and without the EM algorithm.

Keywords: bivariate lomax distribution, censored sample, EM algorithm, missing data, simulation study

Inverse Lomax distribution that is one of the inverted distributions which is very flexible to analyze the situation with non-monotonic failure rate, Singh, et al.¹ If a random variable Y has Lomax distribution, then x=1/y has an Inverse Lomax distribution (ILD) Singh, et al.² The Inverse Lomax Distribution (ILD) is used in random modeling of life components that have a decreasing failure rate. Like other distributions included in the family of generalized Beta distribution, Inverse Lomax Distribution also has application in actuarial sciences and economy and Kleiber and Kotz.³

 Inverse Lomax was implemented on geophysical databases McKenzie, et al. (2011)4 These databases were about sizes of land fires in the state of California Rahman et al.,⁵ carried out research about statistical inference and Prediction on inverse Lomax distribution via Bayesian inferences. Kleiber⁶ tackled Inverse Lomax distribution to possess the Lorenz ordering relationship between ordered statistics.

In this article, maximum likelihood estimates (MLEs) of the parameters of the bivariate Inverse Lomax distribution (MOBIL) of Marshall-Olkin⁷ are obtained based on censored samples. The censoring time (T) is supposed to be independent of the life times X,Y of the two components. For example, this situation happens in medical studies of organs with pairs like kidneys, eyes, lungs, or any other paired organs of an individual like a two-component system that is interdependent. Failure of an individ ual might censor failure of either one of the paired organ or both.

This scheme of censoring is right censoring, There are similar situation in engineering science whenever sub-systems are considered having two components with life times (X,Y) being independent of the life time (T) of the entire system. However, failure of the main system may censor failure of either one component or both. Maximum likelihood estimators of the parameters for the case of univariate right censoring were derived by Hanagal (a,b).^8,9 In addition, censoring might happen in different ways. Patients may not follow up during the study. Some patients might decide to move somewhere else. Thereupon, the experimenter may not follow them again, or the patients may not continue to cooperate because of bad side effects of the treatment. These cases are called withdrawal from the study. Useful information is presented by patients with censored data. Therefore, they should not be omitted from the analysis. Due to lack of data about real processes, the data in this study are derived from the BVL of Marshall-Olkin⁷ with Matlab software. Subsequently, the parameters are assessed using EM algorithm.

Let X be a random variable with the following CDF as follows;

$F_{IL}\left(x\right)={(1+\frac{\theta }{x})}^{-\alpha } ,x>0 , \alpha ,\theta >0 \left(1\right)$

The distribution of this form is said to be a Inverse Lomax distribution with parameters α and θ, will be denoted by IL(α,θ). The PDF and the hazard function of Inverse Lomax distribution with parameter (α,θ) will be

$f_{IL}(x;\alpha ,k)=\frac{\theta \alpha }{x^2}{(1+\frac{\theta }{x})}^{-\alpha -1} ,x>0 \left(2\right)$

Suppose U₁,U₂ and U₃ are three independent random variables such that $U_i~ IL({\alpha }_i,\theta )$ for i = 1, 2 and 3 it is assumed that ${\alpha }_1, {\alpha }_2, {\alpha }_3, \theta >0$. Define $X=\max (u_1,u_3)$ and $Y=\max (u_2,u_3)$.

Then we say that the bivariate vector $(X,Y)$ has a bivariate Inverse Lomax distribution Marshall-Olkin.⁷

The joint cumulative function of X and Y is defined as

$F_{X,Y}(x,y)={(1+\frac{\theta }{x})}^{-{\alpha }_1}{(1+\frac{\theta }{y})}^{-{\alpha }_2}{(1+\frac{\theta }{\min (x, y)})}^{-{\alpha }_3}$$=F_{IL}(x;{\alpha }_1,\theta )F_{IL}(y;{\alpha }_2,\theta )F_{IL}(\min (x, y).;{\alpha }_3,\theta ) \left(3\right)$

$F_{X,Y}(x,y)=\{\begin{array}{c}{F}_1(x,y) 0<x<y<\infty \\ F_2(x,y) 0<y<x<\infty \\ F_0\left(x\right) 0<x=y<\infty \end{array} \left(4\right)$

where

$F_1(x,y)=F(y;{\alpha }_2+{\alpha }_3,\theta )F(x;{\alpha }_1,\theta )= {(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)}{(1+\frac{\theta }{y})}^{-{\alpha }_2}$

$F_2(x,y)=F(x;{\alpha }_1+{\alpha }_3,\theta )F(y;{\alpha }_2,\theta )={(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)}{(1+\frac{\theta }{x})}^{-{\alpha }_1}$

$F_0\left(x\right)=F(x;\alpha ,\theta )={(1+\frac{\theta }{x})}^{-\alpha }$

and $\alpha ={\alpha }_1+{\alpha }_2+{\alpha }_3$

The joint probability density function $f_{X,Y}(x,y)$ of X and Y takes the form

$f_{X,Y}(x,y)=\{\begin{array}{c}{f}_1(x,y) if x<y\\ f_2(x,y) if y<x\\ f_0\left(x\right) if x=y\end{array} \left(5\right)$

where

$f_1(x,y)=\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 y^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{y})}^{-{\alpha }_2-1}$

$f_2(x,y)=\frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{x^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{x})}^{-{\alpha }_1-1}$

and

$f_0\left(x\right)=\frac{{\theta }^2{\alpha }_3}{x^2}{(1+\frac{\theta }{x})}^{-\alpha -1}$

The density functions of $X|\left\{(x,y)|x>y\right\}, Y|\left\{(x,y)|y>x\right\}$ and Z=min (X,Y) are given as follows:

$f_{X|\left\{(x,y)|x>y\right\}}\left(x\right)=\frac{\theta ({\alpha }_1+{\alpha }_3)}{x^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{y})}^{{\alpha }_1+{\alpha }_3} , x>y$

$f_{Y|\left\{(x,y)|y>x\right\}}\left(y\right)=\frac{\theta ({\alpha }_2+{\alpha }_3)}{y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{x})}^{{\alpha }_2+{\alpha }_3} , y<x$

$f_Z\left(z\right)=\frac{\theta \alpha }{z^2}{(1+\frac{\theta }{z})}^{-\left(\alpha \right)-1} , Z=\min (x,y)$

This article aims at deriving an estimation method for the parameters of a bivariate Inverse Lomax distribution of Marshall-Olkin by the EM algorithm. In Section 2, the EM algorithm is presented, while in Section 3, the parameters are estimated by applying the EM algorithm. Finally in Section 4, we present the results of a simulation study.

The expectation-maximization (EM) algorithm was introduced by Dempster et al.¹⁰ The algorithm is considered a repeated procedure used to find the maximum likelihood estimates of data that are missing, incomplete, unobserved or censored data. Because it is easy to implement, the impact of the EM algorithm has had great effect, not only as a tool for computation but also as a method of solving complicated statistical problems.

 The basic idea behind the method is to transform a set of incomplete data into a complete data problem for which the required maximization is computationally more tractable and stable numerically. Each repetition raises the likelihood, which is finally converged to a local maximum. The complete data set x can be viewed as consisting of vectors $(t,t^{*})$, where t is the observed incomplete data, and $t^{\text{*}}$ is the missing data.

The EM Algorithm method has been applied by several authors, for example, see Qin, et al.,¹¹ Rudolf,¹² Ning, et al.,¹³ Zanini, et al.,¹⁴ Açıkgöz,¹⁵ Kalabatsos,¹⁶ Attia, et al.¹⁷ and Hanagal and Ahmadi.¹⁸

The iterations

The objective is to draw inferences about the parameter vectors We will useto denote the likelihood function where t is the vector of observed data. Let $t^{\text{*}}$ represent the vector of missing data. Starting with a guessed value for the parameter $\underset{\_}{\alpha }$, carry out the following iterations:

• Replace the missing data $t^{\text{*}}$ by their expectation given the guessed value of the parameter vector and the observed data. Let this conditional expectation be denoted by $\widetilde{t^{\text{*}}}$.
• Maximize $L(\underset{\_}{\alpha },t^{\text{*}}\text{|}t)$ with respect to $\underset{\_}{\alpha }$ replacing the missing data $t^{\text{*}}$ by their expected values. This is equivalent to maximizing $L\left(\underset{\_}{\alpha },E\left[T^{\text{*}}\text{|}t\right]\right)$.
• Re-estimate the missing values $t^{\text{*}}$ using their conditional expectation based on the updated $\underset{\_}{\alpha }$.
• Re-estimate $\underset{\_}{\alpha }$ and continue until the difference between a new iterated value and the previous iterated value is less than 0.00001.

Execution of the algorithm

• Consider the logarithm of the likelihood log $L(\underset{\_}{\alpha },t^{\text{*}}\text{|}t)$.
• Let the conditional expectation of log $L(\underset{\_}{\alpha },t^{\text{*}}\text{|}t)$ with respect to $t^{\text{*}}|({\underset{\_}{\alpha }}^{\left(k\right)},t)$ at k^th step be noted by $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$, where ${\underset{\_}{\alpha }}^{\left(k\right)}$ is the current guess of $\underset{\_}{\alpha }$.Then the EM-steps are as follows:

1. E-step: Calculate $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$, that is, that is, the expectation of the log-likelihood with respect to the conditional distribution of the missing data, given the observed data and the current guess of $\underset{\_}{\alpha }$.
2. M-step: Maximize $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$ with respect to $\underset{\_}{\alpha }$ and set the result equal to ${\underset{\_}{\alpha }}^{(\text{k}+1)}$, the new value of the parameter vector. Since ${\underset{\_}{\alpha }}^{(\text{k}+1)}$ maximizes $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$ the M- Step results in $Q\left({\underset{\_}{\alpha }}^{(k+1)}|{\underset{\_}{\alpha }}^{\left(k\right)}\right)\ge Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$ for all $\underset{\_}{\alpha }\in \Omega$, implying that ${\underset{\_}{\alpha }}^{(\text{k}+1)}$ is a solution to equation,$\frac{\partial Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial \underset{\_}{\alpha }}=0.$

The two steps are repeated iteratively until the difference between two successive iterations is less than 0.00001. This iterative procedure leads to a monotonic increase of log $L\left(\underset{\_}{\alpha } ,E \left[T^{*}\right|t]\right)$:

$\log L\left({\underset{\_}{\alpha }}^{(k+1)} ,E \left[T^{*}\right|t]\right)\ge \log L\left(\underset{\_}{\alpha } ,E \left[T^{*}\right|t]\right) for k = 1,2,\mathrm{...} \left(8\right)$

Since the likelihood increases in each step, the EM algorithm converges generally to a local maximum. When there is no closed form solution of the M-step, a numerical algorithm as, for example, the Newton-Raphson procedure, may be used for iteratively computing ${\underset{\_}{\alpha }}^{\left(k\right)}$. In fact, in this paper the Newton-Raphson procedure is used to obtain maximum likelihood estimates of $\underset{\_}{\alpha }$ at the (k +1)^th iteration as follows:

${\underset{\_}{\alpha }}^{(k+1)}={\underset{\_}{\alpha }}^{\left(k\right)}-{{\left(\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial \underset{\_}{\alpha }\partial \acute{\underset{\_}{\alpha }}}\right)}^{-1}}_{\underset{\_}{\alpha }={\underset{\_}{\alpha }}^{\left(k\right)}}{\left(\frac{\partial Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial \underset{\_}{\alpha }}\right)}_{\underset{\_}{\alpha }={\underset{\_}{\alpha }}^{\left(k\right)}} \left(9\right)$

The iterative procedure is carried out until the difference ${\underset{\_}{\alpha }}^{(k+1)}-{\underset{\_}{\alpha }}^{\left(k\right)}<0.0001$, For more detail see Hanagal and Ahmadi.¹⁸

 The density function of (X,Y) is given by

$f_{X,Y}(x,y)=\{\begin{array}{c}\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 y^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{y})}^{-{\alpha }_2-1} for 0<x<y<\infty \\ \frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{x^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{x})}^{-{\alpha }_1-1} for 0<y<x<\infty \\ \frac{{\theta }^2{\alpha }_3}{x^2}{(1+\frac{\theta }{x})}^{-\alpha -1} for x=y\end{array}\left(10\right)$

For the bivariate life time distribution, we use the univariate censoring plan presented by Hanagal(a,b)^8,9 as the persons do not join the study at the same time and withdrawal or death of a person or ending the study will censor life times of both components. The time of censoring is independent of life times of both components, This is the standard univariate right.

Suppose that there are n independent pairs of components, for example, paired kidneys, lungs, eyes, ears in an individual under study and i^th pair of the components have life times $(X_i,Y_i)$ and censoring time (T_i ). The life times associated with i-th pair of the components are given by

$(X_i,Y_i)=\{\begin{array}{c}(X_i,Y_i) if\max (X_i,Y_i)<T_i\\ (X_i,T_{iy}) if X_i<T_{iy}<Y_i\\ (T_{ix},Y_i) if Y_i<T_{ix}<X_i\\ (T_{ixy},T_{ixy}) if\min (X_i,Y_i)>T_{ixy}\end{array} \left(11\right)$

for i = 1,2,...,n

where T_ix,T_iy and T_ixy represent the unobserved random variables x<y, x>y, x=y respectively.

The likelihood of the sample of size n after discarding factors which do not contain any of the parameters of interest is given as follows

$L(\underset{\_}{\alpha },t_{ix},t_{iy},t_{ixy})=\prod _{j=i}^6\prod _{i=1}^{n_j}{f}_{X,Y}^j(x_i,y_i) \left(12\right)$

where

$\begin{array}{c}\begin{array}{c}{f^1}_{X,Y}(x,y)=\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 y^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{y})}^{-{\alpha }_2-1} , 0<x<y<t\\ {f^2}_{X,Y}(x,y)=\frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{x^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{x})}^{-{\alpha }_1-1} , 0<x<y<t\\ {f^3}_{X,Y}(x,y)=\frac{{\theta }^2{\alpha }_3}{x^2}{(1+\frac{\theta }{x})}^{-\alpha -1} , 0< x=y<t\end{array}\\ {f^4}_{X,Y}(x,y)=\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 {t_y}^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{t_y})}^{-{\alpha }_2-1} ,0<x<t<t_y=y \\ \begin{array}{c}{f^5}_{X,Y}(x,y)=\frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{{t_x}^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{t_x})}^{-{\alpha }_1-1} , 0<y<t<t_x=x \\ {f^6}_{X,Y}(x,y)=\frac{{\theta }^2{\alpha }_3}{{t_{xy}}^2}{(1+\frac{\theta }{t_{xy}})}^{-\alpha -1} , 0<t<t_{xy}=\min (x,y)\end{array}\end{array}\}\left(13\right)$

where $\sum _{j=1}^6{n}_j=n$,$n_1, n_2, n_3, n_4, n_5 and n_6$ be the numbers of realizations falling in the range corresponding to ${f^j}_{X,Y}(x,y)$ and j=1,2,...,6, respectively.${f^1}_{X,Y}(x,y)$, ${f^2}_{X,Y}(x,y)$, ${f^4}_{X,Y}(x,y)$ and ${f^5}_{X,Y}(x,y)$ are density functions with respect to the Lebesque measure on R2, while , ${f^3}_{X,Y}(x,y)$,and ${f^6}_{X,Y}(x,y)$ are density functions with respect to the Lebesque measure on R.

Let the range of variability corresponding to ${f^j}_{X,Y}(x,y)$ be denoted by $A_j$, j=1, 2, …., 6.and $A_1,\dots ,A_6$ are disjoint sets and letting $E_1=A_1\cup {A_4}^{}$,$E_2=A_2\cup A_5,E_3=A_3\cup A_6,E_4=E_1\cup E_2$,$E_5=E_1\cup A_2\cup A_3$ and $E_6=E_2\cup A_1$, the log-likelihood function $\ln L_c(\underset{\_}{\alpha },t_{ix},t_{iy},t_{ixy})$ can be written as follows:

$\ln L=(n_1+n_4)(\ln ({\alpha }_1+{\alpha }_3)+\ln {\alpha }_2)+(n_2+n_5)(\ln ({\alpha }_2+{\alpha }_3)+\ln {\alpha }_1)+(n_3+n_6)\ln {\alpha }_3$

$-2\sum _{i\in E_5}\ln \left(x_i\right)-2\sum _{i\in E_6}\ln \left(y_i\right)-({\alpha }_1+{\alpha }_3+1)\sum _{i\in E_1}\ln (1+\frac{\theta }{x_i})-(\alpha +1)\sum _{i\in A_3}\ln (1+\frac{\theta }{x_i})$

$-({\alpha }_1+1)\sum _{i\in A_2}\ln (1+\frac{\theta }{x_i})-({\alpha }_2+{\alpha }_3+1)\sum _{i\in E_2}\ln (1+\frac{\theta }{y_i})-({\alpha }_2+1)\sum _{i\in A_6}\ln (1+\frac{\theta }{y_i})$

$-2\sum _{i\in A_4}\ln \left(t_{yi}\right)-2\sum _{i\in A_5}\ln \left(t_{xi}\right) -({\alpha }_1+1)\sum _{i\in A_5}\ln (1+\frac{\theta }{t_{xi}})-({\alpha }_2+1)\sum _{i\in A_4}\ln (1+\frac{\theta }{t_{yi}})$

$-2\sum _{i\in A_6}\ln \left(t_{xyi}\right)-(\alpha +1)\sum _{i\in A_6}\ln (1+\frac{\theta }{t_{xyi}})+n\ln \left(\theta \right)+2(n_1+n_4+n_2+n_5)\ln \left(\theta \right) \left(14\right)$

Estimation parameter with EM algorithm

The E-step and M- step are obtained as follow:

E- step:

The unobserved random variables, T_ix,T_iy and T_ixy follow the distributions as stated in (7). The conditional distributions of $T_x|\left\{(t_x,t,y)|(t_x>t>y)\right\}$, $T_y|\left\{(t_y,t,x)|(t_y>t>x)\right\}$ and $T_{xy}|\left\{(t_{xy},t)|(t_{xy}>t)\right\}$ are given as follows:

$\begin{array}{c}{f}_{T_x|(t_x,t,y)}\left(t_x\right)=\frac{\theta ({\alpha }_1+{\alpha }_3)}{{t_x}^2}{(1+\frac{\theta }{t_x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{t})}^{{\alpha }_1+{\alpha }_3} , y<t<t_x \\ f_{T_y|(t_y,t,x)}\left(t_y\right)=\frac{\theta ({\alpha }_2+{\alpha }_3)}{{t_y}^2}{(1+\frac{\theta }{t_y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{t})}^{{\alpha }_2+{\alpha }_3}, x<t<t_y\\ f_{T_{xy}|(t_{xy},t)}\left(t_{xy}\right)=\frac{\theta \alpha }{{t_{xy}}^2}{(1+\frac{\theta }{t_{xy}})}^{-\left(\alpha \right)-1}{(1+\frac{\theta }{t})}^{\left(\alpha \right)} ,t>t_{xy}\end{array}\}\left(15\right)$

The values of the first moments of the conditional unobserved random variables $(1+\frac{T_x}{k})|\left\{(t_x,t,y)|(t_x>t>y)\right\} ,(1+\frac{T_y}{k})|\left\{(t_y,t,x)|(t_y>t>x)\right\} and$
 (1+Txy/k)|{(txy,t)|(txy>t)} are as follows.

$\begin{array}{c}E(1+T_x/k)|\left\{(t_x,t,y)|(t_x>t>y)\right\}=\frac{-({\alpha }_1+{\alpha }_3)}{{\alpha }_1+{\alpha }_3-1} (1+\frac{\theta }{t})={\gamma }_1(1+\frac{\theta }{t}) \left(say\right)\\ E(1+T_y/k)|\left\{(t_y,t,x)|(t_y>t>x)\right\}=\frac{-({\alpha }_2+{\alpha }_3)}{{\alpha }_2+{\alpha }_3-1}(1+\frac{\theta }{t})={\gamma }_2(1+\frac{\theta }{t}) \left(say\right)\\ E(1+T_{xy}/k)|\left\{(t_{xy},t)|(t_{xy}>t)\right\}=\frac{-\alpha }{\alpha -1}(1+\frac{t}{k}) ={\gamma }_3(1+\frac{\theta }{t}) \left(say\right)\end{array}\}\left(16\right)$

The conditional expectation $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$ is obtained as follows:

$Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)=(n_1+n_4)(\ln ({\alpha }_1+{\alpha }_3)+\ln {\alpha }_2)+(n_2+n_5)(\ln ({\alpha }_2+{\alpha }_3)+\ln {\alpha }_1)+(n_3+n_6)\ln {\alpha }_3$ $-2\sum _{i\in E_5}\ln \left(x_i\right)-2\sum _{i\in E_6}\ln \left(y_i\right)-({\alpha }_1+{\alpha }_3+1)\sum _{i\in E_1}{\dot{x}}_i-(\alpha +1)\sum _{i\in A_3}{\dot{x}}_i$ $-({\alpha }_1+1)\sum _{i\in A_2}{\dot{x}}_i-({\alpha }_2+{\alpha }_3+1)\sum _{i\in E_2}{\dot{y}}_i-({\alpha }_2+1)\sum _{i\in A_6}{\dot{y}}_i$ $+2\sum _{i\in A_4}\ln (\frac{{\gamma }_2^{\left(k\right)}-1}{\theta }+\frac{{\gamma }_2^{\left(k\right)}}{t_i})+2\sum _{i\in A_5}\ln (\frac{{\gamma }_1^{\left(k\right)}-1}{\theta }+\frac{{\gamma }_1^{\left(k\right)}}{t_i}) -({\alpha }_1+1)\sum _{i\in A_5}({\gamma }_1^{\left(k\right)}+{\dot{t}}_i)$ $-({\alpha }_2+1)\sum _{i\in A_4}({\gamma }_2^{\left(k\right)}+{\dot{t}}_i) +2\sum _{i\in A_6}\ln (\frac{{\gamma }_3^{\left(k\right)}-1}{\theta }+\frac{{\gamma }_3^{\left(k\right)}}{t_i})-(\alpha +1)\sum _{i\in A_6}({\gamma }_3^{\left(k\right)}+{\dot{t}}_i)$ $+ n\ln \left(\theta \right)+2(n_1+n_4+n_2+n_5)\ln \left(\theta \right)$ (17)

suppose is θ known and ${{\gamma }_1}^{\left(k\right)}, {{\gamma }_2}^{\left(k\right)}, {{\gamma }_3}^{\left(k\right)}$ are expressed in terms of ${{\alpha }_1}^{\left(k\right)}, {{\alpha }_2}^{\left(k\right)}, {{\alpha }_3}^{\left(k\right)}$ at K^th iteration, and

${\dot{x}}_i=\ln (1+\frac{\theta }{x_i}), {\dot{y}}_i=\ln (1+\frac{\theta }{y_i}) and {\dot{t}}_i=\ln (1+\frac{\theta }{t_i}).$

M-step:

The following equations are obtained by equating the partial derivatives of $Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)$ with respect to ${\alpha }_1,{\alpha }_{2}and {\alpha }_3$ to zero:

$\frac{n_2+n_5}{{\alpha }_1}+\frac{n_1+n_4}{{\alpha }_1+{\alpha }_3}-\sum _{i\in D_1}{\dot{x}}_i-\sum _{i\in A_5}({\gamma }_1^{\left(k\right)}+{\dot{t}}_i)-\sum _{i\in A_6}({\gamma }_3^{\left(k\right)}+{\dot{t}}_i)=0 \left(18\right)$

$\frac{n_1+n_4}{{\alpha }_2}+\frac{n_2+n_5}{{\alpha }_2+{\alpha }_3}-\sum _{i\in A_3}{\dot{x}}_i-\sum _{i\in D_2}{\dot{y}}_i-\sum _{i\in A_4}({\gamma }_2^{\left(k\right)}+{\dot{t}}_i)-\sum _{i\in A_6}({\gamma }_3^{\left(k\right)}+{\dot{t}}_i)=0 \left(19\right)$

$\frac{n_1+n_4}{{\alpha }_1+{\alpha }_3}+\frac{n_2+n_5}{{\alpha }_2+{\alpha }_3}+\frac{n_3+n_6}{{\alpha }_3}-\sum _{i\in A_3}{\dot{x}}_i-\sum _{i\in E_1}{\dot{x}}_i-\sum _{i\in E_2}{\dot{y}}_i-\sum _{i\in A_6}({\gamma }_3^{\left(k\right)}+{\dot{t}}_i)=0 \left(20\right)$

Where $D_1=E_1\cup A_2\cup A_3$ and $D_2=E_2\cup A_6.$

The above likelihood equations are solved for the maximum likelihood estimates $({{\hat{\alpha }}_1}^{\left(k\right)}, {{\hat{\alpha }}_2}^{\left(k\right)}, {{\hat{\alpha }}_3}^{\left(k\right)})$ using the Newton-Raphson procedure. Below the observed elements of the symmetric information matrix at the current guess $\underset{\_}{\alpha }={\underset{\_}{\alpha }}^{\left(k\right)}$ are given which is required in the Newton-Raphson iterative procedure.

$\begin{array}{l}-\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_1^2}=\frac{(n_2+n_5)}{{{\alpha }_1}^2}+\frac{(n_1+n_4)}{{({\alpha }_1+{\alpha }_3)}^2}\\ -\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_1\partial {\alpha }_2}=0\\ -\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_1\partial {\alpha }_3}=\frac{(n_1+n_4)}{{({\alpha }_1+{\alpha }_3)}^2}\\ -\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_2^2}=\frac{(n_1+n_4)}{{\alpha }_2^2}+\frac{(n_2+n_5)}{{({\alpha }_2+{\alpha }_3)}^2}\\ -\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_2\partial {\alpha }_3}=\frac{(n_2+n_5)}{{({\alpha }_2+{\alpha }_3)}^2}\\ -\frac{{\partial }^{2}Q\left(\underset{\_}{\alpha }|{\underset{\_}{\alpha }}^{\left(k\right)}\right)}{\partial {\alpha }_3^2}=\frac{(n_1+n_4)}{{({\alpha }_1+{\alpha }_3)}^2}+\frac{(n_2+n_5)}{{({\alpha }_2+{\alpha }_3)}^2}+\frac{(n_3+n_6)}{{\alpha }_3^2}\end{array}\}\left(21\right)$

As mentioned above the iterative process is carried out until the following condition is met by two subsequent solutions:

$L({\underset{\_}{\alpha }}^{(k+1)},E[t_{ix},t_{iy},t_{ixy}|x_i,y_i,t_i])-L({\underset{\_}{\alpha }}^{\left(k\right)},E[t_{ix},t_{iy},t_{ixy}|x_i,y_i,t_i])<0.0001 \left(22\right)$

Estimation parameter without EM algorithm

The likelihood of the sample of size n after discarding factors which do not contain any of the parameters of interest is given as follows

$L(\underset{\_}{\alpha },t)=\prod _{j=i}^6\prod _{i=1}^{n_j}{f}_{X,Y}^j(x_i,y_i) \left(23\right)$

where

$\begin{array}{l}{f^1}_{X,Y}(x,y)=\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 y^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{y})}^{-{\alpha }_2-1} , 0<x<y<t\\ {f^2}_{X,Y}(x,y)=\frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{x^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{x})}^{-{\alpha }_1-1} , 0<x<y<t\\ {f^3}_{X,Y}(x,y)=\frac{{\theta }^2{\alpha }_3}{x^2}{(1+\frac{\theta }{x})}^{-\alpha -1} , 0< x=y<t\\ {f^4}_{X,Y}(x,y)=\frac{({\alpha }_1+{\alpha }_3){\alpha }_2 {\theta }^2}{x^2 t^2}{(1+\frac{\theta }{x})}^{-({\alpha }_1+{\alpha }_3)-1}{(1+\frac{\theta }{t})}^{-{\alpha }_2-1} ,0<x<t<y\\ {f^5}_{X,Y}(x,y)=\frac{({\alpha }_2+{\alpha }_3){\alpha }_1{\theta }^2}{t^2 y^2}{(1+\frac{\theta }{y})}^{-({\alpha }_2+{\alpha }_3)-1}{(1+\frac{\theta }{t})}^{-{\alpha }_1-1} , 0<y<t<x\\ {f^6}_{X,Y}(x,y)=\frac{{\theta }^2{\alpha }_3}{t^2}{(1+\frac{\theta }{t})}^{-\alpha -1} , 0<t<\min (x,y)\end{array}\}\left(24\right)$

The log-likelihood function $\ln L_c(\underset{\_}{\alpha },t)$ can be written as follows:

$\ln L=(n_1+n_4)(\ln ({\alpha }_1+{\alpha }_3)+\ln {\alpha }_2)+(n_2+n_5)(\ln ({\alpha }_2+{\alpha }_3)+\ln {\alpha }_1)+(n_3+n_6)\ln {\alpha }_3$ $-2\sum _{i\in E_5}\ln \left(x_i\right)-2\sum _{i\in E_6}\ln \left(y_i\right)-({\alpha }_1+{\alpha }_3+1)\sum _{i\in E_1}{\dot{x}}_i-(\alpha +1)\sum _{i\in A_3}{\dot{x}}_i$ $-({\alpha }_1+1)\sum _{i\in A_2}{\dot{x}}_i-({\alpha }_2+{\alpha }_3+1)\sum _{i\in E_2}{\dot{y}}_i-({\alpha }_2+1)\sum _{i\in A_6}{\dot{y}}_i$ $-2\sum _{i\in A_4}\ln \left(t_i\right)-2\sum _{i\in A_5}\ln \left(t_i\right) -({\alpha }_1+1)\sum _{i\in A_5}{\dot{t}}_i-({\alpha }_2+1)\sum _{i\in A_4}{\dot{t}}_i$ $-2\sum _{i\in A_6}\ln \left(t_i\right)-(\alpha +1)\sum _{i\in A_6}{\dot{t}}_i+n\ln \left(\theta \right)+2(n_1+n_4+n_2+n_5)\ln \left(\theta \right) \left(25\right)$

The following likelihood equations are obtained by equating the partial derivatives of $\ln \left(L\right)$ with respect to ${\alpha }_1,{\alpha }_{2}and {\alpha }_3$ to zero:

$\frac{n_2+n_5}{{\alpha }_1}+\frac{n_1+n_4}{{\alpha }_1+{\alpha }_3}-\sum _{i\in D_1}{\dot{x}}_i-\sum _{i\in A_5\cup A_6}{\dot{t}}_i=0 \left(26\right)$

$\frac{n_1+n_4}{{\alpha }_2}+\frac{n_2+n_5}{{\alpha }_2+{\alpha }_3}-\sum _{i\in A_3}{\dot{x}}_i-\sum _{i\in D_2}{\dot{y}}_i-\sum _{i\in A_4\cup A_6}{\dot{t}}_i=0 \left(27\right)$

$\frac{n_1+n_4}{{\alpha }_1+{\alpha }_3}+\frac{n_2+n_5}{{\alpha }_2+{\alpha }_3}+\frac{n_3+n_6}{{\alpha }_3}-\sum _{i\in A_3}{\dot{x}}_i-\sum _{i\in E_1}{\dot{x}}_i-\sum _{i\in E_2}{\dot{y}}_i-\sum _{i\in A_6}{\dot{t}}_i=0 \left(28\right)$

The sample data are generated based on following algorithms:

Step 1: Generate ${\text{u}}_{\text{i}}$ using the Lomax distribution white parameters ${\text{α}}_1,{\text{α}}_2{\text{and α}}_3$.

Step 2: Take $X=\min (u_1,u_3)$ and $Y=\min (u_2,u_3)$ and, therefore (X,Y) follows a bivariate Inverse Lomax distribution of Marshall-Olkin type.

Step 3: Generate t_i using the Inverse Lomax distribution with $\text{β},\text{θ}$, where t_is are the censoring times.

For two cases with respect to the ${\alpha }_{i}s$, we generated 1000 sets of samples. Each set consisted of three samples with sizes n = 20, 35, 50 and 100. The corresponding maximum likelihood estimates are displayed in Table 1 together with the empirical standard deviation. The estimates denoted by $ML{E}_{em}$ and $S{E}_{em}$ are obtained by using the EM algorithm, while the estimates denoted by MLE and SE are obtained without the EM algorithm.

Table 1 and Table 2 shows that the estimates obtained by using the EM algorithm gave a smaller empirical standard error than those obtained without the EM algorithm.

Moreover, the estimates MLE_em are close to the true parameter values and the standard errors (SE_em) decrease as the sample size increases. The estimates for both methods are obtained by taking the mean of the 1000 maximum likelihood estimates and the mean of the 1000 standard errors from the 1000 samples of size n = 20, 35,50 and 100.

Also note that with a large sample size estimates of the values of both methods (with and without EM algorithm) from each other as well as the standard error, which means that when the sizes of large samples may give both methods convergent results and then fit whichever is applicable.

None.

None.

1. Singh SK, Singh U,  Kumar D. Bayes estimators of the reliability function and parameters of inverted exponential distribution using informative and non-informative priors. Journal of Statistical computation and simulation. 2012;83(12):2258–2269.
2. Singh A, Singh U, Singh SK. On hybrid censored inverse lomax distribution: application to the survival data. Statistica. 2016;76(2)185–203.
3. Kleiber C, Kotz S. Statistical size distributions in economics and actuarial sciences. New Jersey: John Wiley & Sons, Inc., Hoboken; 2003:1–340.
4. McKenzie D, Miller C, Falk DA. The Landscape ecology of fire. New York: Springer; 2011:1–312.
5. Rahman J, Aslam M, Ali S. Estimation and prediction of inverse lomax model via bayesian approach. Caspian Journal of Applied Sciences Research. 2013;2(3):43–56.
6. Kleiber C. Lorenz ordering of order statistics from log-logistic and related distributions. Journal of Statistical Planning and Inference. 2004;120(1-2):13–19.
7. Marshall AW, Olkin I. A multivariate exponential distribution. Journal of the American Statistical Association. 1967;62:30–44.
8. Hanagal DD. Some inference results in bivariate exponential distributions based on censored samples. Communications in Statistics-Theory & Methods. 1992a;21(5):1273–1295.
9. Hanagal DD. Some inference results in modified Freund’s bivariate exponential distribution. Biometrical Journal. 1992b;34:745–756.
10. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via de EM algorithm (with discussion). Journal of the Royal Statistical Society. B 1977;39(1):1–38.
11. Qin J, Ning J, Liu H, et al. Maximum likelihood estimations and EM algorithms with length-biased data. J Am Stat Assoc. 2011;106(496):1434–1449.
12. Rudolf KK. Robust estimation by expectation maximization algorithm. J Geod. 2013;87:107–116.
13. Ning J, Chen Y, Piao J. Maximum likelihood estimation and EM algorithm of Copas-like selection model for publication bias correction. Biostatistics. 2017;18(3):495–504.
14. Zanini P, Said S, Cavalcante CC, et al. Stochastic EM Algorithm for Mixture Estimation on Manifolds. IEEE. 2017.
15. İnci Açıkgöz. Parameter estimation with profile likelihood method and penalized EM algorithm in normal mixture distributions. Journal of Statistics and Management Systems. 2018;21(7):1211–1228.
16. Karabatsos G. Fast Search and Estimation of Bayesian Nonparametric Mixture Models Using a Classification Annealing EM Algorithm. 2019.
17. Attia AF, Shaban SA, Amer YM. Parameter Estimation for the Bivariate Lomax Distribution Based on Censored Samples. Applied Mathematical Sciences. 2014;8(35):1711–1721.
18. Hanagal DD, Ahmadi KA.  Parameter Estimation for the Bivariate Exponential Distribution by the EM Algorithm Based on Censored Samples. Economic Quality Control. 2008;23(2):257–266.