In this paper a new extension of the Lindley distribution is presented using the Weibull link function introduced and studied by Tahir et al.,¹ to develop a Weibull-Lindley distribution. We derive and discuss the mathematical and Statistical properties of the subject distribution along with its reliability analysis and inference for the parameters. Finally, the Weibull-Lindley distribution has been used to model four lifetime datasets and the results show that the proposed generalization performs better than the other known extensions of the Lindley distribution considered for the study.

Keywords: Lindley distribution, Weibull-Lindley distribution, mathematical properties, reliability function, parameter estimation, applications.

The Lindley distribution introduced by Lindley et al.,² in the context of Bayesian analysis as a counter example of fiducial statistics, is defined by its probability density function (PDF) and cumulative distribution function (CDF) as

$G\left(x\right)=1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}$ (1.1)

And

$g\left(x\right)=\frac{{\theta }^2}{\theta +1}\left(1+x\right){{\displaystyle e}}^{-\theta x}$ (1.2)

respectively. For $x>0,\theta >0,$ where $\theta$ is the scale parameter of the Lindley distribution.

Details of this distribution, its mathematical and statistical properties, estimation of its parameter and application including the superiority of Lindley distribution over exponential distribution has been done by Ghitany et al.,³ We have so many generalized families of distributions proposed by different researchers that are used in extending other distributions to produce compound distributions with better performance. These are several ways of adding one or more parameters to a distribution function which makes the resulting distribution richer and more flexible for modeling data. A brief summary of some of these methods or families of distribution include the beta generalized family (Beta-G) by Eugene et al.,⁴ the Kumaraswamy-G by Cordeiro et al.⁵  Transmuted family of distributions by Shaw et al.,⁶ Gamma-G (type 1) by Zografos et al.,⁷ McDonald-G by Alexander et al.,⁸ Gamma-G (type 2) by Risti et al.,⁹ Gamma-G (type 3) by Torabi et al.,¹⁰ Log-gamma-G by Amini et al.,⁹ Exponentiated T-X by Alzaghal et al.,¹² Exponentiated-G (EG) by Cordeiro et al.,¹³ Logistic-G by Torabi et al.,¹⁴ Gamma-X by Alzaatreh et al.,¹⁵ Logistic-X by Tahir et al.,¹⁶ Weibull-X by Alzaatreh et al.,¹⁷ Weibull-G by Bourguignon et al.,¹⁸ a new Weibull-G family by Tahir et al.,¹ a Lomax-G family by Cordeiro et al.,¹⁹ a new generalized Weibull-G family by Cordeiro et al.,²⁰ and Beta Marshall-Olkin family of distributions by Alizadeh et al.,²¹ and some other families of the distributions.

Hence, there are also some generalizations of the Lindley distribution recently proposed in the literature such as the transmuted Lindley distribution by Merovci et al.,²³ the exponentiated Power Lindley distribution by Ashour et al.,²⁴ Generalized Lindley distribution by Nadarajah et al.,²⁴ Transmuted Generalized Lindley distribution by Elgarhy et al.,²⁵ Extended Power Lindley distribution by Alkarni et al.,²⁶ a two-parameter Lindley distribution by Shanker et al.,²⁷ the Lomax-Lindley distribution by Yahaya et al.,²⁸ Transmuted Two-Parameter Lindley distribution by Al-khazaleh et al.,²⁹ and a three-parameter Lindley distribution by Shanker et al.,³⁰ The aim of this article is to introduce a new continuous distribution called Weibull-Lindley distribution (WLnD) from the proposed family by Tahir et al.,¹ The remaining parts of this article are presented in sections as follows: We defined the new distribution and give its plots in section 2.1. Section 2.2 derived some properties of the new distribution. Section 2.3 proposes some reliability functions of the new distribution. The order statistics for the new distribution are also given in section 2.4. The maximum likelihood estimates (MLEs) of the unknown model parameters of the new distribution are obtained in section 2.5. In section 3 we carryout application of the proposed model with others to four lifetime datasets. Lastly, in section 4, we give the summary of our work and concluding remarks.

Construction of Weibull-Lindley distribution (WLnD)

In the next section, we have defined the cdf and pdf of the Weibull-Lindley distribution (WLnD) using the method proposed by Tahir et al.¹ According to Tahir et al.,¹  the formula or Weibull link function for deriving the cdf and pdf of any Weibull-based continuous distribution is defined as:

$F(x)={\displaystyle \underset{0}{\overset{-\log \left[G(x)\right]}{\int }}\alpha \beta t^{\beta -1}{{\displaystyle e}}^{-\alpha t^{\beta }}dt}={{\displaystyle e}}^{-\alpha {\left\{-\log \left[G(x)\right]\right\}}^{\beta }}$  (2.1.1)

And

$f(x)=\alpha \beta \frac{g(x)}{G(x)}{\left\{-\log \left[G(x)\right]\right\}}^{\beta -1}{{\displaystyle e}}^{-\alpha {\left\{-\log \left[G(x)\right]\right\}}^{\beta }}$  (2.1.2)

respectively, where $g\left(x\right)$  and $G\left(x\right)$  are the pdf and cdf of any continuous distribution to be generalized respectively and $\alpha$  and $\beta$  are the two additional new parameters responsible for the shape of the distribution.

Using equation (1.1) and (1.2) in (2.1.1) and (2.1.2) and simplifying, we obtain the cdf and pdf of the Weibull-Lindley random variable X as:

$F\left(x;\theta ,\alpha ,\beta \right)={{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}$  (2.1.3)

And

$f(x;\theta ,\alpha ,\beta )=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}$  (2.1.4)

respectively.

For $x>0;\theta ,\alpha ,\beta >0$ ; where $\theta$  is a scale parameter and $\alpha$  and $\beta$  is a shape parameters of the Weibull-Lindley distribution.

The following is a graphical representation of the pdf and cdf of the Weibull-Lindley distribution. Given some values of the parameters $\alpha ,\beta \&\theta$ , we provide some possible graphs for the pdf and the cdf of the WLnD as shown in Figure 1&2 below: Figure 1 indicates that the WLnD is a skewed distribution and can take various forms. This means that distribution can be very useful for datasets that are skewed.

From the above cdf plot, the cdf increases when X increases, and approaches 1 when X becomes large, as expected.

In this section, we defined and discuss some properties of the WLnD distribution.

The Quantile function

This function is derived by inverting the cdf of any given continuous probability distribution. It is used for obtaining some moments like skewness and kurtosis as well as the median and for generation of random variables from the distribution in question. Hyndman et al.,²¹ defined the quantile function for any distribution in the form  Q(u) = $F^{-1}$ (u)  where Q(u) is the quantile function of F(x) for 0 < u <1

Taking F(x) to be the cdf of the Weibull-Lindley distribution and inverting it as above will give us the Quantile function as follows:

$F\left(x\right)={{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}=u$  (2.2.1)

Simplifying equation (2.2.1) above, we obtain:

  $Q\left(u\right)=X_q=-1-\frac{1}{\theta }-\frac{1}{\theta }W\left(-\frac{\left(\theta +1\right)\left(1-\exp \left\{-{\left(-\frac{\ln u}{\alpha }\right)}^{\frac{1}{\alpha }}\right\}\right)}{{{\displaystyle e}}^{\theta +1}}\right)$  (2.2.2)

By using (2.2.2) above, the median of X from the WLnD is simply obtained by setting u=0.5 while random numbers can be generated from WLnD by setting $X=Q\left(u\right)$ , where u is a uniform variate on the unit interval (0,1) and $W(.)$ represents the negative branch of the Lambert function.

Skewness and kurtosis

The quantile based measures of skewness and kurtosis will employed due to non-existence of the classical measures in some cases. The Bowley’s measure of skewness based on quartiles by Kenney et al.,³² is given as;

< $SK=\frac{Q\left(\frac{3}{4}\right)-2Q\left(\frac{1}{2}\right)+Q\left(\frac{1}{4}\right)}{Q\left(\frac{3}{4}\right)-Q\left(\frac{1}{4}\right)}$  (2.2.3)

while the Risti et al.,9 kurtosis based on octiles is given by;

$KT=\frac{Q\left(\frac{7}{8}\right)-Q\left(\frac{5}{8}\right)-Q\left(\frac{3}{8}\right)+\left(\frac{1}{8}\right)}{Q\left(\frac{6}{8}\right)-Q\left(\frac{1}{4}\right)}$  (2.2.4)

where Q(.) is any quartile or octile of interest.

Moments

Moments of a random variable are very important in distribution theory because they are used to study some of the most important features and characteristics of a random variable such as mean, variance, skewness and kurtosis.

Let X denote a continuous random variable, the n^th moment of X is given by;

${{\displaystyle \mu }}_n^{\text{'}}=E\left[{{\displaystyle X}}^n\right]={\displaystyle \underset{0}{\overset{\infty }{\int }}{{\displaystyle x}}^{n}f(x)dx}$ (2.2.5)

Considering f(x) to be the pdf of the Weibull-Lindley distribution as given in equation (2.1.4)

${{\displaystyle \mu }}_n^{\text{'}}=E\left[{{\displaystyle X}}^n\right]={\displaystyle \underset{0}{\overset{1}{\int }}{{\displaystyle x}}^{n}f(x)dx}$

Recall that from equation (2.1.4),

$f(x;\theta ,\alpha ,\beta )=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}$  (2.2.6)

Let

$A={{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}$

Then, using a power series expansion for A, we can write A as:

$A={{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}={\displaystyle \sum _{i=0}^{\infty }\frac{{\left(-1\right)}^i{\alpha }^i}{i!}}{\left(-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right)}^{\beta i}$

Substituting for the expansion above in equation (2.2.6), we have;

$f(x)=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}{\displaystyle \sum _{i=0}^{\infty }\frac{{\left(-1\right)}^i{\alpha }^i}{i!}}{\left(-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right)}^{\beta i}$

$f\left(x\right)=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}{\displaystyle \sum _{i=0}^{\infty }\frac{{\left(-1\right)}^i{\alpha }^i}{i!}}{\left(-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right)}^{\beta \left(i+1\right)-1}$  (JKJK) (2.2.7)

Also, let

$B={\left(-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right)}^{\beta \left(i+1\right)-1}$

Now, considering the following formula from Tahir et al.,¹ which holds for B for i≥1, and then we can write B as follows:

${\left(-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right)}^{\beta \left(i+1\right)-1}={\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{j+k+l}\left(\beta \left(i+1\right)\right)}{\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}{\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}^l$  (2.2.8)

Where for (for j≥0) P_j,0=1 and (for k=1,2,…..)

$P_{j,k}=k^{-1}{\displaystyle \sum _{m=1}^k{\left(-1\right)}^m\frac{\left[m\left(j+1\right)-k\right]}{\left(m+1\right)}{P}_{j,k-m}}$  (2.2.9)

Combining equation (2.2.8) and (2.2.9) and inserting the above power series in equation (2.2.7) and simplifying, we have: 

$f\left(x\right)=\frac{\alpha \beta {\theta }^2}{\left(\theta +1\right)}{\displaystyle \sum _{i=0}^{\infty }\frac{{\left(-1\right)}^i{\alpha }^i}{i!}}{\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{j+k+l}\left(\beta \left(i+1\right)\right)}{\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}^{l-1}$

$f\left(x\right)=\frac{\alpha \beta {\theta }^2}{\left(\theta +1\right)}{\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{i+j+k+l}{\alpha }^{i+1}\left(\beta \left(i+1\right)\right)}{i!\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}^{l-1}}$ (2.2.10)

Now,if l is a positive non-integer, we can expand the last term in (2.2.10) as:

${\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}^{l-1}={\displaystyle \sum _{m=0}^{\infty }{\left(-1\right)}^m}\left(\begin{array}{c}l-1\\ m\end{array}\right){\left[\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}^m$  (2.2.11)

Therefore, f(x) becomes:

$f\left(x\right)=\frac{\alpha \beta {\theta }^2}{\left(\theta +1\right)}{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{i+j+k+l+m}{\alpha }^{i+1}\left(\beta \left(i+1\right)\right)}{i!\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}l-1\\ m\end{array}\right)\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}\left(1+x\right){{\displaystyle e}}^{-\theta x\left(1+m\right)}}}{\left[1+\frac{\theta x}{\theta +1}\right]}^m$  (2.2.12)

Using power series expansion on the last term in equation (2.2.12), we have

${\left(1+\frac{\theta x}{\theta +1}\right)}^m={\displaystyle \sum _{p=0}^{\infty }\left(\begin{array}{c}m\\ p\end{array}\right)}{\left(\frac{\theta x}{\theta +1}\right)}^p={\displaystyle \sum _{p=0}^{\infty }\left(\begin{array}{c}m\\ p\end{array}\right)}{\left(\frac{\theta }{\theta +1}\right)}^p{x}^p$  (2.2.13)

Now, substituting equation (2.2.13), the power series expansion in equation (2.2.12) above, one gets:

$f\left(x\right)=\frac{\alpha \beta {\theta }^2}{\left(\theta +1\right)}{\displaystyle \sum _{p=0}^{\infty }{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{i+j+k+l+m}{\alpha }^{i+1}\left(\beta \left(i+1\right)\right)}{i!\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}l-1\\ m\end{array}\right)\left(\begin{array}{c}m\\ p\end{array}\right){\left(\frac{\theta }{\theta +1}\right)}^p\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}\left(1+x\right)x^p{{\displaystyle e}}^{-\theta x\left(1+m\right)}}}}$

This implies that:

$f(x)={\eta }_{i,j,k,l,m,p}\left(1+x\right)x^p{{\displaystyle e}}^{-\theta x\left(1+m\right)}$  (2.2.14)

Where

${\eta }_{i,j,k,l,m,p}=\frac{\alpha \beta {\theta }^2}{\left(\theta +1\right)}{\displaystyle \sum _{p=0}^{\infty }{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{k,l=0}^{\infty }{\displaystyle \sum _{j=0}^k\frac{{\left(-1\right)}^{i+j+k+l+m}{\alpha }^{i+1}\left(\beta \left(i+1\right)\right)}{i!\left(\beta \left(i+1\right)-1-j\right)}\left(\begin{array}{c}k-\left(\beta \left(i+1\right)-1\right)\\ k\end{array}\right)}}\left(\begin{array}{c}l-1\\ m\end{array}\right)\left(\begin{array}{c}m\\ p\end{array}\right){\left(\frac{\theta }{\theta +1}\right)}^p\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}\left(\beta \left(i+1\right)-1\right)+k\\ l\end{array}\right)P_{j,k}}}}$

Hence,

${{\displaystyle \mu }}_n^{\text{'}}=E\left[{{\displaystyle X}}^n\right]={\displaystyle \underset{0}{\overset{\infty }{\int }}{{\displaystyle x}}^{n}f(x)dx}={\displaystyle \underset{0}{\overset{\infty }{\int }}{\eta }_{i,j,k,l,m,p}\left(1+x\right)x^{n+p}{{\displaystyle e}}^{-\theta x\left(1+m\right)}dx}$

${{\displaystyle \mu }}_n^{\text{'}}=E\left[{{\displaystyle X}}^n\right]={\displaystyle \underset{0}{\overset{\infty }{\int }}{{\displaystyle x}}^{n}f(x)dx}={\eta }_{i,j,k,l,m,p}\left[{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^{n+p}{e}^{-\theta x\left(1+m\right)}dx+{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^{n+p+1}{e}^{-\theta x\left(1+m\right)}dx}}\right]$  (2.2.15)

Also, using integration by substitution method in equation (2.2.15); we obtain the following:

 Let $u=\theta x\left(1+m\right)\Rightarrow x=\frac{u}{\theta \left(1+m\right)}$ ; $\frac{du}{dx}=\theta \left(1+m\right)$  and $dx=\frac{du}{\theta \left(1+m\right)}$

Substituting for $u$ , $x$ and $dx$  in equation (2.2.15) and simplifying; we have:

${{\displaystyle \mu }}_n^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left(\frac{u}{\theta \left(1+m\right)}\right)}^{n+p}{e}^{-u}\left(\frac{du}{\theta \left(1+m\right)}\right)+{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left(\frac{u}{\theta \left(1+m\right)}\right)}^{n+p+1}{e}^{-u}\left(\frac{du}{\theta \left(1+m\right)}\right)}}\right]$

${{\displaystyle \mu }}_n^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[\frac{1}{{\left(\theta \left(1+m\right)\right)}^{n+p+1}}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}^{n+p}{e}^{-u}du+\frac{1}{{\left(\theta \left(1+m\right)\right)}^{n+p+2}}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}^{n+p+1}{e}^{-u}du}}\right]$

${{\displaystyle \mu }}_n^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[\frac{1}{{\left(\theta \left(1+m\right)\right)}^{n+p+1}}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}^{n+p+1-1}{e}^{-u}du+\frac{1}{{\left(\theta \left(1+m\right)\right)}^{n+p+2}}{\displaystyle \underset{0}{\overset{\infty }{\int }}{u}^{n+p+2-1}{e}^{-u}du}}\right]$  (2.2.16)

Again recall that ${\displaystyle \underset{0}{\overset{\infty }{\int }}{t}^{k-1}{e}^{-t}dt=\Gamma \left(k\right)}$ and that ${\displaystyle \underset{0}{\overset{\infty }{\int }}{t}^k{e}^{-t}dt={\displaystyle \underset{0}{\overset{\infty }{\int }}{t}^{k+1-1}{e}^{-t}dt=\Gamma \left(k+1\right)}}$

Thus we obtain the n^th ordinary moment of X for the Weibull-Lindley distribution as follows:

${{\displaystyle \mu }}_n^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(n+p+1\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+1}}+\frac{\Gamma \left(n+p+2\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+2}}\right]$  (2.2.17)

The Mean

The mean of the WLnD can be obtained from the n^th moment of the distribution when n=1 as follows:

${{\displaystyle \mu }}_n^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(n+p+1\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+1}}+\frac{\Gamma \left(n+p+2\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+2}}\right]$

$E\left[X\right]={{\displaystyle \mu }}_1^{\text{'}}={\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(p+2\right)}{{\left(\theta \left(1+m\right)\right)}^{p+2}}+\frac{\Gamma \left(p+3\right)}{{\left(\theta \left(1+m\right)\right)}^{p+3}}\right]$  (2.2.17)

The Variance

The n^th central moment or moment about the mean of X, say ${\mu }_n$ , can be obtained as

${{\displaystyle \mu }}_n=E{\left[X-{\mu }_1^{\text{'}}\right]}^n={\displaystyle \sum _{i=0}^n{(-1)}^i\left(\begin{array}{c}n\\ i\end{array}\right){\mu }_1^{\text{'}i}{\mu }_{n-i}^{\text{'}}}$  (2.2.18)

The variance of X for WLnD is obtained from the n^th central moment when n=2, that is, the variance of X is the n^th central moment of order two (n=2) and is given as follows:

$Var(X)=E[X^2]-{\left\{E[X]\right\}}^2$ (2.2.19)

$Var(X)={\mu }_2{}^{\text{'}}-{\left\{{\mu }_1{}^{\text{'}}\right\}}^2$

$Var(X)={\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(p+3\right)}{{\left(\theta \left(1+m\right)\right)}^{p+3}}+\frac{\Gamma \left(p+4\right)}{{\left(\theta \left(1+m\right)\right)}^{p+4}}\right]-{\left\{{\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(p+2\right)}{{\left(\theta \left(1+m\right)\right)}^{p+2}}+\frac{\Gamma \left(p+3\right)}{{\left(\theta \left(1+m\right)\right)}^{p+3}}\right]\right\}}^2$  (2.2.20)

The coefficients variation, skewness and kurtosis measures can also be calculated from the non-central moments using some well-known relationships.

Moment generating function

The mgf of a random variable X can be obtained by

${{\displaystyle M}}_x(t)=E\left[{{\displaystyle e}}^{tx}\right]={\displaystyle \underset{0}{\overset{\infty }{\int }}{{\displaystyle e}}^{tx}f(x)dx}$  (2.2.21)

Using power series expansion in equation (2.2.21) and simplifying the integral therefore we have;

$M_x(t)={\displaystyle \sum _{n=0}^{\infty }\frac{t^n}{n!}{\mu }_n^{\text{'}}}={\displaystyle \sum _{n=0}^{\infty }\frac{t^n}{n!}\left\{{\eta }_{i,j,k,l,m,p}\left[\frac{\Gamma \left(n+p+1\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+1}}+\frac{\Gamma \left(n+p+2\right)}{{\left(\theta \left(1+m\right)\right)}^{n+p+2}}\right]\right\}}$  (2.2.22)

where n and t are constants, t is a real number and ${\mu }_n^{\text{'}}$  denotes the n^th ordinary moment of X .

Characteristics function

The characteristics function of a random variable X is given by;

${\phi }_x(t)=E\left[e^{itx}\right]=E\left[\cos (tx)+i\sin (tx)\right]=E\left[\cos (tx)\right]+E\left[i\sin (tx)\right]$  (2.2.23)

Simple algebra and power series expansion proves that

${{\displaystyle \varphi }}_x(t)={\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{t}^{2n}}{\left(2n\right)!}{{\displaystyle \mu }}_{2n}^{\text{'}}}+i{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{t}^{2n+1}}{\left(2n+1\right)!}{{\displaystyle \mu }}_{2n+1}^{\text{'}}}$  (2.2.24)

Where ${\mu }_{2n}^{\text{'}}$  and ${\mu }_{2n+1}^{\text{'}}$ are the moments of X for n=2n and n=2n+1 respectively and can be obtained from ${\mu }_n^{\text{'}}$  in equation (2.2.17)

Some reliability functions

In this section, we present some reliability functions associated with WLnD including the survival and hazard functions.

The Survival function

The survival function describes the likelihood that a system or an individual will not fail after a given time. It tells us about the probability of success or survival of a given product or component. Mathematically, the survival function is given by:

$S\left(x\right)=1-F\left(x\right)$  (2.3.1)

Taking F(x) to be the cdf of the Weibull-Lindley distribution, substituting and simplifying (2.3.1) above, we get the survival function of the WLnD as:

$S\left(x\right)=1-{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}$  (2.3.2)

Below is a plot of the survival function at chosen parameter values in Figure 3. The figure above revealed that the probability of survival for any random variable following a Weibull-Lindley distribution decreases as the values of the random variable increases, that is, as time goes on, probability of life decreases. This implies that the Weibull-Lindley distribution can be used to model random variables whose survival rate decreases as their age grows.

The Hazard function

Hazard function as the name implies is also called risk function, it gives us the probability that a component will fail or die for an interval of time. The hazard function is defined mathematically as;

$h\left(x\right)=\frac{f\left(x\right)}{1-F\left(x\right)}=\frac{f\left(x\right)}{S\left(x\right)}$  (2.3.3)

Taking f(x) and F(x) to be the pdf and cdf of the proposed Weibull-Lindley distribution given previously, we obtain the hazard function as:

$h\left(x\right)=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\left(1-{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}\right)}$ (2.3.4)

The following is a plot of the hazard function at chosen parameter values in Figure 4.

Figure 4 above shows the behavior of hazard function of the WLnD and it means that the probability of failure for any Weibull-Lindley random variable increases as the time or age of the variable grows or increases, that is, as time goes, the probability of failure or death increases and becomes constant after some times.

Order statistics

Order statistics are used widely over the years for solving a huge set of problems such as in robust statistical estimation and detection of outliers, characterization of probability distributions and goodness of fit tests, entropy estimation, analyses of censored samples, reliability analysis, quality control and strength of materials. Suppose $X_1,X_2,\mathrm{......},X_n$  is a random sample from a distribution with pdf, f(x), and let $X_{1:n}<X_{2:n}<\mathrm{......}<X_{a:n}$  denote the corresponding order statistic obtained from this sample. The pdf, $f_{a:n}\left(x\right)$  of the $a^{th}$  order statistic can be defined as;

${{\displaystyle f}}_{a:n}\left(x\right)=\frac{n!}{\left(a-1\right)!\left(n-a\right)!}{\displaystyle \sum _{k=0}^{n-a}{{\displaystyle (-1)}}^k\left(\begin{array}{c}n-a\\ k\end{array}\right)}f\left(x\right){{\displaystyle F\left(x\right)}}^{k+a-1}$  (2.4.1)

where f(x) and F(x) are the pdf and cdf of the Weibull-Lindley distribution respectively.

Using (2.1.3) and (2.1.4), the pdf of the $a^{th}$  order statistics $X_{a:n}$ , can be expressed from (2.4.1) as;

${{\displaystyle f}}_{a:n}(x)=\frac{n!}{(a-1)!(n-a)!}{\displaystyle \sum _{k=0}^{n-a}{{\displaystyle (-1)}}^k\left(\begin{array}{c}n-a\\ k\end{array}\right)}\left[\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]{{\displaystyle e}}^{\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}}\right]*{{\displaystyle \left[{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}\right]}}^{a+k-1}$  (2.4.2)

Hence, the pdf of the minimum order statistic $X_{(1)}$  and maximum order statistic $X_{(n)}$  of the WLnD are given by;

${{\displaystyle f}}_{1:n}(x)=n{\displaystyle \sum _{k=0}^{n-1}{{\displaystyle (-1)}}^k\left(\begin{array}{c}n-1\\ k\end{array}\right)}\left[\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]{{\displaystyle e}}^{\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}}\right]*{{\displaystyle \left[{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}\right]}}^k$  (2.4.3)

And

${{\displaystyle f}}_{n:n}(x)=n\left[\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]{{\displaystyle e}}^{\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}}\right]{{\displaystyle \left[{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}\right]}}^{n-1}$  (2.4.4)

,p>respectively.

Parameter estimation via maximum likelihood

Let $X_1,X_2,\mathrm{....},X_n$ > be a sample of size ‘n’ independently and identically distributed random variables from the WLnD with unknown parameters α, β and Ө defined previously. The pdf of the WLnD is given from (2.1.3) as

$f\left(x\right)=\frac{\alpha \beta {\theta }^2\left(1+x\right){{\displaystyle e}}^{-\theta x}{\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta -1}{{\displaystyle e}}^{-\alpha {\left\{-\log \left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]\right\}}^{\beta }}}{\left(\theta +1\right)\left[1-\left[1+\frac{\theta x}{\theta +1}\right]{{\displaystyle e}}^{-\theta x}\right]}$

The likelihood function is given by;

$L\left(\underset{\_}{X}|\alpha ,\beta ,\theta \right)=\frac{{{\displaystyle \left(\alpha {\beta }^{\alpha }{\theta }^2\right)}}^n{\displaystyle \prod _{i=1}^n\left\{\left(1+x_i\right){{\displaystyle e}}^{-\theta x_i}\right\}}{\displaystyle \prod _{i=1}^n{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta -1}}}{{\left(\theta +1\right)}^n{\displaystyle \prod _{i=1}^n\left\{1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right\}{{\displaystyle e}}^{\alpha {\displaystyle \sum _{i=1}^n{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta }}}}}$  (2.5.1)

Let the log-likelihood function, $l=\log L\left(\underset{\_}{X}|\alpha ,\beta ,\theta \right)$ , therefore

$l=n\log \alpha +n\log \beta +2n\log \theta -n\log \left(\theta +1\right)+{\displaystyle \sum _{i=1}^n\log \left(1+x_i\right)}-\theta {\displaystyle \sum _{i=1}^n{x}_i}+\left(\beta -1\right){\displaystyle \sum _{i=1}^n\log \left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}$

$-{\displaystyle \sum _{i=1}^n\log \left(\left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right)}-\alpha {\displaystyle \sum _{i=1}^n{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta }}$  (2.5.2)

Differentiating $l$  partially with respect to α, β and Ө respectively gives;

$\frac{\partial l}{\partial \alpha }=\frac{n}{\alpha }-{\displaystyle \sum _{i=1}^n{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta }}$  (2.5.3)

$\frac{\partial l}{\partial \beta }=\frac{n}{\beta }+{\displaystyle \sum _{i=1}^n\log \left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}-\alpha {\displaystyle \sum _{i=1}^n{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta }\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}$  (2.5.4)

$\frac{\partial l}{\partial \theta }=\frac{2n}{\theta }-\frac{n}{\theta +1}-{\displaystyle \sum _{i=1}^n{x}_i}+\left(\beta -1\right){\displaystyle \sum _{i=1}^n\left\{\frac{x_i{{\displaystyle e}}^{-\theta x_i}\left[\frac{1}{{\left(\theta +1\right)}^2}-\frac{\theta x_i+\theta +1}{\theta +1}\right]}{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}\left\{\left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}\right\}}-{\displaystyle \sum _{i=1}^n\left\{\frac{x_i{{\displaystyle e}}^{-\theta x_i}\left[\frac{1}{{\left(\theta +1\right)}^2}-\frac{\theta x_i+\theta +1}{\theta +1}\right]}{\left(\left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right)}\right\}}$

$+\alpha \beta {\displaystyle \sum _{i=1}^n\left\{\frac{x_i{{\displaystyle e}}^{-\theta x_i}\left[\frac{1}{{\left(\theta +1\right)}^2}-\frac{\theta x_i+\theta +1}{\theta +1}\right]{\left\{-\log \left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]\right\}}^{\beta -1}}{\left[1-\left[1+\frac{\theta x_i}{\theta +1}\right]{{\displaystyle e}}^{-\theta x_i}\right]}\right\}}$  (2.5.5)

Equating equations (2.5.3), (2.5.4) and (2.5.5) to zero and solving for the solution of the non-linear system of equations will give us the maximum likelihood estimates of parameters $\alpha ,\beta \&\theta$  respectively. However, the above equations cannot be solved manually due to their complexity unless numerically with the help of statistically inclined computer programs like Python, R, SAS, etc., when data sets are available.

This section presents four datasets, their descriptive statistics, graphics and applications to some selected extensions of the Lindley distribution including the classical Lindley distribution. We have compared the performance of the Weibull-Lindley distribution (WLlD) to some families of Lindley distribution such as Lomax-Lindley distribution (LLlD), Two-parameter Lindley distribution (TPLlD), Transmuted Lindley distribution (TLlD) and the Lindley distribution (LlD).

 Data sets and their nature

In this section, four different datasets and their summary are presented for fitting the above listed distributions. The available data sets I, II, III, and IV and their respective summary statistics are provided in Table 1–5 respectively as follows;

Dataset I: This dataset represents the lifetime’s data relating to relief times (in minutes) of 20 patients receiving an analgesic and reported by Gross et al.,³⁴ and has been used by Shanker et al.,³⁵ Table 1.

Dataset II: This data represent the survival times in weeks for male rats Lawless et al.,³⁶ (Table 2).

Data set III: This data set represents the strength of 1.5cm glass fibers initially collected by members of staff at the UK national laboratory. It has been used by Bourguignon et al.,¹⁸ Afify et al.,³⁷  Barreto Souza et al.,³⁸ Oguntunde et al.,³⁹ Ieren  et al.,⁴⁰ as well as Smith et al.,⁴¹ Its summary is given as follows: (Table 3).

Dataset IV: This dataset represents 59 observations of the monthly actual taxes revenue in Egypt (in 1,000 million Egyptian pounds) between January 2006 and November 2010. The data has been previously used by Owoloko et al.⁴² The descriptive statistics for this data are as follows:

We also provide some histograms and densities for the three data sets as shown in Figure 5–8 below respectively (Table 4).

From the summary statistics of the four data sets, we found that data sets I and IV are positively skewed, while II is a bit negatively skewed or approximately normal and III is negatively skewed. Also, data sets I, III and IV have higher kurtosis while II have low level or degree of peakness.

Analysis of data

These four different datasets presented above were used to fit all the above listed Lindley distributions by applying the formulas of the test statistics in section 2 in order to get the best fitted model and the results are presented as follow in the four tables for each dataset below: (Table 5).

From Table 5, the values of the parameter MLEs and the corresponding values of ƖƖ, AIC, A^*, W^* and K-S for each model show that the Weibull-Lindley distribution (WLlD) has better performance compared to the other four models namely: Lomax-Lindley distribution (LLlD), Two-parameter Lindley distribution (TPLlD), Transmuted Lindley distribution (TLlD) and the Lindley distribution (LlD) and hence becomes the best fitted distribution based the data set I (Table 6).

Again the results in Table 6 above shows that the Weibull-Lindley distribution (WLlD) fits the second dataset better than the other four models (LLlD, TPLlD, TLlD and LlD) because the values of the statistics; AIC, A^*, W^* and K-S are smaller for the WLlD than the other models and therefore it is considered as the best fitted distribution based the data set II (Table 7).

The results from Table 7 also agrees with the previous results that the WLlD is more flexible compared to the three other models this also agrees with the fact that generalizing any continuous distribution provides a compound distribution with at least better fit than the classical distribution (i.e Lindley) irrespective of the nature of the data used provide it is asymmetry Table 8.

Lastly, our results in Table 8 provides the same results as obtained in the above previous tables with the Weibull-Lindley distribution performing better than the other three distributions considered in this study.

The following figures displayed the histogram and estimated densities of the fitted models for the four real life data sets used in this study.

From the estimated density plots in Figures 9 we can observe that though there is no big difference between the performance of the other four models, it is very clear that the performance of the Weibull-Lindley distribution (WLlD) remains the best and consistent irrespective of the nature the datasets as compared to the Lomax-Lindley distribution (LLlD), Two-parameter Lindley distribution (TPLlD), Transmuted Lindley distribution (TLlD) and the Lindley distribution (LlD).

Furthermore, the performance of the Weibull-Lindley could be attributed to the fact that the Weibull-Lindley distribution is heavy-tailed and highly skewed to the right with excellent flexibility which allows it to take various shapes depending on the parameter values and it also exhibit some degree of kurtosis all of which are features of the four datasets used in this research, hence, the Weibull-Lindley distribution will be more appropriate for lifetime datasets which are positively skewed with a higher degree of peakness as well as those that are approximately normal with observations above zero.

Hence, having demonstrated earlier in Tables 5–8, we have a similar conclusion based on figure 3-5 that the Weibull-Lindley distribution has a better fit for the four data sets considered in this study.