Inference for three-parameter M-Wright distributions with applications

TL;DR

This paper develops point and interval estimators for three-parameter M-Wright distributions, including methods for uncertainty quantification, and demonstrates their effectiveness on synthetic and real data.

Contribution

It introduces novel estimation procedures for the three-parameter M-Wright family, including uncertainty quantification and asymptotic covariance analysis.

Findings

Interval estimator for scale outperforms previous methods when location is zero.

Asymptotic covariance structure enables correlation estimation between parameters.

Proposed methods are validated on synthetic, age, and height data.

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 \cite{cah12} 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.

Figures

Equations

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

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

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

$$\mathbf{E}\left(e^{-\beta X}\right)={\varphi }_X(\beta )=\sum _{j=1}^{\infty }\frac{{\left(-\beta \right)}^j}{j!\Gamma (1+\alpha j)}$$

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

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

$${\mu }_x=\frac{1}{\alpha \Gamma (\alpha )},\text{and}{\sigma }_x^2=\frac{1}{\alpha \Gamma (2\alpha )}-\frac{1}{{\left(\alpha \Gamma (\alpha )\right)}^2},$$

$$\frac{{\sigma }_x}{{\mu }_x}=\sqrt{\frac{2\alpha \Gamma (\alpha )\Gamma (1+\alpha )}{\Gamma (1+2\alpha )}-1}=\{\begin{array}{ll}1& \alpha =0,\\ \sqrt{\pi /2-1}& \alpha =1/2,\\ 0& \alpha =1.\end{array}$$

$$M_{\alpha ,\rho ,\mu }(x)=\frac{1}{\rho }{M}_{\alpha }\left(\frac{x-\mu }{\rho }\right),x>\mu ,$$

$$\mathbf{E}{X}^{\kappa }=\frac{{\rho }^{\kappa }\Gamma (1+\kappa )}{\Gamma (1+\alpha \kappa )},{\varphi }_X(\beta )=E_{\alpha }(-\beta \rho ),{\varphi }_{\rho }(\beta )={\beta }^{1-1/\alpha }{e}^{-x\beta },\text{and}X\stackrel{d}{=}\rho S^{-\alpha }.$$

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

$${\mu }_{X^{{}^{\prime }}}=\log (\rho )+\gamma (\alpha -1),\text{and}{\sigma }_{X^{{}^{\prime }}}^2=\frac{{\pi }^2}{6}\left(1-{\alpha }^2\right),$$

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

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

$$\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,$$

$${\sigma }_{\hat{\alpha }\hat{\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 }_{\hat{\rho }\hat{\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},$$

$${\mu }_3^{{}^{\prime }}=2({\alpha }^3-1)\zeta (3)\text{and}{\mu }_4^{{}^{\prime }}=\frac{{\pi }^4({\alpha }^4-10{\alpha }^2+9)}{60},$$

$$\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{{\scriptstyle d}}{{\longrightarrow}}\textsl{N}\left[\bm{0},\bf{\Sigma}\right]$$

$$\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),$$

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

$$\text{g}({\mu }_{X^{{}^{\prime }}},{\sigma }_{X^{{}^{\prime }}}^2)={\left(\sqrt{1-\frac{6{\sigma }_{X^{{}^{\prime }}}^2}{{\pi }^2}},\exp \left({\mu }_{X^{{}^{\prime }}}+\gamma (1-\alpha )\right)\right)}^{\text{T}}$$

$$\dot{\text{g}}({\mu }_{X^{{}^{\prime }}},{\sigma }_{X^{{}^{\prime }}}^2)=$$

$$\left(\begin{array}[]{cc}0&\exp\left(\mu_{X^{{}^{\prime}}}+\gamma(1-\alpha)\right)\\ &\\ -3\Big{/}\left(\pi^{2}\sqrt{1-\frac{6\sigma_{X^{{}^{\prime}}}^{2}}{\pi^{2}}}\;\right)&\;\;\Big{(}3\gamma\exp(\mu_{X^{{}^{\prime}}}+\gamma)\sqrt{1-\frac{6\sigma_{X^{{}^{\prime}}}^{2}}{\pi^{2}}}\;\Big{)}/\left(\pi^{2}\sqrt{1-\frac{6\sigma_{X^{{}^{\prime}}}^{2}}{\pi^{2}}}\;\right)\\ \end{array}\right).\qed$$

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

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

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

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

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

$$P(\hat{\mu }>\mu +q_{\nu }\rho )$$

Full text

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:

[TABLE]

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

[TABLE]

Moreover,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

giving the mean and variance as

[TABLE]

correspondingly. The coefficient of variation is straightforward to calculate as

[TABLE]

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:

[TABLE]

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.

Case 2.1: $\bm{\mu=0}$

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

*where *

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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{{\scriptstyle iid}}{{=}}\;X$. The $(1-\nu)100\%$ confidence intervals for $\alpha$ and $\rho$ can be approximated as

[TABLE]

and

[TABLE]

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\leq 1000$. Note that the large-sample formula is faster to calculate than the resampling-based method especially for large sample sizes.

Case 2.2: $\bm{\mu\neq 0}$

Consider the location-scale structure

[TABLE]

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

[TABLE]

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

[TABLE]

which suggests that

[TABLE]

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<<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.

3 Symmetric M-Wright distribution

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

[TABLE]

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,

[TABLE]

where ‘$ind$’ means independent.

Case 3.1: $\bm{\mu=0}$

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Case 3.2: $\bm{\mu\neq 0}$

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

[TABLE]

and

[TABLE]

*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}$ $\left(ARE\left(\overline{X},\widetilde{X}\right)\right)$ is

[TABLE]

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{{\scriptstyle iid}}{{=}}\;X$ in (32). From Proposition 3, the approximate mean-based $(1-\nu)\%$ confidence interval for $\mu$ is

[TABLE]

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

[TABLE]

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$.

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$.

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$\_$Data$\_$

Dinov$\_$020108$\_$HeightsWeights). 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.

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.

6 Acknowledgment

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.

7 References

• Beghin and Orsingher (2010) Beghin, L., Orsingher, E., 2010. Poisson type processes governed by fractional and higher-order recursive diffferential equations. Electronic Journ Proby (15), 684-709.
• Cahoy (2012) Cahoy, D.O., 2012. Moment estimators for the two-parameter M-Wright distribution. Computational Statistics 27(3), 487-497.
• Cahoy (2012b) Cahoy, D.O., 2012. Estimation and simulation for the M-Wright function Communications in Statistics - Theory and Methods 41(8), 1466-1477.
• Cahoy and Polito (2012) Cahoy, D.O., Polito, F., 2012. Simulation and estimation for the fractional Yule process. Methodology and Computing in Applied Probability 14(2), 383-403.
• Hyndman and Fan (1996) Hyndman, R. J., Fan, Y. ,1996. Sample quantiles in statistical packages, American Statistician 50, 361-365.
• Mainardi et. al (2010) Mainardi, F., Mura, A., Pagnini, G., 2010. The M-Wright function in time-fractional diffusion processes: a tutorial survey. Int’l J of Diff’l Equations, vol. (2010), Article ID 104505, 29 pages, doi:10.1155/2010/104505.
• Mura et. al (2008) Mura, A., Taqqu, M.S., Mainardi, F., 2008. Non-Markovian diffusion equations and processes: analysis and simulations, Physica A, 387, 5033–5064.
• Pagnini and Scalas (2014) Pagnini, G., Scalas, E., 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 et. al (2005) Piryatinska, A., Saichev, A.I., Woyczynski, W.A., 2005. Models of anomalous diffusion:the subdiffusive case, Physica A: Statistical Physics 349, 375-424.
• Zolotarev (1986) Zolotarev, V.M. (1986) One-dimensional Stable Distributions: Translations of Mathematical Monographs. American Mathematical Society, vol 65, Printed in United States of America.