A Bivariate Power Lindley Survival Distribution

Abstract

We introduce and investigate the properties of new families of univariate and bivariate distributions based on the survival function of the Lindley distribution. The univariate distribution, to reflect the nature of its construction, is called a power Lindley survival distribution. The basic distributional properties of this model are described. Maximum likelihood estimates of the parameters of the distribution are studied and the corresponding information matrix is identified. A bivariate power Lindley survival distribution is introduced using the technique of conditional specification. The corresponding joint density and marginal and conditional densities are derived. The product moments of the distribution are obtained, together with bounds on the range of correlations that can be exhibited by the model. Estimation of the parameters of the model is implemented by maximizing the corresponding pseudo-likelihood function. The asymptotic variance–covariance matrix of these estimates is investigated. A simulation study is performed to illustrate the performance of these parameter estimates. Finally some examples of model fitting using real-world data sets are described.

1. Introduction

The Lindley (L) distribution, which was first introduced by Lindley [1], has been found to be useful in many research areas. In particular, it has been applied in survival analysis in many scientific arenas, including engineering, health, economics, etc. Ghitany et al. [2] studied the statistical properties of this distribution and showed that it was preferred to the exponential distribution in several applications. Several generalizations of the L model have been proposed by a variety of authors who have studied the properties of the extended models and have investigated their suitability in various application areas. Ghitany et al. [3] proposed a class of weighted L distributions. Ramos and Louzada [4] introduced a generalized weighted L model. Ristic and Balakrishnan [5] investigated the gamma–Lindley model. MirMostafaee et al. [6] studied the beta L distribution. Bakouch et al. [7] introduced an extended L distribution. Ghitany et al. [3] developed a two-parameter weighted L distribution. Nadarajah et al. [8] described a generalized L (GL) distribution. Ghitany et al. [9] introduced the power L distribution, Ashour and Eltehiwy [10] proposed the exponentiated power L distribution, Asgharzadeh et al. [11] introduced a Weibull Lindley distribution, Khokhar et al. [12] studied the Zografos Balakrishnan Power Lindley Distribution, Algarni [13] introduced a new generalized Lindley distribution, and Chhetri et al. [14] in the Cubic Rank Transmuted Lindley Distribution, among others.

This is a representative listing of variations on the Lindley theme, but it is not intended to be a complete list. The density function of the L distribution is of the form

(1)$f_L\left(x\right)={\displaystyle \frac{{\theta }^2}{\theta +1}}(1+x)e^{-\theta x},x>0,$

where the parameter $\theta$ is positive. If a random variable X has density (1) we write $X\sim L\left(\theta \right).$ The corresponding distribution function is

(2)$F_L\left(x\right)=1-{\displaystyle \frac{\theta +1+\theta x}{\theta +1}}{e}^{-\theta x},x>0,$

The survival function ($S_L\left(x\right)$) and hazard or failure rate function ($r_L\left(x\right)$) of this distribution are, respectively,

$S_L\left(x\right)={\displaystyle \frac{\theta +1+\theta x}{\theta +1}}{e}^{-\theta x}\mathrm{and}{r}_L\left(x\right)={\displaystyle \frac{{\theta }^2(x+1)}{\theta +1+\theta x}}.$

Lindley noted that $r_L\left(0\right)={\displaystyle \frac{{\theta }^2}{\theta +1}}$ and that $r_L\left(x\right)$ is an increasing function of x and of $\theta$. Moreover, ${\displaystyle \frac{{\theta }^2}{\theta +1}}<r\left(x\right)<\theta$. Since the failure function increases in X, the distribution is $IFR$ and the validity of the following chain of implication is well established. $IFR\Rightarrow IFRA\Rightarrow NBU\Rightarrow NBUE$ where IFR, IFRA, NBU, NBUE denote an increasing failure rate, an increasing failure rate average, new being better than used, and new being better than used in expectation, respectively (see Barlow and Proschan, [15]).

The first objective of this research project is to provide a flexible extension of the univariate L distribution that can adapt to a range of asymmetry and kurtosis features of data. The flexible model is based on the power of the survival function of the classical L model. The second objective is to extend discussion to the bivariate case, because there are few bivariate extensions of even the basic L model. The distributional properties of the proposed bivariate power L model are investigated and a method is proposed for estimation of its parameters.

The paper evolves as follows: In Section 2, we deliver the power Lindley survival model. In Section 3, we perform inference in the power Lindley model. In Section 4, we deliver the bivariate power Lindley survival model. In Section 5, estimation of the parameters is accomplished using pseudo-likelihood. In Section 6, a simulation study is carried out. In Section 7 two applications are made to real data sets.

2. The Power Lindley Survival Model

In this paper, we will introduce a new extension of the L model called the power L survival (PLS) distribution. The distribution function of this model is of the form

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

in which both $\theta$ and $\alpha$ are positive parameters. The corresponding density function is

(4)$f_{PLS}\left(x\right)={\displaystyle \frac{\alpha {\theta }^2}{{(\theta +1)}^{\alpha }}}(1+x){(\theta +1+\theta x)}^{\alpha -1}{e}^{-\alpha \theta x},x>0.$

If a random variable X has this density, we will write $X\sim PLS(\theta ,\alpha ).$ Figure 1 shows the form of the PLS density for various values of the parameters, varying the value of the parameter $\theta$ for the cases $\left(a\right)$ $\alpha =0.5$ and $\left(b\right)$ $\alpha =2.0$.

Observe that $PLS(\theta ,1)=L\left(\theta \right),$ and that the PLS distribution is more flexible than the L distribution and, as a consequence of observations of Lindley [1], the PLS distribution can be expected to outperform the exponential model in certain applications. The PLS distribution is thus a flexible alternative to the exponential and L models because of the introduction of the shape parameter $\alpha$.

The survival function and the hazard function of the PLS distribution are, respectively,

$S_{PLS}\left(x\right)={\left({\displaystyle \frac{\theta +1+\theta x}{\theta +1}}\right)}^{\alpha }{e}^{-\alpha \theta x}=S_L^{\alpha }\left(x\right)\mathrm{and}{r}_{PLS}\left(x\right)={\displaystyle \frac{\alpha {\theta }^2(x+1)}{\theta +1+\theta x}}=\alpha r_L\left(x\right).$

The fact that the survival function of the PLS distribution is the $\alpha$ power of the Lindley survival function is the reason for the name selected for the PLS distribution. Note that, as a consequence of the relationship beteen the L and PLS survival functions, the corresponding hazard rates are proportional.

One may verify that

• $r_{PLS}\left(x\right)$ is an increasing function in x and of $\theta ,$ for each $\alpha$.

• $r_{PLS}\left(0\right)={\displaystyle \frac{\alpha {\theta }^2}{\theta +1}}$ and ${lim}_{x\to \infty }{r}_{PLS}\left(x\right)=\alpha \theta ,$ so that ${\displaystyle \frac{\alpha {\theta }^2}{\theta +1}}<r_{PLS}\left(x\right)<\alpha \theta$.

• For $\alpha >1,$$r_{PLS}\left(x\right)>r_L\left(x\right)$ and for $\alpha <1,$ $r_{PLS}\left(x\right)<r_L\left(x\right)$.

• Since the PLS distribution has an increasing failure rate, the following chain of implications hold for it: $IFR\Rightarrow IFRA\Rightarrow NBU\Rightarrow NBUE$

A situation in which a PLS distribution will be encountered is one in which we are dealing with the lifetime of a system consisting of n independent L distributed (i.e., $L\left(\theta \right)$) components connected in series. The system lifetime is, therefore, a random variable with a $PLS(\theta ,n)$ distribution. Allowing the parameter n in this model to be a positive real number instead of an integer, we are led to the full PLS model.

The PLS distribution has finite moments of all orders. The expressions for its moments involve the incomplete gamma function denoted by $F_G(x:\alpha ,\beta )$ and defined by

$F_G(x;\alpha ,\beta )={\int }_0^x{\displaystyle \frac{u^{\alpha -1}{e}^{-(u/\beta )}}{\Gamma \left(\alpha \right){\beta }^{\alpha }}}du.$

Thus, we have

(5)$\mu =E\left(X\right)={\int }_0^{\infty }[1-F_{PLS}\left(x\right)]dx={\displaystyle \frac{e^{\alpha (\theta +1)}\Gamma \left(\alpha \right)\left[1-F_G(\theta +1;\alpha +1,\alpha )\right]}{{\alpha }^{\alpha }\theta {(\theta +1)}^{\alpha }}}.$

For $r=2,3,\dots$ the r’th moment is expressible as

(6)$\begin{array}{cc}\hfill E\left(X^r\right)& ={\displaystyle \frac{\alpha e^{\alpha (\theta +1)}}{{(\theta +1)}^{\alpha }}}{\sum }_{j=0}^r{\sum }_{k=0}^{r-j}\left({\displaystyle \genfrac{}{}{0pt}{}{r}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r-j}{k}}\right){\displaystyle \frac{{\theta }^{k+1-r}\Gamma (\alpha +j)[1-F_G(\theta +1;\alpha +j,\alpha )]}{{\alpha }^{\alpha +j}}}+\\ & {\displaystyle \frac{\alpha e^{\alpha (\theta +1)}}{{(\theta +1)}^{\alpha }}}{\sum }_{j=0}^{r+1}{\sum }_{k=0}^{r+1-j}\left({\displaystyle \genfrac{}{}{0pt}{}{r+1}{j}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r-j+1}{k}}\right){\displaystyle \frac{{\theta }^{k-r}\Gamma (\alpha +j)[1-F_G(\theta +1;\alpha +j,\alpha )]}{{\alpha }^{\alpha +j}}}\end{array}$

From these formulas, expressions can be obtained for the variance, skewness and kutosis of the distribution.

A small-scale simulation study has been conducted to evaluate the range of possible values for the coefficients of asymmetry and kurtosis for $\alpha \in (0.05,1500].$ Calculations were performed using the function integrate in R software [16]. The values of $\theta$ considered were between $0.05$ and 300. The results indicated that the range of values for the coefficients of asymmetry and kurtosis for the model were given by $\sqrt{{\beta }_1}\in [0.4710,2.8813]$ and ${\beta }_2\in [2.0236,11.4217]$. These intervals include the corresponding ranges for the Lindley model which were, respectively, $(\sqrt{2},2)$ and $(6,9).$ This confirms the observation that has already been made that the $PLS$ model is more flexible than the L model.

An expression is available for the moment-generating function of the PLS distribution; thus,

(7)$M_X\left(t\right)={\displaystyle \frac{e^{\frac{\alpha \theta -t}{\theta }(\theta +1)}\Gamma \left(\alpha \right)}{{\alpha }^{\alpha -1}{(\theta +1)}^{\alpha }}}\left[F_G\left(\theta +1;\alpha +1,{\displaystyle \frac{\alpha \theta -t}{\theta }}\right)-F_G\left(\theta +1;\alpha ,{\displaystyle \frac{\alpha \theta -t}{\theta }}\right)\right]$

for $t<\alpha \theta .$

Quantile Function

The p’th quantile ($0<p<1$) of the PLS distribution can be expressed as

$Z_p=F_L^{-1}(1-{(1-p)}^{{\displaystyle \frac{1}{\alpha }}}),$

where $F_L^{-1}(·)$ is the Lindley quantile function. We can thus use the same technique as that used by Jodrá [17] to generate variables with L distributions to instead generate PLS variables. An alternative expression for the PLS quantiles involving a special function is available in the form

$Z_p=-1-{\displaystyle \frac{1}{\theta }}-{\displaystyle \frac{1}{\theta }}{W}_{-1}(-(\theta +1)e^{-(\theta +1)}{(1-p)}^{1/\alpha }),$

where $W_{-1}(·)$ is the Lambert function with branch $-1.$

Thus, to generate a random variable, X, with a $PLS(\theta ,\alpha )$ distribution, we can generate U, a $uniform(0,1)$ variable, and set

$X=-1-{\displaystyle \frac{1}{\theta }}-{\displaystyle \frac{1}{\theta }}{W}_{-1}(-(\theta +1)e^{-(\theta +1)}{(1-U)}^{1/\alpha }).$

3. Inference

For a random sample of size n from the $PLS(\theta ,\alpha )$ distribution, the log-likelihood function of the parameter vector $(\theta ,\alpha )$, omitting the constant term, is given by

(8)$\ell (\theta ,\alpha )=nln\left({\displaystyle \frac{\alpha {\theta }^2}{{(\theta +1)}^{\alpha }}}\right)+(\alpha -1)\sum _{i=1}^{n}ln(\theta +1+\theta x_i)-\alpha \theta \sum _{i=1}^n{x}_i$

The score functions, defined as the derivatives of the log-likelihood function with respect to the parameters, are

(9)$U\left(\alpha \right)={\displaystyle \frac{n}{\alpha }}-nln(\theta +1)+\sum _{i=1}^{n}ln(\theta +1+\theta x_i)-\theta \sum _{i=1}^n{x}_i$

(10)$U\left(\theta \right)={\displaystyle \frac{2n}{\theta }}-{\displaystyle \frac{n\alpha }{\theta +1}}+(\alpha -1)\sum _{i=1}^n{\displaystyle \frac{x_i+1}{\theta +1+\theta x_i}}-\alpha \sum _{i=1}^n{x}_i$

Equating these expressions to zero, we have the likelihood equations whose solutions provide us with the maximum likelihood estimates of $\theta$ and $\alpha .$

The observed information matrix, which is minus the matrix of the second derivatives of the log-likelihood with respect to the parameters, has elements of the form

(11)$j_{\theta \theta }={\displaystyle \frac{2n}{{\theta }^2}}-{\displaystyle \frac{n\alpha }{{(\theta +1)}^2}}+(\alpha -1)\sum _{i=1}^n{\displaystyle \frac{{(x_i+1)}^2}{{(\theta +1+\theta x_i)}^2}},$

(12)$j_{\theta \alpha }={\displaystyle \frac{n}{\theta +1}}-\sum _{i=1}^n{\displaystyle \frac{x_i+1}{\theta +1+\theta x_i}}+\sum _{i=1}^n{x}_i$

(13)$j_{\alpha \alpha }={\displaystyle \frac{n}{{\alpha }^2}}.$

The elements of the expected information matrix (or Fisher information), defined as the expected values of the elements of the observed information matrix, are given by

(14)$i_{\theta \theta }={\displaystyle \frac{2}{{\theta }^2}}-{\displaystyle \frac{\alpha }{{(\theta +1)}^2}}+{\displaystyle \frac{\alpha (\alpha -1)e^{\alpha \theta }}{{\theta }^2{(\theta +1)}^{\alpha }}}\sum _{j=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha -3}{j}}\right){\displaystyle \frac{\Gamma (j+4)\left[1-F_G(1;j+4,\alpha \theta )\right]}{{\alpha }^{j+4}}}$

(15)$i_{\alpha \theta }=\mu +{\displaystyle \frac{1}{\theta +1}}-{\displaystyle \frac{\alpha e^{\alpha \theta }}{\theta {(\theta +1)}^{\alpha }}}\sum _{j=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha -2}{j}}\right){\displaystyle \frac{\Gamma (j+3)\left[1-F_G(1;j+3,\alpha \theta )\right]}{{\alpha }^{j+3}}}$

(16)$i_{\alpha \alpha }={\displaystyle \frac{1}{{\alpha }^2}}.$

A bivariate normal approximation for the joint distribution of the maximum likelihood estimates, $\widehat{\theta }$ and $\widehat{\alpha }$, can then be used to construct confidence intervals for $\theta$ and $\alpha .$

4. A Bivariate Model

For the construction of a bivariate PLS (BPLS) model, we will make use of the approach discussed by Arnold at al. [18] based on conditional distributions.

According to Arnold et al. [18], a two-dimensional random vector $(X_1,X_2)$ has a distribution that is conditionally specified, if the conditional distribution of $X_1$ given that $X_2=x_2$ for each $x_2$ is a member of a specified parametric family of distributions and also that the conditional distribution of $X_2$ given that $X_1=x_1$ for each $x_1$ is a member of a possibly different specified parametric family of distributions.

Suppose now that the joint BPLS distribution function $F_{BPLS}(x_1,x_2)$ of the random vector $(X_1,X_2)$ is such that the conditional distributions of $X_1$ given $X_2=x_2$ and the conditional distributions of $X_2$ given $X_1=x_1$ are all members of the PLS family of distributions which are absolutely continuous with respect to the Lebesgue measure. We denote this by writing

(17)$X_1|X_2=x_2\sim PL{S}_1({\theta }_1,\omega \left(x_2\right))$

and

(18)$X_2|X_1=x_1\sim PL{S}_2({\theta }_2,\tau \left(x_1\right)),$

where $\omega$ and $\tau$ are positive dependence functions which are to be determined.

In such a case, we can recognize that, for fixed choices of ${\theta }_1$ and ${\theta }_2$, we have conditionals in given one-parameter exponential families and we can identify the corresponding joint density using a result from Arnold and Strauss [19]. We can make the following argument. If $f_{X_1}\left(x_1\right)$ and $f_{X_2}\left(x_2\right)$ are marginal densities for a joint PLS density $f_{BPLS}(x_1,x_2)$ with conditional densities given by (17) and (18), then it follows that

(19)$\begin{array}{cc}\hfill f_{BPLS}(x_1,x_2)& =\tau \left(x_1\right)f_{X_1}\left(x_1\right)f_L(x_2;{\theta }_2){\{1-F_L(x_2;{\theta }_2)\}}^{\tau \left(x_1\right)-1}\hfill \\ \hfill & =\omega \left(x_2\right)f_{X_2}\left(x_2\right)f_L(x_1;{\theta }_1){\{1-F_L(x_1;{\theta }_1)\}}^{\omega \left(x_2\right)-1},\hfill \end{array}$

since both expressions are representations of the joint density. For this equality to hold, it must be the case that (see Arnold et al. [20] or Arnold and Strauss [21])

(20)$\omega \left(x_2\right)={\alpha }_1-{\alpha }_{12}ln[1-F_L(x_2;{\theta }_2)]$

and

(21)$\tau \left(x_1\right)={\alpha }_2-{\alpha }_{12}ln[1-F_L(x_1;{\theta }_1)],$

with ${\alpha }_1,{\alpha }_2,$ positive real constants and ${\alpha }_{12}\ge 0$.

Then, using the theorems that appear in the work of Arnold and Strauss [22] and Arnold et al. ([18], chapter 4), we obtain

(22)$\begin{array}{cc}\hfill f_{BPLS}(x_1,x_2)& {\displaystyle =k\left(\underline{\alpha }\right)\prod _{j=1}^2\left[\frac{{\theta }_j^2{({\theta }_j+1+{\theta }_j{x}_j)}^{{\alpha }_j-1}(1+x_j)}{{({\theta }_j+1)}^{{\alpha }_j}}\right]e^{{\displaystyle -\sum _{j=1}^2{\alpha }_j{\theta }_j{x}_j}}}\hfill \\ \hfill & \times e^{-{\alpha }_{12}\left[{\displaystyle \prod _{j=1}^2\left\{ln\left(\frac{{\theta }_j+1+{\theta }_j{x}_j}{{\theta }_j+1}\right)-{\theta }_j{x}_j)\right\}}\right]}\hfill \end{array}$

where the parameter vector $\underline{\alpha }={({\theta }_1,{\theta }_2,{\alpha }_1,{\alpha }_2,{\alpha }_{12})}^{\prime }$ and $k\left(\underline{\alpha }\right)$ is a normalizing constant.

The corresponding conditional densities are

(23)$\begin{array}{cc}\hfill f_{X_1|X_2}\left(x_1\right|x_2)& ={\displaystyle \frac{{\theta }_1^2\omega \left(x_2\right){({\theta }_1+1+{\theta }_1{x}_1)}^{\omega \left(x_2\right)-1}(1+x_1)}{{({\theta }_1+1)}^{\omega \left(x_2\right)}}}{e}^{-{\theta }_1\omega \left(x_2\right)x_1}\hfill \end{array}$

where

$\omega \left(x_2\right)={\alpha }_1-{\alpha }_{12}ln\left({\displaystyle \frac{{\theta }_2+1+{\theta }_2{x}_2}{{\theta }_2+1}}\right)+{\alpha }_{12}{\theta }_2{x}_2$

and

(24)$\begin{array}{cc}\hfill f_{X_2|X_1}\left(x_2\right|x_1)& ={\displaystyle \frac{{\theta }_2^2\tau \left(x_1\right){({\theta }_2+1+{\theta }_2{x}_2)}^{\tau \left(x_1\right)-1}(1+x_2)}{{({\theta }_2+1)}^{\tau \left(x_1\right)}}}{e}^{-{\theta }_2\tau \left(x_1\right)x_2}\hfill \end{array}$

where

$\tau \left(x_1\right)={\alpha }_2-{\alpha }_{12}ln\left({\displaystyle \frac{{\theta }_1+1+{\theta }_1{x}_1}{{\theta }_1+1}}\right)+{\alpha }_{12}{\theta }_1{x}_1$

The marginal densities are given by

(25)$\begin{array}{cc}\hfill f_{X_1}\left(x_1\right)& =k\left(\underline{\alpha }\right){\displaystyle \frac{{\theta }_1^2{({\theta }_1+1+{\theta }_1{x}_1)}^{{\alpha }_1-1}(1+x_1)e^{-{\alpha }_1{\theta }_1{x}_1}}{{({\theta }_1+1)}^{{\alpha }_1}\left[{\alpha }_2-{\alpha }_{12}\left(ln({\theta }_1+1+{\theta }_1{x}_1)-ln({\theta }_1+1)\right)+{\alpha }_{12}{\theta }_1{x}_1\right]}}\hfill \end{array}$

and

(26)$\begin{array}{cc}\hfill f_{X_2}\left(x_2\right)& =k\left(\underline{\alpha }\right){\displaystyle \frac{{\theta }_2^2{({\theta }_2+1+{\theta }_2{x}_2)}^{{\alpha }_2-1}(1+x_2)e^{-{\alpha }_2{\theta }_2{x}_2}}{{({\theta }_2+1)}^{{\alpha }_2}\left[{\alpha }_1-{\alpha }_{12}\left(ln({\theta }_2+1+{\theta }_2{x}_2)-ln({\theta }_2+1)\right)+{\alpha }_{12}{\theta }_2{x}_2\right]}}\hfill \end{array}$

The case of independence corresponds to setting ${\alpha }_{12}=0.$

The conditional distributions corresponding to the conditional densities in (23) and (24) are of the forms

(27)$F_{X_1|X_2}\left(x_1\right|x_2)=1-{\left({\displaystyle \frac{{\theta }_1+1+{\theta }_1{x}_1}{{\theta }_1+1}}\right)}^{\omega \left(x_2\right)}{e}^{-{\theta }_1\omega \left(x_2\right)x_1}$

and

(28)$F_{X_2|X_1}\left(x_2\right|x_1)=1-{\left({\displaystyle \frac{{\theta }_2+1+{\theta }_2{x}_2}{{\theta }_2+1}}\right)}^{\tau \left(x_1\right)}{e}^{-{\theta }_2\tau \left(x_1\right)x_2}.$

Contour graphs of the BPLS density are displayed in Figure 2 for two representative parameter vectors. In both cases, it is concluded that the distribution is unimodal.

Expressions for the (non-linear) regression functions can be written involving the Lambert function. Thus,

(29)$E\left(X_1\right|X_2=x_2)=-{\displaystyle \frac{{\theta }_1+1}{{\theta }_1}}-{\displaystyle \frac{1}{{\theta }_1}}{\int }_0^1{W}_{-1}\left({\displaystyle \frac{{\theta }_1+1}{e^{{\theta }_1+1}}}(y-1)\right)y^{\omega \left(x_2\right)-1}dy$

and

(30)$E\left(X_2\right|X_1=x_1)=-{\displaystyle \frac{{\theta }_2+1}{{\theta }_2}}-{\displaystyle \frac{1}{{\theta }_2}}{\int }_0^1{W}_{-1}\left({\displaystyle \frac{{\theta }_2+1}{e^{{\theta }_2+1}}}(y-1)\right)y^{\tau \left(x_1\right)-1}dy$

The variances (${\sigma }_1^2$ and ${\sigma }_2^2$), the covariance (${\sigma }_{12}$) and the correlation (${\rho }_{12}$) can be written in terms of the Lambert function and the exponential integral function (defined below).

${\sigma }_j^2={\displaystyle \frac{{({\theta }_j+1)}^2}{{\alpha }_{12}{\theta }_j^2}}{e}^{{\displaystyle \frac{{\alpha }_1{\alpha }_2}{{\alpha }_{12}}}}\left[-Ei\left(0\right)+{\displaystyle \frac{2}{{\theta }_j+1}}{h}_1\left(y\right)+{\displaystyle \frac{1}{{({\theta }_j+1)}^2}}{h}_2\left(y\right)\right]{\displaystyle \frac{e^{\frac{2{\alpha }_1{\alpha }_2}{{\alpha }_{12}}}}{{\alpha }_{12}^2{\theta }_j^2}}{\left[({\theta }_j+1)(-Ei\left(0\right))+h_1\left(y\right)\right]}^2$

where

$-Ei\left(x\right)={\int }_x^{\infty }{\displaystyle \frac{e^{-z}}{z}}dz,h_1\left(y\right)={\int }_0^{\infty }{\displaystyle \frac{W_{-1}\left(\frac{{\theta }_j+1}{e^{{\theta }_j+1}}\left(e^{({\alpha }_{12}y-{\alpha }_1{\alpha }_2)/\left({\alpha }_j{\alpha }_{12}\right)}-1\right)\right)e^{-y}}{y}}dy$

and

$h_2\left(y\right)={\int }_0^{\infty }{\displaystyle \frac{W_{-1}^2\left(\frac{{\theta }_j+1}{e^{{\theta }_j+1}}\left(e^{({\alpha }_{12}y-{\alpha }_1{\alpha }_2)/\left({\alpha }_j{\alpha }_{12}\right)}-1\right)\right)e^{-y}}{y}}dy.$

(31)$\begin{array}{cc}\hfill {\sigma }_{12}& {\displaystyle =\frac{1}{{\theta }_1{\theta }_2}{\int }_0^1{\int }_0^{1}g(y_1,y_2)\left[{\displaystyle \prod _{j=1}^2({\theta }_j+1)+\prod _{j=1}^2{W}_{-1}\left(\frac{{\theta }_j+1}{e^{{\theta }_j+1}}(y_j-1)\right)}\right]\prod _{j=1}^2{y}_j^{{\alpha }_j-1}d{y}_j+}\\ & {\displaystyle \frac{1}{{\theta }_1{\theta }_2}{\int }_0^1{\int }_0^{1}g(y_1,y_2)\left[{\displaystyle \sum _{j=1}^2\prod _{j^{\prime }=1,j^{\prime }\ne j}^2({\theta }_j+1)W_{-1}\left(\frac{{\theta }_{j^{\prime }}+1}{e^{{\theta }_{j^{\prime }}+1}}(y_{j^{\prime }}-1)\right)}\right]\prod _{j=1}^2{y}_j^{{\alpha }_j-1}d{y}_j}\end{array}$

and ${\rho }_{12}={\displaystyle \frac{{\sigma }_{12}}{{\sigma }_1{\sigma }_2}},$ where

$g(y_1,y_2)=y_1^{-{\alpha }_{12}ln\left(y_2\right)}-{\displaystyle \frac{1}{({\alpha }_1-{\alpha }_{12}ln\left(y_2\right))({\alpha }_2-{\alpha }_{12}ln\left(y_1\right))}}.$

The values of the correlation coefficient were calculated for a variety of parametric configurations and are displayed in Table 1 and Table 2. The range of correlations encountered was [−0.8667, 0.9814], which is broader than the range of correlations in many well-known bivariate survival models.

The transformation $Y_1=-ln\left(1-F_L\left(X_1\right)\right)\mathrm{and}{Y}_2=-ln\left(1-F_L\left(X_2\right)\right)$ yields the bivariate exponential conditionals model discussed in detail by Arnold et al. [18].

(32)$f_{Y_1{Y}_2}(y_1,y_2)=\mathbf{k}({\alpha }_1,{\alpha }_2,{\alpha }_{12})exp\left(-{\alpha }_1{y}_1-{\alpha }_2{y}_2-{\alpha }_{12}{y}_1{y}_2\right)$

Note that the exponential conditionals distribution exhibits only a very limited range of negative correlations.

Consistent asymptotically normal moment-based estimates of ${\alpha }_1$, ${\alpha }_2$ and ${\alpha }_{12}$ as functions of the $y_i$s are available.

${\tilde{\alpha }}_1={\displaystyle \frac{\tilde{\gamma }}{{\overline{y}}_1(\tilde{\gamma }+I(\tilde{\gamma }-1))}},{\tilde{\alpha }}_2={\displaystyle \frac{\tilde{\gamma }}{{\overline{y}}_2(\tilde{\gamma }+I(\tilde{\gamma }-1))}}\mathrm{and}{\tilde{\alpha }}_{12}={\displaystyle \frac{\tilde{\gamma }(\tilde{\gamma }-1)}{{\overline{y}}_1{\overline{y}}_2(\tilde{\gamma }+I(\tilde{\gamma }-1))}},$

where $\tilde{\gamma }={\displaystyle \frac{I}{1+{\rho }_{Y_1{Y}_2}I}}$ with ${\rho }_{Y_1{Y}_2}=cor(y_1,y_2)$ and $I=cv\left(y_1\right)cv\left(y_2\right),$ where cor is the usual Pearson correlation between $Y_1\mathrm{and}{Y}_2$ and $cv\left(y\right)=\sqrt{S_y^2}/\overline{y}.$ If we, in addition, have available consistent estimates of ${\theta }_1$ and ${\theta }_2$, we can make use of these exponential conditionals estimates in conjunction with an estimated inverse transformation to develop estimates for the parameters of the BPLS distribution.

A Bivariate Lindley Conditionals Model

If we wish to identify a bivariate Lindley survival model with Lindley conditionals (BLCs), it is tempting to merely set ${\alpha }_1={\alpha }_2=1$ in the BPLS model. This will indeed result in a valid bivariate survival model, but it will still have power Lindley conditional distributions. To identify the model with Lindley conditionals, we must return to definition of the Lindley distribution and recognize that it corresponds to a one-parameter exponential family of distributions with parameter $\theta$. Applying the result obtained by Arnold and Strauss [22], we can identify the BLC model as one with conditionals

$X_1|X_2=x_2\sim L({\theta }_1+{\theta }_{12}{x}_2),X_2|X_1=x_1\sim L({\theta }_2+{\theta }_{12}{x}_1)$

and joint density given by

(33)$f(x_1,x_2;\underline{\theta })\propto (1+x_1)(1+x_2)exp\left\{-{\theta }_1{x}_1-{\theta }_2{x}_2-{\theta }_{12}{x}_1{x}_2\right\},x_1,x_2>0.$

The marginal and conditional densities for this model are given as follows:

$f\left(x_1\right)=k\left(\theta \right){\displaystyle \frac{(x_1+1)(1+{\theta }_2+{\theta }_{12}{x}_1)}{{({\theta }_2+{\theta }_{12}{x}_1)}^2}}exp(-{\theta }_1{x}_1),x_1>0,$

$f\left(x_2\right)=k\left(\theta \right){\displaystyle \frac{(x_2+1)(1+{\theta }_1+{\theta }_{12}{x}_2)}{{({\theta }_1+{\theta }_{12}{x}_2)}^2}}exp(-{\theta }_2{x}_2),x_2>0.$

$f\left(x_1\right|x_2)={\displaystyle \frac{(x_1+1){({\theta }_1+{\theta }_{12}{x}_2)}^2}{(1+{\theta }_1+{\theta }_{12}{x}_2)}}exp(-({\theta }_1+{\theta }_{12}{x}_2)x_1),x_1>0,$

$f\left(x_2\right|x_1)={\displaystyle \frac{(x_2+1){({\theta }_2+{\theta }_{12}{x}_1)}^2}{(1+{\theta }_2+{\theta }_{12}{x}_1)}}exp(-({\theta }_2+{\theta }_{12}{x}_1)x_2),x_2>0.$

Note that the BLC density can be written in the form

(34)$f(x_1,x_2;\underline{\theta })\propto p(x_1,x_2)exp\left\{-{\theta }_1{x}_1-{\theta }_2{x}_2-{\theta }_{12}{x}_1{x}_2\right\},x_1,x_2>0.$

Thus, it can recognized as a weighted bivariate exponential conditionals density with weight function $p(x_1,x_2).$ If we wish to decide whether a given data set will be best fitted by a BPLS or by a BLC model, we cannot simply test the hypothesis $H_0:{\alpha }_1={\alpha }_2=1.$ The models are not nested and the decision regarding which will best fit the data will likely involve comparing the maximized pseudo-likelihoods of the two models.

5. Estimation

In order to estimate the parameters of a BPLS model based on a sample from that distribution, we will employ the method of pseudo-likelihood which utilizes the conditional likelihoods. The reason for this is that the nature of the normalizing constant in the BPLS density makes use of the usual maximum likelihood (ML) approach somewhat difficult. The pseudo-likelihood approach involve maximization of the logarithm of the product of the conditional densities, which avoids the necessity of dealing with the normalizing constant, $k\left(\underline{\alpha }\right)$. In this case, although $k\left(\underline{\alpha }\right)$ is troublesome, ML estimation could be accomplished by a computer-intensive search involving repeated evaluation of the expression for the normalizing constant which, though available, involves the exponential integral function. Thus,

$k\left(\underline{\alpha }\right)={\displaystyle \frac{c{e}^{-1/c}}{-Ei(1/c)}},$

where

$c={\displaystyle \frac{{\alpha }_{12}}{{\alpha }_1{\alpha }_2}}.$

In this situation, pseudo-likelihood is an attractive alternative.

Pseudo-Likelihood Estimation for the BPLS Distribution

An early reference for the technique known as pseudo-likelihood estimation is Besag [23]. In the bivariate case, the method involves replacing the log-likelihood function by the product of the two conditional likelihood functions and seeking parameter values that will maximize this conditional objective function. A convenient reference for discussion of pseudo-likelihood estimates and their properties is Arnold and Strauss [22]. Such estimates are typically consistent and asymptotically normal, though somewhat less efficient than the more elusive ML estimates.

The pseudo likelihood corresponding to a single observation $(X_1,X_2)$ from a BPLS distribution is defined to be

(35)$L_p\left(\underline{\beta }\right)=f_{X_1|X_2}\left(x_1\right|x_2)f_{X_2|X_1}\left(x_2\right|x_1).$

The pseudo-likelihood function corresponding to sample of size n from the BPLS distribution is then given by

$L_p^{\left(n\right)}\left(\underline{\beta }\right)=\prod _{i=1}^n{f}_{X_1|X_2}\left(x_{i1}\right|x_{i2})f_{X_2|X_1}\left(x_{i2}\right|x_{i1}).$

The pseudo-likelihood estimate of the parameter vector $\underline{\beta }$ is that value of $\underline{\beta }$ which maximizes $L_p^{\left(n\right)}\left(\underline{\beta }\right)$. Equivalently, it is the value that maximizes the log-pseudo-likelihood, denoted by ${\ell }_p^{\left(n\right)}\left(\underline{\beta }\right)$. For a sample from the BPLS distribution, we have

(36)$\begin{array}{cc}\hfill {\ell }_p\left(\underline{\beta }\right)& =2nlog\left({\theta }_1{\theta }_2\right)+{\sum }_{i=1}^n(ln\left(b_2\left(x_{i1}\right)\right)+ln\left(b_1\left(x_{i2}\right)\right))+{\sum }_{i=1}^{n}ln\left((1+x_{i1})(1+x_{i2})\right)\\ & +{\sum }_{i=1}^n((b_1\left(x_{i2}\right)-1)ln({\theta }_1+1+{\theta }_1{x}_{i1})+(b_2\left(x_{i1}\right)-1)ln({\theta }_2+1+{\theta }_2{x}_{i2}))\\ & -{\sum }_{i=1}^n(b_1\left(x_{i2}\right)ln({\theta }_1+1)+b_2\left(x_{i1}\right)ln({\theta }_2+1))-{\sum }_{i=1}^n({\theta }_2{x}_{i2}{b}_2\left(x_{i1}\right)+{\theta }_1{x}_{i1}{b}_1\left(x_{i2}\right)).\end{array}$

where $b_1\left(x_2\right)=\omega \left(x_2\right)$ and $b_2\left(x_1\right)=\tau \left(x_1\right).$

The pseudo-score functions are defined to be the partial derivatives of the log-pseudo-likelihood with respect to each of the parameters in the model, denoted by $U_p\left(\underline{\mathit{\beta }}\right)={(U_p\left({\theta }_1\right),U_p\left({\theta }_2\right),U_p\left({\alpha }_1\right),U_p\left({\alpha }_2\right),U_p\left({\alpha }_{12}\right))}^{\prime }.$

In the case of the BPLS model, the elements of the pseudo-score vector are given as follows:

$\begin{array}{cc}\hfill U_p\left({\theta }_j\right)& ={\displaystyle \frac{2n}{{\theta }_j}}-{\sum }_{i=1}^n{\displaystyle \frac{x_{ij}+1}{{\theta }_j+1+{\theta }_j{x}_{ij}}}+{\sum }_{i=1}^n{\displaystyle \frac{x_{ij}{b}_j\left(x_{i{j}^{\prime }}\right)\left(1-{({\theta }_j+1)}^2-{\theta }_j({\theta }_j+1)x_{ij}\right)}{({\theta }_j+1)({\theta }_j+1+{\theta }_j{x}_{ij})}}-\\ & {\displaystyle \frac{{\alpha }_{12}}{{\theta }_j+1}}{\sum }_{i=1}^n{x}_{ij}{\displaystyle \frac{1-{({\theta }_j+1)}^2-{\theta }_j({\theta }_j+1)x_{ij}}{{\theta }_j+1+{\theta }_j{x}_{ij}}}\left[{\displaystyle \frac{1}{b_{j^{\prime }}\left(x_{ij}\right)}}+ln\left({\displaystyle \frac{{\theta }_{j^{\prime }}+1+{\theta }_{j^{\prime }}{x}_{i{j}^{\prime }}}{{\theta }_{j^{\prime }}+1}}\right)-{\theta }_{j^{\prime }}{x}_{i{j}^{\prime }}\right]\end{array}$

for $j=1,2$ and $j^{\prime }=1,2$ with $j^{\prime }\ne j,$

$U_p\left({\alpha }_j\right)=\sum _{i=1}^n\left[{\displaystyle \frac{1}{b_j\left(x_{i{j}^{\prime }}\right)}}+ln\left({\displaystyle \frac{{\theta }_j+1+{\theta }_j{x}_{ij}}{{\theta }_j+1}}\right)-{\theta }_j{x}_{ij}\right]$

for $j=1,2$ and $j^{\prime }=1,2$ with $j^{\prime }\ne j$ and

$\begin{array}{cc}\hfill U_p\left({\alpha }_{12}\right)& =-{\sum }_{i=1}^n\left[{\displaystyle \frac{1}{b_2\left(x_{i1}\right)}}+ln\left({\displaystyle \frac{{\theta }_2+1+{\theta }_2{x}_{i2}}{{\theta }_2+1}}\right)-{\theta }_2{x}_{i2}\right]\left[ln\left({\displaystyle \frac{{\theta }_1+1+{\theta }_1{x}_{i1}}{{\theta }_1+1}}\right)-{\theta }_1{x}_{i1}\right]\\ & -{\sum }_{i=1}^n\left[{\displaystyle \frac{1}{b_1\left(x_{i2}\right)}}+ln\left({\displaystyle \frac{{\theta }_1+1+{\theta }_1{x}_{i1}}{{\theta }_1+1}}\right)-{\theta }_1{x}_{i1}\right]\left[ln\left({\displaystyle \frac{{\theta }_2+1+{\theta }_2{x}_{i2}}{{\theta }_2+1}}\right)-{\theta }_2{x}_{i2}\right].\end{array}$

The estimating equations consist of the elements of the pseudo-score vector set equal to zero. The solution to these equations is the vector of maximum pseudo-likelihood (MPL) estimates. The solutions are typically obtained by iterative numerical means, such as Newton–Raphson or quasi-Newton.

The asymptotic distribution MPL estimates for the BPLS model can be identified using the results of Arnold and Strauss [21].

Thus, the pseudo-likelihood estimators $\widehat{\underline{\mathit{\beta }}}$ of $\mathit{\beta }$ are consistent and asymptotically normally distributed with an asymptotic covariance matrix given by

${\Sigma }_p=J^{-1}\left(\underline{\mathit{\beta }}\right)K\left(\underline{\mathit{\beta }}\right)J^{-1}\left(\underline{\mathit{\beta }}\right)$

(see Arnold and Strauss, 1991 [19]), where for $l,m=1,2$

$K_{lm}\left(\underline{\mathit{\beta }}\right)=\mathbb{E}\left[\left\{{\displaystyle \frac{\partial {\ell }_p\left(\underline{\mathit{\beta }}\right)}{\partial {\beta }_l}}\right\}{\left\{{\displaystyle \frac{\partial {\ell }_p\left(\underline{\mathit{\beta }}\right)}{\partial {\beta }_m}}\right\}}^{\prime }\right]\mathrm{and}{J}_{lm}\left(\underline{\mathit{\beta }}\right)=-\mathbb{E}\left[{\displaystyle \frac{{\partial }^2{\ell }_p\left(\underline{\mathit{\beta }}\right)}{\partial {\beta }_l\partial {\beta }_m}}\right].$

Typically MPL estimates are less efficient than ML estimates (see Tibaldi et al., [24]), but the relative ease of computation offsets the loss of efficiency.

As a consistent estimator of the asymptotic covariance matrix of the MPL estimates, we can use the sandwich estimator proposed by Cheng and Riu [25]. This estimator can be obtained as follows.

Let $U_{pi}\left(\underline{\mathit{\beta }}\right)={\displaystyle \frac{\partial {\ell }_{pi}\left(\underline{\mathit{\beta }}\right)}{\partial \underline{\mathit{\beta }}}},$ be the vector of pseudo-scores for the ith observation, and then define

${\widehat{J}}_n\left(\underline{\mathit{\beta }}\right)=-{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\displaystyle \frac{\partial U_{pi}\left(\underline{\mathit{\beta }}\right)}{\partial {\underline{\mathit{\beta }}}^{\prime }}}{|}_{\widetilde{\underline{\mathit{\beta }}}}$, which is the sum over all of the observations of the matrix of second derivatives of

${\ell }_p\left(\underline{\mathit{\beta }}\right)$ evaluated at the pseudo-estimator $\widetilde{\underline{\mathit{\beta }}}.$

Define

${\widehat{K}}_n\left(\underline{\mathit{\beta }}\right)={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{U}_{pi}\left(\underline{\mathit{\beta }}\right)U_{pi}{\left(\underline{\mathit{\beta }}\right)}^{\prime }{|}_{\widetilde{\underline{\mathit{\beta }}}},$ then a consistent estimator of the asymptotic covariance matrix is given as follows:

$\widehat{\Sigma }\left(\widetilde{\underline{\mathit{\beta }}}\right)={\displaystyle \frac{1}{n}}{\widehat{J}}_n^{-1}\left(\widetilde{\underline{\mathit{\beta }}}\right){\widehat{K}}_n\left(\widetilde{\underline{\mathit{\beta }}}\right){\widehat{J}}_n^{-1^{\prime }}\left(\widetilde{\underline{\mathit{\beta }}}\right).$

Setting $w_{ij}={\displaystyle \frac{{\theta }_j+1+{\theta }_j{x}_{ij}}{{\theta }_j+1}},$ $q_{ij}=x_{ij}{\displaystyle \frac{1-{({\theta }_j+1)}^2-{\theta }_j({\theta }_j+1)x_{ij}}{({\theta }_j+1)({\theta }_j+1+{\theta }_j{x}_{ij})}}$ and $z_{ij}={\displaystyle \frac{(1+2{\theta }_j)x_{ij}^2+2({\theta }_j+1)x_{ij}}{{({\theta }_j+1)}^2{({\theta }_j+1+{\theta }_j{x}_{ij})}^2}},$ the elements of the matrix of second derivatives of the log-pseudo-likelihood function with respect to the parameters

$H\left(\underline{\mathit{\beta }}\right)=\left\{h_{{\beta }_j{\beta }_{j^{\prime }}}\right\}$ of the BPLS distribution are given as follows:

$\begin{array}{cc}\hfill h_{{\theta }_j{\theta }_j}& =-{\displaystyle \frac{2n}{{\theta }_j^2}}+{\alpha }_{12}{\sum }_{i=1}^n{\displaystyle \frac{b_{j^{\prime }}\left(x_{ij}\right)z_{ij}-{\alpha }_{12}{q}_{ij}^2}{b_{j^{\prime }}^2\left(x_{ij}\right)}}+{\sum }_{i=1}^n{\displaystyle \frac{{(1+x_{ij})}^2}{{({\theta }_j+1+{\theta }_j{x}_{ij})}^2}}-{\sum }_{i=1}^n{b}_j\left(x_{i{j}^{\prime }}\right)z_{ij}\\ & +{\alpha }_{12}{\sum }_{i=1}^n{z}_{ij}\left[ln\left(w_{i{j}^{\prime }}\right)-{\theta }_{j^{\prime }}{x}_{i{j}^{\prime }}\right]\end{array}$

$\begin{array}{c}\hfill h_{{\theta }_j{\theta }_{j^{\prime }}}=-{\alpha }_{12}\sum _{i=1}^n{x}_{i{j}^{\prime }}{q}_{i{j}^{\prime }}{q}_{ij}-{\alpha }_{12}\sum _{i=1}^n{x}_{ij}{q}_{ij}{q}_{i{j}^{\prime }},\end{array}$

$\begin{array}{c}\hfill h_{{\alpha }_j{\theta }_j}=\sum _{i=1}^n{x}_{ij}{q}_{ij},\end{array}\begin{array}{c}\hfill h_{{\alpha }_{j^{\prime }}{\theta }_j}=-{\alpha }_{12}\sum _{i=1}^n{\displaystyle \frac{x_{ij}{q}_{ij}}{b_{j^{\prime }}^2\left(x_{ij}\right)}},\end{array}$

$\begin{array}{c}\hfill h_{{\alpha }_{12}{\theta }_j}=-\sum _{i=1}^n{x}_{ij}{q}_{ij}\left[2(ln\left(w_{i{j}^{\prime }}\right)-{\theta }_{j^{\prime }}{x}_{i{j}^{\prime }})+{\displaystyle \frac{b_{j^{\prime }}\left(x_{ij}\right)+{\alpha }_{12}(ln\left(w_{ij}\right)-{\theta }_j{x}_{ij})}{b_{j^{\prime }}^2\left(x_{ij}\right)}}\right],\end{array}$

$\begin{array}{c}\hfill h_{{\alpha }_j{\alpha }_j}=-\sum _{i=1}^n{\displaystyle \frac{1}{b_j^2\left(x_{i{j}^{\prime }}\right)}},\end{array}\begin{array}{c}\hfill h_{{\alpha }_{j^{\prime }}{\alpha }_j}=0,\end{array}\begin{array}{cc}\hfill h_{{\alpha }_{12}{\alpha }_j}& =\sum _{i=1}^n{\displaystyle \frac{ln\left(w_{i{j}^{\prime }}\right)-{\theta }_{j^{\prime }}{x}_{i{j}^{\prime }}}{b_j^2\left(x_{i{j}^{\prime }}\right)}}\hfill \end{array}$

and

$h_{{\alpha }_{12}{\alpha }_{12}}=-\sum _{i=1}^n\sum _{j=1}^2{\left[{\displaystyle \frac{ln\left(w_{ij}\right)-{\theta }_j{x}_{ij}}{b_{j^{\prime }}\left(x_{ij}\right)}}\right]}^2.$

We then have

${\widehat{J}}_n\left(\underline{\mathit{\beta }}\right)=-{\displaystyle \frac{1}{n}}H\left(\underline{\mathit{\beta }}\right){|}_{\widetilde{\underline{\mathit{\beta }}}},$

from which we may obtain the estimated covariance matrix of Cheng and Riu [25] for the BPLS model.

6. Numerical Results

6.1. Univariate Simulation

A simulation study was carried out to investigate the behavior of the ML estimate of the parameters $\theta$ and $\alpha$ in the PLS distribution. A set of 5000 samples from the PLS distribution were simulated for each combination of parameter values $\theta =0.5,2.5,7.5$ and $\alpha =0.5,2.5,4.5$ with sample sizes $n=30,50,150$ and 300 respectively. For fixed values of the parameters $\theta$ and $\alpha$, and utilizing a sample of n$U(0,1)$ variables, we generate X’s of the following form

$X=-1-{\displaystyle \frac{1}{{\theta }_0}}-{\displaystyle \frac{1}{{\theta }_0}}{W}_{-1}(-({\theta }_0+1)e^{-({\theta }_0+1)}{(1-U)}^{1/{\alpha }_0}),$

which have the $PLS({\theta }_0,{\alpha }_0)$ distribution. For this random sample we construct the corresponding log-likelihood function $\ell ({\theta }_0,{\alpha }_0;X)$ which is used to obtain the ML estimates of the parameters using the optim function in R. This procedure was repeated $m=5000$ times and using the estimates obtained, we computed values of the bias ($\left|Bias\right|$) and root of the mean square error ($\sqrt{MSE}$) of the estimates in the forms

$|\mathrm{Bias}\left({\theta }_j\right)|={\displaystyle \frac{1}{m}}\left|\sum _{k=1}^m\left({\displaystyle \frac{{\widehat{\theta }}_j^{\left(k\right)}-{\theta }_j}{{\theta }_j}}\right)\right|;\mathrm{RMSE}\left({\theta }_j\right)=\sqrt{{\displaystyle \frac{1}{m}}\sum _{k=1}^m{\left({\widehat{\theta }}_j^{\left(k\right)}-{\theta }_j\right)}^2,}$

respectively, where ${\widehat{\theta }}_j^{\left(k\right)}$ is the estimator of ${\theta }_j$ for the jth sample and m the iterations number.

The average relative bias and average root mean squared error of the ML estimates of the parameters are presented in Table 3 and Table 4.

As is to be expected, since the estimates are consistent, the relative bias and the root mean squared errors of the estimates decrease in almost all cases as the sample size increases. The results in the Tables confirm the expected good performance of ML estimation for this model.

6.2. Simulation Study for the BPLS Distribution

A simulation study, similar to that discussed in the previous subsection, was carried out for the BPLS distribution, but in this case employing MPL estimates of the model parameters, $\underline{\beta }=({\theta }_1,{\theta }_2,{\alpha }_1,{\alpha }_2,{\alpha }_{12}),$ That is to say, realizations of the bivariate exponential conditionals distribution (discussed in detail by Arnold et al. [18], in particular; see equation 32), were simulated, and using the following inverse transformation $Y_1=-ln\left(1-F_L(X_1;{\theta }_1)\right)\mathrm{and}{Y}_2=-ln\left(1-F_L(X_2;{\theta }_2)\right)$, we obtained a sample from the $BPLS({\theta }_1,{\theta }_2,{\alpha }_1,{\alpha }_2,{\alpha }_{12}).$ The corresponding pseudo-likelihood function was then maximized numerically using the optim function provided in the R-project.

The relative bias and root mean squared errors of the estimates were investigated in a small study involving 1000 simulated random samples of each of the sizes $n=30,60,90$ and 200 for three different parametric scenarios: scenario 1 with $\underline{\beta }=(0.50,0.75,1.25,2.75,1.0)$, scenario 2 with $\underline{\beta }=(1.50,0.25,0.75,1.50,0.50)$, and scenario 3 with $\underline{\beta }=(0.25,0.50,2.50,2.0,1.50).$ The results of this simulation study are presented in Table 5.

The results in the table confirm that the absolute value of the bias and the root mean squared error of the MPL estimators of the model parameters decrease as the sample sizes increase, and that the estimators are approximately unbiased. The performance of the MPL estimators for moderate sample sizes is deemed to be satisfactory, and the performance for larger sample sizes is more than satisfactory, as is to be expected because of the consistent asymptotic normality of such estimates, as discussed in Section 5.

7. Applications of the Models

7.1. Application 1

To illustrate the utility of the PLS model, we analyze a real data set dealing with the active repair times (X) for an airborne communication transceiver. The data are taken from Jorgensen [26].

Descriptive statistics for the data set are presented in Table 6. The quantities $\sqrt{b_1}$ and $b_2$ indicate the sample asymmetry and kurtosis coefficients, respectively.

Table 6 reveals a positively skewed distribution for the variable (X) denoting active repair times. Moreover, the asymmetry and kurtosis coefficients for X are quite far from what would be expected with the L distribution.

To model the active repair times, we use the L, GL and PLS models. To compare the quality of the fit of each of the models, we use the AIC (Akaike, [27])) criterion, namely AIC $=-2\widehat{\ell }(·)+2p,$ where p denotes the number of parameters in the model. We also consider the BIC (Schwarz, [28])) criterion, namely BIC $=-2\widehat{\ell }(·)+ln\left(n\right)p,$ and the modified AIC criterion, typically called the consistent AIC (CAIC, see Bozdogan, [29])), namely CAIC=$-2\widehat{\ell }(·)+(1+ln\left(n\right))p,$ where again, p is the number of parameters for the model being considered. The best model is considered to be the one with the smallest AIC (or BIC, or CAIC).

Maximum likelihood estimators and estimated standard errors (in parenthesis) for the L, GL and PLS models were computed by maximizing likelihood using the function optim in R. The results are presented in Table 7, together with the values of the AIC, BIC and CAIC criteria for each of the models. To obtain initial values for $\theta$ and $\alpha$, the LindleyR and fitdisrplus libraries of R-project were used, with which an estimate for the $\theta$ parameter of the L distribution could be obtained. With this value, we use Equation (9), the score of the $\alpha$ parameter of the PL equal to zero (score equation for $\alpha$ of the PLS), and from this equation we obtain an initial value for the $\alpha$ parameter. From the table, we determine that model PLS exhibits the best fit to the data set, according to the AIC, BIC or CAIC criteria. The graphs in Figure 3a–c reveal that the PLS model’s fit is quite good.

Figure 3a–c show, respectively, the empirical cdf for the active repair times (solid line), the qq-plot for the PLS model, and the qq-plot with envelopes for the PLS model. In (a), the dotted line corresponds to the cdf for the PLS model, calculated with the estimates of the parameters of the PLS model.

Additionally, the Kolmogorov–Smirnov test statistic ($KS$) for the fit to the PLS model yields the value $KS=0.150$, with a corresponding p-value $=0.759$, giving support to the hypothesis that the active repair time variable has a PLS distribution.

Figure 4a,b shows the qq-plot for the L and GL distributions, calculated with the estimates of the parameters in each model; it can be observed that both models do not fit well with the data set studied.

Finally, we consider testing the hypothesis of no difference between the PLS and the L distributions for the data set under study, which corresponds to testing the hypotheses

$H_0:\alpha =1\mathrm{versus}{H}_1:\alpha \ne 1,$

using the statistic

$\Lambda ={\displaystyle \frac{{\ell }_L\left(\widehat{\mathit{\theta }}\right)}{{\ell }_{PLS}\left(\widehat{\mathit{\theta }}\right)}}.$

For the available data, we find

$-2log(\Lambda )=6.416$

which is greater than the 5% chi square critical value, ${\chi }_{1,95\%}^2=3.8414.$ Hence, the PLS model seems to be a useful alternative to be used for modeling active repair time data.

7.2. Application 2

In this example, a bivariate data set will be considered. The suitability of the BPLS model will be assessed by comparing its comportment with that of three competing bivariate models. The data set contains 51 data points corresponding to the 51 largest cities in the United States. For each city measurements are provided for the average precipitation per day in millimeters ($X_1$) and the average daily maximum temperature in Celcius degrees ($X_2$). The data are available from the US National Center for Climatalogical Data (NCDC) at https://www.ncdc.noaa.gov (accessed on 25 August 2022).

To analyse this data set, we consider three models: The Morgenstern type bivariate Lindley (BLM) distribution (see Vaidyanathan et al. [30]); the bivariate exponential of Gumbel (BEG) distribution (see Gumbel [31]); and the BPLS distribution. To compare the fits of the models, we use the AIC, BIC and CAIC criteria ML estimators and estimated standard errors (in parenthesis), as three of the bivariate models were computed using the function optim in R. For estimation of the parameters, we first fitted the marginal densities for each coordinate and combined these parametric estimates using the method explained in Section 4 to obtain the values used as initial values in the iterative process of estimation using the full bivariate data set.

For the BPLS model, estimates were obtained by maximizing the pseudo-likelihood.

Parameter estimates for model parameters are presented in Table 8, together with the AIC and BIC criteria. Contours for the BLM, BEG and the BPLS models are presented in Figure 5, which clearly indicates a better fit for the BPLS model.

For the BPLS model, we may consider

$H_0:({\alpha }_1,{\alpha }_2,{\alpha }_3)=(1,1,0)\mathrm{versus}{H}_1:({\alpha }_1,{\alpha }_2,{\alpha }_3)\ne (1,1,0)$

Note that, under $H_0$, we have a bivariate L distribution model with independent marginals for $(X_1,X_2).$

An appropriate test involves use of a Wald-type statistic whose distribution follows from the asymptotic normality of the pseudo-maximum likelihood estimator $\widehat{\underline{\mathit{\beta }}}$. This statistic can be defined as

$W_n={(C\widehat{\underline{\mathit{\beta }}}-b)}^{\prime }{\left(\widehat{G}\left(\widehat{\underline{\mathit{\beta }}}\right)\right)}^{-1}(C\widehat{\underline{\mathit{\beta }}}-b),$

where $\widehat{G}\left(\widehat{\underline{\mathit{\beta }}}\right)$ is a submatrix of $\widehat{\Sigma }\left(\widetilde{\underline{\mathit{\beta }}}\right)$ corresponding to the vector $({\alpha }_1,{\alpha }_2,{\alpha }_3),$Π$C=\left({\mathbf{0}}_{3\times 2}{\mathbf{I}}_3\right)$ and $b=(1,1,0)$, which, under the null hypothesis, follows the chi square distribution with three degrees of freedom.

With the given data, we obtain $W_n=110.9424$ with a $p_{value}<<0.05,$ indicating that the BPLS model is significantly better at the 5% level.

To evaluate the goodness-of-fit of the BPLS model to the given data set, we use the multivariate Kolmogorov–Smirnov test of goodness-of-fit proposed by Justel et al. [32]). For the case of a bivariate distribution, the bivariate Kolmogorov–Smirnov (BKS) statistic is of the form

$d_n=\underset{(x_1,x_2)\in {\mathbb{R}}^2}{sup}|F_n(x_1,x_2)-F(x_1,x_2)|,$

where $F_n$ is the empirical distribution function of the sample and F is some specified distribution function. When $F,$ is not known, the bivariate Kolmogorov–Smirnov statistic is defined by

$dn\left(F\right)=max\{D^1,D^2\},$

where

$D^1=\underset{(x_1,x_2)\in {\mathbb{R}}^2}{sup}|G_n(y_1,y_2)-y_1\times y_2|$

using the transformation $y_1=F_{X_1}\left(x_1\right)$ and $y_2=F_{X_1|X_2}\left(x_1\right|x_2)$ and

$D^2=\underset{(x_1,x_2)\in {\mathbb{R}}^2}{sup}|G_n(y_2,y_1)-y_2\times y_1|$

using the transformation $y_2=F_{X_2}\left(x_2\right)$ and $y_1=F_{X_2|X_1}\left(x_2\right|x_1),$ where $G_n$ is the empirical distribution function of the sample.

In our case involving the BPLS model, we find

$d_n\left(BPHL\right)=max\{0.1322,0.1787\}=0.1787,$

which is less than the $5\%$ critical value provided in Table 1 of Justel et al. [32] (the critical value is not provided for $n=51$, but by referring to the given values for $n=50$ and $n=60$ we can justify our claim). We may conclude that the fit of the data to the BPLS model is acceptable.

Author Contributions

Conceptualization, G.M.-F., B.C.A. and H.W.G.; software, G.M.-F.; methodology, G.M.-F., B.C.A. and H.W.G.; formal analysis, G.M.-F. and B.C.A.; investigation, G.M.-F., B.C.A. and H.W.G.; writing—original draft preparation, G.M.-F., B.C.A. and H.W.G.; writing—review and editing, G.M.-F., B.C.A. and H.W.G.; validation, G.M.-F., B.C.A. and H.W.G.; resources, H.W.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The authors are very grateful to the anonymous referees for their efforts in improving the manuscript. G. Martínez-Flórez-G acknowledge the support given by Universidad de Córdoba, Montería, Colombia.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Graphs of the PLS density for [Forumla omitted. See PDF.]0.5 (line of points and dashes), 1.25, 1.75 (broken line), 2.5 (continuous line) with (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.].

Figure 2. Contours for the BPLS model (a) BPLS (0.15, 0.25, 2.5, 2.0, 1.50) and (b) BPLS (0.25, 0.75, 1.75, 1.25, 0.75).

Figure 3. Graphs for the PLS model (a) cdf (b) qq-plot and (c) envelope.

Figure 4. q-q plot for (a) L model and (b) GL Model.

Figure 5. Contours for (a) BLM distribution, (b) BEG distribution and (c) BPLS distribution.

Table 1

Correlation coefficients for the BPLS model.

		${\theta }_1=0.50$ and ${\theta }_2=0.75$			${\theta }_1=0.75$ and ${\theta }_2=1.50$			${\theta }_1=1.30$ and ${\theta }_2=1.70$
${\alpha }_2$	${\alpha }_1$	${\alpha }_{12}=0.25$	${\alpha }_{12}=1.0$	${\alpha }_{12}=1.75$	${\alpha }_{12}=0.25$	${\alpha }_{12}=1.0$	${\alpha }_{12}=1.75$	${\alpha }_{12}=0.25$	${\alpha }_{12}=1.0$	${\alpha }_{12}=1.75$
	0.25	−0.0453	−0.0279	−0.0215	−0.1542	−0.0958	−0.0743	−0.3329	−0.2079	−0.1622
0.1	0.75	−0.0429	−0.0263	−0.0190	−0.1633	−0.1008	−0.0737	−0.3871	−0.2401	−0.1773
	1.25	−0.0389	−0.0233	−0.0157	−0.1722	−0.1036	−0.0709	−0.5090	−0.3039	−0.2112
	0.25	−0.0457	−0.0275	−0.0200	−0.1564	−0.0950	−0.0701	−0.3423	−0.2099	−0.1562
0.25	0.75	−0.0375	−0.0207	−0.0124	−0.1433	−0.0803	−0.0494	−0.3443	−0.1952	−0.1225
	1.25	−0.0290	−0.0135	−0.0050	−0.1283	−0.0611	−0.0245	−0.3791	−0.1825	−0.0775
	0.25	−0.0419	−0.0243	−0.0163	−0.1442	−0.0848	−0.0577	−0.3246	−0.1930	−0.1331
0.50	0.75	−0.0261	−0.0104	−0.0017	−0.1003	−0.0417	−0.0091	−0.2471	−0.1056	−0.0264
	1.25	−0.0118	0.0025	0.0113	−0.0526	0.0086	0.0469	−0.1577	0.0211	0.1349
	0.25	−0.0374	−0.0206	−0.0124	−0.1293	−0.0727	−0.0448	−0.2998	−0.1705	−0.1068
0.75	0.75	−0.0153	−0.0007	0.0081	−0.0593	−0.0047	0.0285	−0.1505	−0.0151	0.0685
	1.25	0.0039	0.0171	0.0259	0.0161	0.0727	0.1114	0.0467	0.2164	0.3369
	0.25	−0.0248	−0.0100	−0.0016	−0.0871	−0.0366	−0.0077	−0.2246	−0.0962	−0.0227
1.5	0.75	0.0125	0.0251	0.0338	0.0473	0.0952	0.1290	0.1302	0.2625	0.3587
	1.25	0.0437	0.0550	0.0637	0.1924	0.2426	0.2822	0.6350	0.8083	0.9552

Table 2

Correlation coefficients for the BPLS model.

		${\mathit{\theta }}_1=3.50$ and ${\mathit{\theta }}_2=5.75$			${\mathit{\theta }}_1=2.75$ and ${\mathit{\theta }}_2=6.50$			${\mathit{\theta }}_1=5.0$ and ${\mathit{\theta }}_2=2.75$
${\mathit{\alpha }}_2$	${\mathit{\alpha }}_1$	${\mathit{\alpha }}_{12}=2.0$	${\mathit{\alpha }}_{12}=2.75$	${\mathit{\alpha }}_{12}=3.50$	${\mathit{\alpha }}_{12}=2.0$	${\mathit{\alpha }}_{12}=2.75$	${\mathit{\alpha }}_{12}=3.50$	${\mathit{\alpha }}_{12}=2.0$	${\mathit{\alpha }}_{12}=2.75$	${\mathit{\alpha }}_{12}=3.50$
	2.25	0.2874	0.3062	0.3220	0.3485	0.3700	0.3878	0.3646	0.3887	0.4088
1.0	3.0	0.3265	0.3423	0.3557	0.3948	0.4125	0.4274	0.4063	0.4263	0.4433
	3.75	0.3547	0.3684	0.3802	0.4274	0.4426	0.4555	0.4367	0.4540	0.4689
	4.5	0.3763	0.3885	0.3991	0.4520	0.4653	0.4767	0.4603	0.4757	0.4890
	2.25	0.3796	0.3919	0.4027	0.4291	0.4427	0.4546	0.4946	0.5104	0.5241
2.0	3.0	0.4175	0.4278	0.4369	0.4714	0.4827	0.4926	0.5334	0.5465	0.5580
	3.75	0.4453	0.4542	0.4621	0.5021	0.5118	0.5204	0.5627	0.5739	0.5839
	4.5	0.4670	0.4748	0.4819	0.5258	0.5344	0.5420	0.5860	0.5959	0.6048
	2.25	0.4480	0.4566	0.4643	0.4939	0.5033	0.5117	0.5878	0.5987	0.6083
3.5	3.0	0.4865	0.4936	0.5000	0.5361	0.5438	0.5507	0.6263	0.6352	0.6432
	3.75	0.5151	0.5211	0.5267	0.5671	0.5737	0.5796	0.6559	0.6635	0.6704
	4.5	0.5376	0.5429	0.5478	0.5914	0.5971	0.6024	0.6798	0.6864	0.6925
	2.25	0.4988	0.5051	0.5108	0.5438	0.5506	0.5569	0.6540	0.6618	0.6689
5.5	3.0	0.5384	0.5435	0.5483	0.5868	0.5924	0.5975	0.6934	0.6998	0.7057
	3.75	0.5680	0.5724	0.5764	0.6188	0.6235	0.6278	0.7241	0.7295	0.7345
	4.5	0.5915	0.5953	0.5989	0.6440	0.6481	0.6519	0.7492	0.7539	0.7583

Table 3

Simulations for $\widehat{\theta }$ for the PLS model with 5000 iterations for $\theta$ = 0.5, 2.5, 7.5 and $\alpha$ = 0.5, 2.5, 4.5 with sample sizes n = 30, 50, 100 and 300, respectively.

		$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=100$		$\mathit{n}=300$
$\mathit{\theta }$	$\mathit{\alpha }$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$
	0.5	0.0634	0.3384	0.0632	0.3019	0.0594	0.3463	0.0246	0.1677
0.5	2.5	0.0207	0.1659	0.0116	0.1124	0.0209	0.1359	0.0005	0.0349
	4.5	0.0076	0.0768	0.0040	0.0618	0.0018	0.0404	0.0014	0.0290
	0.5	0.2408	0.7052	0.1018	0.4233	0.0413	0.3438	0.0114	0.1619
2.5	2.5	0.0266	0.7298	0.0110	0.6882	0.0052	0.3668	0.0045	0.4449
	4.5	0.0072	0.3723	0.0071	0.3507	0.0006	0.2592	0.0002	0.1827
	0.5	0.0245	0.1888	0.0161	0.1902	0.0052	0.1445	0.0003	0.0407
7.5	2.5	0.0230	0.2440	0.0118	0.2012	0.0106	0.2336	0.0021	0.1091
	4.5	0.0254	0.4401	0.0093	0.3756	0.0013	0.2473	0.0061	0.2347

Table 4

Simulations for $\widehat{\alpha }$ for the PLS model with 5000 iterations for $\theta$ = 0.5, 2.5, 7.5 and $\alpha$ = 0.5, 2.5, 4.5 with sample sizes n = 30, 50, 150, and 300, respectively.

		$\mathit{n}=30$		$\mathit{n}=50$		$\mathit{n}=100$		$\mathit{n}=300$
$\mathit{\theta }$	$\mathit{\alpha }$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$	$\left|\mathit{B}\mathit{i}\mathit{a}\mathit{s}\right|$	$\sqrt{\mathit{M}\mathit{S}\mathit{E}}$
	0.5	0.1646	0.4287	0.1086	0.3726	0.0651	0.2775	0.0333	0.1945
0.5	2.5	0.0222	0.3653	0.0292	0.3712	0.0237	0.3727	0.0170	0.2220
	4.5	0.0322	0.3272	0.0305	0.3103	0.0124	0.4404	0.0100	0.2139
	0.5	0.2217	0.5397	0.0873	0.3250	0.0298	0.1284	0.0079	0.0594
2.5	2.5	0.3146	0.8976	0.2523	0.8676	0.0868	0.4347	0.0513	0.5333
	4.5	0.1305	0.5626	0.0935	0.6519	0.0323	0.5542	0.0194	0.3495
	0.5	0.0200	0.1006	0.0109	0.0774	0.0037	0.0524	0.0014	0.0291
7.5	2.5	0.0889	0.4399	0.0494	0.3323	0.0267	0.2325	0.0083	0.1319
	4.5	0.1174	0.5963	0.0799	0.4741	0.0346	0.3322	0.0210	0.2233

Table 5

Simulations for the $BPLS\left(\underline{\mathit{\beta }}\right)$ model with 1000 iterations and with sample sizes n = 30, 60, 90 and 200, respectively.

Scenario 1
n	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$
30	0.1236	0.6551	0.3269	1.0653	1.0899	2.1701	0.6024	2.6989	1.6199	2.5336
60	0.0862	0.5256	0.251	0.8143	0.8629	1.8099	0.5326	2.3808	1.2446	2.1412
90	0.0622	0.4936	0.1920	0.6754	0.6280	1.4775	0.4337	2.0347	0.9237	1.8002
200	0.0498	0.3112	0.1328	0.5381	0.3910	1.1129	0.3562	1.7097	0.5017	1.3194
Scenario 2
n	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$
30	1.0233	0.7638	0.1750	0.4802	1.0195	1.6341	0.4081	1.7194	2.1485	2.2890
60	0.9859	0.5943	0.1240	0.3482	0.9494	1.3375	0.2200	1.3114	2.0599	1.8930
90	0.9111	0.4675	0.1092	0.2816	0.9412	1.1967	0.0978	1.1076	2.0110	1.6285
200	0.8156	0.4537	0.0769	0.1576	0.8658	0.9036	0.0213	0.8472	1.8631	1.2091
Scenario 3
n	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$	$\left|Bias\right|$	$\sqrt{MSE}$
30	0.2784	0.5608	0.1469	0.6254	0.4837	2.0776	0.7651	2.2199	0.7959	2.1892
60	0.2730	0.5033	0.1065	0.5472	0.1677	1.7553	0.6526	1.8444	0.3805	1.7982
90	0.2350	0.4009	0.0843	0.4638	0.1510	1.6180	0.6017	1.6769	0.3212	1.5325
200	0.2342	0.2880	0.0536	0.3308	0.0525	1.2179	0.3898	1.2856	0.1142	1.0140

Table 6

Descriptive statistics for the active repair times (X).

n	mean	variance	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$
40	4.0125	26.6801	2.6158	10.0226

Table 7

Parameter estimates and estimated standard errors for L, GL and PLS distributions.

Parameter	L	GL	PLS
$\theta$	0.4241(0.0306)	0.3588(0.0588)	67.9071(5.4333)
$\alpha$	–	0.7459(0.1676)	0.0036(0.0006)
AIC	199.582	199.821	195.161
BIC	201.271	203.199	198.545
CAIC	202.271	205.199	200.545

Table 8

Parameter estimates and estimated standard errors for BLM, BEG and BPLS distributions.

Parameter	BLM	BEG	BPLS
${\theta }_1$	0.0266(0.0026)		0.0082(0.0003)
${\theta }_2$	0.1880(0.0196)		0.0752(0.0010)
${\alpha }_1$		0.0120(0.0017)	5.0075(1.5483)
${\alpha }_2$		0.0807(0.0150)	3.0053(1.2481)
${\alpha }_{12}$	0.9756(0.4637)	0.9936(0.4999)	4.8559(2.6724)
AIC	854.5633	887.3073	843.1525
BIC	860.3583	893.1028	852.8116
CAIC	863.3583	896.1028	857.8116

References

1. Lindley, D.V. Fiducial Distributions and Bayes’ Theorem. J. R. Stat. Soc. Ser. B; 1958; 20, pp. 102-107. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1958.tb00278.x]

2. Ghitany, M.E.; Atieh, B.; Nadarajah, S. Lindley distribution and its application. Math. Comput. Simulat.; 2008; 78, pp. 493-506. [DOI: https://dx.doi.org/10.1016/j.matcom.2007.06.007]

3. Ghitany, M.E.; Al-qallaf, F.; Al-Mutairi, D.K.; Hussain, H.A. A two parameter weighted Lindley distribution and its applications to survival data. Math. Comput. Simulat.; 2011; 81, pp. 1190-1201. [DOI: https://dx.doi.org/10.1016/j.matcom.2010.11.005]

4. Ramos, P.L.; Louzada, F. The generalized weighted Lindley distribution: Properties, estimation, and applications. Cogent. Math.; 2016; 3, 1256022. [DOI: https://dx.doi.org/10.1080/23311835.2016.1256022]

5. Ristić, M.M.; Balakrishnan, N. The gamma exponentiated exponential distribution. J. Stat. Comput. Simul.; 2012; 82, pp. 1191-1206. [DOI: https://dx.doi.org/10.1080/00949655.2011.574633]

6. MirMostafaee, S.M.T.K.; Mahdizadeh, M.; Nadarajah, S. The beta Lindley distribution. J. Data Sci.; 2015; 13, pp. 603-626. [DOI: https://dx.doi.org/10.6339/JDS.201507_13(3).0010]

7. Bakouch, H.S.; Al-Zahrani, B.M.; Al-Shomrani, A.A.; Marchi, V.A.A.; Louzada, F. An extended Lindley distribution. J. Korean Stat. Soc.; 2012; 41, pp. 75-85. [DOI: https://dx.doi.org/10.1016/j.jkss.2011.06.002]

8. Nadarajah, S.; Bakouch, H.S.; Tahmasbi, R.A. Generalized Lindley distribution. Sankhya B; 2011; 73, pp. 331-359. [DOI: https://dx.doi.org/10.1007/s13571-011-0025-9]

9. Ghitany, M.E.; Al-Mutairi, D.K.; Balakrishnan, N.; Al-Enezi, L.J. Power Lindley distribution and associated inference. Comput. Stat. Data Anal.; 2013; 64, pp. 20-33. [DOI: https://dx.doi.org/10.1016/j.csda.2013.02.026]

10. Ashour, S.K.; Eltehiwy, M.A. Exponentiated power Lindley distribution. J. Adv. Res.; 2015; 6, pp. 895-905. [DOI: https://dx.doi.org/10.1016/j.jare.2014.08.005]

11. Asgharzadeh, A.; Nadarajah, S.; Sharafi, F. Weibull Lindley distribution. REVSTAT; 2018; 16, pp. 87-113.

12. Khokhar, J.; Khalil, R.; Shahid, N. Zografos Balakrishnan Power Lindley Distriution. J. Data Sci.; 2020; 18, pp. 279-298.

13. Algarni, A. On a new generalized lindley distribution: Properties, estimation and applications. PLoS ONE; 2021; 16, e0244328. [DOI: https://dx.doi.org/10.1371/journal.pone.0244328] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33626075]

14. Chhetri, S.; Mdziniso, N.; Ball, C. Cubic Rank Transmuted Lindley Distribution with Applications. Rev. Colomb. Estadística-Appl. Stat.; 2022; 45, pp. 65-83. [DOI: https://dx.doi.org/10.15446/rce.v45n1.93548]

15. Barlow, R.E.; Proschan, F. Statistical Theory of Reliability and Life Testing. Probability Models; To Begin With: Silver Springs, MD, USA, New York, NY, USA, 1981.

16. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021; Available online: https://www.R-project.org/ (accessed on 15 January 2024).

17. Jodrá, P. Computer generation of random variables with Lindley or Poisson–Lindley distribution via the Lambert W function. Math. Comput. Simulat.; 2010; 81, pp. 851-859. [DOI: https://dx.doi.org/10.1016/j.matcom.2010.09.006]

18. Arnold, B.C.; Castillo, E.; Sarabia, J.M. Conditional Specification of Statistical Models; Springer Series in Statistics; Springer: Berlin/Heidelberg, Germany, 1999.

19. Arnold, B.C.; Strauss, D.J. Bivariate distributions with conditionals in prescribed exponential families. J. R. Stat. Soc. Ser. B.; 1991; 53, pp. 365-375. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1991.tb01829.x]

20. Arnold, B.C.; Martínez-Flórez, G.; Gómez, H.W. Bivariate Proportional Hazard Models: Structure and Inference. Symmetry; 2022; 14, 2073. [DOI: https://dx.doi.org/10.3390/sym14102073]

21. Arnold, B.C.; Strauss, D.J. Bivariate distributions with exponential conditionals. J. Am. Statist. Assoc.; 1988; 83, pp. 522-527. [DOI: https://dx.doi.org/10.1080/01621459.1988.10478627]

22. Arnold, B.C.; Strauss, D.J. Pseudolikelihood estimation: Some examples. Sankhya Ser. B; 1991; 53, pp. 233-243.

23. Besag, J. Statistical Analysis of Non-Lattice Data. J. R. Stat. Soc. Ser. D; 1975; 24, pp. 179-195. [DOI: https://dx.doi.org/10.2307/2987782]

24. Tibaldi, F.; Molenberghs, G.; Burzykowsky, T.; Geys, H. Pseudolikelohood estimation for a marginal multivariate survival model. Stat. Med.; 2004; 23, pp. 947-963. [DOI: https://dx.doi.org/10.1002/sim.1664]

25. Cheng, C.; Riu, J. On estimating linear relationships when both variables are subject to heteroscedastic measurement errors. Technometrics; 2006; 48, pp. 511-519. [DOI: https://dx.doi.org/10.1198/004017006000000237]

26. Jorgensen, B. Statistical Properties of the Generalized Inverse Gaussian Distribution; Springer: New York, NY, USA, 1982.

27. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control.; 1974; 19, pp. 716-723. [DOI: https://dx.doi.org/10.1109/TAC.1974.1100705]

28. Schwarz, G. Estimating the dimension of a model. Ann. Stat.; 1978; 6, pp. 461-464. [DOI: https://dx.doi.org/10.1214/aos/1176344136]

29. Bozdogan, H. Model selection and Akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika; 1987; 52, pp. 345-370. [DOI: https://dx.doi.org/10.1007/BF02294361]

30. Vaidyanathan, V.S.; Varghese, A.S. Morgenstern type bivariate Lindley distribution. Stat. Optim. Inf. Comput.; 2016; 4, pp. 132-146. [DOI: https://dx.doi.org/10.19139/soic.v4i2.183]

31. Gumbel, E.J. Bivariate Exponential Distributions. J. Am. Statist. Assoc.; 1960; 55, pp. 698-707. [DOI: https://dx.doi.org/10.1080/01621459.1960.10483368]

32. Justel, A.; Peña, D.; Zamar, Z. A multivariate Kolmogorov-Smirnov test of goodness of fit. Stat. Probab. Lett.; 1997; 35, pp. 251-259. [DOI: https://dx.doi.org/10.1016/S0167-7152(97)00020-5]