A Lindley-Type Distribution for Modeling High-Kurtosis Data

Abstract

This article proposes a heavy-tailed distribution for modeling positive data. The proposal arises with the ratio of independent random variables, specifically, a Lindley distribution divided by a beta distribution. This leads to a three-parameter extension of the Lindley distribution capable of modeling high levels of kurtosis. The main structural properties of the proposed distribution are derived. The skewness and kurtosis behavior of the distribution are described. Parameter estimation is discussed under consideration of the moment and maximum likelihood methods. Finally, in order to avoid the parameter non-identifiability problem, a two-parameter version of the proposed distribution is derived. The usefulness of this special case is illustrated by fitting data in two real scenarios.

1. Introduction

The Lindley (L) distribution [1,2] is one of the alternatives commonly used to model the lifetime of mechanical units. This distribution has various attractive properties, for example: (1) Although it has a single parameter, its probability density function (pdf) can present both monotonic decreasing and unimodal shapes; (2) like other classical distributions, such as Weibull and gamma, the hazard rate function (hrf) of the L distribution exhibits an increasing shape. However, unlike the Weibull and gamma hrf’s, the hrf of the L distribution is different from 0 at time $t=0$ and converges to a finite quantity (a function of the parameter) as $t\to \infty$; and (3) many of its properties have closed forms, which makes it easy to implement computationally.

A random variable Z has a L distribution with shape parameter $\theta >0$, denoted as $Z\sim \mathrm{L}\left(\theta \right)$, if its pdf is given by

$$\begin{array}{c}\hfill f_Z(z;\theta )=\frac{{\theta }^2}{(1+\theta )}(1+z)e^{-\theta z},z>0.\end{array}$$

(1)

Ghitany et al. [3] carry out a detailed study of the properties of the L distribution, among which it is shown that Equation (1) corresponds to the pdf of a mixture distribution with exponential and gamma components. Although the L distribution can model an important variety of data whose histograms exhibit decreasing or unimodal behavior, performance can be poor when the data exhibit high levels of skewness and/or kurtosis. In this sense, it is possible to find in the literature various extensions of the L distribution capable of capturing high levels of skewness and/or kurtosis. Good examples of this are the weighted-Lindley [4], generalized Lindley [5], power-Lindley [6], quasi-Lindley [7], and two-parameter Lindley [8,9] distributions, which are more flexible than the L distribution, especially in terms of skewness. The three-parameter extensions proposed by Zakerzadeh and Dolati [10], Shanker et al. [11], and Ekhosuehi and Opone [12] are more flexible than the L distribution in both skewness and kurtosis.

An extension of the L distribution oriented especially to model data with high kurtosis levels is the three-parameter Lindley slash (LS) distribution [13]. A random variable Y follows a LS distribution with scale parameter $\sigma >0$, shape parameter $\theta >0$ and kurtosis parameter $\alpha >0$, denoted as $X\sim \mathrm{LS}(\sigma ,\theta ,\alpha )$, if it can be represented as

$$\begin{array}{c}\hfill Y=\sigma \frac{Z}{U^{1/\alpha }},\end{array}$$

(2)

where $Z\sim \mathrm{L}\left(\theta \right)$ and $U\sim \mathrm{Uniform}(0,1)$ are independent.

From Equation (2), it follows that the LS distribution (when $\sigma =1$) converges to the L distribution as $\alpha \to \infty$. The pdf of the LS distribution is given by

$$\begin{array}{c}\hfill f_Y(y;\sigma ,\theta ,\alpha )=\frac{\alpha {\theta }^2}{\sigma (1+\theta )}{\int }_0^1\left(1+\frac{yu}{\sigma }\right)e^{-\frac{\theta yu}{\sigma }}{u}^{\alpha }du,y>0.\end{array}$$

The LS distribution has heavier tails than the L distribution, so it can be used to fit positive data with outliers.

Iriarte and Rojas [14], assuming in Equation (2) that that Z has a power-Lindley distribution, propose a new three-parameter heavy-tailed extension of the L distribution, the slashed-power Lindley distribution. In the literature, it is possible to find a large amount of information on slash-type distributions, which arise assuming certain distributions for the random variables involved in the Equation (2). For some details of slash-type distributions for positive data, see [15,16,17,18,19,20], among others.

In this article, we introduce a new heavy-tailed extension of the L distribution, the extended slash Lindley (ESL) distribution. Based on the results of Rojas et al. [21], the proposed distribution arises from Equation (2) considering $\sigma =\alpha =1$ and that U has a beta distribution with mean $\alpha /(\alpha +\beta )$, $\alpha ,\beta >0$. Thus, Equation (2) gives rise to a new distribution with one shape parameter controlling unimodality and two shape parameters controlling kurtosis.

The fact that the ESL distribution has two parameters controlling kurtosis allows this distribution to be able to capture extremely high levels of kurtosis. However, we observe that this same fact generates a problem of non-identifiability of the parameters. In this sense, we propose a two-parameter special case of the ESL distribution provided with one parameter controlling unimodality and one parameter controlling kurtosis. We emphasize that this special case has heavier tails than the L distribution, and even heavier than the LS distribution, so it can be considered as a viable alternative to model positive data with outliers.

The article is organized as follows. In Section 2, we propose the ESL distribution and derive some of its main structural properties. In Section 3, we discuss parameter estimation considering the moment and maximum likelihood methods. In addition, we carry out a simulation study to evaluate the behavior of the moment and maximum likelihood estimators under the two-parameter ESL distribution. Section 4 presents an application with real data in order to illustrate the usefulness of the proposal. Some final comments are considered in Section 5.

2. ESL Distribution

In this section, we propose the ESL distribution and derive its main structural properties.

2.1. Stochastic Representation

Definition 1.

The random variable X follows the extended slash Lindley distribution, denoted as $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$, if it can be represented as $X=Y/U,$ where $Y\sim \mathrm{L}\left(\theta \right)$ and $U\sim \mathrm{Beta}(\alpha ,\beta )$ are independent random variables.

From Definition 1, it is observed that: (1) If $\beta =1$, then U has a power distribution [22] on the interval $(0,1)$. Consequently, X has a nonscaled $(\sigma =1)$ LS distribution; (2) If $\beta =1$ and $\alpha \to \infty$, then the distribution of X converges to the L distribution.

The following proposition provides the pdf of the ESL distribution.

Proposition 1.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, the pdf of X is given by

$$\begin{array}{c}\hfill f_X(x;\theta ,\alpha ,\beta )=\frac{{\theta }^2}{(\theta +1)B(\alpha ,\beta )}{\int }_0^1(1+xw)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta xw}dw,x>0,\end{array}$$

(3)

where $\theta ,\alpha ,\beta >0$ are shape parameters and $B(a,b)={\int }_0^1{u}^{a-1}{(1-u)}^{b-1}du$ is the beta function.

Proof. 

From Definition 1, we use the Jacobian method to compute the pdf of X. Specifically, considering the transformations $X=Y{U}^{-1}$ and $W=U$, we have $Y=XW$ and $U=W$. Thus, the determinant of the Jacobian of the transformation is

$$\begin{array}{c}\hfill J=\left|\begin{array}{cc}\frac{\partial Y}{\partial X}\hfill & \frac{\partial Y}{\partial W}\hfill \\ \frac{\partial U}{\partial X}\hfill & \frac{\partial U}{\partial W}\hfill \end{array}\right|=w.\end{array}$$

Hence, the joint pdf of X and W is given by

$$\begin{array}{c}\hfill f_{X,W}(x,w)=f_Y(xw;\theta )f_U(w;\alpha ,\beta )w,x>0,0<w<1,\end{array}$$

(4)

where $f_Y(·;\theta )$ denotes the pdf of the L distribution with shape parameter $\theta >0$ and $f_U(·;\beta ,\alpha )$ denotes the pdf of the beta distribution with shape parameters $\alpha ,\beta >0$. Then, the marginal distribution of X is obtained directly by integrating Equation (4), which leads to Equation (3). □

The following result shows that the pdf of the ESL distribution can be written in closed form in terms of the gamma, beta and confluent hypergeometric functions. In this sense, we recall the integral representations of these functions given in Andrews [23]. For $\mathrm{R}\left(a\right)>0$ and $\mathrm{R}\left(b\right)>0$, the gamma and beta functions are given by $\Gamma \left(a\right)={\int }_0^{\infty }{u}^{a-1}{e}^{-u}du$ and $\mathrm{B}(a,b)={\int }_0^1{u}^{a-1}{(1-u)}^{b-1}du$, respectively. For $\mathrm{R}\left(b\right)>\mathrm{R}\left(a\right)>0$, the confluent hypergeometric function of the first kind is given by $\mathrm{M}(a,b,z)=\Gamma \left(b\right)/[\Gamma \left(a\right)\Gamma (b-a)]{\int }_0^1{e}^{zu}{u}^{a-1}{(1-u)}^{b-a-1}du.$ It should be noted that the gamma, beta and confluent hypergeometric functions are implemented in some computer algebra systems (such as MATLAB) and programming languages (such as R [24]). From these functions, we establish the following.

Corollary 1.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then,

$$\begin{array}{c}\hfill f_X(x;\theta ,\alpha ,\beta )=\frac{{\theta }^2\Gamma \left(\beta \right)}{(1+\theta )\mathrm{B}(\alpha ,\beta )}\left[r_1(x;\theta ,\alpha ,\beta )+x{r}_2(x;\theta ,\alpha ,\beta )\right],\end{array}$$

(5)

where

$$\begin{array}{c}\hfill r_j(x;\theta ,\alpha ,\beta )=\frac{\Gamma (\alpha +j)\mathrm{M}(\alpha +j,\alpha +\beta +j,-\theta x)}{\Gamma (\alpha +\beta +j)},j=1,2.\end{array}$$

Proof. 

We note that Equation (3) can be written as

$$\begin{array}{ccc}\hfill f_X(x;\theta ,\alpha ,\beta )& =& \frac{{\theta }^2}{(1+\theta )\mathrm{B}(\alpha ,\beta )}\left[{\int }_0^1{w}^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta xw}dw\right.\hfill \\ & & \left.+x{\int }_0^1{w}^{\alpha +1}{(1-w)}^{\beta -1}{e}^{-\theta xw}dw\right],\hfill \end{array}$$

and the result is obtained by recognizing the confluent hypergeometric function in the previous integrals. □

From Proposition 1 and Corollary 1, we see that

$$\begin{array}{ccc}\hfill f_X(x;\theta ,\alpha ,\beta =1)& =& \frac{\alpha {\theta }^2}{(1+\theta )}{\int }_0^1\left(1+\frac{xu}{\sigma }\right)e^{-\frac{\theta xu}{\sigma }}{u}^{\alpha }du\hfill \\ & =& \frac{\alpha {\theta }^2}{(1+\theta )}\left[r_1(x;\theta ,\alpha ,1)+x{r}_2(x;\theta ,\alpha ,1)\right],\hfill \end{array}$$

which corresponds to the pdf of the non-scaled version of the LS distribution. Note that the closed form above for the LS pdf is, to our knowledge, not known.

On the other hand, we observe that

$$\begin{array}{c}\hfill \underset{\alpha \to \infty }{lim}{f}_X(x;\theta ,\alpha ,\beta =1)=\frac{{\theta }^2}{(1+\theta )}(1+x)e^{-\theta x},\end{array}$$

which corresponds to the pdf of the L distribution.

Figure 1 shows some curves for the ESL pdf considering different values of $\theta$, $\alpha$ and $\beta$. We use the function kummerM [25] of the R programming language to compute the confluent hypergeometric function. We provide the R code in Appendix A.

In Figure 1, it can be seen that:

• The ESL pdf exhibits a unimodal shape when $\theta$ is equal to 0.5 (see top panels) and a monotonic decreasing shape when $\theta$ is equal to 5 (see bottom panels).
• If $\theta$ and $\beta$ are fixed (see left panels), it is observed that the unimodal and monotonic shapes of the pdf are preserved with each choice of $\alpha$. In this case, the different values of $\alpha$ define a set of pdfs with different weights in the tails. The closer to 0 the value of $\alpha$ is, the heavier the tail of the ESL pdf.
• For fixed values of $\theta$ and $\alpha$ (see right panels), something similar is observed. The unimodal and monotonic shapes of the pdf are preserved even when the parameter $\beta$ varies. Here, different values for $\beta$ also define a set of pdf’s with different weights in the tails. The larger the value of $\beta$, the heavier the tail of the ESL pdf.

The above observations suggest that the unimodal and monotone shapes of the ESL pdf depend exclusively on $\theta$, while $\alpha$ and $\beta$ control the kurtosis of the distribution. More details on the behavior of the kurtosis of the ESL distribution are given in Section 2.3.

2.2. Reliability Analysis

The reliability function (rf) and the hazard rate function (hrf) play a central role in the treatment of lifetime data in reliability studies. If X is a random variable representing the failure time of mechanical units, the rf of X, defined as $R_X\left(x\right):=\mathbb{P}(X>x)$, $x>0$, indicates the probability that mechanical units survive beyond the time x. On the other hand, the hrf of X, defined as $h_X\left(x\right):=f_X\left(x\right)/R_X\left(x\right)$ (where $f_X(·)$ represents the pdf of X), measures the propensity of a mechanical units to fail or die depending on the age it has reached.

In what follows, we concentrate on deriving analytic expressions and some graphical representations for the rf and the hrf of the ESL distribution. For this, we first calculate the cdf of the ESL distribution.

Proposition 2.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, the cdf of X is given by

$$\begin{array}{c}\hfill F_X(x;\theta ,\alpha ,\beta )=F_Y(x;\theta )-\frac{{\theta }^2}{(1+\theta )B(\alpha ,\beta )}E\left(x;\theta ,\alpha ,\beta \right),\end{array}$$

where $F_Y(·;·)$ is the cdf of the L distribution and

$$\begin{array}{c}\hfill E(x;\theta ,\alpha ,\beta )={\int }_0^x(1+v)B\left(\frac{v}{x};\alpha ,\beta \right)e^{-\theta v}dv.\end{array}$$

such that $B\left(z;a,b\right)={\int }_0^z{u}^{a-1}{(1-u)}^{b-1}du$ is the incomplete beta function.

Proof. 

Considering the change of variable $v=xw$, the pdf of X given in the Equation (3) can be written as $f_X(x;\theta ,\alpha ,\beta )={\theta }^2/\left[(1+\theta )B(\alpha ,\beta )\right]{\int }_0^x(1+v)v^{\alpha }{\left(1-v/x\right)}^{\beta -1}dv$. Then, the cdf of X is given by

$$\begin{array}{c}\hfill F_X(x;\theta ,\alpha ,\beta )=\frac{{\theta }^2}{(1+\theta )B(\alpha ,\beta )}{\int }_0^x\frac{1}{u^{\alpha +1}}{\int }_0^u(1+v)v^{\alpha }{\left(1-\frac{v}{u}\right)}^{\beta -1}{e}^{-\theta v}dvdu.\end{array}$$

Now, exchanging the order of integration in the previous expression, we obtain that

$$\begin{array}{ccc}\hfill F_X(x;\theta ,\alpha ,\beta )& =& \frac{{\theta }^2}{(1+\theta )B(\alpha ,\beta )}{\int }_0^x\frac{1+v}{v}{e}^{-\theta v}{\int }_v^x{\left(\frac{v}{u}\right)}^{\alpha +1}{\left(1-\frac{v}{u}\right)}^{\beta -1}dudv\hfill \\ & =& 1-\frac{(1+\theta +\theta x)e^{-\theta x}}{(1+\theta )}-\frac{{\theta }^2}{(1+\theta )B(\alpha ,\beta )}{\int }_0^x(1+v)B\left(\frac{v}{x};\alpha ,\beta \right)e^{-\theta v}dv.\hfill \end{array}$$

Finally, the result is obtained by establishing that $E(x;\theta ,\alpha ,\beta )={\int }_0^x(1+v)B(v/x;\alpha ,\beta )$ $e^{-\theta v}dv$ and recognizing the cdf of the L distribution in the previous expression. □

Corollary 2.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, the rf of X is given by

$$\begin{array}{c}\hfill R_X(x;\theta ,\alpha ,\beta )=R_Y(x;\theta )+\frac{{\theta }^2}{(1+\theta )B(\alpha ,\beta )}E\left(x;\theta ,\alpha ,\beta \right),\end{array}$$

and the hrf of X is given by

$$\begin{array}{c}\hfill h_X(x;\theta ,\alpha ,\beta )=\frac{{\theta }^2\Gamma \left(\beta \right)[r_1(x;\theta ,\alpha ,\beta )+x{r}_2(x;\theta ,\alpha ,\beta )]}{(1+\theta )B(\alpha ,\beta )R_Y(x;\theta )+{\theta }^{2}E(x;\theta ,\alpha ,\beta )},\end{array}$$

where $R_Y(x;\theta )=[(1+\theta +\theta x)/(1+\theta )]e^{-\theta x}$ is the rf of the L distribution.

From the previous result, we see that

$$\begin{array}{ccc}\hfill h_X(x;\theta ,\alpha ,\beta =1)& =& \frac{{\theta }^2[r_1(x;\theta ,\alpha ,1)+x{r}_2(x;\theta ,\alpha ,1)]}{\alpha (1+\theta )R_Y(x;\theta )+{\theta }^{2}E(x;\theta ,\alpha ,1)}\mathrm{and}\hfill \\ \hfill \underset{\alpha \to \infty }{lim}{h}_X(x;\theta ,\alpha ,\beta =1)& =& \frac{{\theta }^2(1+x)}{1+\theta +\theta x},\hfill \end{array}$$

which correspond to the hrfs of the SL and L distributions, respectively.

We provide an R code for the computation of the cdf given in Proposition 2. Using this code, together with the R code for the ESL pdf presented in the previous section, the ESL hrf can be easily computed. Figure 2 presents some ESL hrf curves considering different choices for $\theta$, $\alpha$ and $\beta$. From the figure, the following can be noted:

• The ESL hrf can present unimodal or monotonic decreasing shapes.
• In the left panels ($\beta =1$), the unimodal shape of the hrf of the ESL distribution approaches the increasing shape of the hrf of the L distribution when $\alpha$ assumes a sufficiently large value.
• In the right panels, it is observed that the hrf of the ESL distribution is unimodal when $\theta =0.5$ and decreasing when $\theta =5$. It is also observed that the unimodal or increasing forms are preserved for the different choices of the parameters $\alpha$ and $\beta$. Thus, the unimodal and monotonic decreasing shape seems to depend on $\theta$.

2.3. Moments and Related Measurements

Proposition 3.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, for $r=1,2,\dots$ and $\alpha >r$, it follows that rth raw moment is given by

$$\begin{array}{c}\hfill {\mu }_r=E\left(X^r\right)=\frac{r!(\theta +r+1)}{{\theta }^r(\theta +1)}\prod _{i=1}^r\frac{\alpha +\beta -i}{\alpha -i}.\end{array}$$

Proof. 

Using the representation given in Definition 1, we have that

$$\begin{array}{c}\hfill {\mu }_r=E\left(X^r\right)=E\left({\left(\frac{Y}{U}\right)}^r\right)=E\left(Y^r\right)E\left(U^{-r}\right).\end{array}$$

where it follows that $E\left(U^{-r}\right)=\prod _{i=1}^n\frac{\alpha +\beta -i}{\alpha -i}$, $\alpha >r$ and $E\left(Y^r\right)=\frac{r!(\theta +r+1)}{{\theta }^r(\theta +1)}$ is the rth raw moment of the L distribution. □

Corollary 3.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, it follows that

$$\begin{array}{c}\hfill E\left(X\right)=\frac{(\theta +2)(\alpha +\beta -1)}{\theta (\theta +1)(\alpha -1)},\alpha >1and\end{array}$$

$$Var\left(X\right)=\frac{\alpha +\beta -1}{\theta (\theta +1)(\alpha -1)}\left(\frac{2(\theta +3)(\alpha +\beta -2)}{\theta (\alpha -2)}-\theta -2\right),\alpha >2.$$

Corollary 4.

Let $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, Fisher’s skewness $\left(\sqrt{{\beta }_1}\right)$ and kurtosis $\left({\beta }_2\right)$ coefficients are given by

$$\begin{array}{c}\hfill \sqrt{{\beta }_1}=\frac{2{(\alpha -2)}^{1/2}A}{{(\alpha +\beta -1)}^{1/2}(\alpha -3)B^{3/2}},\alpha >3,and{\beta }_2=\frac{3(\alpha -2)C}{(\alpha -3)(\alpha -4)(\alpha +\beta -1)B^2},\alpha >4,\end{array}$$

where

$$\begin{array}{ccc}\hfill A& =& 3{(\theta +1)}^2(\theta +4){(\alpha -1)}^2(\alpha +\beta -2)(\alpha +\beta -3)\hfill \\ & & -3(\theta +1)(\theta +2)(\theta +3)(\alpha -1)(\alpha -3)(\alpha +\beta -1)(\alpha +\beta -2)\hfill \\ & & +{(\theta +2)}^3(\alpha -2)(\alpha -3){(\alpha +\beta -1)}^2,\hfill \\ & & \\ \hfill B& =& 2(\theta +1)(\theta +3)(\alpha -1)(\alpha +\beta -2)-{(\theta +2)}^2(\alpha -2)(\alpha +\beta -1)and\hfill \\ & & \\ \hfill C& =& 8{(\theta +1)}^3(\theta +5){(\alpha -1)}^3(\alpha +\beta -2)(\alpha +\beta -3)(\alpha +\beta -4)\hfill \\ & & -8{(\theta +1)}^2(\theta +2)(\theta +4){(\alpha -1)}^2(\alpha -4)(\alpha +\beta -1)(\alpha +\beta -2)(\alpha +\beta -3)\hfill \\ & & +4(\theta +1){(\theta +2)}^2(\theta +3)(\alpha -1)(\alpha -3)(\alpha -4){(\alpha +\beta -1)}^2(\alpha +\beta -2)\hfill \\ & & -{(\theta +2)}^4(\alpha -2)(\alpha -3)(\alpha -4){(\alpha +\beta -1)}^3.\hfill \end{array}$$

Figure 3 depicts plots of the Fisher’s skewness and kurtosis coefficients of the ESL distribution. These plots show that the coefficients decrease as $\theta$ and $\alpha$ decrease and as $\beta$ increases. Although the coefficients depend on the three parameters, it can be seen that $\theta$ has a slight impact on the skewness and kurtosis, while $\alpha$ and $\beta$ have a great impact, especially on kurtosis.

3. Parameter Estimation

This section discusses parameter estimation for the ESL distribution considering the moment and maximum likelihood methods. In addition, we carry out a simulation study that illustrates the behavior of the estimators under different sample sizes.

3.1. Moment Estimation

Proposition 4.

Let $X_1,\dots ,X_n$ be a random sample of the random variable $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, for $\alpha >3$, the moment estimators ${\widehat{\theta }}_{\mathrm{M}}$, ${\widehat{\alpha }}_{\mathrm{M}}$ and ${\widehat{\beta }}_{\mathrm{M}}$ for θ, α and β are given by

$$\begin{array}{c}\hfill \left[\frac{\overline{X^2}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}+2)}{2\overline{X}({\widehat{\theta }}_{\mathrm{M}}+3)}-1\right]\left[\alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)-2\right]-\beta \left({\widehat{\theta }}_{\mathrm{M}}\right)=0,\\ \hfill {\widehat{\alpha }}_{\mathrm{M}}=\alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)\mathrm{and}{\widehat{\beta }}_{\mathrm{M}}=\beta \left({\widehat{\theta }}_{\mathrm{M}}\right),\end{array}$$

(6)

where

$$\begin{array}{ccc}\hfill \alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)& =& \frac{6[\overline{X}({\widehat{\theta }}_{\mathrm{M}}^2+7{\widehat{\theta }}_{\mathrm{M}}+12)+\overline{X^2}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}^2+6{\widehat{\theta }}_{\mathrm{M}}+8)-\overline{X^3}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}^2+6{\widehat{\theta }}_{\mathrm{M}}+9)]}{{\widehat{\theta }}_{\mathrm{M}}[3\overline{X^2}({\widehat{\theta }}_{\mathrm{M}}^2+6{\widehat{\theta }}_{\mathrm{M}}+8)-2\overline{X^3}({\widehat{\theta }}_{\mathrm{M}}^2+6{\widehat{\theta }}_{\mathrm{M}}+9)]},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \beta \left({\widehat{\theta }}_{\mathrm{M}}\right)& =& \left[\frac{\overline{X^3}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}+3)}{3\overline{X^2}({\widehat{\theta }}_{\mathrm{M}}+4)}-1\right]\left[\alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)-3\right].\hfill \end{array}$$

(8)

Proof. 

The first three equations in the method of moments are given by ${\mu }_j=\overline{X^j}$, with $j=1,2,3$, where ${\mu }_j$ is as in Proposition 3. Thus, combining these equations, we obtain that

$$\begin{array}{ccc}\hfill \beta \left({\widehat{\theta }}_{\mathrm{M}}\right)& =& \left[\frac{\overline{X^2}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}+2)}{2\overline{X}({\widehat{\theta }}_{\mathrm{M}}+3)}-1\right]\left[\alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)-2\right],\hfill \\ \hfill \beta \left({\widehat{\theta }}_{\mathrm{M}}\right)& =& \left[\frac{\overline{X^3}{\widehat{\theta }}_{\mathrm{M}}({\widehat{\theta }}_{\mathrm{M}}+3)}{3\overline{X^2}({\widehat{\theta }}_{\mathrm{M}}+4)}-1\right]\left[\alpha \left({\widehat{\theta }}_{\mathrm{M}}\right)-3\right],\hfill \end{array}$$

where the last of these equations corresponds to the result given in Equation (8) and the algebraic manipulation of the first gives rise to Equation (6). Finally, Equation (7) is obtained by equating the previous equations under a suitable algebra. □

3.2. Maximum Likelihood Estimation

Suppose $X_1,\dots ,X_n$ is a random sample of the random variable $X\sim \mathrm{ESL}(\theta ,\alpha ,\beta )$. Then, from Equation (3), we have that the log-likelihood function can be written as

$$\begin{array}{ccc}\hfill \ell (\theta ,\alpha ,\beta ;x_i)& =& log\prod _{i=1}^n{f}_X(x_i;\theta ,\alpha ,\beta )\hfill \\ & =& 2nlog\left(\theta \right)-nlog(\theta +1)-log\left[B(\alpha ,\beta )\right]+\sum _{i=1}^{n}log\left[K\left(x_i\right)\right],\hfill \end{array}$$

(9)

where

$$\begin{array}{c}\hfill K\left(x_i\right)={\int }_0^1(1+x_{i}w)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw.\end{array}$$

The maximum likelihood (ML) estimators ${\widehat{\theta }}_{\mathrm{ML}}$, ${\widehat{\alpha }}_{\mathrm{ML}}$ and ${\widehat{\beta }}_{\mathrm{ML}}$ for $\theta$, $\alpha$ and $\beta$ can be obtained by taking the partial derivatives of Equation (9) and solving the corresponding system of equations, which is given by the equations

$$\begin{array}{c}\hfill \frac{2n}{\theta }+\frac{n}{\theta +1}+\sum _{i=1}^n\frac{K_1\left(x_i\right)}{K\left(x_i\right)}=0,\end{array}$$

(10)

$$\begin{array}{c}\hfill n\mathsf{\Psi }\left(\alpha \right)-n\mathsf{\Psi }(\alpha +\beta )+\sum _{i=1}^n\frac{K_2\left(x_i\right)}{K\left(x_i\right)}=0,\end{array}$$

(11)

$$\begin{array}{c}\hfill n\mathsf{\Psi }\left(\beta \right)-n\mathsf{\Psi }(\alpha +\beta )+\sum _{i=1}^n\frac{K_3\left(x_i\right)}{K\left(x_i\right)}=0,\end{array}$$

(12)

where

$$\begin{array}{ccc}\hfill K_1\left(x_i\right)& =& -x_i{\int }_0^1(1+x_{i}w)w^{\alpha +1}{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_2\left(x_i\right)& =& {\int }_0^1(1+x_{i}w)log\left(w\right)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_3\left(x_i\right)& =& {\int }_0^1(1+x_{i}w)log(1-w)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \end{array}$$

such that $\mathsf{\Psi }(·)$ is the digamma function.

The asymptotic distribution (under regularity conditions) of the ML estimator of $\delta ={(\theta ,\alpha ,\beta )}^{\prime }$ is $N_3(\delta ,I{\left(\delta \right)}^{-1})$, where $I{\left(\delta \right)}^{-1}$ is the expected information matrix. Taking into account the structure of Equation (9), we observe that it is not easy to derive the analytical expression of this matrix. Thus, we consider an approximation from the observed information matrix, where the elements of this matrix are computed as minus the second partial derivatives of the log-likelihood function with respect to each parameter (assessed in the ML estimates). The observed information matrix is given by

$$\begin{array}{c}\hfill I\left(\delta \right)=\left(\begin{array}{ccc}{I}_{\theta \theta }& I_{\alpha \theta }& I_{\beta \theta }\\ & I_{\alpha \alpha }& I_{\beta \alpha }\\ & & I_{\beta \beta }\end{array}\right),\end{array}$$

where

$$\begin{array}{ccc}\hfill I_{\theta \theta }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial {\theta }^2}=-\sum _{i=1}^n\frac{K_{11}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_1^2\left(x_i\right)}{K^2\left(x_i\right)}+\frac{2n}{{\theta }^2}-\frac{n}{{(\theta +1)}^2},\hfill \\ \hfill I_{\alpha \alpha }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial {\alpha }^2}=-\sum _{i=1}^n\frac{K_{22}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_2^2\left(x_i\right)}{K^2\left(x_i\right)}+n{\mathsf{\Psi }}_1\left(\alpha \right)-n{\mathsf{\Psi }}_1(\alpha +\beta ),\hfill \\ \hfill I_{\beta \beta }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial {\beta }^2}=-\sum _{i=1}^n\frac{K_{33}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_3^2\left(x_i\right)}{K^2\left(x_i\right)}+n{\mathsf{\Psi }}_1\left(\beta \right)-n{\mathsf{\Psi }}_1(\alpha +\beta )\},\hfill \\ \hfill I_{\alpha \theta }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial \alpha \partial \theta }=-\sum _{i=1}^n\frac{K_{12}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_1\left(x_i\right)K_2\left(x_i\right)}{K^2\left(x_i\right)},\hfill \\ \hfill I_{\beta \theta }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial \beta \partial \theta }=-\sum _{i=1}^n\frac{K_{13}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_1\left(x_i\right)K_3\left(x_i\right)}{K^2\left(x_i\right)},\hfill \\ \hfill I_{\beta \alpha }& =& -\frac{{\partial }^2\ell (\theta ,\alpha ,\beta ;x_i)}{\partial \beta \partial \alpha }=-\sum _{i=1}^n\frac{K_{23}\left(x_i\right)}{K\left(x_i\right)}+\sum _{i=1}^n\frac{K_2\left(x_i\right)K_3\left(x_i\right)}{K^2\left(x_i\right)}+{\mathsf{\Psi }}_1(\alpha +\beta ),\hfill \end{array}$$

such that

$$\begin{array}{ccc}\hfill K_{11}\left(x_i\right)& =& x_i^2{\int }_0^1(1+x_{i}w)w^{\alpha +2}{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_{22}\left(x_i\right)& =& {\int }_0^1(1+x_{i}w){log}^2\left(w\right)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill K_{33}\left(x_i\right)& =& {\int }_0^1(1+x_{i}w){log}^2(1-w)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_{12}\left(x_i\right)& =& -x_i{\int }_0^1(1+x_{i}w)log\left(w\right)w^{\alpha +1}{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_{13}\left(x_i\right)& =& -x_i{\int }_0^1(1+x_{i}w)log(1-w)w^{\alpha +1}{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw,\hfill \\ \hfill K_{23}\left(x_i\right)& =& {\int }_0^1(1+x_{i}w)log\left(w\right)log(1-w)w^{\alpha }{(1-w)}^{\beta -1}{e}^{-\theta x_{i}w}dw.\hfill \end{array}$$

and ${\mathsf{\Psi }}_1(·)$ is the trigamma function.

Then, approximate $100(1-\phi )$% confidence intervals for $\theta$, $\alpha$ and $\beta$ can be determined by $\widehat{\theta }\pm z_{\phi /2}{s}_{\widehat{\theta }}$, $\widehat{\alpha }\pm z_{\phi /2}{s}_{\widehat{\alpha }}$ and $\widehat{\beta }\pm z_{\phi /2}{s}_{\widehat{\beta }}$, respectively, where $z_{\phi /2}$ is the upper $(\phi /2)$th percentile of the standard normal distribution, and $s_{\widehat{\theta }}$, $s_{\widehat{\alpha }}$ and $s_{\widehat{\beta }}$ are the diagonal elements of the matrix ${\left[I\left(\delta \right)\right]}^{-1}$ (assessed in the ML estimates).

3.3. Practical Considerations

From Equation (6), we observe that the moment estimator ${\widehat{\theta }}_{\mathrm{M}}$ of $\theta$ does not have an explicit analytical form, so it is necessary to solve Equation (6) numerically to obtain the estimate. We use the function uniroot of the R programming language to solve this problem.

From the system of Equations (10)–(12), we see that the ML estimators ${\widehat{\theta }}_{\mathrm{ML}}$, ${\widehat{\alpha }}_{\mathrm{ML}}$ and ${\widehat{\beta }}_{\mathrm{ML}}$ for $\theta$, $\alpha$ and $\beta$ do not have closed forms. Therefore, the ML estimates must be obtained by numerically solving the system of Equations (10)–(12), which can be quite a complicated problem. In this case, we prefer to obtain the ML estimates by solving the optimization problem ${max}_{(\theta ,\alpha ,\beta )}\ell (\theta ,\alpha ,\beta ;x_i)$, subject to $\theta ,\alpha ,\beta >0$, where $\ell (\theta ,\alpha ,\beta ;x_i)$ is as in Equation (9). We solve this problem using the function optim of the R programming language under consideration of the L-BFGS-B algorithm [26]. This algorithm requires the declaration of initial values in the parameter space to initialize the iterative process. Through simulation experiments, we observe that the estimates obtained are strongly influenced by the initial value considered. The first alternative that we considered to overcome this drawback was to use as initial values the moment estimates obtained from the Proposition 4. However, these estimates present an important bias, so they cannot be considered as good initial values. Finally, we note that this drawback occurs due to the existence of two shape parameters that control the kurtosis of the distribution, $\alpha$ and $\beta$, which leads to a problem of non-identifiability of the parameters. This problem disappears if there is only one shape parameter that controls the kurtosis and, for this, different solutions can be adopted. In particular, based on the kurtosis behavior described in Section 2.3, we assume $\beta =1+100/\alpha$, obtaining a more parsimonious distribution with extremely heavy tails. In this choice for $\beta$, the following can be verified:

• The L distribution can be obtained as a limit case of the two-parameter ESL distribution once $\beta \to 1$ as $\alpha \to \infty$.
• Taking into account the kurtosis of the ESL distribution increases as $\beta$ grows, the two-parameter ESL distribution has very heavy tails since $\beta \to \infty$ as $\alpha \downarrow 0$.
• In the choice $\beta =1+100/\alpha$, the number 100 allows us to obtain a high value for $\beta$ when $\alpha$ assumes a value considered small. For example, we get $\beta =1+100/2=51$ when $\alpha =2$. Thus, the representation of the ESL random variable given in Proposition (1) assumes a Beta$(\alpha =2,\beta =51)$ distribution for the variable U in the denominator, thus defining a distribution with extremely heavy tails.

For practical purposes, in accordance with the above, we recommend the reader to use the two-parameter version of the ESL distribution.

3.4. Simulation Study

In this section, we carry out a simulation study to evaluate the behavior of the moment and ML estimators assuming that $\beta =1+100/\alpha$. The study is developed considering the following stages:

• We choose the values 0.5 and 2 for the parameter $\theta$ and the values 3, 4 and 5 for the parameter $\alpha$. In this way, we establish three scenarios where the pdf of the ESL distribution is unimodal (when $\theta =0.5$) and three other scenarios where the pdf is decreasing (when $\theta =2$).
• Under the consideration of the choices above, we generate 1000 pseudo-random samples of size $n=100$, 200, 300, 400 and 500 from the two-parameter ESL distribution. Samples are generated as follows:

  (a)

  Generate Y having a L$\left(\theta \right)$ distribution.

  • Generate $S\sim Exponential\left(\theta \right)$.
  • Generate $W\sim Gamma(2,\theta )$.
  • Generate $V\sim Bernoulli\left(\frac{\theta }{1+\theta }\right)$.
  • Compute $Y=VS+(1-V)W$.

  (b)

  Generate X having a ESL$(\theta ,\alpha )$ distribution.

  • Generate $U\sim Beta(\alpha ,1+100/\alpha )$.
  • Compute $X=Y{U}^{-1}$.

  Step (a) is based on the results proposed by Ghitany et al. [3] for the L distribution, while step (b) is based on the representation of the ESL random variable given in Definition 1. An R code is provided in Appendix A.

• For each sample generated, we obtain the ML estimates by maximizing Equation (9) using the optim function of the R programming language.

Table 1. The average estimate (AE), the average of the estimated bias (AB) and the standard deviation (SD) of the estimates obtained in each scenario and sample size considered in the study.

Parameter		${\widehat{\mathit{\theta }}}_{\mathbf{M}}$			${\widehat{\mathit{\alpha }}}_{\mathbf{M}}$			${\widehat{\mathit{\theta }}}_{\mathbf{ML}}$			${\widehat{\mathit{\alpha }}}_{\mathbf{ML}}$
$\theta$	$\alpha$	AE	SD	AB	AE	SD	AB	AE	SD	AB	AE	SD	AB
$n=100$
0.5	3	0.351	0.147	−0.149	3.965	1.103	0.965	0.495	0.189	−0.005	3.250	0.771	0.250
	4	0.413	0.190	−0.087	5.035	1.602	1.035	0.506	0.201	0.006	4.354	1.206	0.354
	5	0.451	0.216	−0.049	6.197	2.301	1.197	0.515	0.220	0.015	5.490	1.676	0.490
2	3	1.370	0.556	−0.630	4.041	1.123	1.041	2.110	0.897	0.110	3.299	0.895	0.299
	4	1.662	0.835	−0.338	5.213	1.897	1.213	2.210	1.236	0.210	4.422	1.564	0.422
	5	1.829	0.976	−0.171	6.420	2.691	1.420	2.229	1.321	0.229	5.686	2.496	0.686
$n=200$
0.5	3	0.381	0.128	−0.119	3.637	0.681	0.637	0.499	0.134	−0.001	3.113	0.482	0.113
	4	0.438	0.167	−0.062	4.623	1.025	0.623	0.501	0.138	0.001	4.160	0.736	0.160
	5	0.465	0.172	−0.035	5.687	1.409	0.687	0.505	0.149	0.005	5.261	1.100	0.261
2	3	1.508	0.526	−0.492	3.699	0.756	0.699	2.043	0.688	0.043	3.131	0.595	0.131
	4	1.778	0.744	−0.222	4.684	1.160	0.684	2.070	0.755	0.070	4.207	0.940	0.207
	5	1.885	0.837	−0.115	5.840	1.771	0.840	2.084	0.795	0.084	5.382	1.550	0.382
$n=300$
0.5	3	0.396	0.121	−0.104	3.521	0.576	0.521	0.499	0.108	−0.001	3.079	0.383	0.079
	4	0.460	0.166	−0.040	4.432	0.821	0.432	0.501	0.116	0.001	4.077	0.581	0.077
	5	0.478	0.157	−0.022	5.467	1.105	0.467	0.504	0.115	0.004	5.160	0.814	0.160
2	3	1.603	0.512	−0.397	3.529	0.632	0.329	2.024	0.549	0.024	3.089	0.479	0.089
	4	1.836	0.678	−0.164	4.488	0.919	0.488	2.063	0.625	0.063	4.098	0.743	0.098
	5	1.910	0.749	−0.090	5.581	1.270	0.581	2.075	0.667	0.075	5.213	1.118	0.213
$n=400$
0.5	3	0.410	0.114	−0.090	3.424	0.508	0.424	0.500	0.095	0.000	3.045	0.329	0.045
	4	0.460	0.144	−0.040	4.378	0.730	0.378	0.501	0.095	0.001	4.078	0.481	0.078
	5	0.483	0.153	−0.017	5.359	0.946	0.359	0.502	0.104	0.002	5.112	0.699	0.112
2	3	1.653	0.486	−0.347	3.451	0.563	0.451	2.017	0.478	0.017	3.053	0.412	0.053
	4	1.858	0.652	−0.142	4.406	0.794	0.406	2.030	0.516	0.030	4.083	0.625	0.083
	5	1.911	0.643	−0.089	5.470	1.042	0.470	2.033	0.524	0.033	5.170	0.875	0.170
$n=500$
0.5	3	0.418	0.112	−0.082	3.384	0.477	0.384	0.500	0.088	0.000	3.024	0.293	0.024
	4	0.468	0.142	−0.032	4.334	0.716	0.334	0.500	0.091	0.000	4.012	0.473	0.012
	5	0.492	0.138	−0.008	5.251	0.849	0.251	0.501	0.090	0.001	5.050	0.605	0.050
2	3	1.663	0.464	−0.337	3.434	0.520	0.434	2.015	0.426	0.015	3.051	0.364	0.051
	4	1.869	0.601	−0.131	4.361	0.742	0.361	2.028	0.472	0.028	4.080	0.572	0.080
	5	1.962	0.637	−0.038	5.349	0.964	0.349	2.031	0.463	0.031	5.122	0.776	0.122

reports the average estimate (AE), the average of the estimated bias (AB), and the standard deviation (SD) of the estimates obtained in each scenario and each sample size considered. Note that the AEs approach the true values of the parameters and the SDs become small as the sample size increases, an expected result since the ML and moment estimators are consistent.

Comparing the results obtained with each method, we can see that the ML method provides estimates with less bias. It is also observed that the ML estimates for $\alpha$ are more efficient in all the sample sizes considered. Finally, we observe that the SDs of the moment estimates of $\theta$ are smaller than the corresponding SDs obtained with the ML method when the sample size is less than or equal to 200. However, in this case, the AEs provided by the moment method are much more distant from the real value of $\theta$ than the AEs provided by the ML method. For sample sizes greater than 200, the ML estimates for $\theta$ are less biased and more efficient.

4. Real Data Analysis

In this section, by fitting data featuring outliers, we compare the performance of the two-parameter ESL distribution with that of the L, power-Lindley (PL) [6], and LS distributions. The pdfs of the L and LS distributions can be consulted in Section 1, while the pdf of the PL distribution is given by $f(x;\theta ,\alpha )=[\left(\alpha {\theta }^2\right)/(\theta +1)](1+x^{\alpha })x^{\alpha -1}{e}^{-\theta x^{\alpha }}$, with $x>0$ and $\theta ,\alpha >0$. Note that the PL distribution is one of the most popular extensions to the L distribution, also note that this distribution has the same parameter dimension as the two-parameter ESL distribution.

We evaluate the comparative performance of the fitted distributions using the Akaike Information Criterion (AIC) [27] and the Bayesian Information Criterion (BIC) [28]. On the other hand, we consider the goodness-of-fit tests proposed by Chen and Balakrishnan [29] to assess the quality of fit of the distributions. Specifically, based on the modified Cramer-von Mises ($W^{*}$) and Anderson–Darling ($A^{*}$) statistics, we test the hypothesis $H_0:x_1,\dots ,x_n$ is an observed random sample from a continuous distribution $F(x;\delta )$, where $F(·;·)$ is known, but $\delta$ must be estimated efficiently. Here, $H_0$ is rejected at significance level equal to 0.05 if $W^{*}>0.126$ and $A^{*}>0.752$. Under a significance level equal to 0.10, $H_0$ is rejected when $W^{*}$ and $A^{*}$ are greater than 0.104 and 0.631, respectively.

For the analyzed data, we evidenced the presence of outliers by elaborating the traditional box-plot and the adjusted box-plot [30] for skewed data. Note that the second plot includes a robust measure of skewness in the determination of the whiskers, which helps to avoid misreporting of outliers. We use the adjbox function [31] available in the R programming language for the elaboration of this plot.

4.1. Time Failures Related to Pascal Programming

The data considered in this application were collected from a single-user workstation at the Center for Software Reliability and represents the time failures related to Pascal programming. Some descriptive statistics for these data are as follows: sample size, 104; average, 147.8; standard deviation, 245.095; Fisher’s skewness coefficient, 3; Fisher’s kurtosis coefficient, 14.579. Figure 4 presents the traditional box-plot and the adjusted box-plot for the time failure data in which the existence of outliers can be observed.

Previous studies with these data can be found in Lyu [32] and Astorga et al. [33].

Computing initially the moment estimators under the ESL distribution, we have the following estimates: ${\widehat{\theta }}_{\mathrm{M}}=0.274$ and ${\widehat{\alpha }}_{\mathrm{M}}=2.730$. Using the moment estimates as initial values, the ML estimates are computed and presented in

Table 2. The ML estimates (with standard errors in parentheses), the maximum values of the corresponding log-likelihood functions, the values associated with the information criteria (AIC and BIC) and with the modified statistics ($W^{*}$ and $A^{*}$) for each fitted distribution to the time failure data.

Parameter	L	PL	LS	ESL
$\sigma$	-	-	$2.735\times {10}^4$	-
			($1.046\times {10}^4$)
$\theta$	0.013	0.215	$1.169\times {10}^3$	1.819
	(0.001)	(0.034)	($5.483\times {10}^2$)	(0.004)
$\alpha$	-	0.480	0.799	1.284
		(0.030)	0.148	(0.003)
ℓ	−703.450	−600.236	−602.341	−599.780
AIC	1408.900	1204.472	1210.683	1203.561
BIC	1411.545	1209.761	1218.616	1208.850
W${}^{*}$	0.353	0.160	0.119	0.084
A${}^{*}$	2.087	0.978	0.808	0.566

with the standard errors in parentheses. In addition, we also report the maximum value (ℓ) of the corresponding log-likelihood function for each fitted distribution. Note that the largest ℓ value is reported for the ESL distribution.

Table 2 also reports for each fitted distribution the values associated with the information criteria and with the modified statistics. In the table, based on the $A^{*}$ statistic, considering a level of significance equal to 0.05, it is observed that only the ESL distribution adequately fits the time failure data. An analogous observation can be considered from the values of both statistics under a level of significance equal to 0.10. On the other hand, also in Table 2, it is observed that the ESL distribution is the one with the lowest AIC and BIC values, suggesting that this distribution should be selected to model the time failure data.

Figure 5 depicts plots of the fitted L, PL, LS and ESL distributions using the ML estimates. Notice that the fitted ESL distribution has heavier tails.

4.2. State Personal Income Data

We consider a set of annual observations on state personal income (total, nominal, in millions of dollars) recorded in the 48 continental United States. Stock and Watson [34] analyze these data in a study on the demand for tobacco. Some descriptive statistics for this data are as follows: sample size, 104; average, 99.879; standard deviation, 120.541; Fisher’s skewness coefficient, 2.817; Fisher’s kurtosis coefficient, 13.275. We provide this dataset in Appendix B. Figure 6 presents the traditional box-plot and the adjusted box-plot for the state personal income data in which the existence of outliers can be observed.

With these data, we obtain the following moment estimates: ${\widehat{\theta }}_M=0.207$ and ${\widehat{\alpha }}_M=3.702$.

Table 3. The ML estimates (with standard errors in parentheses), the maximum values of the corresponding log-likelihood functions, the values associated with the information criteria (AIC and BIC) and with the modified statistics ($W^{*}$ and $A^{*}$) for each fitted distribution to the state personal income data.

Parameter	L	PL	LS	ESL
$\sigma$	-	-	$4.260\times {10}^4$	-
			($7.528\times {10}^3$)
$\theta$	0.019	0.083	682.380	0.312
	(0.001)	(0.019)	192.712	(0.089)
$\alpha$	-	0.705	2.614	3.033
		(0.046)	1.147	(0.484)
ℓ	−553.170	−536.240	−536.596	−534.902
AIC	1108.341	1076.480	1079.193	1073.804
BIC	1110.906	1081.609	1086.886	1078.933
$W^{*}$	0.140	0.084	0.074	0.051
$A^{*}$	1.003	0.653	0.579	0.432

reports the ML estimates for the parameters of each distribution fitted to the state personal income data. In addition, the maximum values of the corresponding log-likelihood functions and the values associated with the information criteria and the modified statistics are reported. In this table, based on the statistics modified under a level of significance equal to 0.05, it can be seen that the L distribution is the only distribution that does not properly fit the state personal income data. On the other hand, the ESL distribution is the one with the lowest AIC and BIC values, which suggests that this distribution should be selected to fit the state personal income data.

Figure 7 depicts plots of the fitted L, PL, LS and ESL distributions using the ML estimates. Notice that the fitted ESL distribution has heavier tails.

5. Concluding Remarks

This article proposes a heavy-tailed extension of the Lindley distribution, the extended slash Lindley (ESL) distribution. The proposal arises as the ratio of independent random variables, specifically, a Lindley distribution divided by a beta distribution. In this way, a three-parameter distribution with heavier tails than the classical Lindley distribution is obtained.

The main structural properties such as the pdf, cdf, rf and hrf are derived. The pdf is written in closed form in terms of the confluent hypergeometric function, which facilitates its computational implementation. Both the pdf and the hrf can exhibit monotonically decreasing or unimodal shapes. In addition, closed expressions are derived for the raw moments of the ESL distribution, which are then used to describe the behavior of the Fisher’s skewness and kurtosis coefficients. We note that two of the parameters of the ESL distribution ($\alpha$ and $\beta$) control the kurtosis of the distribution, while the remaining parameter ($\theta$) controls the unimodal or monotonic shape of the pdf.

Although the presence of two shape parameters controlling the kurtosis allows the ESL distribution to be able to capture high levels of kurtosis, we observe that this generates a problem of non-identifiability of the parameters, which was evidenced by simulation experiments. This problem disappears when there is a single parameter that controls the kurtosis, so several solutions to this problem can be adopted. In particular, we propose the special case obtained by considering $\beta =1+100/\alpha$, which results in a two-parameter version capable of capturing high levels of kurtosis. We verify that the pdf and the hrf of this special case, as in the three-parameter version, can have unimodal and monotonic decreasing shapes. We also note that the classical Lindley distribution can be derived as a limit case of this special case.

We discuss parameter estimation considering moment and ML methods. It is noted that both methods require the use of numerical procedures. For the two-parameter ESL distribution, through simulation experiments, we observe that the ML method provides estimates with less bias than the estimates obtained with the moment method.

Finally, when fitting data that present a high level of kurtosis, we observe that the two-parameter ESL distribution can present a better performance than the classical Lindley distribution and some of its known extensions.