A New Class of Exponentiated Beta-Skew-Laplace Distribution

Abstract

This paper introduce two new families of distributions that allow fitting unimodal, bimodal or trimodal data sets. Statistical properties such as distribution function, moments, moment generating function and stochastic representation of these new families are studied in details. The problem of estimating parameters is addressed by considering the maximum likelihood method and Fisher information matrices are derived. A small Monte Carlo simulation study is conducted to examine the performance of the obtained estimators. The methodology developed is illustrated with three real data applications.

Key words
Asymmetry; beta-skew-Laplace distribution; information matrix; maximum likelihood estimation; multimodality

Introduction

Skew-symmetric distributions for modeling skew and bimodal data sets have been proposed by many researchers. Some of these proposals are obtained from generalizations of the standard normal distribution by incorporating additional parameters. Elal-Olivero (2010) for example, introduces the bimodal-normal (BN) distribution and defines the alpha-skew-normal (ASN) density of families which has sufficient flexibility to fit unimodal as well as bimodal type data. Sobhan et al. (2016) extended the ASN model by adding a new parameter to the distribution. This extension, which is called alpha-beta-skew-normal (ABSN) allows to fit data sets of up to four modes and the admissible intervals for asymmetry and kurtosis parameters are wider than those of the skew-normal (SN) and ASN distributions. Bolfarine et al. (2018) proposed a new family of distributions using methods analogous to those of Kim (2005) and Arnold et al. (2009) with the main characteristic that their normalization constant has a closed and simple form, in addition, the Fisher information matrix is non-singular, which guarantees the asymptotic properties of the maximum likelihood estimators of the parameters when large samples are available. Martínez-Flórez et al. (2018) present a flexible-normal model, which extends the SN model and is suitable for fitting censored symmetric and asymmetric uni-bimodal data. One of the main advantages of this model is that involves fewer parameters to be estimated than the mixture of normal distributions.

Other proposals of bimodal distributions have been considered in the statistical literature. Ma & Genton (2004) proposed a flexible class of skew-symmetric distributions for which the probability density function has the form of a product of a symmetric density and a skewing function. The new class of distribution is able to capture skewness, heavy tails and multimodality systematically. Gómez et al. (2009) consider an extension to the skew‐normal model through the inclusion of an additional parameter which can lead to both uni- and bi-modal distributions. Elal-Olivero et al. (2009) consider an extension of the family of skew-elliptical distributions studied by Azzalini (1985) by using Bayesian approach. The new class is referred to as bimodal skew-elliptical (BSE) distributions and the model is able to adopt both uni- and bimodal shapes.

Although some proposals in statistical literature has been developed to deal with the problem of uni-bimodal, symmetric or asymmetric data analysis, there are few works in the statistical literature able of fitting multimodal data (three or more modes), among these works stand out the works of Sobhan et al. (2016) and Ma & Genton (2004). Therefore, it is important to introduce new distribution functions that manage to model data sets that have multiple modes and can fit to different degrees of asymmetry and kurtosis.

The main objective of this work is to introduce two new distributions that allow fitting symmetric data with up to three modes and it is obtained by considering the Laplace distribution. The Laplace distribution has been considered previously by Shams-Harandi & Alamatsaz (2013) to introduce a new class of skew distributions of Elal-Olivero’s type. This class of distributions has been shown to be flexible enough to support both unimodal and bimodal shapes with increasing and U-shaped hazard function for the truncated case at the origin.

Probability distributions such as those introduced by Elal-Olivero et al. (2009), which is obtained by incorporating an additional parameter to the normal distribution have shown to be very useful in the statistical literature, since it allows to expand the range of possibilities to fit data sets, especially with more than a mode, in this work, we also consider the alpha-power family of distributions introduced by Durrans (1992), which has shown to be a useful alternative to the asymmetric distribution of Azzalini (1985) since it has a range of kurtosis higher, see Pewsey et al. (2012). The two new distributions introduced in this paper are capable of adjusting data with two or more modes and a higher degree of kurtosis compared to some competing models and In addition, data with a high degree of asymmetry can also be fitted.

The rest of the article is organized as follows: Section “The Beta-Skew-Laplace Distribution“ introduces the BSL distribution which we will denote it by BSL$\left(\beta \right)$. We discuss some of its important features and properties and consider its truncated situations as well. Further, we consider a location-scale family and obtain moments of this new family. We provide a stochastic representation and a method of generating data from this distribution. Maximum likelihood estimation of parameters of the location-scale BSL$\left(\beta \right)$ distribution is also considered. In Section “Exponentiated Beta-Skew-Laplace Distribution“ we extend the BSL distribution to the exponentiated case and discuss some of its important features and properties. We consider a location-scale family and the extension to censored data is also presented. Section “Simulation Study“ presents a small Monte Carlo simulation study. Section “Illustrations“ describes fitting of the proposed model to some real data sets and compares it with several rival models.

THE BETA-SKEW-LAPLACE DISTRIBUTION

In order to introduce the new class of skew-Laplace distribution, in this section we define the Bimodal-Laplace distribution of order $k$. First, we recall that a random variable $Y$ is said to have a Laplace distribution with location parameter $\mu$ and scale parameter $\sigma$ if its probability density function (pdf) is given by

$$f(y;\mu ,\sigma )=\frac{1}{2\sigma }exp(-\frac{|y-\mu |}{\sigma }),y\in ℝ,$$

which is denoted by L$(\mu ,\sigma )$. We consider the standard Laplace distribution for simplicity, i.e., the case $\mu =0$ and $\sigma =1$, and we define the bimodal-Laplace (BL) distribution of order $k$ as follows:

Definition 1. The random variable $X$ is said to have a beta-skew-Laplace distribution with parameter $\beta$, denoted by $$X\sim \text{B}SL(\beta )$$, if $X$ has pdf given by

$$f(x;\beta )=\frac{({(1-\beta x^3)}^2+1)}{2(1+360{\beta }^2)}\frac{1}{2}exp(-\left|x\right|),x\in ℝ,$$ (1)

where $\beta \in ℝ$ is a parameter of asymmetry.

Notice that the function (1) for given $\beta$ is a proper pdf since $f(x;\beta )\ge 0$ for all ${(x,\beta )}^{\top }\in {ℝ}^2$, and by the moments of standard Laplace distribution, we can see that

$${\int }_{ℝ}[{(1-\beta x^3)}^2+1]\frac{1}{2}exp(-\left|x\right|),dx=2+720{\beta }^2,\forall \beta \in ℝ.$$

Figure 1 presents the shape of the density function of the BSL distribution for some selected values of the parameter $\beta$. From Figure 1 can be seen that the proposed model is very flexible and the parameter $\beta$ have significant effects on the possible number of modes of the distribution.

Figure 1
Probability density function of $𝐁𝐒𝐋\left(𝛃\right)$ for some selected values of $𝛃$.

The following properties of the BSL distribution are deducted immediately from the definition:

1. If $\beta =0$, then $$X\sim \text{L}(0,1)$$.

2. If $\beta \to \pm \infty$, then $$X\stackrel{d}{\to }\text{B}L(3)$$.

3. $$-X\sim \text{B}SL(-\beta )$$.

Proposition 1. The density function of the BSL$\left(\beta \right)$ distribution has at most three modes.

Proof. Differentiating (1) with respect to $x$ we obtain

$$\begin{array}{rl}{f}^{\prime }(x;\beta )& =\frac{\exp (-|x|)}{4(1+360{\beta }^2)}[-{\beta }^2{x}^5|x|+6{\beta }^2{x}^5+2\beta x^2|x|-6\beta x^2-2\frac{x}{|x|}],x\ne 0\\ & =\frac{1}{4(1+360{\beta }^2)}\{\begin{array}{ll}(-{\beta }^2{x}^6+6{\beta }^2{x}^5+2\beta x^3-6\beta x^2-2)\exp (-x),& \text{if }x>0\\ \\ ({\beta }^2{x}^6+6{\beta }^2{x}^5-2\beta x^3-6\beta x^2+2)\exp (x),& \text{if }x<0\end{array}\end{array}$$ (2)

Clearly, each expression in (2) has at most six zeros. Thus, the function $f(x;\beta )$ can have at most three modes. ◻

Proposition 2. Let $$X\sim \text{B}SL(\beta )$$, then for $k=1,2,3\dots$

$$𝔼\left[X^{2k}\right]=\frac{1}{2(1+360{\beta }^2)}(2\Gamma (2k+1)+{\beta }^2\Gamma (2k+7)),$$ (3)

and

$$𝔼\left[X^{2k-1}\right]=-\frac{\beta \Gamma (2k+3)}{2(1+360{\beta }^2)}.$$ (4)

where $\Gamma (\cdot )$ is the Gamma function.

Proof. Knowing that, for $k=1,2,3,\dots$, $𝔼\left[Z^{2k}\right]=\Gamma (2k+1)$ and $𝔼\left[Z^{2k-1}\right]=0$, when $$Z\sim \text{L}(0,1)$$, we obtain

$$\begin{array}{cc}\hfill 𝔼\left[X^{2k}\right]& =\frac{1}{2(1+360{\beta }^2)}(2𝔼[Z^{2k}]-2\beta 𝔼[Z^{2(k+2)-1}]+{\beta }^2𝔼[Z^{2(k+3)}])\hfill \\ \hfill & =\frac{1}{2(1+360{\beta }^2)}(2\Gamma (2k+1)+{\beta }^2\Gamma (2k+7)).\hfill \end{array}$$

on the other hand

$$\begin{array}{cc}\hfill 𝔼\left[X^{2k-1}\right]& =\frac{1}{2(1+360{\beta }^2)}(2𝔼[Z^{2k-1}]-2\beta 𝔼[Z^{2(k+1)}]+{\beta }^2𝔼[Z^{2(k+3)-1}])\hfill \\ \hfill & =-\frac{\beta \Gamma (2k+3)}{2(1+360{\beta }^2)}\mathrm{.◻}\hfill \end{array}$$

Proposition 3. Let $$X\sim \text{B}SL(\beta )$$, then the mean and the variance of $X$, and the indices of skewness and kurtosis, which are denoted by $\sqrt{{\beta }_1}$ and ${\beta }_2$, respectively, are given by

1. $𝔼\left[X\right]=-{\displaystyle \frac{24\beta }{1+360{\beta }^2}}$.

2. Var$\left[X\right]={\displaystyle \frac{2(1+10152{\beta }^2+3628800{\beta }^4)}{1+720{\beta }^2+129600{\beta }^4}}$.

3. $\sqrt{{\beta }_1}={\displaystyle \frac{-576\beta (1-1662{\beta }^2-745200{\beta }^3)}{[2(1+10152{\beta }^2+3628800{\beta }^4){]}^{3/2}}}$.

4. $${\beta }_2={\displaystyle \frac{24(1+216{\beta }^2(343+24{\beta }^2(15997+2028600{\beta }^2+680400029{\beta }^4)))}{[2(1+10152{\beta }^2+3628800{\beta }^4){]}^2}}$$.

Figure 2 presents the behavior of the skewness and kurtosis indices as a function of the $\beta$ parameter. Notice that, the length of the admissible intervals for the kurtosis parameter of the BSL distribution are larger than the corresponding intervals of the alpha-skew-normal (ASN) model of Elal-Olivero (2010), the skew-normal (SN) model by Azzalini (1985), the power-normal (PN) model, see Pewsey et al. (2012); and the alpha-power skew normal (SNAP) distribution by Martínez-Flórez et al. (2014) which are (1.700,3.7489), (3.000,3.869), (1.717,4.3556) and (4.5328,5.4386), respectively.

Figure 2
The graphs of skewness (on the left) and kurtosis (on the right) of $𝐁𝐒𝐋\left(𝛃\right)$ as function of $𝛃$.

Proposition 4. If $ℱ(x;\beta )$ is the cumulative distribution function of $$X\sim \text{B}SL(\beta )$$ then

$$ℱ(x;\beta )=\{\begin{array}{cc}{\displaystyle \frac{2exp\left(x\right)+2\beta \Gamma (4,-x)+{\beta }^2\Gamma (7,-x)}{4(1+360{\beta }^2)}},\hfill & \mathrm{\text{if}}x<0\mathrm{\text{,}}\hfill \\ \hfill \\ 1-{\displaystyle \frac{2exp(-x)-2\beta \Gamma (4,x)+{\beta }^2\Gamma (7,x)}{4(1+360{\beta }^2)}},\hfill & \mathrm{\text{if}}x\ge 0\mathrm{\text{,}}\hfill \end{array}$$

where $\Gamma (n,x)$ is the upper incomplete gamma function defined as:

$$\Gamma (n,x)={\int }_x^{+\infty }{t}^{n-1}{e}^{-t}dt;\mathrm{\text{for}}x\in ℝ,n\in \mathbb{N}$$

Proof. For $x<0$, we obtain

$$\begin{array}{cc}\hfill ℱ(x;\beta )& =P(X\le x)=\frac{1}{4(1+360{\beta }^2)}{\int }_{-\infty }^x(2-2\beta t^3+{\beta }^2{t}^6)exp\left(t\right)dt\hfill \\ \hfill & =\frac{1}{4(1+360{\beta }^2)}[2exp\left(x\right)+2\beta {\int }_{-x}^{+\infty }{t}^{3}exp(-t)dt+{\beta }^2{\int }_{-x}^{+\infty }{t}^{6}exp(-t)dt]\hfill \\ \hfill & =\frac{1}{4(1+360{\beta }^2)}[2exp\left(x\right)+2\beta \Gamma (4,-x)+{\beta }^2\Gamma (7,-x)]\hfill \end{array}$$

On the other hand, for $x\ge 0$, we obtain

$$\begin{array}{cc}\hfill ℱ\left(x\right)& =1-P(X>x)=1-\frac{1}{4(1+360{\beta }^2)}{\int }_x^{+\infty }(2-2\beta t^3+{\beta }^2{t}^6)exp(-t)dt\hfill \\ \hfill & =1-\frac{1}{4(1+360{\beta }^2)}[2{\int }_x^{+\infty }exp(-t)dt-2\beta {\int }_x^{+\infty }{t}^{3}exp(-t)dt+{\beta }^2{\int }_x^{+\infty }{t}^{6}exp(-t)dt]\hfill \\ \hfill & =1-\frac{1}{4(1+360{\beta }^2)}[2exp(-x)-2\beta \Gamma (4,x)+{\beta }^2\Gamma (7,x)]◻\hfill \end{array}$$

Remark 1. We consider a truncated BSL(β) distribution to obtain a new and useful lifetime distribution. A random variable T has a truncated beta-skew Laplace distribution (at zero), denoted by TBSL(β), if its pdf is given by

$$f\left(t\right)=\frac{1+{(1-\beta t^3)}^2}{2(1+6\beta +360{\beta }^2)}exp(-t),t>0,$$ (5)

We consider a truncated Note that, if β=0, (5) reduces to the standard exponential distribution.

The survival function and hazard rate of a random variable $T$ following a TBSL(β) distribution is given by

$$S_T\left(t\right)=P(T>t)=\frac{2exp(-t)-2\beta \Gamma (4,t)+{\beta }^2\Gamma (7,t)}{2(1+6\beta +360{\beta }^2)},t>0,$$

and

$$h_T\left(t\right)=\frac{2-2\beta t^3+{\beta }^2{t}^6}{2exp(-t)-2\beta \Gamma (4,t)+{\beta }^2\Gamma (7,t)}exp(-t),t>0,$$

respectively. Some shapes of the hazard function are presented in Figure 3.

Figure 3
Illustrations of the hazard rate of $𝐁𝐒𝐋\left(𝛃\right)$ for selected values of $𝛃$.

Proposition 5. If $M_X\left(t\right)$ is the moment generating function (MGF) of $$X\sim \text{B}SL(\beta )$$ then

$$M_X(t)=\frac{[(M_1(t)+M_1(-t))-6\beta (M_4(t)-M_4(-t))+360{\beta }^2(M_7(t)+M_7(-t))]}{2(1+360{\beta }^2)}$$

where $M_i\left(t\right)$ is the MGF of a random variable Gamma$(i,1)$ given by:

$$M_i\left(t\right)={(1-t)}^{-i}$$

Proof. From Equation (1), the MGF of the BSL distribution is expressed as

$$\begin{array}{cc}\hfill M_X\left(t\right)& =\frac{{\int }_{-\infty }^{+\infty }{e}^{tx}(2-2\beta x^3+{\beta }^2{x}^6)exp(-\left|x\right|)dx}{4(1+360{\beta }^2)}\hfill \\ \hfill & =\frac{{\int }_0^{+\infty }(2+2\beta x^3+{\beta }^2{x}^6)e^{-(1+t)x}dx+{\int }_0^{+\infty }(2-2\beta x^3+{\beta }^2{x}^6)exp(-(1-t)x)dx}{4(1+360{\beta }^2)}\hfill \\ \hfill & =\frac{(M_1\left(t\right)+M_1(-t))-6\beta (M_4\left(t\right)-M_4(-t))+360{\beta }^2\beta (M_7\left(t\right)+M_7(-t))}{2(1+360{\beta }^2)}\hfill \end{array}$$

where $M_i\left(t\right)={(1-t)}^{-i}.$ ◻

Remark 2. Let $$X\sim \text{B}SL(\beta )$$. The beta-skew-Laplace family of distributions with location and scale parameters can be defined as the distribution of $Z=\mu +\sigma X$, which has pdf given by

$$f_{BSL}(z;\mu ,\sigma ,\beta )=\frac{{(1-\beta {\left(\frac{z-\mu }{\sigma }\right)}^3)}^2+1}{4\sigma (1+360{\beta }^2)}exp(-\frac{|z-\mu |}{\sigma }),z\in ℝ,$$ (6)

and it is denoted by $$Z\sim \text{B}SL(\mu ,\sigma ,\beta )$$.

The distribution function associated to the density function (6) is given by

$${ℱ}_{BSL}(z;\mu ,\sigma ,\beta )=\{\begin{array}{cc}{\displaystyle \frac{2exp\left(\frac{z-\mu }{\sigma }\right)+2\beta \Gamma (4,-\frac{z-\mu }{\sigma })+{\beta }^2\Gamma (7,-\frac{z-\mu }{\sigma })}{4(1+360{\beta }^2)}},\hfill & \mathrm{\text{if}}z<\mu \mathrm{\text{,}}\hfill \\ \hfill \\ 1-{\displaystyle \frac{2exp(-\frac{z-\mu }{\sigma })-2\beta \Gamma (4,\frac{z-\mu }{\sigma })+{\beta }^2\Gamma (7,\frac{z-\mu }{\sigma })}{4(1+360{\beta }^2)}},\hfill & \mathrm{\text{if}}z\ge \mu \mathrm{\text{,}}\hfill \end{array}$$ (7)

and the moment of order $r$ by

$$𝔼\left[Z^r\right]=\sum _{i=0}^r\left(\genfrac{}{}{0ex}{}{r}{i}\right){\sigma }^i{\mu }^{r-i}𝔼\left[X^i\right],r=1,2,\dots$$

where $𝔼\left[X^i\right]$ is the $i$-th moment of $$X\sim \text{B}SL(\beta )$$.

Maximum likelihood estimation

We consider a random sample $$\mathbf{z}=(z_1,\dots ,z_n{)}^{\mathrm{\top }}$$ of size $n$ from BSL$(\mu ,\sigma ,\beta )$ distribution. The log-likelihood distribution is given by

$$\ell (\mu ,\sigma ,\beta )=-nlog\sigma -nlog4(1+360{\beta }^2)+\sum _{i=1}^{n}log(1+{(1-\beta x_i^3)}^2)-\sum _{i=1}^n\left|x_i\right|$$ (8)

where $x_i=(z_i-\mu )/\sigma$. Note that $\ell (\mu ,\sigma ,\beta )$ is a continuous function in each parameter, however, it is not differentiable at $z_i=\mu$, $i=1,\dots ,n$. We assume $z_i\ne \mu$, $i=1,\dots ,n$, then, the score functions can be obtained by taking the derivative of the likelihood function with respect to parameters $\mu$, $\sigma$ and $\beta$. By equating the score functions to zero and solving the system equations, we provide the maximum likelihood estimators (MLEs) of $\mu$, $\sigma$ and $\beta$, which do not have a closed form. The MLEs can be obtained by numerical method such as Newton-Rapshon type procedure. We use the maxLike function of the statistical software R Development Core Team (2018) for maximizing the likelihood function $\ell (\mu ,\sigma ,\beta )$ which use the Newton-Rapshon method. The Fisher information matrix, which is denoted by $$I_{\mathbf{y}}(\mathrm{\Theta })$$, where $\Theta =(\mu ,\sigma ,\beta )$, has entries given by $i_{{\theta }_k{\theta }_{k\prime }}=-𝔼\left[{\partial }^2\ell /\partial {\theta }_k\partial {\theta }_{k\prime }\right]$, for $k,k\prime =1,2,3$, with $({\theta }_1,{\theta }_2,{\theta }_3)=(\mu ,\sigma ,\beta )$.

We can observe that the BSL$(\mu ,\sigma ,\beta )$ family does not satisfy the regularity conditions, because the Laplace distribution, as particular case of the BSL family density, is not differentiable at $\mu$. However, under a weaker assumption that the density is absolutely continuous, which is the case for the Laplace density (see Kotz et al. (2001)), we can obtain the observed Fisher information matrix $$J_{\mathbf{y}}(\hat{\mathrm{\Theta }})$$, with entries $j_{{\theta }_k{\theta }_{k\prime }}=-{\partial }^2\ell /\partial {\theta }_k\partial {\theta }_{k\prime }$, for $k,k\prime =1,2,3$. Thus, we can find a lower bound for the standard errors (SE) of the MLEs as the square root of the diagonal elements of the observed Fisher information matrix.

EXPONENTIATED BETA-SKEW-LAPLACE DISTRIBUTION

Durrans (1992) introduced an asymmetric version of the normal distribution known as exponentiated normal or power-normal model. Its pdf is given by

$$f_{PN}(x;\alpha )=\alpha \varphi \left(x\right)(\Phi \left(x\right){)}^{\alpha -1},x\in ℝ,\alpha \in {ℝ}^{+},$$ (9)

The power-normal is denoted shortly by PN$\left(\alpha \right)$ and the $\alpha$ parameter controls the asymmetry and kurtosis of the distribution. The PN model is an alternative to the skew-normal model of Azzalini (1985) to fit asymmetric data and it was studied widely by other authors, see for example Pewsey et al. (2012) and Gupta & Gupta (2008). Based on the general structure of the model in (11), which is known as the alpha-power (AP) family of distributions with pdf given by

$$f_{AP}(x;\alpha )=\alpha f\left(x\right)(ℱ\left(x\right){)}^{\alpha -1},x\in ℝ,\alpha \in {ℝ}^{+},$$(10)

where $ℱ$ is an absolutely continuous cdf, and $f=dℱ/dx$ is the respective pdf, in this section we introduce an extension of the BSL distribution to exponentiated families (10) by considering the $\alpha$ parameter.

Definition 2. The random variable $X$ is said to have an exponentiated beta-skew-Laplace distribution, which we will denote by $$X\sim \text{E}BSL(\mu .\sigma ,\alpha ,\beta )$$, if $X$ has pdf given by

$$f_{EBSL}(x;\mu ,\sigma ,\beta )=\alpha f_{BSL}(x;\mu ,\sigma ,\beta )({ℱ}_{BSL}(x;\mu ,\sigma ,\beta ){)}^{\alpha -1},x\in ℝ,$$ (11)

where $f_{BSL}(\cdot )$ and ${ℱ}_{BSL}(\cdot )$ are the pdf and the cdf of the BSL distribution, respectively.

The graph in Figure 4 shows the behavior of EBSL distribution for some selected values of $\alpha$ and $\beta$ for the standard case, that is, $\mu =0$ and $\sigma =1$ which is denoted by EBSL$(\alpha ,\beta )$. It is observed from the graph that, the distribution has bimodal and trimodal behavior.

Figure 4
Probability density function of EBSL$(𝛂,𝛃)$ for selected values of $𝛂$ and $𝛃$.

Properties of the EBSL distribution

Let $$X\sim \text{E}BSL(\mu .\sigma ,\alpha ,\beta )$$,

1. If $\alpha =1$, then $$X\sim \text{B}SL(\mu ,\sigma ,\beta )$$.

2. If $\beta =0$, the pdf (11) reduces to the following pdf

   $$f_{EL}(x;\mu ,\sigma ,\alpha )=\alpha f_L(x;\mu ,\sigma )({ℱ}_L(x;\mu ,\sigma ){)}^{\alpha -1},x\in ℝ,$$(12)
   where $f_L(\cdot )$ and ${ℱ}_L(\cdot )$ are the pdf and cdf of the Laplace distribution. In sequel Equation (12) is referred as the exponentiated Laplace (EL) distribution, $$X\sim \text{E}L(\mu ,\sigma ,\alpha )$$.

3. If $\alpha =1$ and $\beta =0$, then \(X\sim{\textrm L}(\mu,\sigma)\).

4. If $\beta \to \pm \infty$, then $$X\stackrel{d}{\to }\text{B}L(\mu ,\sigma ,3)$$, that is, the bimodal-symmetric Laplace of order 3 distribution of scale-location follows.

Lemma 1. The density function (11) has at most three modes.

Lemma 2. Let $$Y\sim \text{E}BSL(\alpha ,\beta )$$, then the $k$-th moment of $Y$ is given by

$$𝔼[Y^k]=\frac{\alpha }{4(1+360{\beta }^2)}{\int }_{-\infty }^{+\infty }(2y^k-2\beta y^{k+3}+{\beta }^2{y}^{k+6})exp(-\left|y\right|)[{ℱ}_{BSL}(y;\alpha ,\beta ){]}^{\alpha -1}dy.$$ (13)

Maximum likelihood estimation

We consider a random sample $$\mathbf{z}=(z_1,\dots ,z_n{)}^{\mathrm{\top }}$$ of size $n$ from EBSL$(\mu ,\sigma ,\alpha ,\beta )$ distribution. The log-likelihood distribution is given by

$$\begin{array}{rl}\ell (\mu ,\sigma ,\alpha ,\beta )& amp;=n\log \alpha -n\log \sigma -n\log 4(1+360{\beta }^2)-\sum _{i=1}^n|x_i|\\ & amp;+\sum _{i=1}^n\log \{1+{(1-\beta x_i^3)}^2\}+(\alpha -1)\sum _{i=1}^n\log {\mathcal{F}}_{BSL}(x_i;\mu ,\sigma ,\beta )\end{array}$$

where $x_i=(z_i-\mu )/\sigma$. The function $\ell (\mu ,\sigma ,\alpha ,\beta )$ is continuous in each parameter, but it is not differentiable at $z_i=\mu$, $i=1,\dots ,n$. To obtain the MLEs of $\mu$, $\sigma$, $\alpha$ and $\beta$ and their corresponding standard errors, we proceed similarly as in BSL model.

Extension to Censored Data

Based on the goodness of the EBSL distribution to fit asymmetric uni-multimodal data, in this section we introduce the censored EBSL model which we will be denote by CEBSL.

Definition 3. Suppose that the random variable $X$ follows EBSL distribution and consider a random sample $\mathbf{\text{𝐗}}=(X_1,\dots ,X_n)$, where only the $X_i$ values greater than a constant $k$ are recorded. In addition, for values $X_i\le k$ only the value of $k$ is recorded. Therefore, for $i=1,2,\dots ,n$, the observed values $X_i^o$, can be written as

$X_i^o=\{\begin{array}{cc}k,\hfill & \mathrm{\text{if}}{X}_i\le k,\hfill \\ X_i,\hfill & \mathrm{\text{if}}{X}_i>k.\hfill \end{array}$

The resulting sample is said to be a censored EBSL and we say that $X$ is a censored random variable EBSL. We will use the notation $X\sim \mathrm{\text{CEBSL}}(\mu ,\sigma ,\alpha ,\beta )$.

From Definition 3 it follows that $P(X_i^o=k)=P(X_i\le k)=({ℱ}_{BSL}((k-\mu )/\sigma ){)}^{\alpha }$ and for the observations $X_i^o=X_i$ the distribution of $X_i^o$ is the same of $X_i$, i.e., $X_i^o\sim \mathrm{\text{EBSL}}(\mu ,\sigma ,\alpha ,\beta )$. For convenience, we choose to work with the case of left-censored data, however, the followings results can be extended to other types of censorship.

Properties of the CEBSL Model

Let $X\sim \mathrm{\text{CEBSL}}(\mu ,\sigma ,\alpha ,\beta )$,

1. If $\alpha =1$, then $X\sim \mathrm{\text{CBSL}}(\mu ,\sigma ,\beta )$, where CBSL indicates the censored beta-skew-Laplace model.

2. If $\beta =0$, then $X\sim \mathrm{\text{CEL}}(\mu ,\sigma ,\alpha ,\nu )$, where CEL indicates the censored exponentiated Laplace model.

3. If $\alpha =1$ and $\beta =0$, then $X\sim \mathrm{\text{CL}}(\mu ,\sigma ,\nu )$, that is, the censored Laplace model follows.

The estimates of the parameters of the model can be obtained by using maximum likelihood method, where the log-likelihood function is given by

$$\begin{array}{rl}\ell (\theta \theta ;\mathbf{X})& amp;=\alpha {\sum }_0\log {\mathcal{F}}_{BSL}(y_k;\alpha ,\beta )+n_1\log \alpha -n_1\log \sigma -n_1\log 4(1+360{\beta }^2)\\ & amp;+{\sum }_1\log \{1+{(1-\beta y_i^3)}^2\}-{\sum }_1|y_i|+(\alpha -1){\sum }_1\log {\mathcal{F}}_{BSL}(y_i;\mu ,\sigma ,\beta )\end{array}$$

where $y_k=(k-\mu )/\sigma$, $y_i=(x_i-\mu )/\sigma$, and ${\sum }_1$ and ${\sum }_0$ are the sum over censored individuals and uncensored individuals, respectively, and $n_1$ is the number of uncensored individuals.

In order to identify atypical observations and/or model misspecification, we analyzed the transformation of the martingale residual, $r_{M{T}_i}$, proposed by Barros et al. (2010). These residuals are defined by

$$r_{MTi}=\text{s}gn(r_{Mi})\sqrt{-2[r_{Mi}+{\rho }_i\log ({\rho }_i-r_{M_i})]},i=1,\dots ,n$$

where $r_{Mi}={\rho }_i+logS(x_i;\widehat{𝛉})$ is the martingal residual proposed by Ortega et al. (2003), where ${\rho }_i=0,1$ indicates whether the $i$th observation is censored or not, respectively, sgn$\left(r_{Mi}\right)$ denotes the sign of $r_{Mi}$ and $S(y_i;\widehat{𝛉})=P_{\widehat{𝛉}}(Y_i>y_i)$ represents the survival function evaluated at $y_i$, where $\widehat{𝛉}$ are the MLE for $𝛉$.

SIMULATION STUDY

In this section, a small Monte Carlo simulation is conducted to evaluate the maximum likelihood estimation of the BSL$(\mu ,\sigma ,\beta )$ distribution parameters. We considered samples sizes of $n=100$, $200$, $500$ and $1,000$; the true values of the parameter were $\beta =-1.5,$ $-0.15$, $0$, $0.15$ and $2.0$ while, without loss of generality we considered values for the location and scale parameter set to $\mu =0$ and $\sigma =1$, respectively. We generated $2,000$ artificial samples and to examine the performance of the MLEs, for each sample size and for each estimate, which we denote by ${\hat{\theta }}_j$, we computed the bias, denoted by $𝔼\left[{\hat{\theta }}_j\right]-{\theta }_j$ and the mean square error (MSE), defined by $𝔼\left[{({\hat{\theta }}_j-\theta )}^2\right]$, for $j=1,2,3$. Table I presents the bias and the MSE of the MLEs of the parameters $\mu$, $\sigma$ and $\beta$.

Table I
Bias and MSE of the MLE of unknonw parameter for the BSL distribution.

𝛽	$n$	Estimates of $\mu$			Estimates of $\sigma$			Estimates of $\beta$
		Bias	MSE		Bias	MSE		Bias	MSE
−1.50	100	−0.0816	0.2139		−0.0121	0.0935		0.0375	0.0827
	200	−0.0533	0.1522		−0.0101	0.0717		0.0109	0.0561
	500	−0.0149	0.0813		−0.0084	0.0309		0.0085	0.0406
	1,000	−0.0081	0.0412		−0.0055	0.0092		0.0024	0.0205
−0.15	100	−0.0567	0.1893		−0.0053	0.0848		0.0158	0.0619
	200	−0.0381	0.1066		−0.0051	0.0629		0.0048	0.0302
	500	−0.0051	0.0450		−0.0038	0.0114		0.0038	0.0265
	1,000	−0.0029	0.0283		−0.0027	0.0097		0.0009	0.0194
0	100	0.0316	0.1401		−0.0858	0.0759		0.0071	0.0230
	200	0.0118	0.0574		−0.0460	0.0503		0.0034	0.0096
	500	0.0007	0.0121		−0.0177	0.0111		0.0012	0.0039
	1.000	0.0019	0.0045		−0.0085	0.0041		0.0005	0.0010
0.15	100	0.0212	0.1716		0.0011	0.0917		−0.0058	0.0736
	200	0.0091	0.0916		0.0007	0.0678		−0.0011	0.0313
	500	0.0042	0.0434		0.0002	0.0213		−0.0004	0.0114
	1,000	0.0034	0.0209		0.0001	0.0086		−0.0003	0.0093
2.00	100	0.0973	0.2315		0.0132	0.1055		−0.0410	0.0882
	200	0.0575	0.1589		0.0097	0.0721		−0.0179	0.0441
	500	0.0238	0.0721		0.0082	0.0294		−0.0075	0.0209
	1,000	0.0097	0.0294		0.0035	0.0108		−0.0017	0.0109

We can see from Table I that MSE of the MLEs are decreasing in $n$ and become zero as $n$ goes to infinity. Also, we can see that the absolute value of the bias of the MLEs are decreasing and tends to zero, and so the estimates are asymptotically unbiased. For $\beta >0$ ($\beta <0$) the MLE $\hat{\beta }$ underestimates (overestimates) the parameter $\beta$. The MLE $\hat{\mu }$ underestimates (overestimates) the parameter $\mu$ when $\beta <0$ ($\beta >0$). For $\beta >0$ ($\beta <0$) the MLE $\hat{\sigma }$ overestimates (underestimates) the parameter $\sigma$.

ILLUSTRATIONS

In this section, we apply the proposed models to some real data sets. We use the statistical software R version 3.5.3 with the package maxLike for maximizing the corresponding likelihood functions (R Development Core Team 2018). For comparing purposes of various models, the AIC (Akaike 1974), BIC (Schwarz 1978) and CAIC (Bozdogan 1987) information criteria were used.

Illustration 1: Old faithful geyser data

The second data set concerns to measurements about wait times between the eruptions (in minutes) of the Old Faithful geyser in Yellowstone National Park, Wyoming, U.S. More information of these data can be seen in Azzalini & Bowman (1990) who takes a look at this set of data. In addition to the EBSL and BSL models, we also fit the ASL and the ABSL models. The results of the fit of these models can be seen in Table II. The standard errors of the estimators were obtained by using the observed information matrix, and again, to compare the fitted models, we use the AIC, BIC and CAIC criteria. According to any of these criteria, the best model is EBSL. Figure 5 shows the histogram of the data together with the different fitted models.

Table II
Parameter estimates (S.E) of the fitted models to the old faithful geyser data set.

Estimates	Distribution
	ASL	BSL	ABSL	EBSL
$\hat{\mu }$	69.001(0.682)	68.345(0.740)	68.996(0.623)	64.000(0.529)
$\hat{\sigma }$	4.295(0.160)	1.942(0.046)	2.130(0.089)	2.166(0.053)
$\hat{\alpha }$	-4.009(1.150)	–	-1.216(0.428)	1.772(0.123)
$\hat{\beta }$	–	-0.169(0.022)	-0.175(0.041)	0.169(0.021)
$\hat{\ell }$	-1057.236	-1054.877	-1046.016	-1038.073
AIC	2120.472	2115.754	2100.032	2084.146
BIC	2131.289	2126.571	2112.206	2098.569
CAIC	2139.451	2134.733	2125.338	2109.452

Illustration 2: Stellar abundances data

The third data set is related to measurements for 68 solar-type stars4, which are available in the package astrodatR (Feigelson 2014) of the software R Development Core Team (2018) under the name Stellar abundances. These data were previously analyzed by Mattos et al. (2018) using the Scale Mixture of Skew Normal Censored Regression (SMSNCR) models. We take only the response variable: log N(Be), which represents the log of the abundance of beryllium scaled to Sun’s abundance (i.e. the Sun has $logN\left(Be\right)=0.0$).

In astronomical research, a previously identified sample of objects (stars, galaxies, quasars, $X$-ray sources, etc.) is observed at some new wavebands. According to Feigelson (2014), due to limited sensitivities, some objects may be undetected, leading to upper limits in their derived luminosities. For this dataset we have 12 left-censored data points, i.e. 12 undetected beryllium measurement, that represents 19.35% of observations. Table III presents the ML estimates for the parameters of the censored alpha-skew-Laplace (CASL), censored beta-skew-Laplace (CBSL), censored alpha-beta-skew-Laplace (CABSL) and censored exponentiated beta-skew-Laplace (CEBSL) models, together with their corresponding standard errors. Table III also compares the fit of the four models using the model selection criteria (AIC, BIC and CAIC). Note that the CEBSL model has better fit than the CASL, CBSL and CABSL models. The plots of $r_{M{T}_i}$ with generated confidence envelopes are presented in Figure 6. From this figure, we can see clearly that the CBSL, CABSL and CEBSL models fit better the data than the CASL model, since, in that cases, there are not observations which lie outside the envelopes.

Table III
Parameter estimates (S.E) of the fitted models to the Stellar abundances data set.

Estimates	Distribution
	CASL	CBSL	CABSL	CEBSL
$\hat{\mu }$	0.7493(0.0003)	0.4495(0.0028)	0.5520(0.0066)	1.3891(0.0003)
$\hat{\sigma }$	0.1552(0.0003)	0.0873(0.0001)	0.0792 (0.0001)	0.0068(0.0000)
$\hat{\alpha }$	-1.4180(0.0786)	–	-0.5165(0.2665)	0.0146(0.0000)
$\hat{\beta }$	–	-0.1974 0.0043	-0.1708(0.0098)	0.2682(1.2871)
$\hat{\ell }$	-32.330	-38.577	-38.749	-23.493
AIC	70.660	83.154	85.499	54.986
BIC	77.319	89.813	94.378	63.864
CAIC	80.319	92.813	98.377	67.864

Figure 5
Old faithful geyser data: Density function of the fitted models.

Figure 6
Stellar abundances data. Envelopes of transformed martingale residuals for CASL, CBSL, CABSL and CEBSL models.

CONCLUSIONS

In this paper, two new class of unimodal, bimodal as well as trimodal skew distribution is proposed. The main statistical properties and the problem of estimation of the parameters are studied in details by using maximum likelihood method. The models extend the Laplace distribution to trimodal symmetric and asymmetric case. An extension to censored data is also considered. Furthermore, we have shown that such distributions are more flexible than certain rival models and they fit better to some real data sets.

REFERENCES

• AKAIKE H. 1974. A new look at statistical model identification. IEEE Trans Automat Contr 19(6): 716-723.
• ARNOLD BC, GÓMEZ HW & SALINAS HS. 2009. On multiple constraint skewed models. Statistics 43(3): 279-293.
• AZZALINI A. 1985. A class of distributions which includes the normal ones. Scand J Stat 12(2): 171-178.
• AZZALINI A & BOWMAN AW. 1990. A Look at some data on the old faithful geyser. J R Stat Soc Ser C Appl Stat 39(3): 357-365.
• BARROS M, GALEA M, GONZÁLEZ M & LEIVA V. 2010. Influence diagnostics in the tobit censored response model. Stat Methods Appl 19(3): 379-397. doi:10.1007/ s10260-010-0135-y.
• BOLFARINE H, MARTÍNEZ-FLÓREZ G & SALINAS HS. 2018. Bimodal symmetric-asymmetric power-normal families. Commun Stat Theory Methods 47(2): 259-279.
• BOZDOGAN H. 1987. Model selection and akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika 52(3): 345-370.
• DURRANS SR. 1992. Distributions of fractional order statistics in hydrology. Water Resour Res 28(6): 1649-1655.
• ELAL-OLIVERO D. 2010. Alpha-skew-normal distribution. Proyecciones (Antofagasta) 29(12): 224-240.
• ELAL-OLIVERO D, GÓMEZ HW & QUINTANA FA. 2009. Bayesian modeling using a class of bimodal skew-elliptical distributions. J Stat Plan Inference 139(4): 1484-1492.
• FEIGELSON ED. 2014. astrodatR: Astronomical Data. URL https://cran.r-project.org/web/packages/astrodatR/ R package v. 0.1.
  » https://cran.r-project.org/web/packages/astrodatR/
• GÓMEZ HW, ELAL-OLIVERO D, SALINAS H & BOLFARINE H. 2009. Bimodal extension based on the skew-normal distribution with application to pollen data. Environmetrics 22(1): 50-62.
• GUPTA RD & GUPTA RC. 2008. Analyzing skewed data by power-normal model. Test 17(1): 197-210.
• KIM HJ. 2005. On a class of two-piece skew-normal distribution. Statistics 39(6): 537-553.
• KOTZ S, KOZUBOWSKI T & PODGORKI K. 2001. The Laplace distribution and generalizations: A revisits with applications to communications, economics, Engineering and finance. Springer.
• MA Y & GENTON MG. 2004. Flexible class of skew-symmetric distributions. Scand J Stat 31(3): 459-468.
• MARTÍNEZ-FLÓREZ G, BOLFARINE H & GÓMEZ HW. 2014. Skew-normal alpha-power model. Statistics 48(6): 1414-1428.
• MARTÍNEZ-FLÓREZ G, BOLFARINE H & GÓMEZ HW. 2018. Censored bimodal symmetric - asymmetric families. Stat Interface 11: 237-249.
• MATTOS T, GARAY AM & LACHOS VH. 2018. Likelihood-based inference for censored linear regression models with scale mixtures of skew-normal distributions. J Appl Stat 45(11): 2039-2066.
• ORTEGA EM, BOLFARINE H & PAULA GA. 2003. Influence diagnostics in generalized log-gamma regression models. Comput Stat Data Anal 42(1-2): 165-186.
• PEWSEY A, GÓMEZ HW & BOLFARINE H. 2012. Likelihood-based inference for power distributions. Test 21(4): 775-789.
• R DEVELOPMENT CORE TEAM. 2018. R: A Language and Environment for Statistical Computing. URL http://www.R-project.org ISBN 3-900051-07-0.
  » http://www.R-project.org
• SCHWARZ G. 1978. Estimating the dimension of a model. Ann Stat 6(2): 461-464.
• SHAMS-HARANDI S & ALAMATSAZ MH. 2013. Alpha–skew–Laplace distribution. Stat Probab Lett 83(3): 774-782.
• SOBHAN S, DOOSTPARAST M & JAMALIZADEH A. 2016. The alpha-beta skew normal distribution: properties and applications. Statistics 50(2): 338-349.

Publication Dates

• Publication in this collection
  24 Oct 2022
• Date of issue
  2022

History

• Received
  22 Dec 2019
• Accepted
  2 June 2020