An Extension of the Akash Distribution: Properties, Inference and Application

1. Introduction

The slash distribution is an extension of the normal distribution. Its representation is the quotient between two independent random variables, one normal and the other a power of the uniform distribution. In this way we say that S has a slash distribution if:

$$\begin{array}{ccc}S& =& U_1/U_2,\end{array}$$

(1)

where $U_1\sim N(0,1)$, $U_2\sim Beta(q,1),$$U_1$ is independent of $U_2$ and $q>0$, its representation can be seen in Johnson et al. [1]. This distribution has heavier tails than the normal distribution, that is, it has greater kurtosis. Properties of this family are discussed in Rogers and Tukey [2] and Mosteller and Tukey [3]. The maximum likelihood estimation of the location and scale parameters are discussed in Kafadar [4]. Wang and Genton [5] provide a multivariate version of the slash distribution and a multivariate skew version. Gomez et al. [6] and Gómez and Venegas [7] extend the slash distribution using the family of univariate and multivariate elliptic distributions. This methodology to increase the weight of the queues has also been used in distributions with positive support, for example, Gómez et al. [8] in the Birnbaum-Saunders distribution, Olmos et al. [9,10] in the half-normal and generalized half-normal distributions, Astorga et al. [11] in the Muth power distribution and Rivera et al. [12] in the Rayleigh distribution, among others.

Based on the work of Rivera et al. [12], the scale mixture of Rayleigh (SMR) distribution is introduced. We say that $Y\sim SMR(\theta ,q)$ with $\theta >0$ and $q>0$. Then the probability density function (pdf) of Y is

$$\begin{array}{c}\hfill f_Y(y;\theta ,q)={\displaystyle \frac{qy}{2\theta {\left(\frac{y^2}{2\theta }+1\right)}^{\frac{q}{2}+1}}},y>0.\end{array}$$

(2)

Also, a necessary distribution in the development of this paper is the gamma distribution, whose pdf is given by

$$g(t;a,b)=\frac{b^a}{\Gamma \left(a\right)}{t}^{a-1}{e}^{-bt},$$

(3)

where $a,b,t>0$. Its corresponding cumulative distribution function (cdf) is denoted by:

$$G(z;a,b)={\int }_0^{z}g(t;a,b)dt$$

(4)

Shanker [13] introduced the Akash distribution and applied it to real lifetime data sets from medical science and engineering. Thus, we say that a random variable Y has an Akash distribution (AK) with shape parameter $\theta$ if its pdf is given by

$$f_Y(y;\theta )=\frac{{\theta }^3}{{\theta }^2+2}(1+y^2)exp(-\theta y),$$

(5)

where $\theta ,y>0$ and we denote it by $Y\sim AK\left(\theta \right).$ The parameter $\theta$ is a shape parameter, and if we add a scale parameter the pdf is given by

$$f_Y(y;\sigma ,\theta )=\frac{{\theta }^3}{\sigma ({\theta }^2+2)}(1+y^2/{\sigma }^2)exp(-\theta y/\sigma ),$$

(6)

where $\sigma >0$ is a scale parameter, $\theta >0$ is a shape parameter and we denote it by $Y\sim AK(\sigma ,\theta ).$

Extensions of the AK distribution are carried out by Shanker and Shukla [14,15], among others. Both extensions consider adding a parameter and we will compare them with the new distribution. The two-parameter Akash distribution (TPAD) introduced by Shanker and Shukla [14] has the following pdf:

$$f_Y(y;\theta ,\alpha )=\frac{{\theta }^3}{\alpha {\theta }^2+2}(\alpha +y^2)exp(-\theta y),$$

(7)

where $\theta ,\alpha ,y>0$ and we denote it by $Y\sim TPAD(\theta ,\alpha ).$

The power Akash distribution (PAD), introduced by Shanker and Shukla [15], has the following pdf:

$$f_Y(y;\theta ,\alpha )=\frac{\alpha {\theta }^3}{{\theta }^2+2}(1+\alpha y^{2\alpha })y^{\alpha -1}exp(-\theta y^{\alpha }),$$

(8)

where $\theta ,\alpha ,y>0$ and we denote it by $Y\sim PAD(\theta ,\alpha ).$

The main objective of this paper is to introduce an extension of the AK distribution given in (6), making use of the slash methodology, in order to obtain a new distribution with greater kurtosis to be able to accomcodate outliers.

The paper evolves as follows: In Section 2 we deliver the new distribution and present its properties. In Section 3 we perform inference using the method of moments and maximum likelihood via the EM algorithm, a simulation study is also carried out. In Section 4 we apply the distribution to a real data set and compare it with other extensions of the AK distribution. In Section 5 we provide some conclusions.

2. New density and its properties

In this section we introduce the representation, density and properties of the new distribution.

2.1. Representation

The representation of this new distribution is given by

$$X=\frac{Y}{Z},$$

(9)

where $Y\sim AK\left(\theta \right)$, $Z\sim Beta(q,1)$, Y and Z are independent random variables with $\theta ,q>0$. We name the distribution of X slash AK (SAK) and denote it by $X\sim SAK(\theta ,q)$.

2.2. Density function

The following Proposition shows the pdf of the SAK distribution is generated using the representation given in (9).

Proposition 1.

Let $X\sim SAK(\theta ,q)$. Then, the pdf of X is given by

$$f_X(x;\theta ,q)=\frac{q^2\Gamma \left(q\right)x^{-(q+1)}}{({\theta }^2+2){\theta }^q}\left\{{\theta }^{2}G(\theta x;q+1,1)+(q+1)(q+2)G(\theta x;q+3,1)\right\},$$

(10)

where $\theta ,q,x>0$ and G is the cdf of the gamma distribution given in $\left(4\right)$.

Proof. 

Using the representation given in (9) and procedures based on the Jacobian method, we get the result. □

Observation 1.

Table 1 and Figure 1 show that as the value of the parameter q diminishes, the weight of the right tail increases.

In particular, Table 1 compares $P(X>x)$ in the AK and SAK distributions for different values of x.

2.3. Properties

The following Proposition gives the cdf in closed form. It depends on G, which is the cdf of the gamma distribution given in $\left(4\right)$.

Proposition 2.

Let $X\sim SAK(\theta ,q)$. Then, the cdf of X is given by

$$F_X(x;\theta ,q)=\frac{({\theta }^2+2G(\theta x;3,1)){\left(\theta x\right)}^q-{\theta }^{3}q\Gamma \left(q\right)G(\theta x;q,1)-\Gamma (q+3)G(\theta x;q+3,1)}{({\theta }^2+2){\left(\theta x\right)}^q},$$

(11)

where $\theta ,q,x>0$ and G is given in $\left(4\right)$.

Proof. 

The result follows from a direct application of the definition of a cdf. □

2.3.1. Reliability analysis

The reliability function $r\left(t\right)=1-F\left(t\right)$ and the hazard function $h\left(t\right)=\frac{f\left(t\right)}{r\left(t\right)}$ of the SAK distribution are given in the following corollary.

Corollary 1.

Let $T\sim SAK(\theta ,q)$. Then, the $r\left(t\right)$ and $h\left(t\right)$ of T are given by

1. 

$r\left(t\right)=1-\frac{({\theta }^2+2G(\theta t;3,1)){\left(\theta t\right)}^q-{\theta }^{3}q\Gamma \left(q\right)G(\theta t;q,1)-\Gamma (q+3)G(\theta t;q+3,1)}{({\theta }^2+2){\left(\theta t\right)}^q},$

2. 

$h\left(t\right)=\frac{q^2\Gamma \left(q\right)\left({\theta }^{2}G(\theta t;q+1,1)+(q+1)(q+2)G(\theta t;q+3,1)\right)}{t\left(2(1-G(\theta t;3,1)){\left(\theta t\right)}^q+{\theta }^{3}qGamma\left(q\right)G(\theta t;q,1)-\Gamma (q+3)G(\theta t;q+3,1)\right)},$

where $\theta ,q>0$.

In Figure 2, we introduce the Hazard function of the SAK distribution for several values of $\sigma$ and q.

2.3.2. Right tail of the SAK distribution

According to Rolski et al. [16] a distribution has a heavy right tail if

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; sup}\left(-\frac{logr\left(t\right)}{t}\right)=0.\end{array}$$

The following result shows that the SAK distribution is heavy-tailed.

Proposition 3.

The distribution of the random variable $T\sim SAK(\theta ,q)$ is heavy-tailed.

Proof. 

Applying L’Hospital’s rule twice we have,

$$\begin{array}{ccc}\hfill \underset{t\to \infty }{lim\; sup}\left(-\frac{logr\left(t\right)}{t}\right)& =& \underset{t\to \infty }{lim\; sup}\left(\frac{f_T(t;\sigma ,q)}{1-F_T(t;\sigma ,q)}\right)\hfill \\ & =& \underset{t\to \infty }{lim\; sup}\left(\frac{q+1}{t}-\frac{{\theta }^{3}g(\theta t;q+1,1)+(q+1)(q+2)\theta g(\theta t;q+3,1)}{{\theta }^{2}G(\theta t,q+1,1)+(q+1)(q+2)G(\theta t,q+3,1)}\right)\hfill \\ & =& 0.\hfill \end{array}$$

□

The following Proposition shows that the SAK distribution can be represented as a scale mixture between the AK and Beta distributions.

Proposition 4.

If $X|Z=z\sim AK(z^{-1},\theta )$ and $Z\sim Beta(q,1)$ then $X\sim SAK(\theta ,q)$.

Proof. 

The marginal pdf of X is given by

$$\begin{array}{c}\hfill f_X(x;\theta ,q)={\int }_0^1{f}_{X|Z}\left(x\right)f_Z\left(z\right)dz=\frac{{\theta }^3}{{\theta }^2+2}{\int }_0^{1}z(1+z^2{x}^2)exp(-\theta zx)q{z}^{q-1}dz,\end{array}$$

and using $\left(3\right)$, $\left(4\right)$ and $\left(5\right)$ , this result is obtained. □

The following result shows that when the parameter q tends to infinity, the AK distribution is obtained.

Proposition 5.

Let $X\sim SAK(\theta ,q)$. If $q\to \infty$ then X converges in law to a random variable $Y\sim AK\left(\theta \right)$.

Proof. 

Using its representation $X=\frac{Y}{Z}$ we analyze the convergence of this quotient, where $Y\sim AK\left(\theta \right)$ and $Z\sim$Beta$(q,1)$. In the Beta($q,1$) distribution we have, $\mathbb{E}\left[Z\right]=\frac{q}{1+q}$ and $Var\left[Z\right]=\frac{q}{(q+2){(q+1)}^2}$. Then, applying Chebychev’s inequality for Z, we have $\forall ϵ>0$

$$P\left[|Z-\mathbb{E}[Z\left]\right|>ϵ\right]\le \frac{Var\left(Z\right)}{{ϵ}^2}=\frac{q}{(q+2){(q+1)}^2{ϵ}^2}.$$

(12)

If $q\to \infty$ then the right hand side of (12) tends to zero, i.e. $W=Z-\mathbb{E}\left[Z\right]$ converges in probability to 0. Also $\mathbb{E}\left[Z\right]=\frac{q}{1+q}⟶1,q\to \infty$, then we have,

$$Z=W+\mathbb{E}\left[Z\right]\stackrel{\mathcal{P}}{⟶}1,q\to \infty .$$

Since $Y\sim AK\left(\theta \right)$, applying Slutsky’s Lemma to $X=\frac{Y}{Z}$, we have

$$X\stackrel{\mathcal{L}}{⟶}Y\sim AK\left(\theta \right),q\to \infty .$$

Thus, for increasing values of q, X converges in law to a $AK\left(\theta \right)$ distribution. □

2.3.3. Moments

In this subsection we obtain the moments of the SAK distribution. To achieve this aim, the next lemma will be useful.

Lemma 1.

Let $Y\sim AK(\sigma ,\theta )$ with $\sigma ,\theta >0$. For $r>0$, $\mathbb{E}\left[Y^r\right]$ exists if and only if $q>r$ and in this case

$$\begin{array}{ccc}\hfill \mathbb{E}\left[Y^r\right]& =& {\displaystyle \frac{{\sigma }^r\left(r!{\theta }^2+(r+2)!\right)}{{\theta }^r({\theta }^2+2)}}.\hfill \end{array}$$

(13)

Proof. 

The r-th moment of the random variable $V\sim AK\left(\theta \right)$ is given by Shanker [13], which is $\mathbb{E}\left(V^r\right)={\displaystyle \frac{r!{\theta }^2+(r+2)!}{{\theta }^r({\theta }^2+2)}}$, then calculating the r-th moment of the random variable $Y=\sigma V$, where $\sigma$ is a parameter of scale, the result is obtained. □

The moments of a SAK distribution are given in the following Proposition 6,

Proposition 6.

Let $X\sim SAK(\theta ,q)$ with σ and $q>0$. For $r>0$, $\mathbb{E}\left[X^r\right]$ exists if and only if $q>r$ and in this case

$$\begin{array}{ccc}\hfill {\mu }_r=\mathbb{E}\left[X^r\right]& =& {\displaystyle \frac{q\left(r!{\theta }^2+(r+2)!\right)}{{\theta }^r({\theta }^2+2)(q-r)}}.\hfill \end{array}$$

(14)

Proof. 

Using the representation given in the Proposition 4 and by Lemma 1, we get

$$\begin{array}{ccc}\hfill {\mu }_r=\mathbb{E}\left[X^r\right]& =& \mathbb{E}\left[\mathbb{E}\left(X^r|Z\right)\right]=\mathbb{E}\left[{\displaystyle \frac{Z^{-r}\left(r!{\theta }^2+(r+2)!\right)}{{\theta }^r({\theta }^2+2)}}\right]\hfill \\ & =& {\displaystyle \frac{r!{\theta }^2+(r+2)!}{{\theta }^r({\theta }^2+2)}}\mathbb{E}\left[Z^{-r}\right]={\displaystyle \frac{r!{\theta }^2+(r+2)!}{{\theta }^r({\theta }^2+2)}}{\int }_0^{1}q{z}^{q-r-1}dz.\hfill \end{array}$$

Solving the integral gives the result. □

From Proposition 6, the explicit expression of the noncentral moments, ${\mu }_r=\mathbb{E}\left[X^r\right]$, for $r=1,2,3,4$ and the variance of $X\sim SAK(\theta ,q)$, $Var\left(X\right)$ follow.

Corollary 2.

Let $X\sim SAK(\theta ,q)$ with θ and $q>0$. From (14), the following noncentral moments and the variance of X, $Var\left(X\right)$, are obtained

$$\begin{array}{ccc}& & {\mu }_1=\frac{q{\kappa }_6}{\theta {\kappa }_2(q-1)},q>1,{\mu }_2=\frac{2q{\kappa }_{12}}{{\theta }^2{\kappa }_2(q-2)},q>2,\hfill \\ \\ & & {\mu }_3=\frac{6q{\kappa }_{20}}{{\theta }^3{\kappa }_2(q-3)},q>3,{\mu }_4=\frac{24q{\kappa }_{30}}{{\theta }^4{\kappa }_2(q-4)},q>4,\hfill \end{array}$$

$$\begin{array}{ccc}& & Var\left(X\right)=\frac{q\left[2{\kappa }_{12}{\kappa }_2{(q-1)}^2-q{\kappa }_6^2(q-2)\right]}{{\theta }^2{\kappa }_2^2{(q-1)}^2(q-2)},q>2.\hfill \end{array}$$

where ${\kappa }_i={\theta }^2+i.$

Remark 1.

Note that when $q\to \infty$, $Var\left(X\right)\to \frac{{\theta }^4+16{\theta }^2+12}{{\theta }^2{({\theta }^2+2)}^2}$, which is the variance of an $AK\left(\theta \right)$ distribution.

The next Corollary gives the asymmetry coefficient, $\sqrt{{\beta }_1}$, of a $SAK(\theta ,q)$ model.

Corollary 3.

Let $X\sim SAK(\theta ,q)$ with $\theta >0$ and $q>3$. Then the skewness coefficient of X is:

$$\begin{array}{ccc}\hfill \sqrt{{\beta }_1}& =& {\displaystyle \frac{2\sqrt{q-2}\left[3{\kappa }_{20}{\kappa }_2^2{(q-1)}^3(q-2)-3q{\kappa }_2{\kappa }_6{\kappa }_{12}{(q-1)}^2(q-3)+q^2{\kappa }_6^3(q-2)(q-3)\right]}{\sqrt{q}(q-3){\left[2{\kappa }_2{\kappa }_{12}{(q-1)}^2-q(q-2){\kappa }_6^2\right]}^{3/2}}}\hfill \end{array}$$

Proof. 

Recall that

$$\begin{array}{c}\hfill \sqrt{{\beta }_1}={\displaystyle \frac{\mathbb{E}\left[{(X-\mathbb{E}\left(X\right))}^3\right]}{{\left(Var\left(X\right)\right)}^{3/2}}}={\displaystyle \frac{{\mu }_3-3{\mu }_1{\mu }_2+2{\mu }_1^3}{{({\mu }_2-{\mu }_1^2)}^{3/2}}},\end{array}$$

where ${\mu }_1$, ${\mu }_2$ and ${\mu }_3$ were given in Corollary 2. □

Also, the kurtosis coefficient, ${\beta }_2$, of a $SAK(\theta ,q)$ distribution is given in the following Corollary.

Corollary 4.

Let $X\sim SAK(\theta ,q)$ with $\theta >0$ and $q>4$. Then the kurtosis coefficient of X is

$$\begin{array}{ccc}\hfill {\beta }_2& =& {\displaystyle \frac{3(q-2)\left(8{\kappa }_2^3{\kappa }_{30}{q}_1-8q{\kappa }_6{\kappa }_{20}{\kappa }_2^2{q}_2+4q^2{\kappa }_6^2{\kappa }_{12}{\kappa }_2{q}_3-q^3{\kappa }_6^4{q}_4\right)}{q(q-3)(q-4){\left[2{\kappa }_{12}{\kappa }_2{(q-1)}^2-q{\kappa }_6^2(q-2)\right]}^2}}.\hfill \end{array}$$

where $q_1={(q-1)}^4(q-2)(q-3)$, $q_2={(q-1)}^3(q-2)(q-4)$, $q_3={(q-1)}^2(q-3)(q-4)$ and $q_4=(q-2)(q-3)(q-4).$

Proof. 

Recall that

$$\begin{array}{c}\hfill {\beta }_2={\displaystyle \frac{\mathbb{E}\left[{(X-\mathbb{E}\left(X\right))}^4\right]}{{\left(Var\left(X\right)\right)}^2}}={\displaystyle \frac{{\mu }_4-4{\mu }_1{\mu }_3+6{\mu }_1^2{\mu }_2-3{\mu }_1^4}{{({\mu }_2-{\mu }_1^2)}^2}},\end{array}$$

where ${\mu }_1$, ${\mu }_2$, ${\mu }_3$, and ${\mu }_4$ were given in Corollary 2. □

Remark 2.

It can be verified that for $q\to \infty$ skewness and kurtosis coefficients converge to $\frac{2({\theta }^6+30{\theta }^4+36{\theta }^2+24)}{{({\theta }^4+16{\theta }^2+12)}^{3/2}}$ and $\frac{3(3{\theta }^8+128{\theta }^6+408{\theta }^4+576{\theta }^2+240)}{{({\theta }^4+16{\theta }^2+12)}^2}$ respectively, which coincide with the corresponding coefficients for the AK(θ) distribution (see Shanker, 2015).

The results of Table 2 show that the values of the skewness and kurtosis coefficients depend on the parameters $\theta$ and q and that as q decreases, the skewness and kurtosis coefficients increase. On the other hand, as q increases, the skewness and kurtosis coefficients are those of the AK($\theta$) distribution (Proposition 5).

3. Inference

In this section we study the estimation the parameters by the method of moments and ML via the EM algorithm. We also carry out some simulations to study the behavior of the ML estimators.

3.1. Method of moment estimators

Let $X_1,...,X_n$ be a random sample from $X\sim SAK(\theta ,q)$. Consider the first two sample moments, denoted by $\overline{X}=\frac{{\sum }_{i=1}^n{X}_i}{n}$ and $\overline{X^2}=\frac{{\sum }_{i=1}^n{X}_i^2}{n}$, respectively.

Proposition 7.

Given $X_1,\dots ,X_n$ a random sample from $X\sim SAK(\theta ,q)$ with $q>2$, the moment method estimators of θ and q are

$$\begin{array}{ccc}& & {\widehat{q}}_M={\displaystyle \frac{\overline{X}{\widehat{\theta }}_M({\widehat{\theta }}_M^2+2)}{{\widehat{\theta }}_M({\widehat{\theta }}_M^2+2)\overline{X}-{\widehat{\theta }}_M^2-6}},\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & \overline{X^2}{\widehat{\theta }}_M\left[2({\widehat{\theta }}_M^2+6)-{\widehat{\theta }}_M\overline{X}({\widehat{\theta }}_M^2+2)\right]-2\overline{X}({\widehat{\theta }}_M^2+12)=0,\hfill \end{array}$$

(16)

where (16) must be solved numerically to obtain ${\widehat{\theta }}_M$. Then ${\widehat{\theta }}_M$ must be replaced in (15) to get ${\widehat{q}}_M$.

Proof. 

Consider the method of moment equations

$$\begin{array}{ccc}\hfill \mathbb{E}\left[X\right]=\frac{q({\theta }^2+6)}{\theta ({\theta }^2+2)(q-1)}& =& \overline{X}\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill \mathbb{E}\left[X^2\right]=\frac{2q({\theta }^2+12)}{{\theta }^2({\theta }^2+2)(q-2)}& =& \overline{X^2}\hfill \end{array}$$

(18)

Solving the equation (17) for the parameter q we obtain (15). Then the value of ${\widehat{q}}_M$ is substituted into the equation (18) and the equation given in (16) is obtained. □

3.2. ML estimation

Let $X_1,\dots ,X_n$ be a random sample from $X\sim SAK(\theta ,q)$. Then the log-likelihood function is

$$\begin{array}{ccc}\hfill l(\theta ,q)& =& c(\theta ,q)-(q+1)\sum _{i=1}^{n}log\left(x_i\right)+\sum _{i=1}^{n}log\left[{\theta }^{2}G(\theta x_i;q+1,1)+(q+1)(q+2)G(\theta x_i;q+3,1)\right]\hfill \end{array}$$

where $c(\theta ,q)=2nlog\left(q\right)+nlog(\Gamma \left(q\right))-nlog({\theta }^2+2)-nqlog\left(\theta \right).$ Taking partial derivatives in $l(\theta ,q)$ with respect to $\theta$ and q and setting them equal to zero, we get

$$\begin{array}{ccc}& & \sum _{i=1}^n\frac{2\theta G(\theta x_i;q+1,1)+{\theta }^{2}J(x_i,q+1)+(q+1)(q+2)J(x_i,q+3)}{{\theta }^{2}G(\theta x_i;q+1,1)+(q+1)(q+2)G(\theta x_i;q+3,1)}=\frac{2n\theta }{{\theta }^2+2}+\frac{nq}{\theta },\hfill \\ \hfill \\ & & \sum _{i=1}^n\frac{{\theta }^{2}H(x_i;q+1)+(2q+3)G(\theta x_i;q+3,1)+(q+1)(q+2)H(x_i;q+3)}{{\theta }^{2}G(\theta x_i;q+1,1)+(q+1)(q+2)G(\theta x_i;q+3,1)}=\eta (\theta ,q)-\sum _{i=1}^{n}log\left(x_i\right),\hfill \end{array}$$

where $J(x_i,m)=x_{i}g(\theta x_i;m,1)$, $H(x_i;v)={\int }_0^{\theta x_i}log\left(t\right)g(t;v,1)dt-\psi \left(v\right)G(\theta x_i;v,1)$ and $\eta (\theta ,q)=\frac{2n}{q}+n(\psi \left(q\right)-log\left(\theta \right)).$ Solving numerically this system of equations to find the ML estimates may be a difficult task due to the functions it involves. However, an EM algorithm can be implemented (see Dempster et al. [17]) to obtain the ML estimates. The following subsection is dedicated to achieving this goal.

3.3. EM Algorithm

An alternative stochastic representation for the SAK model is given by

$$\begin{array}{cc}\hfill X_i\mid U_i=u_i,Z_i=z_i& \sim G(1+2u_i,\theta z_i),\hfill \\ \hfill U_i& \sim Bern\left(\frac{2}{{\theta }^2+2}\right),\hfill \\ \hfill Z_i& \sim Beta\left(q,1\right).\hfill \end{array}$$

where $U_i$ and $Z_i$, $i=1,\dots ,n,$ represent non-observable variables. This representation can be used for an alternative estimation procedure based on the EM algorithm (Dempster et al. [17]). In this context, the observed data are given by ${\mathit{D}}_o={\mathit{x}}^{\top }$, where ${\mathit{x}}^{\top }=(x_1,\dots ,x_n)$. The vectors ${\mathit{z}}^{\top }=(z_1,\dots ,z_n)$ and ${\mathit{u}}^{\top }=(u_1,\dots ,u_n)$ are the latent variables and the vector ${\mathit{D}}_c={({\mathit{x}}^{\top },{\mathit{z}}^{\top },{\mathit{u}}^{\top })}^{\top }$ are the complete data. Note that the joint distribution of $(X_i,U_i,Z_i)$ is given by

$$\begin{array}{cc}\hfill f(x_i,u_i,z_i)& =f(x_i\mid u_i,z_i)\times f\left(u_i\right)\times f\left(z_i\right)\hfill \\ \hfill & =\frac{{\left(\theta z_i\right)}^{1+2u_i}}{\Gamma (1+2u_i)}{x}_i^{2u_i}{e}^{-\theta z_i{x}_i}\times {\left(\frac{2}{{\theta }^2+2}\right)}^{u_i}{\left(\frac{{\theta }^2}{{\theta }^2+2}\right)}^{1-u_i}\times q{z}_i^{q-1}\hfill \\ \hfill & =\frac{q{\theta }^3{z}_i^{2u_i+q}{2}^{u_i}}{({\theta }^2+2)\Gamma (1+2u_i)}{x}_i^{2u_i}{e}^{-\theta z_i{x}_i}.\hfill \end{array}$$

Therefore, up to a constant that does not depend on the vector of parameters $\psi =(\theta ,q)$, the complete log-likelihood function for the model is given by

$$\begin{array}{c}\hfill {\ell }_c(\psi ;D_c)=n\left[logq+3log\theta -log({\theta }^2+2)\right]+\sum _{i=1}^n\left[qlog{z}_i-\theta x_i{z}_i\right].\end{array}$$

With this, the expected value of ${\ell }_c(\psi ;D_c)$, given the observed data, is

$$\begin{array}{c}\hfill Q(\psi \mid {\psi }^{\left(k\right)})=n\left[logq+3log\theta -log({\theta }^2+2)\right]+\sum _{i=1}^n\left[q{\widehat{\kappa }}_i^{\left(k\right)}-\theta x_i{\widehat{z}}_i^{\left(k\right)}\right],\end{array}$$

where ${\widehat{z}}_i^{(k)}=\mathbb{E}(Z_i\mid x_i,\psi ={\widehat{\psi }}^{(k)})$ and ${\widehat{\kappa }}_i^{(k)}=\mathbb{E}(log{Z}_i\mid x_i,\psi ={\widehat{\psi }}^{(k)})$. Note that

$$f(z_i,u_i\mid x_i)\propto \underset{Z_i\mid u_i,x_i\sim T{G}_{(0,1)}(2u_i+q+1,\theta x_i)}{\underbrace{\frac{{\left(\theta x_i\right)}^{2u_i+q+1}}{\Gamma (2u_i+q+1)}\frac{z_i^{(2u_i+q+1)-1}{e}^{-\theta x_i{z}_i}}{G(1;2u_i+q+1,\theta x_i)}}}\times \underset{U_i\mid x_i\sim Bern\left({\nu }_i\right)}{\underbrace{\frac{\Gamma (2u_i+q+1)}{\Gamma (2u_i+1)}{\left(\frac{2}{{\theta }^2}\right)}^{u_i}G(1;2u_i+q+1,\theta x_i)}},$$

(19)

where ${\nu }_i=\Gamma (q+3)G(\theta x_i;q+3)/[{\theta }^2\Gamma (q+1)G(\theta x_i;q+1)+\Gamma (q+3)G(\theta x_i;q+3)]$, $G(x;a)={\int }_0^x\frac{1}{\Gamma \left(a\right)}{t}^{a-1}{e}^{-t}dt$ is the cdf for the gamma model and ${\mathrm{TG}}_{(0,1)}(a,b)$ denotes de gamma distribution with shape a and rate b truncated in the interval (0,1). Therefore, using properties of conditional expectations, we have $\mathbb{E}(Z_i\mid x_i)=\mathbb{E}\left(\mathbb{E}(Z_i\mid U_i,x_i)\mid x_i\right)$ and by (19) such expectations are simple to be computed. In a similar manner, we can compute $\mathbb{E}(log{Z}_i\mid x_i)$, obtaining as results

$$\begin{array}{cc}\hfill \mathbb{E}(Z_i\mid x_i)& =\frac{{\nu }_i(q+3)G(\theta x_i,q+4)}{\theta x_{i}G(\theta x_i,q+3)}+\frac{(1-{\nu }_i)(q+1)G(\theta x_i,q+2)}{\theta x_{i}G(\theta x_i,q+1)},\hfill \\ \hfill \mathbb{E}(log{Z}_i\mid x_i)& =\frac{{\nu }_i}{\Gamma (q+3)G(1;q+3,\theta x_i)}{\int }_0^{\theta x_i}log\left(\frac{w_i}{\theta x_i}\right)w_i^{q+2}{e}^{-w_i}d{w}_i\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & +\frac{(1-{\nu }_i)}{\Gamma (q+1)G(1;q+1,\theta x_i)}{\int }_0^{\theta x_i}log\left(\frac{w_i}{\theta x_i}\right)w_i^q{e}^{-w_i}d{w}_i.\hfill \end{array}$$

(21)

Therefore, the kth iteration of the algorithm comprises the following steps:

• E-step: Given ${\widehat{\theta }}^{(k-1)}$ and ${\widehat{q}}^{(k-1)}$, for $i=1,\dots ,n$ compute ${\widehat{z}}_i^{\left(k\right)}$ and ${\widehat{\kappa }}_i^{\left(k\right)}$ using equations (20) and (21), respectively.
• M1-step: Update ${\widehat{q}}^{\left(k\right)}$ as

  $${\widehat{q}}^{\left(k\right)}=\frac{-n}{{\sum }_{i=1}^n{\widehat{\kappa }}_i^{\left(k\right)}}.$$
• M2-step: Update ${\widehat{\theta }}^{\left(k\right)}$ as the solution for the non-linear equation

  $$\frac{3}{\theta }-\frac{2\theta }{{\theta }^2+2}=\frac{1}{n}\sum _{i=1}^n{x}_i{\widehat{z}}_i^{\left(k\right)}.$$

The E, M1 and M2 steps are repeated until convergence is obtained, i.e. until the maximum distance between the estimate obtained in two consecutive iterations is less than a specified value.

In the following subsection we run some simulations to study the behavior of the ML estimators.

3.4. Simulation study

Table 3 shows the empirical bias (bias), the average of the standard errors (SE), the root of the empirical mean squared error (RMSE), and the 95% coverage probability (CP) based on the asymptotic distribution for the ML estimators of the parameters. Table 3 shows that the performance of the estimator improves as n increases.

4. Application

In this section we analyze a real data set showing that the SAK distribution can be more appropriate than other commongly used distributions to model heavy right-tailed data. The data correspond to plasma beta-carotene levels (ng/ml) of 314 patients. This data set contains 14 variables and is available online at “http://Lib.stat.cmu.edu/datasets/Plasma Retinol”. In this study, we consider the variable Betaplasma. The medical interest in this variable comes from the fact that low levels of plasma beta-carotene may be associated with higher risk of developing certain types of cancer. In Table 4 we present some descriptive statistics including the sample skewness, $b_1$, and sample and kurtosis $b_2$. We may observe high kurtosis in this data set.

The moment estimates for the parameters of the SAK distribution are ${\widehat{\theta }}_M=0.025$ and ${\widehat{q}}_M=2.810$. These estimates are useful starting values, required to implement maximum likelihood estimation using numerical methods. Table 5 shows the ML estimates for the parameters of the PAD, SMR and SAK models. For each model we report the value of the log-likelihood. It can be seen that the SAK model presented a larger value of log-likelihood than the other models.

In order to compare the fit of the distributions, we considered the usual Akaike information criterion (AIC), introduced by Akaike [18], and the Bayesian information criterion (BIC), proposed by Schwarz [19]. It is known that AIC=$2k-2loglik$ and BIC=$klogn-2loglik$ where k is the number of parameters in the model, n is the sample size and $loglik$ is the maximized value of the log-likelihood function. Table 6 shows the AIC and BIC for each model, indicating that the SAK distribution leads to a better fit than the other distributions. Figure 3 presents the histogram for the data together with the fitted densities.

We also calculated the quantile residuals (QR). If the model is appropriate for the data, the QR should be a sample of the standard normal distribution (see Dunn and Smyth, [20]). This assumption can be validated with traditional normality tests, such as the Anderson-Darling (AD), Cramér-von Mises (CVM) and Shapiro-Wilkes (SW) tests. Figure 4 shows the qqplot of the quantile residuals of the three fitted distributions. All three tests suggest that the SAK model provides a better fit for this data set.

5. Discussion

This paper presents an extension of the AK distribution based on the slash methodology. Some properties of this new distribution are derived. It is also compared with two other distributions using a real data set. Estimation is done through ML via the EM algorithm. The new SAK distribution is an alternative to fit heavy-tailed right-skewed data. Additional features of the SAK distribution are:

• The distribution has two representations, one based on the quotient of two independent random variables and another based on a scale mixture between the AK and Beta distributions.
• The pdf, cdf and hazard function of the SAK distribution are explicit and are represented by the cdf of the gamma distribution.
• The distribution has a heavy right tail.
• The distribution contains the AK distribution as a limit, that is, when the parameter q tends to infinity in the distribution SAK, the AK distribution is obtained.
• The moments and the coefficients of skewness and kurtosis are explicit.
• In the application, observing the AIC and BIC and the Anderson-Darling, Cramér-von Mises and Shapiro-Wilkes statistical tests, we may conclude that the SAK distribution fits the Betaplasma data set better than the PAD and SMR distributions, which are also extensions of the AK distribution.