The Nadarajah-Haghighi Lindley distribution

Abstract

Abstract: We define a new lifetime model based on compounding the Lindley and Nadarajah-Haghighi distributions. The proposed distribution is very competitive to other lifetime models. Some of its mathematical properties are investigated including generating function, mean residual life, moments, Bonferroni and Lorenz curves and mean deviations. We discuss the estimation of the model parameters by maximum likelihood. We provide a simulation study and two applications to real data for illustrative purposes. We prove empirically that the new distribution yields good fits to both data sets, and it can be a useful alternative for other classical lifetime models.

Key words
Compounding approach; Exponential distribution; Lifetime data; Lindley distribution; Nadarajah-Haghighi distribution

1 - INTRODUCTION

The Lindley distribution was pioneered by(Lindley 1958) in the context of fiducial and Bayesian inference. It is a mixture of the exponential and length-biased exponential distributions. Let $Y$ be a Lindley random variable with parameter $\gamma >0$ having probability density function (pdf)

$g(y)=\frac{{\gamma }^2}{1+\gamma }(1+y)e^{-\gamma y},y>0,$

where the mixing proportion is $\gamma /(1+\gamma )$. The survival function of $Y$ is

$\overline{G}(y)=\frac{(1+\gamma +\gamma y)}{1+\gamma }{e}^{-\gamma y}.$

Various of its statistical properties were discussed in details by(Ghitany, Atieh, and Nadarajah 2008). The authors also showed that the Lindley distribution is quite competitive with the exponential distribution. (Gupta and Singh 2013) studied the parameter estimation of this distribution with hybrid censored data. (Krishna and Kumar 2011) considered the estimation of its model parameters for progressively type II right censored sample and(Mazucheli and Achcar 2011) applied this distribution to competing risks in lifetime data.

In distribution theory context, some generalizations are obtained based on transformations of the Lindley distribution. We refer the reader to(Nadarajah, Bakouch, and Tahmasbi 2011) for the generalized (or exponentiated) Lindley,(Bakouch et. al 2012) for the extended Lindley,(Ghitany et. al 2011) and(Al-Mutairi, Ghitany, and Kundu 2015) for the weighted Lindley,(Ghitany et. al 2013) for the power Lindley and(Ashour and Eltehiwy 2015) for the exponentiated power Lindley distributions.

Another technique that has been considered is the discrete-continuous compounding approach. It is defined by the minimum of $N$ independent and identically continuous random variables, where $N$ is a discrete random variable. (Adamidis and Loukas 1998) pioneered this method by introducing the exponential geometric distribution. We also find in the literature some models obtained from compositions of the Lindley and other discrete distributions. (Sankaran 1970) defined the discrete Poisson-Lindley by combining the Poisson and Lindley distributions.(Zamani and Ismail 2010) presented the negative binomial Lindley distribution. The zero-truncated Poisson-Lindley and Pareto Poisson-Lindley distributions were introduced by(Ghitany, Al-Mutairi, and Nadarajah 2008) and(Asgharzadeh, Bakouch, and Esmaeili 2013), respectively.

The Nadarajah-Haghighi (NH) distribution is a generalization of the exponential distribution first defined by (Nadarajah and Haghighi 2011). Let $Z$ denote a NH random variable with parameters $\alpha >0$ and $\lambda >0$. The pdf and survival function of $Z$

$q(z)=\alpha \lambda (1+\lambda z{)}^{\alpha -1}{e}^{1-(1+\lambda z{)}^{\alpha }}$

and

$\overline{Q}(z)=e^{1-(1+\lambda z{)}^{\alpha }},$

respectively. Some generalizations of the NH distribution have been proposed in recent years, such as the exponentiated Nadarajah-Haghighi (Lemonte (2013)) and gamma Nadarajah-Haghighi (Bourguignon et. al (2015)), among others. By using the discrete-continuous compounding approach, we have the Poisson gamma Nadarajah-Haghighi (Ortega et. al (2015)) and geometric Nadarajah-Haghighi (Marinho (2016)) distributions.

A comprehensive review on compounding method for generating distributions can be found in (Tahir et. al 2016). They pointed out a different compounding approach by taking the minimum between two continuous distributions. In the discrete-continuous compositions, $N$ is a discrete random variable representing the number of identical elements having some continuous distribution. For the continuous-continuous compositions, we suppress the condition to be identically distributed and set $N=2$. Some well-known continuous-continuous compounded models are the additive Weibull ((Xie and Lai 1995), Lemonte, Cordeiro, and Ortega (2014)), exponential-Weibull (Cordeiro et al. 2014) and generalized exponential-exponential (Popovíc, Ristíc, and Cordeiro (2015) ) distributions, among others.

In this paper, we introduce a new continuous-continuous compounded model referred to as the Nadarajah-Haghighi Lindley (NHL) distribution. The new three-parameter distribution is obtained by compounding the Lindley and NH distributions. We assume that $Y$ and $Z$ are independent random variables and define $X=min(Y,Z)$ as a NHL random variable, whose survival function is given by

$\overline{F}(x)=\overline{G}(x)\overline{Q}(x).$

The cumulative distribution function (cdf) of $X$ is

$F(x)=1-\frac{(1+\gamma +\gamma x)}{1+\gamma }{e}^{1-\gamma x-(1+\lambda x{)}^{\alpha }},$ (1)

where $x>0.$ The parametric space of the cdf in(1) can be defined as

$\begin{array}{c}\hfill \Theta =\{(\alpha ,\lambda ,\gamma ):\alpha \ge 0,\gamma \ge 0,\lambda \ge 0,max(\alpha ,\gamma )>0,max(\lambda ,\gamma )>0\}.\end{array}$

The pdf and hazard rate function (hrf) of $X$ are given by

$f(x)=\frac{(1+\gamma +\gamma x)[\gamma +\alpha \lambda (1+\lambda x{)}^{\alpha -1}]-\gamma }{\gamma +1}{e}^{1-\gamma x-(1+\lambda x{)}^{\alpha }}$ (2)

and

$\begin{array}{c}\hfill h(x)=\frac{(1+\gamma +\gamma x)[\gamma +\alpha \lambda (1+\lambda x{)}^{\alpha -1}]-\gamma }{\gamma +\gamma x+1},\end{array}$

respectively. Henceforth, we consider that $X\sim \mathrm{\text{NHL}}(\alpha ,\lambda ,\gamma )$. The proposed distribution contains as special models some well-known distributions. For $\gamma =0$, the NHL reduces to the NH distribution. If $\gamma =0$ and $\alpha =1$, we have the exponential distribution. For $\alpha =0$ or $\lambda =0$, we have the Lindley distribution.

Figure1 displays plots of the pdf of X for some parameter values. The new density presents decreasing and reverse J shaped curve. Figure2 reveals that the NHL distribution can have decreasing, increasing, upside-down bathtub and bathtub-shaped hazard functions. This feature makes the new distribution very competitive to the Weibull, gamma and exponential distributions that exhibit only monotonic hazard rates. According to (Nadarajah, Bakouch, and Tahmasbi 2011) this is a major weakness because most empirical life systems have bathtub shapes for their hrfs.

Figure 1
Plots of the NHL density.

Figure 2
Plots of the NHL hazard function.

The NHL distribution has the following theoretical motivations: (i) It can be useful in engineering and reliability for modeling a system having two sub-systems functioning in series independently at a given time. If $Y$ and $Z$ denote the lifetimes of these independent sub-systems following the Lindley and NH distributions, then the lifetime $X$ of the system has the NHL distribution. (Cordeiro, Ortega, and Lemonte 2014) studied a similar situation for the exponential-Weibull model. (ii) The stochastic representation $X=min\{Y,Z\}$ can arise in several biological and medical applications. (iii) (Nadarajah and Haghighi 2011) mentioned some advantages of the NH model such as the ability to model data with mode fixed at zero and the fact that it can be interpreted as a truncated Weibull distribution. The NHL distribution also accumulates this advantages since it has as special model the NH distribution. (iv) Some of its mathematical properties are easily obtained. Further, the NHL distribution has some practical motivations: (i) It can be used for modeling data with bathtub and unimodal failure rates, which are very common in many applied areas. (ii) The finite sample behavior of the maximum likelihood estimates of the its parameters is adequate. (iii) It is really a very competitive model to well-known distributions such as the Weibull, exponentiated Weibull, NH and Lindley distributions as proved empirically in two applications to real data in sec 7.

The rest of this paper is outlined as follows. In secs 2-4, we obtain a range of mathematical properties of the NHL distribution. In sec 5, we consider the maximum likelihood method to estimate the model parameters. We perform a simulation study in sec 6. Two real data applications are provided in sec 7. Some concluding remarks are offered in sec 8.

2 - MEAN RESIDUAL LIFE

The mean residual life is a relevant characteristic to the design of safe systems in a wide variety of applications in engineering and reliability. Given that a component survives up to time $x>0$, the residual life is defined by

$m(x)=IE(X-x|X>x)=\frac{1}{1-F(x)}\underset{x}{\overset{\infty }{\int }}[1-F(t)]dt.$

It represents the period beyond $x$ until the time of failure. Note that

$m(x)=\frac{e^{x\gamma +(1+\lambda x{)}^{\alpha }}}{(1+\gamma +\gamma x)}\underset{x}{\overset{\infty }{\int }}(1+\gamma +\gamma t)e^{-\gamma t-(1+\lambda t{)}^{\alpha }}dt.$

We consider the integral

$J=\underset{x}{\overset{\infty }{\int }}(1+\gamma +\gamma t)e^{-\gamma t-(1+\lambda t{)}^{\alpha }}dt.$

Setting $u=(1+\lambda t{)}^{\alpha }$, we have $t=(u^{1/\alpha }-1)/\lambda .$ So, the above integral reduces to

$J=\frac{1}{\alpha \lambda }\underset{(1+\lambda x{)}^{\alpha }}{\overset{\infty }{\int }}{u}^{(1-\alpha )/\alpha }[1+\gamma +\gamma {\lambda }^{-1}(u^{1/\alpha }-1)]e^{-u-\gamma {\lambda }^{-1}(u^{1/\alpha }-1)}du.$

By expanding

$e^{-\gamma {\lambda }^{-1}(u^{1/\alpha }-1)}=\sum _{i=0}^{\infty }\frac{(-1{)}^i{\gamma }^i(u^{1/\alpha }-1{)}^i}{{\lambda }^{i}i!},$

and using the binomial theorem for $(u^{1/\alpha }-1{)}^i$, we can write

$\begin{array}{ccc}\hfill m(x)& =\hfill & \hfill \frac{e^{\gamma x+(1+\lambda x{)}^{\alpha }}}{\alpha \lambda (1+\gamma +\gamma x)}\sum _{i=0}^{\infty }{\left(\frac{\gamma }{\lambda }\right)}^i{\Psi }_i(\alpha ,\lambda ,\gamma ,x),\end{array}$ (3)

where

$\begin{array}{ccc}\hfill {\Psi }_i(\alpha ,\lambda ,\gamma ,x)& =\hfill & \hfill \sum _{j=0}^i\frac{(-1{)}^{-j}}{(i-j)!j!}[\frac{\gamma }{\lambda }\Gamma (\frac{j+2}{\alpha },(1+\lambda x{)}^{\alpha })\\ \hfill & +\hfill & \hfill (1+\gamma -\frac{\gamma }{\lambda })\Gamma (\frac{j+1}{\alpha },(1+\lambda x{)}^{\alpha })].\end{array}$

and $\Gamma (a,z)=\Gamma (a)-\gamma (a,z)={\int }_z^{\infty }{t}^{a-1}{e}^{-t}dt$ is the upper incomplete gamma function. Note that (3) can be approximated by setting large values in the upper limit of the sum, say $N$. Hence,

$\begin{array}{ccc}\hfill m(x)& \approx \hfill & \hfill \frac{e^{\gamma x+(1+\lambda x{)}^{\alpha }}}{\alpha \lambda (1+\gamma +\gamma x)}\sum _{i=0}^N{\left(\frac{\gamma }{\lambda }\right)}^i{\Psi }_i(\alpha ,\lambda ,\gamma ,x).\end{array}$ (4)

For illustrative purposes, we provide a numerical study to analyze de convergence of the expansion (4) by comparing its results with those from numerical integration.Table I presents the calculations for $N=5,10,15$ and $20$ and $x=1$. This expansion presents good approximations for the mean residual life. All computations are obtained using the R software. We provide the scripts for these calculations in the Appendix A.

TABLE I
Mean residual life obtained from(4) and by numerical integration.

m(x)	α	λ	γ	Numerical integration	Series expansion
					N=5	N=10
m(1)	1.2	1.5	0.5	0.3970	0.3962	0.3970
	1.5	2.0	1.5	0.1519	0.1081	0.1519
	2.5	1.5	0.5	0.0629	0.0629	0.0629

3 - GENERATING FUNCTION AND MOMENTS

We denote by $M(t)$ the moment generating function (mgf) of $X$. From Equation (2), we obtain

$\begin{array}{cc}\hfill M(t)& =\frac{e}{1+\gamma }{\int }_0^{\infty }\{(\gamma +\gamma x+1)[\gamma +\alpha \lambda (1+\lambda x{)}^{\alpha -1}]-\gamma \}{e}^{x(t-\gamma )-(1+\lambda x{)}^{\alpha }}dx.\hfill \end{array}$

Setting $u=(1+\lambda x{)}^{\alpha }$, we have

$\begin{array}{cc}\hfill M(t)& =\frac{e}{\alpha \lambda (1+\gamma )}{\int }_1^{\infty }{u}^{\frac{1-\alpha }{\alpha }}\{[\gamma +\gamma {\lambda }^{-1}(u^{1/\alpha }-1)+1](\gamma +\alpha \lambda u^{\frac{\alpha -1}{\alpha }})-\gamma \}\hfill \\ \hfill & \times exp\{\frac{(u^{1/\alpha }-1)(t-\gamma )}{\lambda }-u\}du.\hfill \end{array}$

By expanding

$exp\left\{\frac{(u^{1/\alpha }-1)(t-\gamma )}{\lambda }\right\}=\sum _{i=0}^{\infty }\frac{(u^{1/\alpha }-1{)}^i(t-\gamma {)}^i}{{\lambda }^{i}i!},$

using the binomial theorem for $(u^{1/\alpha }-1{)}^i$ and, after some algebra, we obtain (for $\gamma >0$ and $t<\gamma$)

$\begin{array}{cc}\hfill M(t)& =\frac{e}{1+\gamma }\sum _{i=0}^{\infty }{\left(\frac{\gamma -t}{\lambda }\right)}^i{\Psi }_i(\alpha ,\lambda ,\gamma ),\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\Psi }_i(\alpha ,\lambda ,\gamma )=& \sum _{j=0}^i\frac{(-1{)}^j}{j!(i-j)!}[\frac{\gamma }{\lambda }\Gamma (\frac{j+1}{\alpha }+1,1)+\frac{{\gamma }^2}{\alpha {\lambda }^2}\Gamma (\frac{j+2}{\alpha },1)\hfill \\ \hfill & +(1+\gamma -\frac{\gamma }{\lambda })\Gamma (\frac{j+\alpha }{\alpha },1)+\frac{{\gamma }^2(\lambda -1)}{\alpha {\lambda }^2}\Gamma (\frac{j+1}{\alpha },1)].\hfill \end{array}$

The generating function $M(t)$ is useful for computing moments of the NHL distribution by differentiation. Thus, to considering values of the parameters $\alpha ,\lambda$ and $\gamma$ for which the above expansion converges, the $n$th derivative gives

$\begin{array}{c}\hfill M^{(n)}(t)=\frac{(-1{)}^{n}e}{1+\gamma }\sum _{i=n}^{\infty }(i{)}_n\frac{(\gamma -t{)}^{i-n}}{{\lambda }^i}{\Psi }_i(\alpha ,\lambda ,\gamma ),\end{array}$

where $(i{)}_n=i(i-1)\cdots (i-n+1)$ is the falling factorial. Therefore,

$\begin{array}{c}\hfill \mu {\prime }_n=M^{(n)}(0)=\frac{(-1{)}^{n}e}{1+\gamma }\sum _{i=n}^{\infty }(i{)}_n\frac{{\gamma }^{i-n}}{{\lambda }^i}{\Psi }_i(\alpha ,\lambda ,\gamma ).\end{array}$

It is possible to obtain an approximation for $\mu {\prime }_n$ by truncating the previous expansion. Hence, for $N$ large enough,

$\begin{array}{c}\hfill \mu {\prime }_n\approx \frac{(-1{)}^{n}e}{1+\gamma }\sum _{i=n}^N(i{)}_n\frac{{\gamma }^{i-n}}{{\lambda }^i}{\Psi }_i(\alpha ,\lambda ,\gamma ).\end{array}$ (5)

We use the R software to compute $\mu {\prime }_n$ from(5) and compare the results with those obtained by numerical integration. TableII presents the results for some parametric values and $N=5,10,15$ and $20$. In fact, they illustrate that the proposed expansion provides reasonable approximations for the moments. For $n=1$, note that the series converges for all selected parameterizations even for very small $N=5$. More terms may be required when $n$ increases. We provide the scripts for these calculations in Appendix A.

TABLE II
First four moments obtained from(5) and by numerical integration.

$\mu {\prime }_n$	α	λ	γ	Numerical integration	Series expansion
					$N=5$	$N=10$	$N=15$
$\mu {\prime }_1$	1.2	1.5	0.5	0.45973	0.46241	0.45973	0.45973
	1.5	2.0	1.5	0.22280	0.23134	0.22278	0.22280
	2.5	1.5	0.5	0.18676	0.18677	0.18676	0.18676
$\mu {\prime }_2$	1.2	1.5	0.5	0.39537	0.36964	0.39540	0.39537
	1.5	2.0	1.5	0.09018	0.06314	0.09031	0.09018
	2.5	1.5	0.5	0.05757	0.05756	0.05757	0.05757
$\mu {\prime }_3$	1.2	1.5	0.5	0.48550	0.67914	0.48488	0.48550
	1.5	2.0	1.5	0.05092	0.11772	0.05013	0.05092
	2.5	1.5	0.5	0.02317	0.02327	0.02317	0.02317
$\mu {\prime }_4$	1.2	1.5	0.5	0.76476	−0.28626	0.77456	0.76475
	1.5	2.0	1.5	0.03616	−0.08171	0.04024	0.03615
	2.5	1.5	0.5	0.01116	0.01057	0.01116	0.01116

The central moments $({\mu }_n)$ and cumulants $({\kappa }_n)$ of $X$ can be determined from these raw moments using well-known relationships. An alternative expression for the skewness of $X$ can be expressed as (MacGillivray (1986))

$\rho (u)=\rho (u;\alpha ,\beta ,\gamma )=\frac{Q(1-u)+Q(u)-2Q(1/2)}{Q(1-u)-Q(u)},$

where $u\in (0,1)$ and $Q(\cdot )=F^{-1}(\cdot )$ is the quantile function (qf) of $X$. If $X$ follows a symmetric distribution, $[Q(1-u)+Q(u)]/2$ equals the median of $X$ for all $u$, and the numerator of $\rho (u)$ is zero. Thus, the farther $\rho (u)$ is from the horizontal line $\rho (u)=0$, means higher asymmetry. Decrescent form in $\rho (u)$ indicates positive asymmetry. Plots of the MacGillivray skewness for some parameter values are displayed in Figure3. Note that for $\lambda =1.0$ and $\gamma =0.1$, increases in $\alpha$ imply less skewed to the left allowing to model data with a distribution of heavy tails. An opposite result is observed by setting $\alpha =0.5$ and $\gamma =0.1$ and increasing $\lambda$. Less asymmetric distributions arise for higher values of $\lambda$.

Figure 3
Skewness of the NHL distribution for some parameter values.

4 - INCOMPLETE MOMENTS

The $s$th incomplete moment of $X$ is defined by $m_s(y)={\int }_0^y{x}^{n}f(x)dx$. Thus, by inserting (2) in $m_s(y)$, we have

$\begin{array}{cc}\hfill m_s(y)& =\frac{e}{1+\gamma }{\int }_0^y{x}^s\{(\gamma +\gamma x+1)[\gamma +\alpha \lambda (1+\lambda x{)}^{\alpha -1}]-\gamma \}{e}^{-\lambda x-(1+\lambda x{)}^{\alpha }}dx.\hfill \end{array}$

Using the exponential and binomial expansions, we obtain

$\begin{array}{cc}\hfill m_s(y)& =\frac{(-1{)}^{s}e}{\alpha {\lambda }^{1+s}(1+\gamma )}\sum _{i=0}^{\infty }\frac{1}{i!}{\left(\frac{\gamma }{\lambda }\right)}^i{\Psi }_{s,i}(y;\alpha ,\lambda ,\gamma ),\hfill \end{array}$ (6)

where

$\begin{array}{cc}\hfill {\Psi }_{s,i}(y;\alpha ,\lambda ,\gamma )& =\sum _{j=0}^{s+i}\frac{(-1{)}^j(s+i{)}_j}{j!}[\alpha \lambda {\Gamma }_y^{*}(\frac{j+1}{\alpha }+1,)+\frac{{\gamma }^2}{\lambda }{\Gamma }_y^{*}\left(\frac{j+2}{\alpha }\right)\hfill \\ \hfill & +\frac{{\gamma }^2(\lambda -1)}{\lambda }{\Gamma }_y^{*}\left(\frac{j+1}{\alpha }\right)+\alpha (\lambda -\gamma +\lambda \gamma ){\Gamma }_y^{*}(\frac{j}{\alpha }+1)],\hfill \end{array}$

and ${\Gamma }_y^{*}(a)=\gamma (a,(1+\lambda y{)}^{\alpha })-\gamma (a,1)={\int }_1^{(1+\lambda y{)}^{\alpha }}{t}^{a-1}{e}^{-t}dt$.

We compute approximations for $m_s(y)$ by setting values in the upper limit of the sum of equation (6). The results for some parameter values and $y$ are shown in TableIII. This expansion presents good approximations for the incomplete moments. We provide the scripts for these calculations in Appendix A.

TABLE III
Incomplete moments obtained from(6) and by numerical integration.

ms(x)	α	λ	γ	Numerical integration	Series expansion
					$N=2$	$N=5$
m1(0.5)	1.2	1.5	0.5	0.1380	0.1381	0.1380
	1.5	2.0	1.5	0.1559	0.1604	0.1559
	2.5	1.5	0.5	0.1609	0.1610	0.1609
m1(1)	1.2	1.5	0.5	0.3072	0.3095	0.3072
	1.5	2.0	1.5	0.2166	0.2489	0.2161
	2.5	1.5	0.5	0.1866	0.1869	0.1866
m2(0.5)	1.2	1.5	0.5	0.0421	0.0422	0.0421
	1.5	2.0	1.5	0.0421	0.0439	0.0421
	2.5	1.5	0.5	0.0417	0.0417	0.0417
m2(1)	1.2	1.5	0.5	0.1663	0.1681	0.1663
	1.5	2.0	1.5	0.0829	0.1062	0.0825
	2.5	1.5	0.5	0.0574	0.0576	0.0574

Equation(6) is the main result of this sec. In various practical situations, the shape of many distributions can be usefully described by the incomplete moments. For example, the mean deviations about the mean and median and Lorenz and Bonferroni curves are simple applications of the first incomplete moment. Nowadays in the literature, these curves are the most used curves in inequality analysis. Also, note that (6) can be used to approximate $\mu {\prime }_n$ by setting large values in $y$.

For a given probability $\pi$, the Lorenz and Bonferroni curves are defined by $B(\pi )=m_1(q)/(\pi {\mu }_1^{\prime })$ and $L(\pi )=m_1(q)/{\mu }_1^{\prime }$, respectively, where $q=Q(\pi )$ is the qf of $X$ evaluated at $\pi$. Plots of the Lorenz and Bonferroni curves for some parameter values are displayed in Figure4.

Figure 4
(a) Bonferroni curves for the NHL model for some parameter values; (b) Lorenz curves for the NHL model for some parameter values.

Further, the dispersion of $X$ can be measured to some extent by the mean deviations around the mean ${\mu }_1^{\prime }$ and median $m$ given by ${\delta }_1=2{\mu }_1^{\prime }F({\mu }_1^{\prime })-2m_1({\mu }_1^{\prime })$, and ${\delta }_2={\mu }_1^{\prime }-2m_1(m)$, where $m=Q(0.5)$.

5 - MAXIMUM LIKELIHOOD ESTIMATION

Let $x_1,\cdots ,x_n$ be a sample of size $n$ from the NHL$(\alpha ,\lambda ,\gamma )$ distribution and $𝛉=(\alpha ,\lambda ,\gamma {)}^T$ the parameter vector of interest. The log-likelihood function for $𝛉$ based on this sample is given by

$\begin{array}{ccc}\hfill \ell (𝛉)& =\hfill & \hfill n-\gamma \sum _{i=1}^n{x}_i-\sum _{i=1}^n(1+\lambda x_i{)}^{\alpha }-nlog(1+\gamma )\\ \hfill & +\hfill & \hfill \sum _{i=1}^{n}log\{(1+\gamma +\gamma x_i)[\gamma +\alpha \lambda (1+\lambda x_i{)}^{\alpha -1}]-\gamma \}.\end{array}$ (7)

The maximum likelihood estimates (MLEs) of the model parameters can be obtained by maximizing(7). There are several routines for numerical maximization in the R program (optim function), SAS (PROC NLMIXED), Ox (MaxBFGSsub-routine), among others. Alternatively, we can differentiate (7) and solve the resulting nonlinear likelihood equations using the quasi-Newton BFGS and Newton-Raphson algorithms. The score components can be obtained from the authors under request.

Under standard regularity conditions and for large $n$, the distribution of $\widehat{𝛉}$ can be approximated by a trivariate normal $N_3(0,{𝐉(\widehat{𝛉})}^{-1})$ distribution, where $𝐉(𝛉)$ is the observed information matrix. In addition, this approximation can also be used for constructing confidence regions and testing hypotheses on the parameters.

6 - SIMULATION STUDY

In this sec, we provide a Monte Carlo simulation study to evaluate the adequacy of the MLEs of the parameters of the NHL distribution. The simulations are performed by generating observations from eight scenarios with different parameter combinations. Note that the cdf in (1) can only be inverted numerically. However, it is possible generating random numbers by using the qfs of the Lindley and NH distributions.

(Jodrá 2010) used the Lambert $W$ function to generate random numbers from the Lindley distribution. If $Y\sim$ Lindley$(\gamma )$ has cdf $F_Y(y)$, the qf of $Y$ is

$\begin{array}{c}\hfill F_Y^{-1}(u)=-1-\frac{1}{\gamma }-\frac{1}{\gamma }{W}_{-1}[\frac{1+\gamma }{exp(1+\gamma )}(u-1)],0<u<1,\end{array}$

where $W_{-1}$ is the negative branch of the Lambert W function, which can be implemented in R (lamW package) or C (GNU Scientific Library).

On the other hand, if $Z\sim \mathrm{\text{NH}}(\alpha ,\lambda )$ has cdf $F_Z(z)$, the qf of $Z$ is

$\begin{array}{c}\hfill F_Z^{-1}(u)=\frac{1}{\lambda }\{[1-log(1-u){]}^{1/\alpha }-1\},0<u<1.\end{array}$

Then, we can generate an observation $x$ from the NHL distribution as follows:

1. Generate $u_1$ and $u_2$ independently from the uniform $𝒰(0,1)$ distribution;

2. Calculate $y=F_Y^{-1}(u_1)$ (observation from the Lindley distribution);

3. Calculate $z=F_Y^{-1}(u_2)$ (observation from the NH distribution);

4. Finally, calculate $x=min(y,z)$.

The number of Monte Carlo replications is set at $N=10,000$. We use the optim subroutine with analytical derivatives in R for maximizing the log-likelihood function. We compare the performance of the estimators by computing the mean estimates and root mean square errors (RMSEs) from $N$ Monte Carlo replications. We take the sample size as $n\in \{50,100,300,500\}$.

The simulation results are given in Table IV . As expected, under first-order asymptotic theory, the mean estimates of the parameters tend to be closer to the true values and the RMSEs decrease when the sample size $n$ increases for all scenarios. The MLEs of $\lambda$ and $\gamma$ are more accurate than those of $\alpha$. It is noteworthy that the estimates of $\alpha$ have considerable biases, specially in small samples. Nevertheless, we can build bias-corrected point estimator through bootstrap method, see(Davison 1997).

TABLE IV
Monte Carlo results: mean estimates and RMSEs.

Scenario	α	λ	γ	n	Mean			RMSE
					$\hat{\alpha }$	$\hat{\lambda }$	$\hat{\gamma }$	$\hat{\alpha }$	$\hat{\lambda }$	$\hat{\gamma }$
1	0.1	0.2	0.3	50	0.4136	0.1763	0.2710	0.7718	0.0435	0.0735
				100	0.2746	0.1842	0.2865	0.4211	0.0284	0.0520
				300	0.1684	0.1891	0.3006	0.1640	0.0127	0.0228
				500	0.1443	0.1896	0.3048	0.1048	0.0114	0.0188
2	0.5	1.8	1.1	50	0.5808	1.7829	1.1100	0.4232	0.0955	0.0360
				100	0.5474	1.7854	1.1107	0.1645	0.0757	0.0316
				300	0.5266	1.7886	1.1102	0.0418	0.0409	0.0196
				500	0.5212	1.7899	1.1098	0.0311	0.0101	0.0098
3	0.5	2.5	0.1	50	0.5361	2.4629	0.1033	0.0620	0.1472	0.0292
				100	0.5248	2.4749	0.1065	0.0387	0.0954	0.0171
				300	0.5159	2.4844	0.1086	0.0206	0.0411	0.0117
				500	0.5132	2.4863	0.1093	0.0158	0.0323	0.0110
4	0.7	2.8	0.1	50	0.8187	2.5700	0.1540	0.3410	0.5953	0.2359
				100	0.7586	2.6585	0.1259	0.1346	0.4069	0.1371
				300	0.7246	2.7570	0.1121	0.0427	0.1522	0.0578
				500	0.7176	2.7781	0.1101	0.0245	0.0646	0.0228
5	0.7	4	0.5	50	0.8298	3.7395	0.4800	0.5494	0.7593	0.1356
				100	0.7663	3.8218	0.4940	0.1629	0.5583	0.0774
				300	0.7340	3.8921	0.5151	0.0697	0.3781	0.0569
				500	0.7229	3.9379	0.5139	0.0438	0.2558	0.0476
6	0.8	0.4	0.1	50	0.9127	0.3675	0.1323	0.5489	0.0798	0.0737
				100	0.8593	0.3785	0.1209	0.33249	0.04984	0.04455
				300	0.8256	0.3859	0.1143	0.03713	0.02524	0.02618
				500	0.8200	0.3882	0.1115	0.0268	0.0191	0.0204
7	0.9	0.5	0.1	50	1.0241	0.4625	0.1387	0.4564	0.0931	0.0938
				100	0.9832	0.4719	0.1231	0.3026	0.0707	0.0600
				300	0.9294	0.4854	0.1154	0.0416	0.0251	0.0322
				500	0.9238	0.4868	0.1129	0.0401	0.0216	0.0268
8	1.3	0.4	0.1	50	1.5905	0.3703	0.1205	0.9387	0.0926	0.0671
				100	1.4492	0.3801	0.1244	0.7274	0.0602	0.0652
				300	1.3576	0.3860	0.1206	0.2268	0.0326	0.0475
				500	1.3351	0.3895	0.1154	0.0685	0.0129	0.0225

The non-parametric bootstrap procedure for bias correction applied in this paper can be described as follows:

1. Letting $x=x_1,\dots ,x_n$ be a observed random sample, we obtain the MLE $\widehat{𝛉}$ for the NHL parameter vector $𝛉$;

2. From the original data, generate a bootstrap sample $x^{*b}=(x_1^{*b},\dots ,x_n^{*b})$, that is, take a size $n$ random sample from $x$ with replacement;

3. Calculate the MLE of $𝛉$ from the bootstrap sample $x^{*b}$, namely ${\widehat{𝛉}}^{*}$;

4. Repeat the steps (2) and (3) for a very large number $B$ of times, thus obtaining ${\widehat{𝛉}}^{*1},\dots ,{\widehat{𝛉}}^{*𝐁}$;

5. Compute the bias estimate by $B(\widehat{𝛉})=\left[\frac{1}{B}\sum _{b=1}^B{\widehat{𝛉}}^{*𝐛}\right]-𝛉$;

6. The bias corrected estimate of $𝛉$ is given by $\overline{𝛉}=\widehat{𝛉}-B(\widehat{𝛉})$.

In order to evaluate the effectiveness of such correction, we carry out a separate simulation experiment. The numbers of Monte Carlo and bootstrap replications are equal to $1,000$ and $n$ varies in $\{20,30,50,70,100\}$. We study the performance of the bootstrap bias corrected estimates by computing the total relative bias defined in(Cribari-Neto and Soares 2003), which is a measure given by the sum of the absolute values of the individual relative biases. Thus, the total relative bias is an aggregate measure of the biases of the parameters estimates.

The plots in Figure5 show the behavior of the total relative bias for the uncorrected and corrected estimators. They reveal that the corrected estimators are more reliable than the MLEs in small sample sizes. In fact, the corrected estimators outperform the usual MLEs for both scenarios. In addition, it is also observed that the asymptotic property of unbiasedness of both estimators is satisfied, because their biases decrease when the sample size increases.

Figure 5
Total relative bias of the estimator for (a) Scenario 2 and (b) Scenario 5.

7 - APPLICATIONS

In this sec, we fit the NHL distribution to two real data sets. They illustrate the potentiality of this distribution for modeling positive data. The first data set represents the times to reinfection of sexually transmitted diseases (STDs) for eight hundred and seventy seven patients. The data are taken from sec 1.12 of(Klein and Moeschberger 1997). The second data set corresponds to the exceedances of flood peaks (in $m^3/s$) of the Wheaton River near Carcross in Yukon Territory, Canada. The data consist of 72 exceedances for the years $1958-1984$, rounded to one decimal place, see (Choulakian and Stephens 2001). Table5 provides a descriptive summary for both data sets.

TABLE V
Descriptive statistics.

Statistics	Real data sets
	Reinfection times	Exceedances of flood peaks
Mean	369.5268	12.2042
Median	247.0000	9.5000
Mode	5.000	5.0000
Variance	136940.7000	151.2215
Skewness	1.15262	1.4725
Kurtosis	0.51036	2.8895
Maximum	1529.0000	64.0000
Minimum	1.0000	0.1000
$n$	877	72

Note that the first data set has large amplitude and variance. Its measures of central tendency, such as mean, median and mode, are quite distant when compared among them. Besides, both data sets present positive values for the skewness and kurtosis. The exceedances of flood peaks data also present large amplitude and variance, but their values are lower than those for reinfection times.

For modeling these data sets, we fit the NHL distribution and also considered the fits of six related distributions: the ENH introduced by(Lemonte 2013) with pdf

$f(x)=\alpha \beta \lambda \frac{(1+\lambda x{)}^{\alpha -1}exp\{1-(1+\lambda x{)}^{\alpha }\}}{{[1-exp\{1-(1+\lambda x{)}^{\alpha }\}]}^{1-\beta }},x>0,$

where $\alpha >0$ and $\beta >0$ are shape parameters and $\lambda >0$ is a scale parameter; the exponentiated Weibull (EW), whose pdf is

$f(x)=\alpha \beta \lambda x^{\alpha -1}exp(-\lambda x^{\alpha })[1-exp(-\lambda x^{\alpha }){]}^{\beta -1},x>0,$

where $\alpha >0$ and $\beta >0$ are shape parameters and $\lambda >0$ is a scale parameter; the Weibull model arises from the EW model when $\beta =1$; and the NH, Lindley and exponential distributions, which are NHL special models.

We estimate the NHL parameters and the parameters for the above models by maximum likelihood. Based on the results in sec 6, the MLEs obtained for the exceedances of flood peaks data are obtained from the bootstrap bias corrected estimators, with $B=1,000$. The goodness-of-fit statistics considered are: Akaike information criteria (AIC), consistent Akaike information criteria (CAIC), Kolmogorov-Smirnov (KS), Cramér-von Mises ($W^{*}$) and Anderson-Darling ($A^{*}$). The lower are these statistics, the better is the adjustment to the data. The MLEs and goodness-of-fit statistics are calculated using the AdequacyModel script in the R software.

Tables VI and VII list the MLEs (and the corresponding standard errors in parentheses) of the unknown parameters for the fitted models to the first and second data sets, respectively. We note that all distributions present reasonable estimates for the standard errors. Tables VII and IX present the goodness-of-fit statistics for reinfection times and flood peaks, respectively. According to these statistics, the NHL distribution provides a good fit and is quite competitive with the other current distributions.

TABLE VI
MLEs for reinfection data and corresponding standard errors in parentheses.

Distributions	Estimates
NHL($\alpha ,\lambda ,\gamma$)	0.1247	0.0718	0.0031
	(0.0121)	(0.0130)	(0.0001)
ENH($\alpha ,\lambda ,\beta$)	0.8445	0.0030	0.7908
	(0.0347)	(0.0001)	(0.0324)
EW($\alpha ,\lambda ,\beta$)	0.9819	0.0026	0.7447
	(0.0100)	(0.0001)	(0.0293)
Weibull($\alpha ,\lambda$)	0.8464	0.0073
	(0.0206)	(0.0009)
NH($\alpha ,\lambda$)	0.7950	0.0039
	(0.0363)	(0.0003)
Lindley($\gamma$)	0.0053
	(0.0001)
Exp($\lambda$)	0.0026
	(7.568e-05)

TABLE VII
Goodness-of-fit statistics for the models fitted to reinfection times.

	Statistics
Distributions	AIC	CAIC	$W^{*}$	$A^{*}$	KS
NHL	12035.71	12035.74	0.0904	0.9826	0.0332
ENH	12081.29	12081.32	0.4790	3.6459	0.0477
EW	12069.46	12069.49	0.3422	2.7487	0.0576
Weibull	12084.77	12084.78	0.5082	3.8380	0.0552
NH	12121.64	12121.65	0.5782	4.2894	0.0677
Lindley	12828.07	12828.08	0.3532	2.8802	0.2076
Exp	12126.04	12126.04	0.3572	2.8505	0.0844

TABLE VIII
MLEs for exceedances of flood peaks data and corresponding standard errors in parentheses.

Distributions	Estimates
NHL($\alpha ,\lambda ,\gamma$)	0.1693	0.9332	0.0898
	(0.0361)	(0.0190)	(0.0064)
ENH($\alpha ,\lambda ,\beta$)	1.6963	0.0313	0.7302
	(1.1459)	(0.0307)	(0.1353)
EW($\alpha ,\lambda ,\beta$)	1.3392	0.0190	0.5445
	(0.4069)	(0.0303)	(0.2400)
Weibull($\alpha ,\lambda$)	0.9019	0.1095
	(0.0858)	(0.0302)
NH($\alpha ,\lambda$)	0.8464	0.1079
	(0.2493)	(0.0558)
Lindley($\gamma$)	0.1530
	(0.0128)
Exp($\lambda$)	0.0819
	(0.0096)

TABLE IX
Goodness-of-fit statistics for the models fitted to the exceedances of flood peaks.

	Statistics
Distributions	AIC	CAIC	$W^{*}$	$A^{*}$	KS
NHL	505.93	506.28	0.0897	0.5286	0.0979
ENH	507.84	508.20	0.1026	0.6275	0.1066
EW	508.05	508.40	0.1071	0.6504	0.1072
Weibull	506.99	507.17	0.1379	0.7851	0.1052
NH	507.97	508.15	0.1437	0.8149	0.1258
Lindley	530.42	530.48	0.1391	0.8526	0.2409
Exp	506.25	506.31	0.1306	0.7523	0.1422

For both data sets, the NHL distribution yields the best fit under all goodness-of-fit statistics. Figure6 displays the histograms and the estimated densities for three competitive models according to their goodness-of-fit statistics. This graphical inspection also indicates the superiority of the new distribution to the reinfection times and flood peaks data. We also provide the TTT plots and estimated hazard functions of the NHL fitted model for the reinfection times and flood peaks data in Figures7 and8, respectively. They reveal that the NHL hazard function for both fitted models present bathtub shapes, which is in agreement with their corresponding TTT plots. Figure9 displays the plots of the estimated cdfs for the most competitive models and empirical cumulative function to both data sets. These plots illustrate the good adjustment of the NHL distribution. Hence, these results reveal that the proposed distribution can be a very effective alternative to the well-known Weibull, EW and ENH distributions, among others.

Figure 6
Histogram and estimated pdfs of (a) the NHL, ENH and EW models for the reinfection times; and (b) the NHL, ENH and Exp models for the exceedances of flood peaks.

Figure 7
(a) TTT plot; and (b) estimated hazard function for the reinfection times data.

Figure 8
(a) TTT plot; and (b) estimated hazard function for the exceedances of flood peaks data.

Figure 9
Empirical and estimated cdfs of (a) the NHL, ENH and EW models for the reinfection times; and (b) the NHL, ENH and Exp models for the exceedances of flood peaks.

8 - CONCLUSIONS

We introduce the Nadarajah-Haghighi Lindley (NHL) model by compounding the Lindley and Nadarajah-Haghighi distributions. Once we have a composition by taking the minimum of two continuous independent random variables, the proposed distribution might be useful in engineering for modeling the failure time of systems composed of two independent components in series. The NHL distribution has as special cases the Lindley, Nadarajah-Haghighi and exponential distributions. Further, it is a competitive model to the Weibull, exponentiated Weibull and exponentiated Nadarajah-Haghighi distributions, among others. We obtain some structural properties of the proposed distribution, perform the estimation of the parameters by maximum likelihood and provide two applications to real data sets. The new distribution is quite competitive to other classical lifetime models and yields good adjustments in both applications.

AUTHOR CONTRIBUTIONS

The participation of the authors in the production of the manuscript is as follows: first author, characterization of the new distribution, mathematical properties and implementation of computational routines. The second author, application, simulation studies, and computational routines. The third author, review and general correction of the paper.

10 - APPENDIX A Scripts for approximating m_x by the truncated expansion (4)

Scripts for approximating μ^'_n by the truncated expansion (5).

Scripts for approximating m_s(y) by the truncated expansion (6)

REFERENCES

• ADAMIDIS K AND LOUKAS S. 1998. A lifetime distribution with decreasing failure rate. Stat Probabil Lett 39: 35-42
• AL-MUTAIRI DK, GHITANY ME AND KUNDU D. 2015. Inferences on Stress-Strength Reliability from Weighted Lindley Distributions. Commun Stat Theory Methods 44: 4096-4113.
• ASGHARZADEH A, BAKOUCH HS AND ESMAEILI L. 2013. Pareto Poisson–Lindley distribution with applications. J Appl Stat 40: 1717-1734.
• ASHOUR SK AND ELTEHIWY MA. 2015. Exponentiated power Lindley distribution. J Adv Res 6: 895-905.
• BAKOUCH HS, AL-ZAHRANI BM, AL-SHOMRANI AA, MARCHI VA AND LOUZADA F. 2012. An extended Lindley distribution. J Korean Stat Soc 41: 75-85.
• BOURGUIGNON M, LIMA MDCS, LEÃO J, NASCIMENTO ADC, PINHO LGB AND CORDEIRO GM. 2015. A new generalized gamma distribution with applications. Am J Math 34: 309-342.
• CHOULAKIAN V AND STEPHENS M. 2001. Goodness-of-fit for the generalized Pareto distribution. Technometrics 43: 478-484
• CORDEIRO G, ORTEGA E AND LEMONTE AJ. 2014. The exponential–Weibull lifetime distribution. J Stat Comput Sim 47:626-653.
• CRIBARI-NETO F AND SOARES ACN. 2003. Inferencia em modelos heterocedásticos. Rev Bras Econ 57: 319-335.
• DAVISON AC. 1997. Bootstrap methods and their application. Cambridge university press.
• GHITANY M, AL-MUTAIRI D AND NADARAJAH S. 2008a. Zero-truncated Poisson–Lindley distribution and its application. Math Comput Simulat 79: 279-287.
• GHITANY M, ALQALLAF F, AL-MUTAIRI D AND HUSAIN H. 2011. A two-parameter weighted Lindley distribution and its applications to survival data. Math Comput Simulat 81: 1190-1201.
• GHITANY M, ATIEH B AND NADARAJAH S. 2008b. Lindley distribution and its application. Math Comput Simulat 78: 493-506.
• GHITANY ME, AL-MUTAIRI DK, BALAKRISHNAN N AND AL-ENEZI LJ. 2013. Power Lindley Distribution and Associated Inference. Comput Stat Data An 64: 20-33.
• GUPTA PK AND SINGH B. 2013. Parameter estimation of Lindley distribution with hybrid censored data. Int J Syst Assur Eng Manag 4: 378-385.
• JODRÁ P. 2010. Computer generation of random variables with Lindley or Poisson-Lindley distribution via the Lambert W function. Math Comput Simulat 81: 851-859.
• KLEIN J AND MOESCHBERGER M. 1997. Survival Analysis Techniques for Censored and Truncated Data. New York: Springer Verlag.
• KRISHNA H AND KUMAR K. 2011. Reliability estimation in Lindley distribution with progressively type II right censored sample. Math Comput Simulat 82: 281-294.
• LEMONTE AJ. 2013. A new exponential-type distribution with constant, decreasing, increasing, upside-down bathtub and bathtub-shaped failure rate function. Comput Stat Data An 62: 149-170.
• LEMONTE AJ, CORDEIRO GM AND ORTEGA EMM. 2014. On the Additive Weibull Distribution. Commun Stat Theory Methods 43: 2066-2080.
• LINDLEY D. 1958. Fiducial distributions and Bayes’ theorem. J R Stat Soc A 20: 102-107.
• MACGILLIVRAY H. 1986. Skewness and asymmetry: measures and orderings. Anns of Stat 14: 994-1011.
• MARINHO PRD. 2016. Some new families of continuous distributions. Master’s thesis. Doctoral thesis, Universidade Federal de Pernambuco.
• MAZUCHELI J AND ACHCAR J. 2011. The Lindley distribution applied to competing risks lifetime data. Comput Meth Prog Bio 104: 188-192.
• NADARAJAH S, BAKOUCH HS AND TAHMASBI R. 2011. A generalized Lindley distribution. Sankhya B 73: 331-359.
• NADARAJAH S AND HAGHIGHI F. 2011. An extension of the exponential distribution. Statistics 45: 543-558.
• ORTEGA EM, LEMONTE AJ, SILVA GO AND CORDEIRO GM. 2015. New flexible models generated by gamma random variables for lifetime modeling. J Appl Stat 42: 2159-2179.
• POPOVÍC B, RISTÍC M AND CORDEIRO G. 2015. A two parameter distribution obtained by compounding the generalized exponential and exponential distributions. Mediterr J Math 15: 2935-2949.
• SANKARAN M. 1970. The discrete Poisson–Lindley distribution. Biometrics 26: 145-149.
• TAHIR MH, CORDEIRO GM, ALZAATREH A, MANSOOR M AND ZUBAIR M. 2016. The logistic-X family of distributions and its applications. Commun Stat Theory Methods 45: 7326-7349.
• XIE M AND LAI C. 1995. Reliability analysis using an additive Weibull model with bathtub shaped failure rate function. Reliab Eng Syst Safe 52: 87-93.
• ZAMANI H AND ISMAIL N. 2010. Negative binomial-Lindley distribution and its application. J Math Stat 6: 4-9.

Publication Dates

• Publication in this collection
  08 Apr 2019
• Date of issue
  2019

History

• Received
  26 Oct 2017
• Accepted
  9 July 2018