Inference for three-parameter 𝑴 -Wright distributions with applications

Abstract

We propose point estimators for the three-parameter (location, scale, and the fractional parameter) variant distributions generated by a Wright function. We also provide uncertainty quantification procedures for the proposed point estimators under certain conditions. The class of densities includes the three-parameter one-sided and the three-parameter symmetric bimodal $M$ -Wright family of distributions. The one-sided family naturally generalizes the Airy and half-normal models. The symmetric class includes the symmetric Airy and normal or Gaussian densities. The proposed interval estimator for the scale parameter outperformed the estimator derived in Cahoy (2012) when the location parameter is zero. We obtain the asymptotic covariance structure for the scale and fractional parameter estimators, which allows estimation of the correlation. The coverage probabilities of the interval estimators slightly depend on the proposed location parameter estimators. For the symmetric case, the sample mean (or median) is favored than the median (or mean) when the fractional parameter is greater (or lesser) than 0.39106 in terms of their asymptotic relative efficiency. The estimation algorithms were tested using synthetic data and were compared with their bootstrap counterparts. The proposed inference procedures were demonstrated on age and height data.

Keywords

Gaussian skew-normal M-Wright Mittag-Leffler skew-Laplace Airy skew-symmetric Major Leaque Baseball children heights

1. Introduction

The $M$ -Wright function has been increasingly gaining popularity from several areas of study particularly in mathematics, engineering and physics. It is often a probability density function in space which solves time-fractional diffusion processes (see Mura et al., 2008). As a solution, the $M$ -Wright density naturally models the increments or the ‘space’ component of the above processes at any given time. It is also used as a subordinator (as the operational time rather than the physical time) for time-fractional differential equations (Pagnini and Scalas, 2014), for a multi-point probability model of the generalized grey Brownian motion that includes the well-known standard and fractional Brownian motions, and for pure linear birth processes (see Beghin and Orsingher, 2010; Cahoy and Polito, 2012). The single-parameter positive-sided $M$ -Wright function takes the following form:

$M_{\alpha}(x)=\sum\limits_{j=0}^{\infty}\frac{(-x)^{j}}{j!\Gamma[-\alpha j+(1-% \alpha)]}=\frac{1}{\pi}\sum\limits_{j=1}^{\infty}\frac{(-x)^{j-1}}{(j-1)!}% \Gamma(\alpha j)\sin(\pi\alpha j)$ (1)

where $x\in\mathbb{R}^{+}$ , and $0<\alpha<1$ is the fractional parameter. The last equality in the preceding equation follows from the reflection formula for the gamma function

$\Gamma(1-\alpha(j+1))=\frac{\pi}{\Gamma(\alpha(j+1))\sin(\pi\alpha(j+1))}$

and transformation $j+1\to j.$ We have the exponential density $(\alpha=0^{+})$ as a limiting case and the Airy ( $\alpha=1/3$ ) and half-normal ( $\alpha=1/2$ ) (see Mainardi et al., 2010) distributions as special cases where

$M_{1/2}(x)=\frac{1}{\sqrt{\pi}}e^{-x^{2}/4}.$ (2)

Moreover,

$M_{\alpha}(0^{+})=1/\Gamma(1-\alpha),\quad\text{and}\quad M_{1^{-}}(x)=\delta(% x-1),$

where $\delta(\cdot)$ is the generalized Dirac function. The Laplace transform of Eq. (1) is

$\mathbf{E}\left(e^{-\beta X}\right)=\phi_{X}(\beta)=\sum\limits_{j=1}^{\infty}% \frac{\left(-\beta\right)^{j}}{j!\;\Gamma(1+\alpha j)}$ (3)

which is the Mittag-Leffler function. The positive-sided $M$ -Wright random variable has the structural representation

$X\stackrel{d}{=}S^{-\alpha},$ (4)

where $S$ follows an $\alpha^{+}$ -stable distribution (Zolotarev, 1986) with $\phi_{S}(\beta)=\exp(-\beta^{\alpha})$ . The $\kappa$ th moment (see Piryatinska et al., 2005) is known to be

$\mathbf{E}X^{\kappa}=\frac{\Gamma(1+\kappa)}{\Gamma(1+\alpha\kappa)},\qquad% \kappa>-1,$ (5)

giving the mean and variance as

$\mu_{x}=\frac{1}{\alpha\Gamma(\alpha)},\quad\text{and}\quad\sigma_{x}^{2}=% \frac{1}{\alpha\Gamma(2\alpha)}-\frac{1}{\left(\alpha\Gamma(\alpha)\right)^{2}},$ (6)

correspondingly. The coefficient of variation is straightforward to calculate as

$\frac{\sigma_{x}}{\mu_{x}}=\sqrt{\frac{2\alpha\Gamma(\alpha)\Gamma(1+\alpha)}{% \Gamma(1+2\alpha)}-1}=\begin{cases}1&\alpha=0,\\ \sqrt{\pi/2-1}&\alpha=1/2,\\ 0&\alpha=1.\end{cases}$ (7)

The rest of the paper is organized as follows. The one-sided $M$ -Wright density, its properties, and test results are presented in Section 2. The extension to the symmetric case are in Section 3. The applications and concluding remarks are given in Sections 4 and 5, respectively.

2. One-sided M-Wright distribution

The three-parameter one-sided $M$ -Wright density function has the following form:

$M_{\alpha,\rho,\mu}(x)=\frac{1}{\rho}M_{\alpha}\left(\frac{x-\mu}{\rho}\right)% ,\qquad x>\mu,$ (8)

where $\mu\in\mathbb{R}$ and $\rho\in\mathbb{R}^{+}$ are the shift and scale parameters, respectively. Below are some forms of the densities in this family.

Figure 1.

The one-sided $M$ -Wright density for $\alpha=0.1,0.5,0.8,0.9;\rho=1,\;\mu=0$ .

Case 2.1: $\mu=0$

If $X\stackrel{d}{=}M_{\alpha,\rho,0}(x)$ then

$\hskip-21.681pt\mathbf{E}X^{\kappa}=\frac{\rho^{\kappa}\Gamma(1+\kappa)}{% \Gamma(1+\alpha\kappa)},\;\;\phi_{X}(\beta)=E_{\alpha}(-\beta\rho),\;\;\phi_{% \rho}(\beta)=\beta^{1-1/\alpha}e^{-x\beta},\;\;\text{and}\;\;X\stackrel{d}{=}% \rho S^{-\alpha}.$ (9)

Given $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}M_{\alpha,\rho,0}(x),$ and applying the log transformation to the absolute value of the random variable $X$ given in Eq. (9), we obtain

$X^{\prime}\stackrel{d}{=}\log(\rho)-\alpha S^{\prime},$ (10)

where $X^{\prime}=\log(|X|)$ , and $S^{\prime}=\log(S)$ . From Cahoy (2012), the mean and variance are

$\mu_{X^{\prime}}=\log(\rho)+\gamma(\alpha-1),\quad\text{and}\quad\sigma_{X^{% \prime}}^{2}=\frac{\pi^{2}}{6}\left(1-\alpha^{2}\right),$ (11)

respectively, where $\gamma\approx 0.5772156649$ is the Euler’s constant. Moreover, the following point estimators of $\alpha$ and $\rho$ are obtained:

$\widehat{\alpha}=\sqrt{1-\frac{6\hat{\sigma}_{X^{\prime}}^{2}}{\pi^{2}}},\quad% \text{and}\quad\widehat{\rho}=\exp\left(\widehat{\mu}_{X^{\prime}}+\gamma(1-% \widehat{\alpha})\right).$ (12)

Proposition 1. Let $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}M_{\alpha,\rho,0}(x)$ . Then

$\sqrt{n}\left(\begin{array}[]{c}\widehat{\alpha}-\alpha\\ \widehat{\rho}-\rho\end{array}\right)\stackrel{d}{\longrightarrow}\textsl{N}% \left(\bm{0}\;,\;\bm{\Sigma^{\prime}}\right),\qquad n\to\infty,$ (13)

where

$\bm{\Sigma^{\prime}}=\left(\begin{array}[]{cc}\ \sigma_{\widehat{\alpha}% \widehat{\alpha}}&\ \sigma_{\widehat{\alpha}\widehat{\rho}}\\ \ \sigma_{\widehat{\alpha}\widehat{\rho}}&\ \sigma_{\widehat{\rho}\widehat{% \rho}}\\ \end{array}\right),\qquad\sigma_{\widehat{\alpha}\widehat{\alpha}}=\frac{11-% \alpha^{4}}{10\alpha^{2}}-1,$ (14)

$\sigma_{\widehat{\alpha}\widehat{\rho}}=\frac{\rho(10\*\alpha^{2}-11+\*\alpha^% {4})\gamma-\left[60\alpha(\alpha^{3}-1)\zeta(3)\right]/\pi^{2}}{10\*\alpha^{2}},$

$\sigma_{\widehat{\rho}\widehat{\rho}}=\frac{\rho^{2}\left(\;360\alpha(\alpha^{% 3}-1)\gamma\zeta(3)-(\alpha^{2}-1)\pi^{2}(3(11+\alpha^{2})\gamma^{2}+5\alpha^{% 2}\pi^{2})\;\right)}{30\alpha^{2}\pi^{2}},$ (15)

and $\zeta(\cdot)$ is the Riemman zeta function.

Proof Recall the following key results in Cahoy (2012): Let $\mu_{j}^{\prime}=\mathbf{E}\left(X^{\prime}-\mu_{X^{\prime}}\right)^{j},\;j=3,4.$ Then the third and fourth central moments are

$\mu_{3}^{\prime}=2(\alpha^{3}-1)\zeta(3)\quad\text{and}\quad\mu_{4}^{\prime}=% \frac{\pi^{4}(\alpha^{4}-10\alpha^{2}+9)}{60},$ (16)

respectively. In addition, if $\widehat{\mu}_{X^{\prime}}=\overline{X^{\prime}}=\sum_{j=1}^{n}X_{j}^{\prime}% \big{/}n$ and $\widehat{\sigma}_{X^{\prime}}^{2}=\sum_{j=1}^{n}\left(X_{j}^{\prime}-\overline% {X^{\prime}}\right)^{2}\big{/}n$ then it is widely known that

$\sqrt{n}\left(\begin{array}[]{c}\widehat{\mu}_{X^{\prime}}-\mu_{X^{\prime}}\\ \widehat{\sigma}_{X^{\prime}}^{2}-\sigma_{X^{\prime}}^{2}\end{array}\right)% \stackrel{d}{\longrightarrow}\textsl{N}\left[\bm{0},\bf{\Sigma}\right]$ (17)

as $n\to\infty$ , where the variance-covariance matrix $\bf{\Sigma}$ is defined as

$\bf{\Sigma}=\left(\begin{array}[]{cc}\sigma_{X^{\prime}}^{2}&\mu_{3}^{\prime}% \\ \mu_{3}^{\prime}&\mu_{4}^{\prime}-\sigma_{X^{\prime}}^{4}\\ \end{array}\right),$ (18)

$\mu_{3}^{\prime},\mu_{4}^{\prime}$ , and $\sigma_{X^{\prime}}^{2}$ are given in Eqs (11) and (16). Using result Eq. (17) and the multivariate delta method,

$\sqrt{n}\big{(}\textbf{g}(\widehat{\bm{\theta}}_{n})-\textbf{g}(\bm{\theta})% \big{)}\stackrel{d}{\to}\textsl{N}\left(\bm{0},\;\bm{\dot{\textbf{g}}}(\bm{% \theta})^{\text{T}}\bf{\Sigma}\bf{\dot{g}}(\bm{\theta})\right),$ (19)

where $\widehat{\bm{\theta}}_{n}=(\widehat{\mu}_{X^{\prime}},\widehat{\sigma}_{X^{% \prime}}^{2})^{\text{T}},\bf{g}$ is a continuous mapping from $\mathbb{R}^{2}\to\mathbb{R}^{2}$ given as

$\textbf{g}(\mu_{Y^{\prime}},\sigma_{Y^{\prime}}^{2})=\left(\sqrt{1-\frac{6% \sigma_{Y^{\prime}}^{2}}{\pi^{2}}},\exp\left(\mu_{Y^{\prime}}+\gamma(1-\alpha)% \right)\right)^{\text{T}}$

and $\bm{\dot{\textbf{g}}}(\bm{\theta})=\nabla\textbf{g}(\bm{\theta})^{\text{T}}$ is the gradient matrix given by

$\bm{\dot{\textbf{g}}}(\mu_{Y^{\prime}},\sigma_{Y^{\prime}}^{2})=\left(\begin{% array}[]{@{}cc@{}}0&\exp\left(\mu_{Y^{\prime}}+\gamma(1-\alpha)\right)\\ &\\ -3\Bigg{/}\left(\pi^{2}\sqrt{1-{\displaystyle\frac{6\sigma_{Y^{\prime}}^{2}}{% \pi^{2}}}}\right)&\left(3\gamma\exp(\mu_{Y^{\prime}}+\gamma)\sqrt{1-{% \displaystyle\frac{6\sigma_{Y^{\prime}}^{2}}{\pi^{2}}}}\right)\Bigg{/}\left(% \pi^{2}\sqrt{1-{\displaystyle\frac{6\sigma_{Y^{\prime}}^{2}}{\pi^{2}}}}\right)% \\ \end{array}\right).\qed$ (20)

Note that the covariance structure of the scale and fractional parameter estimators given by $\sigma_{\widehat{\alpha}\widehat{\rho}}$ above allows estimation of the correlation.

Corollary 1. Let $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}X$ . The $(1-\nu)$ 100% confidence intervals for $\alpha$ and $\rho$ can be approximated as

$\widehat{\alpha}\;\pm\;z_{\nu/2}\sqrt{\frac{\left[\left(11-\widehat{\alpha}^{4% }\right)/\left(10\widehat{\alpha}^{2}\right)\right]-1}{n}},$ (21)

and

$\widehat{\rho}\;\pm\;z_{\nu/2}\sqrt{\frac{\widehat{\rho}^{2}(360\widehat{% \alpha}(\widehat{\alpha}^{3}-1)\gamma\zeta(3)-(\widehat{\alpha}^{2}-1)\pi^{2}(% 3(11+\widehat{\alpha}^{2})\gamma^{2}+5\widehat{\alpha}^{2}\pi^{2}))}{30n% \widehat{\alpha}^{2}\pi^{2}}},$ (22)

correspondingly, where $z_{\nu/2}$ is the $(1-\nu/2)$ th quantile of the standard normal distribution, and $0<\nu<1$ .

Proof Immediately follows from Proposition 1 and is omitted. ∎

We tested our estimators by simulating the bias ( $100|\widehat{\theta}-\theta)|/\theta$ ), the median absolute deviation (MAD), and the 95% coverage probabilities for the proposed methods and the bootstrap percentile counterparts (with ‘*’) corresponding to several parameter combinations. Table 1 suggests that bias is as large as 5% and as little as 0.21% when $n=10^{4}$ . Reduction in variability is also apparent as the sample size goes large. It can be seen that the smaller the parameter $\alpha$ , the slower the reduction in variability and bias regardless of the sample size. Nevertheless, we conclude that these point estimators are consistent and asymptotically unbiased. Table 2 reveals that the proposed interval estimator of the scale parameter quickly captured (e.g., $n=$ 100 and $\alpha=0.6$ ) the true nominal level than the one in Cahoy (2012) as the sample size goes large. Furthermore, Table 2 illustrates that the large-sample interval estimator outperformed the percentile bootstrap method for estimating $\alpha$ especially when $n\leqslant$ 1000. Note that the large-sample formula is faster to calculate than the resampling-based method especially for large sample sizes.

Table 1

Mean estimates of and dispersions from the true parameters $\alpha$ , and $\rho$

$(\alpha,\rho)$		$n=100$		$n=1000$		$n=10000$
		% Bias	% MAD	% Bias	% MAD	% Bias	% MAD
(0.4, 150)	$\widehat{\alpha}$	32.796	31.795	16.143	19.346	4.706	5.922
	$\widehat{\rho}$	7.899	9.806	3.266	3.854	0.936	1.158
(0.6, 8.77)	$\widehat{\alpha}$	17.443	19.959	5.852	7.176	1.887	2.325
	$\widehat{\rho}$	6.437	7.648	1.955	2.435	0.641	0.809
(0.8, 375)	$\widehat{\alpha}$	7.842	9.505	2.513	3.051	0.763	0.935
	$\widehat{\rho}$	4.423	5.273	1.351	1.706	0.408	0.502
(0.95, 1000)	$\widehat{\alpha}$	2.871	2.686	0.924	1.147	0.294	0.359
	$\widehat{\rho}$	2.329	2.769	0.687	0.838	0.213	0.268

Table 2

Coverage probabilities of 95% interval estimates for different values of $\alpha$ , and $\rho$

$(\alpha,\rho)$		$n=$ 100	$n=$ 1000	$n=$ 10000
(0.4, 150)	$\widehat{\alpha}$	0.917	0.952	0.956
	$\widehat{\rho}$	0.941	0.957	0.949
	$\widehat{\alpha*}$	0.884	0.928	0.950
	$\widehat{\rho*}$	0.942	0.941	0.947
(0.6, 8.77)	$\widehat{\alpha}$	0.950	0.950	0.953
	$\widehat{\rho}$	0.949	0.958	0.951
	$\widehat{\alpha*}$	0.873	0.931	0.944
	$\widehat{\rho*}$	0.943	0.954	0.950
(0.8, 375)	$\widehat{\alpha}$	0.964	0.950	0.954
	$\widehat{\rho}$	0.942	0.958	0.952
	$\widehat{\alpha*}$	0.831	0.925	0.948
	$\widehat{\rho*}$	0.935	0.954	0.948
(0.95, 1000)	$\widehat{\alpha}$	0.960	0.922	0.952
	$\widehat{\rho}$	0.902	0.948	0.950
	$\widehat{\alpha*}$	0.724	0.888	0.931
	$\widehat{\rho*}$	0.922	0.948	0.947

Case 2.2: $\mu\neq 0$

Consider the location-scale structure

$X\stackrel{d}{=}\mu+\rho S^{-\alpha},\quad\text{and}\quad X>\mu.$ (23)

Proposition 2. Let $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}X$ in (23). A $(1-\nu)$ 100% confidence interval for the shift parameter $\mu$ is

$\left(\widehat{\mu}-q_{\nu}\widehat{\rho},\widehat{\mu}\right),$ (24)

where $q_{\nu}>0$ is the $(1-\nu^{1/n})$ th quantile of $M_{\alpha,1,0}$ and $\widehat{\mu}=\min\{X_{i}\}_{i=1}^{n}$ .

Proof Note that

$P(\mu<\widehat{\mu}<\mu+q_{\nu}\rho)=1-\nu,$ (25)

which suggests that

$\displaystyle P(\widehat{\mu}>\mu+q_{\nu}\rho)=\nu$ (26) $\displaystyle=\prod\limits_{i=1}^{n}P(X_{i}>\mu+q_{\nu}\rho)$ (27) $\displaystyle=\prod\limits_{i=1}^{n}P(\mu+\rho S_{i}^{-\alpha}>\mu+q_{\nu}\rho)$ (28) $\displaystyle\implies P(S^{-\alpha}>q_{\nu})=\nu^{1/n}.\qed$ (29)

For reproducibility, we estimate $q_{\nu}$ by generating $10^{6}$ random variates from $M_{\widehat{\alpha},1,0}$ and use the approximately median-unbiased (type 8 of the quantile function in R) estimator to calculate the $(1-\nu^{1/n})$ th quantile as recommended by Hyndman and Fan (1996). Note also that we directly use the point estimators obtained in Case 2.1 after subtracting $\widehat{\mu}$ from the observed data.

Upon testing, Table 3 generally indicates similar observations and conclusions about the estimators of $\alpha$ and $\rho$ as in Table 1. The mean and dispersion of $\widehat{\mu}$ seem to be large when $\nu\approx 1$ . Overall, the proposed point estimators are consistent. In addition, Table 4 shows that the proposed interval estimator for $\mu$ seems to capture the true nominal rate even when the sample size is as small as 100 with $\alpha\ll 1$ . Comparing Tables 3 and 4 with Tables 1 and 2, correspondingly, reveals that the variability induced by the subtraction of the minimum from the data does not seem to seriously affect the performance of the proposed estimators.

Table 3

Mean estimates of and dispersions from the parameters $\alpha$ , $\rho$ , and $\mu$

$(\alpha,\rho,\mu)$		$n=$ 100		$n=$ 1000		$n=$ 10000
		% Bias	% MAD	% Bias	% MAD	% Bias	% MAD
(0.4, 150, $-$ 78)	$\widehat{\alpha}$	32.459	31.174	16.026	19.658	4.919	6.116
	$\widehat{\rho}$	7.706	9.732	3.060	3.559	0.949	1.201
	$\widehat{\mu}$	2.032	2.055	0.220	0.221	0.022	0.023
(0.6, 8.77, 25.2)	$\widehat{\alpha}$	17.258	18.863	6.017	7.427	1.871	2.358
	$\widehat{\rho}$	6.737	8.896	2.012	2.487	0.613	0.761
	$\widehat{\mu}$	0.587	0.599	0.060	0.056	0.006	0.006
(0.8, 375, 375)	$\widehat{\alpha}$	7.597	8.491	2.716	3.407	0.821	1.051
	$\widehat{\rho}$	5.092	5.777	1.374	1.681	0.429	0.550
	$\widehat{\mu}$	3.521	3.404	0.333	0.331	0.036	0.036
(0.95, 1000, 500)	$\widehat{\alpha}$	2.971	2.926	0.940	1.147	0.309	0.397
	$\widehat{\rho}$	13.693	11.921	1.946	1.710	0.297	0.317
	$\widehat{\mu}$	6.936	6.890	2.764	2.690	0.296	0.278

Table 4

Coverage probabilities of 95% interval estimates for different values of $\alpha$ , $\rho$ , and $\mu$

$(\alpha,\rho,\mu)$		$n=$ 100	$n=$ 1000	$n=$ 10000
(0.4, 150, $-$ 78)	$\widehat{\alpha}$	0.921	0.945	0.947
	$\widehat{\rho}$	0.944	0.961	0.956
	$\widehat{\alpha*}$	0.951	0.953	0.947
	$\widehat{\rho*}$	0.952	0.967	0.954
	$\widehat{\mu}$	0.962	0.951	0.945
(0.6, 8.77, 25.2)	$\widehat{\alpha}$	0.955	0.955	0.952
	$\widehat{\rho}$	0.931	0.954	0.958
	$\widehat{\alpha*}$	0.976	0.968	0.963
	$\widehat{\rho*}$	0.940	0.956	0.958
	$\widehat{\mu}$	0.943	0.946	0.949
(0.8, 375, 375)	$\widehat{\alpha}$	0.966	0.940	0.947
	$\widehat{\rho}$	0.865	0.946	0.951
	$\widehat{\alpha*}$	0.980	0.971	0.953
	$\widehat{\rho*}$	0.887	0.945	0.952
	$\widehat{\mu}$	0.911	0.952	0.949
(0.95, 1000, 500)	$\widehat{\alpha}$	0.963	0.947	0.949
	$\widehat{\rho}$	0.945	0.957	0.951
	$\widehat{\alpha*}$	0.996	0.963	0.966
	$\widehat{\rho*}$	0.931	0.959	0.934
	$\widehat{\mu}$	0.907	0.941	0.947

3. Symmetric M-Wright distribution

Replacing $x$ by $|x|$ in Eq. (1) and dividing Eq. (1) by two, the three-parameter symmetric $M$ -Wright density can be written as

$M_{\alpha,\rho,\mu}(x)=\frac{1}{2\rho}M_{\alpha}\left(\frac{|x-\mu|}{\rho}% \right),\qquad x\in\mathbb{R},$ (30)

where $\mu\in\mathbb{R}$ and $\rho\in\mathbb{R}^{+}$ are the location and scale parameters, respectively. The Laplace or double exponential $(\alpha=0^{+})$ is a limiting case while the Gaussian or normal ( $\alpha=1/2$ ) (see Mainardi et al., 2010) distributions are special cases. Moreover,

$X\stackrel{d}{=}\mu+\rho US^{-\alpha},\qquad U\stackrel{ind}{=}(1/2)[\delta(u+% 1)+\delta(u-1)],$ (31)

where ‘ $i n d$ ’ means independent.

Case 3.1: $\mu=0$

The $M$ -Wright function in two variables that is centered at zero satisfies the following transformation:

$\phi_{|X|}(-\beta)=2E_{2\alpha}\left(-\beta^{2}\right).$

When $\alpha=1/2$ , we get the Gaussian density

$\frac{1}{2}M_{1/2,\rho}(-|y|)=\frac{1}{2\sqrt{\pi}\rho^{-1}}\exp\bigg{(}\frac{% -y^{2}}{4\rho^{-2}}\bigg{)}$ (32)

with mean zero and variance $2\rho^{-2}$ . It is easy to show that

$|US^{-\alpha}|\stackrel{d}{=}S^{-\alpha}.$ (33)

The preceding result allows us to estimate the parameters of the two-sided symmetric $M$ -Wright distribution using the properties of its one-sided non-symmetric counterpart. Furthermore, the formula for the integer-order moments of the symmetric two-parameter $M$ -Wright distribution centered at zero can be deduced as

$\mathbf{E}X^{\kappa}=\begin{cases}{\displaystyle\frac{\rho^{\kappa}\Gamma(1+% \kappa)}{\Gamma(1+\alpha\kappa)}}&\text{if}\;j\;\text{is even},\\ 0&\text{if}\;j\;\text{is odd}.\end{cases}$ (34)

For completeness, we reproduce Fig. 3 from (Cahoy, 2012b) to emphasize the flexibility of the symmetric single-parameter $M$ -Wright density.

Figure 2.

The symmetric $M$ -Wright density for $\alpha=$ 0.3, 0.5, 0.8; $\rho=$ 1, $\mu=$ 0.

Figure 3.

The asymptotic relative efficiency of $\overline{X}$ to $\widetilde{X}$ $(\textit{ARE}(\overline{X},$ $\widetilde{X}))$ as a function of $\alpha$ .

Case 3.2: $\mu\neq 0$

Proposition 3. Let $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}X$ in Eq. (31). Then

$\sqrt{n}\left(\overline{X}-\mu\right)\stackrel{d}{\longrightarrow}N\left(0,% \frac{\rho^{2}}{\alpha\Gamma(2\alpha)}\right)$ (35)

and

$\sqrt{n}\left(\widetilde{X}-\mu\right)\stackrel{d}{\longrightarrow}N\left(0,% \rho^{2}\Gamma(1-\alpha)^{2}\right),$ (36)

as $n\to\infty$ where $\widetilde{X}$ is the sample median.

Proof Directly follows from the standard large sample results for mean and median of random samples. ∎

Thus, the asymptotic relative efficiency of $\overline{X}$ to $\widetilde{X}$ $(\textit{ARE}(\overline{X},\widetilde{X}))$ is

$\textit{ARE}(\overline{X},\widetilde{X})=(\alpha\Gamma(2\alpha)\Gamma(1-\alpha% )^{2})^{-1}.$ (37)

Figure 3 displays the asymptotic relative efficiency of $\overline{X}$ to $\widetilde{X}$ as a function of $\alpha$ .

The relative efficiency above equals unity if $\alpha=$ 0.39106. Thus, the sample mean $\overline{X}$ is used for $\alpha>$ 0.39106. Otherwise, the sample median $\widetilde{X}$ is preferred when $\alpha<$ 0.39106 for relatively large samples.

Corollary 2. Let $X_{1},X_{2},\ldots,X_{n}\stackrel{iid}{=}X$ in Eq. (31). From Proposition 3, the approximate mean-based $(1-\nu)$ % confidence interval for $\mu$ is

$\overline{X}\pm z_{\nu/2}\widehat{\rho}\left(\widehat{\alpha}n\Gamma(2\widehat% {\alpha})\right)^{-1/2}$ (38)

while the approximate median-based $(1-\nu)$ % confidence interval for $\mu$ is

$\widetilde{X}\pm z_{\nu/2}\frac{\widehat{\rho}\Gamma(1-\widehat{\alpha})}{% \sqrt{n}}.$ (39)

Proof Directly follows from the central limit theorem and the asymptotic normality of the sample median. ∎

Subtracting $\widehat{\mu}$ from the data and getting the absolute values allow us to use the estimators of $\alpha$ and $\rho$ from the preceding section.

For testing purposes, we used the sample mean as the location parameter estimator as $\alpha$ values are chosen to be at least 0.4. Table 5 suggests negligible increase (due to the variability induced by subtracting the mean from the data) in both bias and MAD for the proposed point estimators of $\alpha$ and $\rho$ in comparison with Table 1 ( $\mu=0$ ) as $n\to\infty$ .

Table 5

Mean estimates of and dispersions from the true parameters $\alpha$ , $\rho$ , and $\mu$

$(\alpha,\rho,\mu)$		$n=$ 100		$n=$ 1000		$n=$ 10000
		% Bias	% MAD	% Bias	% MAD	% Bias	% MAD
(0.4, 150, $-$ 78)	$\widehat{\alpha}$	33.282	31.708	15.402	16.742	4.841	6.076
	$\widehat{\rho}$	7.821	9.771	2.928	3.365	0.958	1.233
	$\widehat{\mu}$	21.865	26.410	6.842	8.610	2.237	2.868
(0.6, 8.77, 25.2)	$\widehat{\alpha}$	17.082	20.509	6.114	7.422	1.959	2.536
	$\widehat{\rho}$	6.272	7.587	1.895	2.311	0.620	0.792
	$\widehat{\mu}$	3.691	4.536	1.136	1.389	0.387	0.494
(0.8, 375, 375)	$\widehat{\alpha}$	8.427	9.125	2.551	3.186	0.837	1.027
	$\widehat{\rho}$	4.479	5.532	1.328	1.636	0.435	0.559
	$\widehat{\mu}$	9.361	11.389	2.980	3.818	1.000	1.306
(0.95, 1000, 500)	$\widehat{\alpha}$	2.841	2.822	0.958	1.192	0.310	0.382
	$\widehat{\rho}$	2.536	2.927	0.695	0.885	0.217	0.271
	$\widehat{\mu}$	16.765	21.059	5.311	6.961	1.657	2.053

We also tested the proposed interval estimators and compared with their bootstrap counterparts (using percentile method). From Table 6, the large-sample interval estimator for $\alpha$ outperformed its bootstrap counterpart especially when $n=100$ .

Table 6

Coverage probabilities of 95% interval estimates for different values of $\alpha$ , $\rho$ , and $\mu$

$(\alpha,\rho,\mu)$		$n=100$	$n=1000$	$n=10000$
(0.4, 150, $-$ 78)	$\widehat{\alpha}$	0.906	0.956	0.951
	$\widehat{\rho}$	0.958	0.960	0.957
	$\widehat{\mu}$	0.940	0.955	0.955
	$\widehat{\alpha*}$	0.859	0.929	0.943
	$\widehat{\rho*}$	0.939	0.957	0.947
	$\widehat{\mu*}$	0.946	0.955	0.958
(0.6, 8.77, 25.2)	$\widehat{\alpha}$	0.956	0.955	0.953
	$\widehat{\rho}$	0.942	0.947	0.955
	$\widehat{\mu}$	0.952	0.956	0.952
	$\widehat{\alpha*}$	0.887	0.938	0.939
	$\widehat{\rho*}$	0.962	0.948	0.955
	$\widehat{\mu*}$	0.95	0.949	0.954
(0.8, 375, 375)	$\widehat{\alpha}$	0.969	0.959	0.945
	$\widehat{\rho}$	0.924	0.954	0.958
	$\widehat{\mu}$	0.940	0.953	0.945
	$\widehat{\alpha*}$	0.862	0.941	0.944
	$\widehat{\rho*}$	0.944	0.949	0.954
	$\widehat{\mu*}$	0.939	0.948	0.943
(0.95, 1000, 500)	$\widehat{\alpha}$	0.976	0.955	0.950
	$\widehat{\rho}$	0.891	0.950	0.953
	$\widehat{\mu}$	0.943	0.954	0.953
	$\widehat{\alpha*}$	0.781	0.908	0.942
	$\widehat{\rho*}$	0.910	0.953	0.951
	$\widehat{\mu*}$	0.953	0.951	0.949

4. Applications

We apply our methods on two real datasets that are available online (used in some researches) using the statistical software R. R codes are also available upon request through dcahoy@latech.edu.

4.1 Ages of major league baseball players

We consider the ages (in years) of 826 Major League Baseball (MLB) players. The data was downloaded from the Statistics Online Computational Resource (SOCR) database (see http://wiki.stat.ucla.edu/socr/index.php/SOCR). The one-sided $M$ -Wright fit to the data yields the point and interval estimates in Table 7. The minimum age of these players tends to be around 25 years old. The confidence interval estimate of the fractional parameter excludes the exponential ( $\alpha=0^{+}$ ) and the Airy ( $\alpha=1/3$ ) distributions but includes the half-normal ( $\alpha=1/2$ ) model. Using the asymptotic bivariate results in Section 2, the correlation between $\widehat{\alpha}$ and $\widehat{\rho}$ can be easily estimated as $-$ 0.989, which indicates a strong inverse linear relationship.

The two-sample Kolmogorov-Smirnov method (using R) was also used to test the fits of 100 simulated data sets (of same size with the observed data) using the parameter estimates. The average $p$ -value (0.841) indicated a reasonably good fit. The succeeding figure demonstrates the $M$ -Wright fit to the SOCR MLB age data with the maximum likelihood fits of gamma(shape=1.2994, rate=0.2605), Weibull(shape=1.2177, scale=5.3071) and lognormal(meanlog = 1.1752, sdlog=1.1292) distributions. By visual inspection, the one-sided $M$ -Wright distribution seems to provide the best fit. The picture also suggests that the one-sided $M$ -Wright had the flexibility to model data populations which have an inflection point (e.g., $\alpha=0.5$ : half-normal) with mode at the origin or minimum and their variants corresponding to $\alpha\approx 0.5.$ It can also be checked that at the origin, the height is $M_{\widehat{\alpha},\widehat{\rho},\widehat{\mu}}(\widehat{\mu})=(\widehat{% \rho}\cdot\Gamma(1-\widehat{\alpha}))^{-1}=0.1352.$

4.2 Human height and weight

The dataset contains 25000 records of human heights (in inches) and can be downloaded from the SOCR website. These data were obtained in 1993 by a Growth Survey of 25000 children from birth to 18 years of age recruited from Maternal and Child Health Centres (MCHC) and schools, and were used to develop Hong Kong’s current growth charts for weight, height, weight-for-age, weight-for-height and body mass index (BMI). Below are the corresponding point and 95% interval estimates for the three parameters. We used the sample mean as the point estimator as $\widehat{\alpha}$ is greater than the cutoff value of 0.39106 above. The interval estimate seems not to favor the double-exponential $(\alpha=0^{+})$ and normal or Gaussian $(\alpha=0.5)$ densities to likely model the distribution of the children’s heights. The estimate of the correlation between $\widehat{\alpha}$ and $\widehat{\rho}$ is $-$ 0.613, which indicates moderate negative association.

The two-sample Kolmogorov-Smirnov method (using R) was again used to test the fits of 100 simulated data sets (of same size with the observed data) using the parameter estimates above. The average $p$ -value (0.586) indicated a reasonably good fit to the data. The following figure demonstrated the fit of the model to the SOCR height data.

Table 7
Estimates for $\mu$ , $\alpha$ , and $\rho$

Parameter	Point estimate	95% Confidence interval
$\mu$	25.020	(24.960, 25.020)
$\alpha$	0.473	(0.338, 0.607)
$\rho$	4.390	(4.094, 4.686)

Table 8

Estimates for $\mu$ , $\alpha$ , and $\rho$

Parameter	Point estimate	95% Confidence interval
$\mu$	67.993	(67.969, 68.017)
$\alpha$	0.481	(0.457, 0.505)
$\rho$	1.352	(1.336, 1.369)

Figure 4.

Model fits to ages of MLB players.

Figure 5.

Symmetric $M$ -Wright fit to 25,000 heights of children from birth to 18 years of age.

5. Concluding remarks

Statistical inference procedures for the three-parameter $M$ -Wright family of distributions were proposed. The point estimators of the location, scale and fractional parameters were proven to be consistent and asymptotically unbiased. The large-sample results allowed quantification of the uncertainty associated with the proposed point estimators. The inference techniques were also demonstrated using real data sets, which indicated the ‘smoothing’ effect of the fractional parameter $\alpha\in(0,1)$ . The proposed location parameter estimators did not seriously affect the properties of the scale and fractional parameter estimates (point and interval). The random number generation algorithms were provided by the structural representations. Improvements of these procedures using robust or Bayesian perspectives and the derivation of the trivariate or joint asymptotic distribution of the location, scale, and fractional estimators would be worth exploring in the future.

Footnotes

Acknowledgments

The authors are grateful to the anonymous reviewers and co-editor-in-chief for their insightful comments and valuable suggestions that significantly improved the article.

References

Beghin

, & Orsingher

(2010). Poisson type processes governed by fractional and higher-order recursive diffferential equations. Electronic Journ Proby, (15), 684-709.

Cahoy

D. O.

(2012). Moment estimators for the two-parameter M-Wright distribution. Computational Statistics, 27(3), 487-497.

Cahoy

D. O.

(2012). Estimation and simulation for the M-Wright function. Communications in Statistics – Theory and Methods, 41(8), 1466-1477.

Cahoy

D. O.

, & Polito

(2012). Simulation and estimation for the fractional Yule process. Methodology and Computing in Applied Probability, 14(2), 383-403.

Hyndman

R. J.

, & Fan

(1996). Sample quantiles in statistical packages. American Statistician, 50, 361-365.

Mainardi

Mura

, & Pagnini

(2010). The M-Wright function in time-fractional diffusion processes: A tutorial survey. Int’l J of Diff’l Equations, 2010, 29. Article ID 104505, doi: 10.1155/2010/104505.

Mura

Taqqu

M. S.

, & Mainardi

(2008). Non-Markovian diffusion equations and processes: Analysis and simulations. Physica A, 387, 5033-5064.

Pagnini

, & Scalas

(2014). Historical notes on the M-Wright/Mainardi function, 2014. Communications in Applied and Industrial Mathematics, 6(1), DOI: 10.1685/journal.caim.496.

Piryatinska

Saichev

A. I.

, & Woyczynski

W. A.

(2005). Models of anomalous diffusion:the subdiffusive case. Physica A: Statistical Physics, 349, 375-424.

10.

Zolotarev

V. M.

(1986). One-dimensional stable distributions: translations of mathematical monographs. American Mathematical Society, 65, Printed in United States of America.

Inference for three-parameter 𝑴 -Wright distributions with applications

Abstract

Keywords

1. Introduction

4.1 Ages of major league baseball players

4.2 Human height and weight

Table 7 Estimates for μ , α , and ρ

Footnotes

Acknowledgments

References

Table 7
Estimates for $\mu$ , $\alpha$ , and $\rho$