A Power Maxwell Distribution with Heavy Tails and Applications

Abstract

In this paper we introduce a distribution which is an extension of the power Maxwell distribution. This new distribution is constructed based on the quotient of two independent random variables, the distributions of which are the power Maxwell distribution and a function of the uniform distribution (0,1) respectively. Thus the result is a distribution with greater kurtosis than the power Maxwell. We study the general density of this distribution, and some properties, moments, asymmetry and kurtosis coefficients. Maximum likelihood and moments estimators are studied. We also develop the expectation–maximization algorithm to make a simulation study and present two applications to real data.

1. Introduction

A distribution related to the normal distribution is the slash distribution. Its stochastic representation is the quotient between two independent random variables: a normal distribution and a function of the uniform distribution. Thus we say that S has a slash distribution if:

$$\begin{array}{ccc}S& =& Z/U^{\frac{1}{q}},\end{array}$$

(1)

where $Z\sim N(0,1)$, $U\sim U(0,1),$Z is independent of U and $q>0$. When $q=1$ we obtain the slash canonical distribution, and when $q\to \infty$ we obtain the standard normal distribution. The density of the canonical slash is:

$$\begin{array}{ccc}p\left(x\right)& =& \left\{\begin{array}{ccc}\frac{\varphi \left(0\right)-\varphi \left(x\right)}{x^2}& & x\ne 0\\ \frac{1}{2}\varphi \left(0\right)& & x=0\end{array}\right.\end{array},$$

(2)

where $\varphi \left(x\right)$ represents the standard normal density defined in Appendix A (see Johnson et al. [1]). This distribution has much heavier tails than the normal distribution; i.e., it has greater kurtosis. Some properties of this family are discussed by Rogers et al. [2] and Mosteller et al. [3]. The maximum likelihood (ML) estimators of scale and location are discussed in Kafadar [4]. One paper Wang et al. [5] introduced a multivariate version of the slash distribution and a skew multivariate version. Reference Gómez et al. [6] extended the slash distribution using the family of univariate and multivariate elliptical distribution. Reference Gómez et al. [7] used this family to extend the Birnbaum–Saunders distribution. Reference Iriarte et al. [8] and Gómez et al. [9] used this methodology to extend the generalized Rayleigh distribution and Gumbel respectively. Reference Olmos et al. [10] also used this methodology to introduce the modified slashed half-normal distribution. Reference Olmos et al. [11] recently introduced the confluent, hypergeometric slashed-Rayleigh distribution using the same methodology.

The Maxwell distribution was first set out by Maxwell [12], and gave the distribution of velocities among the molecules of a gas. Maxwell’s finding was generalized by Boltzmann [13,14,15], to express the distribution of energies among molecules. It has several statistical applications in the areas of physics, chemistry, and physical chemistry, (see Dunbar [16]). For example, in the context of the kinetic molecular theory of gases, a gas contains a large number of particles in rapid motion. Each particle has a different velocity, and each collision between particles changes the velocities of the particles. An understanding of the properties of the gas requires an understanding of the distribution of particle velocities (which is the Maxwell distribution). In addition to these areas, it also has uses in the theory of relativity. Ideal plasmas, which are ionized gases of sufficiently low density, frequently have particle distributions that follow the Maxwell distribution.

Coraddu et al. [17] discussed the physics of nuclear reactions in stellar plasma by checking with special emphasis on the importance of the velocity distribution of ions. Then they claimed that the properties (density and temperature) of the weak-coupled solar plasma were analyzed, showing that the ion velocities should deviate from the Maxwell distribution and could be better described by a "weakly-non-extensive Tsallis" distribution.

Singh et al. [18] introduced the power Maxwell (PM) distribution; i.e., we say that X has a PM distribution if its density function is:

$$f_X(x;\alpha ,\beta )=\frac{4{\alpha }^{3/2}\beta }{\sqrt{\pi }}{x}^{3\beta -1}{e}^{-\alpha x^{2\beta }},x\ge 0,$$

where $\alpha >0$ is a scale parameter and $\beta >0$ is a shape parameter. We denote it as $X\sim PM(\alpha ,\beta )$, where:

$$E\left(X^r\right)={\displaystyle \frac{2}{{\alpha }^{r/2\beta }\sqrt{\pi }}}\Gamma \left({\displaystyle \frac{3\beta +r}{2\beta }}\right)\mathrm{and}{F}_X\left(x\right)={\int }_0^{\alpha x^{2\beta }}{\displaystyle \frac{u^{3/2-1}{e}^{-u}}{\Gamma (3/2)}}du.$$

The main object of this paper is to introduce an extension of the PM distribution using slash methodology. This new distribution presents a greater kurtosis than the PM distribution, so we can use it to model atypical data. In Section 2, we introduce the slash power Maxwell (SPM) distribution, with its stochastic representation, its distribution and its survival and hazard functions. Then we focus on the moments of the distribution, asymmetry and kurtosis coefficients. We devote Section 3 to the study of some properties of the model, such as the mode, the convergence, and the distribution of the stochastic order. We assign Section 4 to inference, wherein we obtain the moments and ML estimators, and establishing the expectation–maximization (EM) algorithm. In Section 5 we present the simulation study, focusing our attention on parameter recovery and criteria comparison. In Section 6, we present two applications to real data, fitting the PM distribution to the datasets; and finally, in Section 7 we present our conclusions.

2. Probability Density Function

2.1. Stochastic Representation

We consider that the random variable Z has a SPM distribution with parameters $\alpha$, $\beta$ and q if it can be represented as

$$Z=\frac{X}{U^{1/q}},$$

(3)

where $X\sim PM(\alpha ,\beta )$ and $U\sim U(0,1)$ are independent, and $q>0$. We denote this by writing $Z\sim SPM(\alpha ,\beta ,q)$, where q is the kurtosis parameter.

2.2. Density Function

The following proposition presents the probability density function (pdf) distribution that can be generated using (3).

Proposition 1.

Let $Z\sim SPM(\alpha ,\beta ,q)$. Then, the pdf of Z is given by

$$f(z;\alpha ,\beta ,q)=\frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }{z}^{q+1}}G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right),z\ge 0,$$

where $\alpha >0$ is the scale parameter, $\beta >0$ is the shape parameter, $q>0$ is the kurtosis parameter, $\Gamma \left(\alpha \right)$ is the gamma function defined in Appendix A and $G(x;a,b)$ is the cumulative distribution function (CDF) of the gamma distribution also defined in Appendix A.

Proof. 

Using the representation given in (3) and computing the Jacobian transformation, we have that

$$\begin{array}{c}\hfill \left.\begin{array}{ccc}\hfill Z& =& X/U^{1/q}\hfill \\ \hfill W& =& U^{1/q}\hfill \end{array}\right\}\Rightarrow \left.\begin{array}{ccc}\hfill X& =& ZW\hfill \\ \hfill U& =& W^q\hfill \end{array}\right\}\Rightarrow J=\left|\begin{array}{cc}{\displaystyle \frac{\partial x}{\partial z}}& {\displaystyle \frac{\partial x}{\partial w}}\\ {\displaystyle \frac{\partial u}{\partial z}}& {\displaystyle \frac{\partial u}{\partial w}}\end{array}\right|=\left|\begin{array}{cc}w& z\\ 0& q{w}^{q-1}\end{array}\right|=q{w}^q.\end{array}$$

Then, $f_{Z,W}(z,w)=\left|J\right|f_{X,U}(zw,w^q)={\displaystyle \frac{4{\alpha }^{3/2}\beta q}{\sqrt{\pi }}}{z}^{3\beta -1}{w}^{3\beta +q-1}exp\left\{-\alpha {\left(zw\right)}^{2\beta }\right\}$, so that by marginalizing with respect to the random variable W, we obtain the density function corresponding to the random variable Z; namely:

$$f_z(z;\alpha ,\beta ,q)={\displaystyle \frac{4{\alpha }^{3/2}\beta q}{\sqrt{\pi }}}{z}^{3\beta -1}{\int }_0^1{w}^{3\beta +q-1}exp\left\{-\alpha {\left(zw\right)}^{2\beta }\right\}dw,$$

so making the variable change $u=\alpha {\left(zw\right)}^{2\beta }$, the result is obtained.□

Figure 1 shows the density of SPM for some values of parameter q. It reveals that the tails become heavier as q becomes smaller.

In

Table 1. Tails comparison for different values for parameters of the SPM and PM (1.5, 1.5) distributions.

Distributions	$\mathit{P}(\mathit{Z}>1.2)$	$\mathit{P}(\mathit{Z}>1.5)$	$\mathit{P}(\mathit{Z}>1.7)$	$\mathit{P}(\mathit{Z}>1.9)$
PM (1.5, 1.5)	$0.159$	$0.018$	$0.002$	$0.000$
SPM (1.5, 1.5, 10)	$0.282$	$0.075$	$0.026$	$0.009$
SPM (1.5, 1.5, 5)	$0.396$	$0.170$	$0.094$	$0.054$
SPM (1.5, 1.5, 3)	$0.509$	$0.293$	$0.203$	$0.146$
SPM (1.5, 1.5, 1)	$0.754$	$0.617$	$0.545$	$0.488$

we perform a brief comparison, illustrating that the tails of the SPM distribution are heavier than those of the PM distribution, showing $P(Z>z)$ for different values of z in the distributions mentioned. It is clear that the SPM distribution may have much heavier tails than the PM distribution.

Proposition 2.

Let $Z\sim SPM(\alpha ,\beta ,q)$. Then the CDF of Z is given by:

$$F_Z\left(z\right)=G(\alpha z^{2\beta },3/2,1)-{\displaystyle \frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }}}G\left(\alpha z^{2\beta },\frac{q+3\beta }{2\beta },1\right).$$

Proof. 

By integrating by parts the kernel of the density

$$I={\int }_0^z{t}^{-q-1}G\left(\alpha t^{2\beta },\frac{q+3\beta }{2\beta },1\right)dt,$$

using $u=G\left(\alpha t^{2\beta },\frac{q+3\beta }{2\beta },1\right)$ and $v=-t^{-q}/q$ we finally obtain

$$I=-{\displaystyle \frac{z^{-q}}{q}}G\left(\alpha z^{2\beta },\frac{q+3\beta }{2\beta },1\right)+{\displaystyle \frac{1}{q}}{\displaystyle \frac{2\beta {\alpha }^{(3\beta +q)/\left(2\beta \right)}}{\Gamma \left({\displaystyle \frac{3\beta +q}{2\beta }}\right)}}{\int }_0^z{t}^{3\beta -1}{e}^{-\alpha t^{2\beta }}dt,$$

and then, by substituting $\mu =\alpha t^{2\beta }$ in the last integral (i.e., $I_2$), we obtain

$$I_2={\displaystyle \frac{{\alpha }^{q/2\beta }\sqrt{\pi }}{2q\Gamma \left({\displaystyle \frac{3\beta +q}{2\beta }}\right)}}G\left(\alpha z^{2\beta },3/2,1\right).$$

Thus the proof is complete.□

Corollary 1.

Let $Z\sim SPM(\alpha ,\beta ,q)$. Then the survival function and the hazard function of Z are given

$$S\left(z\right)=1+\frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }}{G}_1-G_0,h\left(z\right)={\displaystyle \frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)G_1}{z^{q+1}\left[\sqrt{\pi }{\alpha }^{q/2\beta }+2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)G_1-{\alpha }^{q/2\beta }\sqrt{\pi }{G}_0\right]}},$$

where $G_0=G\left(\alpha z^{2\beta };3/2,1\right)$ and $G_1=G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)$.

In Figure 2, we show the CDF, survival and hazard functions. Additionally, we can see that the curve related to the hazard function is unimodal, and that as "q" grows, the curve has longer tails and extends to a greater range.

Proposition 3.

Let $Z\sim SPM(\alpha ,\beta ,q)$. Then the mode of Z can be obtained by solving:

$${\displaystyle \frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }}}\left[2\alpha \beta z^{2\beta -q-2}{g}_1-(q+1)z^{-q-2}{G}_1\right]=0,$$

where $g_1=g\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)$ and $G_1=G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)$. This equation must be solved numerically.

Proof. 

We set ${\displaystyle \frac{d}{dz}}{f}_Z(z;\alpha ,\beta ,q)=0$, we have

$${\displaystyle \frac{2q\Gamma \left({\displaystyle \frac{q+3\beta }{2\beta }}\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }}}\left\{{\displaystyle \frac{\partial }{\partial z}}\left[{\displaystyle \frac{1}{z^{q+1}}}G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)\right]\right\}=0.$$

Considering

$${\displaystyle \frac{\partial }{\partial z}}G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)=2\alpha \beta z^{2\beta -1}g\left(\alpha z^{2\beta };{\displaystyle \frac{q+3\beta }{2\beta }},1\right),$$

the result is obtained.□

2.3. Distribution Relationships

It is easy to see some special cases that are associated with SPM distribution.

• According to a property that we will discuss later in Section 3, if $q\to \infty$ then $Z\stackrel{\mathbb{D}}{\to }X$, where $X\sim PM(\alpha ,\beta$).
• If $\beta =1$ then $Z=Y$ with $\alpha ={\displaystyle \frac{1}{2}}\sigma$, where Y has slash Maxwell (SM) distribution (Iriarte et al. [8]).
• If $q\to \infty$ and $\beta =1$, then $Z\stackrel{\mathbb{D}}{\to }M$, where M has Maxwell distribution.

A diagram showing this relationship is presented in Figure 3.

2.4. Moments

We present a general formula for the rth moment of the SPM distribution.

Proposition 4.

Let $Z\sim SPM(\alpha ,\beta ,q).$ Then, for $r=1,2,\cdots$, it follows that the rth moment of Z can be written as

$${\mu }_r=E\left(Z^r\right)={\displaystyle \frac{2q}{(q-r){\alpha }^{r/2\beta }\sqrt{\pi }}}\Gamma \left({\displaystyle \frac{3\beta +r}{2\beta }}\right),forq>r.$$

Proof. 

Using the stochastic representation for the distribution given in (3), we have that

$$E\left[{\left({\displaystyle \frac{X}{Y^{1/q}}}\right)}^r\right]=E\left(X^r\right)E\left(Y^{-r/q}\right)={\displaystyle \frac{2}{{\alpha }^{r/2\beta }\sqrt{\pi }}}\Gamma \left({\displaystyle \frac{3\beta +r}{2\beta }}\right){\displaystyle \frac{q}{q-r}},$$

where $E\left(Y^{-r/q}\right)=\frac{q}{q-r}$, and $E\left(X^r\right)$ are the moments of the $PM(\alpha ,\beta )$.□

Corollary 2.

If $Z\sim SPM(\alpha ,\beta ,q)$, then it follows that

$$\begin{array}{ccc}\hfill {\mu }_1=E\left(Z\right)={\displaystyle \frac{q}{q-1}}{k}_{1}forq>1& ,& {\mu }_2=E\left(Z^2\right)={\displaystyle \frac{q}{q-2}}{k}_{2}forq>2,\hfill \\ \hfill {\mu }_3=E\left(Z^3\right)={\displaystyle \frac{q}{q-3}}{k}_{3}forq>3& ,& {\mu }_4=E\left(Z^4\right)={\displaystyle \frac{q}{q-4}}{k}_{4}forq>4,\hfill \end{array}$$

$$Var\left(Z\right)={\displaystyle \frac{q}{q-2}}{k}_2-{\left({\displaystyle \frac{q}{q-1}}{k}_1\right)}^{2}forq>2,$$

where $k_r$ are the r-moments of the PM distribution.

Remark 1.

The asymmetry and kurtosis coefficients can be obtained using:

$$\sqrt{{\beta }_1}={\displaystyle \frac{{\mu }_3-3{\mu }_1{\mu }_2+2{\mu }_1^3}{{[{\mu }_2-{\mu }_1^2]}^{3/2}}}and{\beta }_2={\displaystyle \frac{{\mu }_4-4{\mu }_1{\mu }_3+6{\mu }_2{\mu }_1^2-3{\mu }_1^4}{{[{\mu }_2-{\mu }_1^2]}^2}}.$$

Corollary 3.

Let $Z\sim SPM(\alpha ,\beta ,q)$; then the asymmetry coefficient $\left(\sqrt{{\beta }_1}\right)$ and the kurtosis coefficient $\left({\beta }_2\right)$ for $q>3$ and $q>4$ respectively

$$\sqrt{{\beta }_1}={\displaystyle \frac{{(q-2)}^{1/2}[{(q-1)}^3(q-2)k_3-3q{(q-1)}^2(q-3)k_1{k}_2+2q^2(q-2)(q-3)k_1^3]}{\sqrt{q}(q-3){[{(q-1)}^2{k}_2-q(q-2)k_1^2]}^{3/2}}},$$

$${\beta }_2={\displaystyle \frac{(q-2)[{(q-1)}^3(q-2)B_{1,q}+3q^2(q-3)(q-4)k_1^2{B}_{2,q}]}{(q-3)(q-4)q{[{(q-1)}^2{k}_2-q(q-2)k_1^2]}^2}},$$

where $B_{1,q}=(q-1)(q-3)k_4-4q(q-4)k_1{k}_3$, $B_{2,q}=2{(q-1)}^2{k}_2-q(q-2)k_1^2$ and $k_r$ is the rth moment of the $PM$ distribution.

Figure 4a,b show in graphic form the asymmetry and kurtosis coefficients of the PM model for different parameter values. Figure 4c,d show the asymmetry and kurtosis coefficients of the SPM model for different parameter values. It can be seen that both coefficients increase when the parameter q decreases. When observing the Figures for both models, it is evident that the SPM distribution covers a wide range of values for the asymmetry and kurtosis coefficients, depending on the parameter values, implying that the SPM model is sufficiently flexible to model real datasets.

Table 2. Asymmetry and kurtosis values for the SPM distribution.

$\mathit{\beta }$	q	$\sqrt{{\mathit{\beta }}_1}$	${\mathit{\beta }}_2$
$0.5$	$4.1$	$3.365$	188.415
2		$4.716$	479.077
10		$6.506$	767.271
1000		$6.636$	789.663
$0.5$	$4.5$	$2.780$	36.686
2		$3.545$	81.509
10		$5.323$	144.018
1000		$5.466$	149.444
$0.5$	5	2.397	19.256
2		2.670	35.360
10		4.482	70.302
1000		4.648	73.800
$0.5$	7	$1.893$	9.480
2		$1.155$	8.844
10		$3.085$	25.013
1000		$3.381$	27.857
$0.5$	9	$1.763$	8.079
2		$0.595$	5.081
10		$2.475$	16.444
1000		$2.940$	19.755
$0.5$	10	$1.732$	7.791
2		$0.440$	4.360
10		$2.252$	14.224
1000		$2.811$	17.828

reveals that values for the asymmetry and kurtosis coefficients depend on the parameters $\beta$ and q, and that as q decreases and $\beta$ increases, the asymmetry and kurtosis coefficients also increase. On the other hand, as q increases the asymmetry and kurtosis coefficients are those of the PM distribution.

On the other hand, we also have that

$$\underset{q\to \infty }{lim}\sqrt{{\beta }_1}={\displaystyle \frac{2\pi {\Gamma }_3-12{\pi }^{\frac{1}{2}}{\Gamma }_1{\Gamma }_2+16{\left({\Gamma }_1\right)}^3}{{(2{\pi }^{\frac{1}{2}}{\Gamma }_2-4{\left({\Gamma }_1\right)}^2)}^{\frac{3}{2}}}}\mathrm{and}\underset{q\to \infty }{lim}{\beta }_2={\displaystyle \frac{2{\pi }^{\frac{3}{2}}{\Gamma }_4-16\pi {\Gamma }_1{\Gamma }_3+48{\pi }^{\frac{1}{2}}{\Gamma }_2{\left({\Gamma }_1\right)}^2-48{\left({\Gamma }_1\right)}^4}{{(2{\pi }^{\frac{1}{2}}{\Gamma }_2-4{\left({\Gamma }_1\right)}^2)}^2}},$$

are the asymmetry and kurtosis coefficients of the PM distribution, where ${\Gamma }_r=\Gamma \left(\frac{3\beta +r}{2\beta }\right)$.

3. Properties

This section is devoted to studying some properties of the model.

Proposition 5.

Let $Z\sim SPM(\alpha ,\beta ,q)$ and $X\sim PM(\alpha ,\beta )$; then Z converges in distribution to X, as $q\to \infty$.

Proof. 

The proof will use Lehman’s and Slutsky’s theorem (Lehmann [19], p. 70). It is easy to see that $X\stackrel{\mathbb{D}}{\to }X,q\to \infty$, and $U^{1/q}\stackrel{\mathbb{P}}{\to }1,q\to \infty$. Thus, if we define the function $g(X,U^{1/q})=\frac{X}{U^{1/q}}$, then:

$g(X,U^{1/q})\stackrel{\mathbb{D}}{\to }g(X,1)=X$ as $q\to \infty$. Thus, $Z\stackrel{\mathbb{D}}{\to }X$ as $q\to \infty$.□

Proposition 6.

Let $z_{\left(1\right)},\cdots ,z_{\left(n\right)}$ denote the order statistics of a random sample $z_1,\cdots ,z_n$ from $Z\sim SPM(\alpha ,\beta ,q)$. Then the pdf of $Z_{\left(j\right)}$, $f_{Z_{\left(j\right)}}\left(z\right)$, is:

$${\displaystyle \frac{n!}{(j-1)!(n-j)!}}{\displaystyle \frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{\frac{q}{2\beta }}{z}^{q+1}}}{G}_1{\left[G_0-{\displaystyle \frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{\frac{q}{2\beta }}}}{G}_1\right]}^{j-1}{\left[1+{\displaystyle \frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{\frac{q}{2\beta }}}}{G}_1-G_0\right]}^{n-j},$$

where $G_1=G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)$ and $G_0=G\left(\alpha z^{2\beta },3/2,1\right)$.

Proof. 

Let $X_{\left(1\right)},\cdots ,X_{\left(n\right)}$ be the order statistics of a random sample, from a continuous population with CDF $F_X\left(x\right)$ and pdf $f_X\left(x\right)$. Then, the $pdf$ of $X_{\left(j\right)}$ is

$$f_{Z_{\left(j\right)}}\left(z\right)={\displaystyle \frac{n!}{(j-1)!(n-j)!}}{f}_X\left(x\right){\left[F_X\left(x\right)\right]}^{j-1}{[1-F_X\left(x\right)]}^{n-j}.$$

Now, replace $f_X\left(x\right)$ and $F_X\left(x\right)$ with the pdf and CDF of the SPM distribution respectively; then the result is obtained.□

Corollary 4.

Let $z_{\left(1\right)},\cdots ,z_{\left(n\right)}$ denote the order statistics of a random sample $z_1,\cdots ,z_n$ from $Z\sim SPM(\alpha ,\beta ,q)$. Then the pdfs of the minimum $Z_{\left(1\right)}$ and the maximum $Z_{\left(n\right)}$ are respectively:

$$f_{Z_{\left(1\right)}}\left(z\right)=n{\displaystyle \frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/\left(2\beta \right)}{z}^{q+1}}}{G}_1{\left[1+{\displaystyle \frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/\left(2\beta \right)}}}{G}_1-G_0\right]}^{n-1},$$

$$f_{Z_{\left(n\right)}}\left(z\right)=n{\displaystyle \frac{2q\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/\left(2\beta \right)}{z}^{q+1}}}{G}_1{\left[G_0-{\displaystyle \frac{2z^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/\left(2\beta \right)}}}{G}_1\right]}^{n-1},$$

where $G_1=G\left(\alpha z^{2\beta };\frac{q+3\beta }{2\beta },1\right)$ and $G_0=G\left(\alpha z^{2\beta },3/2,1\right)$.

Proposition 7.

Let $Z\sim SPM(\alpha ,\beta ,q)$, and $X\sim PM(\alpha ,\beta )$. with CDF $F_Z$ and $F_X$ respectively. Then $F_Z$ is stochastically greater than $F_X$.

Proof. 

Recall that $F_Z$ is stochastically greater than $F_X\left(x\right)$, if $P(Z>t)\ge P(X>t)$ for every t, which implies $F_Z\left(t\right)\le F_X\left(t\right)$. Now, notice that the second term of the distribution function of Z is always greater than or equal to zero. Then we have:

$$G(\alpha t^{2\beta },3/2,1)\ge G(\alpha t^{2\beta },3/2,1)-{\displaystyle \frac{2t^{-q}\Gamma \left(\frac{q+3\beta }{2\beta }\right)}{\sqrt{\pi }{\alpha }^{q/2\beta }}}G\left(\alpha t^{2\beta },\frac{q+3\beta }{2\beta },1\right),$$

which means that $F_Z\left(t\right)\le F_X\left(t\right)$.□

Table 3. Values of mean, variance, median and mode.

Parameters	Mean	Variance	Median	Mode
$(0.3,1.5,3)$	2.378	2.519	2.071	1.847
$(1.5,1.5,3)$	1.391	0.862	1.211	1.080
$(3,1.5,3)$	1.104	0.543	0.961	0.857
$(1.5,0.5,3)$	1.500	2.750	1.075	0.415
$(1.5,1.5,3)$	1.391	0.862	1.211	1.080
$(1.5,3,3)$	1.428	0.742	1.225	1.099
$(1.5,1.5,2.1)$	1.770	16.439	1.330	1.109
$(0.3,1.5,3)$	1.391	0.862	1.211	1.080
$(0.3,1.5,4)$	1.236	0.336	1.142	1.057
$(0.3,1.5,4.5)$	1.192	0.256	1.119	1.048

shows the quantities of the mean, variance, median and mode when $\alpha$, $\beta$ and q are increasing.

We note that the quantities decrease as long as q increases.

4. Inference

This section is devoted to inference aspects for the SPM model. Parameters are estimated based on the moments and ML methods.

4.1. Moments Estimation

Replacing $E\left(Z\right)$ by the sample mean $\overline{Z}$, gives the equation:

$$\overline{Z}=\frac{q}{q-1}{k}_1,$$

(4)

by which we obtain the moment estimator of q:

$$\tilde{q}={\displaystyle \frac{\overline{Z}}{\overline{Z}-k_1}}.$$

(5)

Therefore, using (5) and replacing $E\left(Z^2\right)$ and $E\left(Z^3\right)$ with the second and third sampling moments, we obtain the following equations:

$$\overline{Z^2}={\displaystyle \frac{\overline{Z}}{2k_1-\overline{Z}}}{k}_2,$$

(6)

$$\overline{Z^3}={\displaystyle \frac{\overline{Z}}{3k_1-2\overline{Z}}}{k}_3,$$

(7)

where $k_r$ is the r-th moment of the PM distribution.

The equation systems generated by (6) and (7) need to be solved numerically leading to the estimators $\tilde{\beta }$ and $\tilde{\alpha }$. The estimator $\tilde{q}$ is obtained from (5), replacing $\alpha$ by $\tilde{\alpha }$ and $\beta$ by $\tilde{\beta }$, which leads to the moment estimators $(\tilde{\alpha },\tilde{\beta },\tilde{q})$ for $(\alpha ,\beta ,q)$.

4.2. ML Estimation

In this section, we present the ML equations for parameters $(\alpha ,\beta ,q)$ of the SPM model. If $Z_1,Z_2,\cdots ,Z_n$ is a random sample from the random variable $Z\sim SPM(\alpha ,\beta ,q)$, the log-likelihood function can be expressed as

$$l(\theta ,\mathbf{z})=nlog\left(2\right)+nlog\left(q\right)-{\displaystyle \frac{n}{2}}log\left(\pi \right)-{\displaystyle \frac{nq}{2\beta }}log\left(\alpha \right)-(q+1)\sum _{i=1}^{n}log\left(z_i\right)+nlog\Gamma \left({\displaystyle \frac{3\beta +q}{2\beta }}\right)+\sum _{i=1}^{n}logH\left(z_i\right),$$

(8)

where $H\left(z_i\right)=G\left(\alpha z_i^{2\beta };\frac{q+3\beta }{2\beta },1\right)$ and $\theta =(\alpha ,\beta ,q)$.

The ML estimators are obtained by maximizing the log-likelihood function given in (8). By deriving the log-likelihood function with respect to each parameter, the following estimating equations are obtained:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial }{\partial \alpha }}l(\theta ,\mathbf{z})& =& {\displaystyle \frac{nq}{2\alpha \beta }}+\sum _{i=1}^n{\displaystyle \frac{H_1\left(z_i\right)}{H\left(z_i\right)}}=0,\hfill \\ \hfill {\displaystyle \frac{\partial }{\partial \beta }}l(\theta ,\mathbf{z})& =& {\displaystyle \frac{nq}{2{\beta }^2}}log\left(\alpha \right)-{\displaystyle \frac{nq}{2{\beta }^2}}{\psi }^{\left(1\right)}\left({\displaystyle \frac{3\beta +q}{2\beta }}\right)+\sum _{i=1}^n{\displaystyle \frac{H_2\left(z_i\right)}{H\left(z_i\right)}}=0,\hfill \\ \hfill {\displaystyle \frac{\partial }{\partial \beta }}l(\theta ,\mathbf{z})& =& {\displaystyle \frac{n}{q}}-{\displaystyle \frac{n}{2\beta }}log\left(\alpha \right)-\sum _{i=1}^{n}log\left(z_i\right)+{\displaystyle \frac{n}{2\beta }}{\psi }^{\left(1\right)}\left({\displaystyle \frac{3\beta +q}{2\beta }}\right)+\sum _{i=1}^n{\displaystyle \frac{H_3\left(z_i\right)}{H\left(z_i\right)}}=0,\hfill \end{array}$$

where $H_1\left(z_i\right)={\displaystyle \frac{\partial }{\partial \alpha }}H\left(z_i\right)$,   $H_2\left(z_i\right)={\displaystyle \frac{\partial }{\partial \beta }}H\left(z_i\right)$,   $H_3\left(z_i\right)={\displaystyle \frac{\partial }{\partial q}}H\left(z_i\right)$ and ${\psi }^{\left(1\right)}$ are the digamma function.  

These equation systems need to be solved numerically using an adequate root finder procedure. The moments estimator could be used to find initial values for the iterative procedure.

The EM algorithm (Dempster et al. [20]) defines an iterative process that allows the likelihood function of a parametric model to be maximized in cases in which some variables of the model are (or in our case, are treated as) "latent" or unknown. The advantage is simple: the EM algorithm is significantly less computationally intensive and far more robust. However, we will use the expectation/conditional maximization (ECM) algorithm (Meng and Rubin [21]), which is an extension of the EM algorithm.

4.3. EM Algorithm

Using the stochastic representation (3), but considering the following parametrization $q^{*}=\frac{q}{\beta }$, we obtain:

$$Z_i=\frac{X_i}{U_i^{1/\beta q^{*}}},$$

(9)

so we have the following hierarchical model

$$Z_i|V_i\sim PM(\alpha V_i^2,\beta )$$

$$V_i\sim \beta (q^{*},1),$$

where $V_i=U_i^{\frac{1}{q^{*}}}$ is the latent variable, and $Z_i$ is the observable variable, for $i=1,\cdots ,n$. Now, consider $\Psi =(\alpha ,\beta ,q^{*})$ the new parameter vector. In this case the complete-data log-likelihood or full log-likelihood function is

$$l_c(\Psi )={\displaystyle \frac{3}{2}}nlog\left(\alpha \right)+nlog\left(\beta \right)+nlog\left(q^{*}\right)+(q^{*}+2)\sum _{i=1}^{n}log\left(v_i\right)+(3\beta -1)\sum _{i=1}^{n}log\left(z_i\right)-\alpha \sum _{i=1}^n{v}_i^2{z}_i^{2\beta },$$

where $l_c(\Psi )$ is the full log-likelihood function. Then $E\left(l_c(\Psi )\right|z_i,\Psi ={\Psi }^{\left(k\right)})$:

$$\begin{array}{ccc}\hfill Q(\Psi ,{\Psi }^{\left(k\right)})& =& {\displaystyle \frac{3n}{2}}log\left(\alpha \right)+nlog\left(\beta \right)+nlog\left(q^{*}\right)+(q^{*}+2)\sum _{i=1}^n{\widehat{log\left(v_i\right)}}^{\left(k\right)}+(3\beta -1)\sum _{i=1}^{n}log\left(z_i\right)\hfill \\ & & -\alpha \sum _{i=1}^n{z}_i^{2\beta }{\widehat{v_i^2}}^{\left(k\right)},\hfill \end{array}$$

where ${\widehat{v_i^2}}^{\left(k\right)}=E\left(V_i^2\right|z_i,\Psi ={\Psi }^{\left(k\right)})$ and ${\widehat{log\left(v_i\right)}}^{\left(k\right)}=E(log\left(V_i\right)|z_i,\Psi ={\Psi }^{\left(k\right)})$. Normally, the next thing to do is find the distribution of $V_i|z_i$

$$f\left(v_i\right|z_i)\propto f(z_i,v_i)\propto v^{q^{*}+2}{e}^{-\alpha v^2{z}^{2\beta }},$$

but instead, we will use the distribution of $W_i=V_i^2|z_i$:

$$f\left(w_i\right|\Psi )\propto w_i^{(q^{*}+3)/2-1}{e}^{-\alpha z^{2\beta }{w}_i}{I}_{(0,1)}\left(w_i\right).$$

Thus, we conclude that W follows a truncated gamma in the interval (0, 1) (shown in the Appendix A), that is, $W_i\sim T{G}_{(0,1)}\left(\frac{q^{*}+3}{2},\alpha z_i^{2\beta }\right)$. Then:

$$\widehat{w_i^{\left(k\right)}}={\displaystyle \frac{G(1;{\displaystyle \frac{q^{*\left(k\right)}+5}{2}},{\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}})}{G(1;{\displaystyle \frac{q^{*\left(k\right)}+3}{2}},{\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}})}}{\displaystyle \frac{q^{*\left(k\right)}+3}{2{\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}}}},$$

(10)

$$\frac{1}{2}{\widehat{log\left[w_i\right]}}^{\left(k\right)}={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\left({\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}}\right)}^{(q^{*\left(k\right)}+3)/2}}{\Gamma \left({\displaystyle \frac{q^{*\left(k\right)}+3}{2}}\right)G(1;{\displaystyle \frac{q^{*\left(k\right)}+3}{2}},{\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}})}}{\int }_0^{1}log\left(w_i\right)w_i^{(q^{*\left(k\right)}+1)/2}{e}^{-{\alpha }^{\left(k\right)}{z}_i^{2{\beta }^{\left(k\right)}}{w}_i}d{w}_i,$$

which correspond to ${\widehat{v_i^2}}^{\left(k\right)}$ and ${\widehat{log\left(v_i\right)}}^{\left(k\right)}$ respectively. The equations result from deriving $Q(\Psi ,{\Psi }^{\left(k\right)})$ with respect to the parameters:

$$\begin{array}{ccc}\hfill {\alpha }^{(k+1)}& =& {\displaystyle \frac{3n}{2{\sum }_{i=1}^n{z}_i^{2{\beta }^{\left(k\right)}}{\widehat{v_i^2}}^{\left(k\right)}}},\hfill \\ \hfill {\beta }^{(k+1)}& =& -n{\left(3\sum _{i=1}^{n}log\left(z_i\right)-{\alpha }^{\left(k\right)}\sum _{i=1}^{n}2z_i^{2{\beta }^{(k+1)}}log\left(z_i\right){\widehat{v_i^2}}^{\left(k\right)}\right)}^{-1},\hfill \\ \hfill q^{*(k+1)}& =& -{\displaystyle \frac{n}{{\sum }_{i=1}^n{\widehat{log\left(v_i\right)}}^{\left(k\right)}}}.\hfill \end{array}$$

(11)

The equation corresponding to ${\beta }^{(k+1)}$ must be solved numerically. Finally, we obtain the following EM algorithm.

Step E: calculate ${\widehat{v_i^2}}^{\left(k\right)}$ and ${\widehat{log\left(v_i\right)}}^{\left(k\right)}$ using (10).

Step M: update ${\widehat{\Psi }}^{(k+1)}$ using (11).  

Once some stop criteria have been completed, calculate $q=q^{*}\beta$ to obtain the estimators of $(\alpha ,\beta ,q)$. For example, one stop criteria could be ${\widehat{\Psi }}^{(k+1)}-{\widehat{\Psi }}^{\left(k\right)}<0.0001$.

5. Simulation

In this section we present the simulation study, which aims to investigate the ML estimation performance for parameters $\alpha$, $\beta$, q under the SPM model; then we focus on comparison, using the AIC (Akaike [22]), BIC (Schwarz [23]), ADR and AD2R (Anderson et al. [24]) criteria. First, we present the algorithm that was used to generate samples for $Z\sim SPM(\mu ,\sigma ,q)$.

• Generate $P_i\sim {\chi }_3^2$ (chi squared with 3 degrees of freedom), $i=1,2,\cdots ,n$.
• Compute $X_i={\left(\frac{\sqrt{P_i}}{\sqrt{2\alpha }}\right)}^{1/\beta }$, $i=1,2,\cdots ,n$.
• Generate $U_i\sim U(0,1)$,   $i=1,2,\cdots ,n$.
• Compute $Z_i=\frac{X_i}{U_i^{1/q}}\sim SPM(\alpha ,\beta ,q)$, $i=1,2,\cdots ,n$,

where ${\chi }_3^2$ and $U(0,1)$ are the chi-squared distribution with three degrees of freedom and the standard uniform distribution respectively.

The algorithm that was used to calculate the ADR and AD2R values is described below:

$$ADR={\displaystyle \frac{n}{2}}-2{\sum }_{i=1}^{n}F\left(x_{\left(i\right)}\right)-{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n(2i-1)log(1-F\left(x_{(n+1-i)}\right)),$$

$$AD2R=-2{\sum }_{i=1}^{n}log(1-F\left(x_{\left(i\right)}\right))+{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\displaystyle \frac{(2i-1)}{(1-F\left(x_{(n+1-i)}\right))}},$$

where $F\left(x\right)$ is the CDF, and $x_{\left(i\right)}$ the i-th statistic of order. To calculate AIC and BIC values:

$$AIC=2k-2l(\widehat{\theta };\mathbf{z}),$$

$$BIC=klog\left(n\right)-2l(\widehat{\theta };\mathbf{z}),$$

where k is the number of parameters, n the sample size and $l(\widehat{\theta };\mathbf{z})$ the value of the log-likelihood function using the ML estimator of $\theta$.

5.1. Parameter Recovery

Using the algorithm, we generate 1000 random samples of sizes $n=50$, $n=100$ and $n=200$ under the SPM model, with different parameter values. In

Table 4. Mean of each maximum likelihood (ML) estimate with standard deviations in parentheses.

n	$(\mathit{\alpha },\mathit{\beta },\mathit{q})$	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{q}}$
50		$1.007\left(0.336\right)$	$0.527\left(0.121\right)$	$1.349\left(1.257\right)$
100	$(1,0.5,1)$	$1.005\left(0.207\right)$	$0.515\left(0.080\right)$	$1.066\left(0.262\right)$
200		$1.001\left(0.135\right)$	$0.507\left(0.054\right)$	$1.029\left(0.147\right)$
50		$0.979\left(0.285\right)$	$0.529\left(0.108\right)$	$4.099\left(3.960\right)$
100	$(1,0.5,2)$	$0.981\left(0.188\right)$	$0.509\left(0.066\right)$	$2.886\left(2.283\right)$
200		$0.987\left(0.130\right)$	$0.503\left(0.043\right)$	$2.318\left(1.053\right)$
50		$0.975\left(0.284\right)$	$0.526\left(0.103\right)$	$6.684\left(5.245\right)$
100	$(1,0.5,3)$	$0.976\left(0.179\right)$	$0.506\left(0.056\right)$	$5.231\left(3.717\right)$
200		$0.982\left(0.128\right)$	$0.503\left(0.040\right)$	$4.237\left(2.540\right)$
50		$0.282\left(0.091\right)$	$1.599\left(0.368\right)$	$4.064\left(3.745\right)$
100	$(0.3,1.5,3.1)$	$0.293\left(0.066\right)$	$1.539\left(0.231\right)$	$3.349\left(0.861\right)$
200		$0.295\left(0.046\right)$	$1.522\left(0.154\right)$	$3.194\left(0.483\right)$
50		$1.56\left(0.559\right)$	$1.581\left(0.364\right)$	$4.119\left(3.951\right)$
100	$(1.5,1.5,3.1)$	$1.536\left(0.367\right)$	$1.549\left(0.239\right)$	$3.314\left(0.819\right)$
200		$1.510\left(0.231\right)$	$1.523\left(0.160\right)$	$3.190\left(0.477\right)$
50		$1.571\left(0.619\right)$	$1.598\left(0.370\right)$	$6.712\left(7.256\right)$
100	$(1.5,1.5,4.1)$	$1.521\left(0.348\right)$	$1.541\left(0.219\right)$	$4.668\left(2.119\right)$
200		$1.504\left(0.209\right)$	$1.519\left(0.143\right)$	$4.277\left(0.820\right)$

we present a summary of the study results. For each sample, an ML estimator was computed numerically using the EM algorithm, where the means and the standard deviations are reported. As expected, as the sample size increases, its standard deviation decreases and therefore the parameter estimate approaches the true simulated value. We conclude that the ML estimates are quite stable.

5.2. Criteria Comparison

We generated 1000 random samples of sizes $n=50$, $n=100$ and $n=200$ under the SPM model using different parameter values. The following tables show the proportions of times that the respective criteria (detailed in the previous section) chose the SPM over the PM model.

Table 5. SPM simulation, criteria comparison.

		AIC		BIC		ADR		AD2R
$\mathit{n}$	$(\mathit{\alpha },\mathit{\beta },\mathit{q})$	SPM	PM	SPM	PM	SPM	PM	SPM	PM
50		0.627	0.373	0.479	0.521	0.804	0.196	0.999	0.001
100	$(0.8,0.7,2.1)$	0.863	0.137	0.772	0.228	0.801	0.199	1.000	0.000
200		0.993	0.007	0.970	0.030	0.798	0.202	1.000	0.000
50		0.643	0.357	0.501	0.499	0.777	0.223	1.000	0.000
100	$(1,1,3.1)$	0.874	0.126	0.758	0.242	0.768	0.232	1.000	0.000
200		0.982	0.018	0.952	0.048	0.794	0.206	1.000	0.000
50		0.432	0.568	0.303	0.697	0.804	0.196	0.998	0.002
100	$(1,1,4.1)$	0.681	0.319	0.539	0.461	0.783	0.217	1.000	0.000
200		0.899	0.101	0.762	0.238	0.795	0.205	0.999	0.001
50		0.993	0.007	0.992	0.008	0.539	0.461	0.999	0.001
100	$(3,4,3.1)$	1.000	0.000	1.000	0.000	0.502	0.498	1.000	0.000
200		1.000	0.000	1.000	0.000	0.453	0.547	0.976	0.024
50		0.880	0.120	0.862	0.138	0.630	0.370	0.997	0.003
100	$(3,4,4.1)$	0.980	0.020	0.980	0.020	0.629	0.371	1.000	0.000
200		0.998	0.002	0.996	0.004	0.587	0.413	0.980	0.020

corresponds to data simulated by the SPM model. Results show that when the data are drawn from the SPM distribution, the four criteria corresponding to the SPM model are usually lower than the criteria for the PM model. We also conclude that for models with heavy tails, ADR and AD2R are better criteria for evaluating model fitness.

6. Applications

We now present two applications to real datasets showing descriptive analysis and ML estimators for the M, PM and SPM models with their AIC, BIC, ADR and AD2R criteria, histograms and Q-Q plots.

6.1. Application 1

We consider a dataset of fund-raising expenses as a percentage of total expenditure for a random sample of 60 charities from the United States. These observations can be found in Devore [25], p. 4. These are:

6.1,12.6,34.7,1.6,18.8,2.2,3.0,2.2,5.6,3.8,2.2,3.1,1.3,1.1,14.1,4.0,21.0,6.1,1.3,20.4

7.5,3.9,10.1,8.1,19.5,5.2,12.0,15.8,10.4,5.2,6.4,10.8,83.1,3.6,6.2,6.3,16.3,12.7,1.3,

0.8,8.8,5.1,3.7,26.3,6.0,48.0,8.2,11.7,7.2,3.9,15.3,16.6,8.8,12.0,4.7,14.7,6.4,17.0,

2.5,16.2

In

Table 6. Descriptive analysis for the charities data.

Mean	S.D.	Median	Interquartile Range	Min.	Max.	Asymmetry	Kurtosis
10.892	12.741	6.800	10.375	0.800	83.100	3.602	19.360

, we can see from the descriptive analysis that this dataset has high kurtosis.

Table 7. ML estimates for the charities data with standard deviations in parentheses.

Parameter	M	PM	SPM
$\alpha$	0.005(0.0006)	0.196(0.047)	0.198(0.059)
$\beta$	—	0.436(0.039)	0.563(0.080)
q	—	—	2.122(0.733)
l	277.564	201.531	199.017
AIC	557.127	407.063	404.034
BIC	559.228	411.251	410.317
ADR	65.062	26.711	25.425
AD2R	2.752 × 10${}^{15}$	2799.457	269.557

shows the ML estimates for the parameters of the M, PM and SPM distributions. For each model, we apply the statistical criteria named at the beginning of this section; we see that all four criteria choose the SPM model over the M and PM models.

Figure 5 depicts the histogram of the dataset, where we can check that the SPM model has a greater reach than the other models.

Finally, Figure 6a,b show the Q-Q plots for the PM and SPM models respectively; Figure 6c shows the empirical plot, which tells us that the SPM model is preferable.

6.2. Application 2

We consider a dataset of the copper content of 24 Bidri products (a traditional Indian handicraft); this study was performed because Bidri crafts are soldered with an alloy containing mainly zinc and some copper (Devore [25], page 33). These are:

2.0, 2.4, 2.5, 2.6, 2.6, 2.7, 2.7, 2.8, 3.0, 3.1, 3.2, 3.3, 3.3, 3.4, 3.6, 3.6, 3.6,

3.7, 4.4, 4.6, 4.7, 4.8, 5.3, 10.1

In

Table 8. Descriptive analysis for the copper data.

Mean	S.D.	Median	Interquartile Range	Min.	Max.	Asymmetry	Kurtosis
3.666	1.612	3.300	1.75	2	10.1	2.760	11.689

, we can see from the descriptive analysis that this dataset has high kurtosis, so it is interesting to see what our model can do here.

Table 9. ML estimates for the copper data with standard deviations in parentheses.

Parameter	M	PM	SPM
$\alpha$	0.094(0.015)	0.097(0.043)	0.006(0.005)
$\beta$	—	0.989(0.130)	2.704(0.517)
q	—	—	3.591(1.083)
l	−42.193	−42.190	−34.567
AIC	86.386	88.380	75.135
BIC	87.564	90.736	78.670
ADR	0.762	0.753	0.106
AD2R	160.036	138.239	2.061

shows the ML estimates for the parameters of the M, PM, and SPM distributions. For each model we apply the statistics criteria named at the beginning of this section, and we see that all four criteria choose the SPM model over the M and PM models.

Figure 7 depicts the histogram of the dataset, where we can check that the SPM model has a greater reach than the other models.

Finally, Figure 8a,b show the Q-Q plots for the PM and SPM models respectively, and Figure 8c shows the empirical plot, which tells us that the SPM model is preferable.

It is interesting to note that the estimate of the shape parameter $\beta$ of the PM model is almost 1; remember that if $\beta =1$, then we obtain the classic M distribution.

7. Conclusions

The slash methodology has often been used recently to extend a variety of well-known distributions, resulting in a flexible distribution that has high kurtosis, and of course heavier tails, as we discuss above. Here we introduced the SPM distribution, where following the same methodology, we extended the PM distribution by using the PM instead of the normal model; see Equation (1). So the PM is a special case of the SPM model. As we saw in Figure 1, the new distribution does indeed present a more flexible kurtosis coefficient than that of the PM model; moreover, as we saw in Table 1, the tails become heavier as parameter q becomes smaller.

We discussed some properties, such as the distribution, survival and hazard functions. We also studied the mode; moments; statistic order; and of course, asymmetry and kurtosis. They all have closed form, which makes computational implementation easier. We also presented two stochastic representations: One is defined as a quotient between two independent random variables: PM in the numerator and the power of a standard uniform distribution in the denominator. The second is a mixture of a PM and a beta distribution, used to establish the EM algorithm.

Parameters could be estimated using the moments method and ML estimators using log-likelihood or EM algorithm methodology, although the EM algorithm is more stable, as we obtained closed expressions.

The results of a simulation study were presented, indicating good parameter recovery. In the simulation study for the criteria, we concluded that the criterion which best explained the fit of the SPM model was AD2R.

In the application we used datasets with high kurtosis, which seems logical considering the object of the work. Four model comparison criteria statistics were used; all four criteria indicate that the SPM model shows the best fit for this dataset. The ADR and AD2R criteria are more reliable for comparing distributions using two real datasets; this can be seen in the application, where these criteria showed a better fit for the SPM model.