Generalized modified slash distribution with applications

Generalized modified slash distribution with applications Jimmy Reyes ∗Inmaculada Barranco-Chamorro†H´ector W. G´omez‡ Abstract In this paper a generalization of the modified slash distribution is introduced. This model is based on the quotient of two independent random variables, whose distributions are a normal and a one-parameter gamma, respectively. The resulting distribution is a new model whose kurtosis is greater than other slash distributions. The probability density function, its properties, moments, and kurtosis coefficient are obtained. Inference based on moment and maximum likelihood methods is carried out. The multivariate version is also introduced. Applications to two datasets are considered. They illustrate the applicability of our results, along with the fact that the new model can provide a better fit to symmetric data with heavy tails than other slash extensions previously introduced in literature. KEY WORDS: Generalized modified slash distribution, Kurtosis, Maximum likelihood, Moments methods. ∗Departamento de Matem´aticas, Facultad de Ciencias B´asicas, Universidad de Antofagasta, Antofagasta, Chile. e-mail: [email protected] †Departamento de Estad´ıstica e Investigaci´on Operativa, Universidad de Sevilla, Sevilla, Espa˜na, e-mail: c[email protected] ‡Departamento de Matem´aticas, Facultad de Ciencias B´asicas, Universidad de Antofagasta, Antofagasta, Chile. e-mail: [email protected] 1 1 Introduction In real-world data, it is quite common to find symmetrical and unimodal histograms with heavy tails that do not fit well to a normal distribution. Slash models are a good option to deal with this kind of situations, in which departures of Gaussianity are a serious problem for the data analyst. This is one the main reason why slash distributions have received a great deal of attention during the last decades. In this context, we face the problem of improving slash models by introducing a generalisation able to model more kurtosis than other slash’s previously proposed in literature. In slash models, the emphasis is on kurtosis, because as Moors [20] pointed out the presence of heavy tails produces high kurtosis. Next, we briefly describe the most relevant precedents of our proposal. The canonic slash distribution was introduced by Rogers and Tukey [25]. This model was defined as the ratio of a N(0,1) normal and a uniform U(0,1) distribution independent. It was proposed as a model for bell shaped data with heavier tails than the corresponding normal distribution. Their theoretical properties can be seen in Rogers and Tukey [25], Mosteller and Tukey [21], and Johnson et al.[17]. Maximum likelihood estimation, with location and scale parameters, can be found in Kafadar [18]. Wang and Genton [28] described multivariate symmetrical and skew-multivariate extensions of the slash-distribution. G´omez et al.[12] and G´omez and Venegas [13] extended the slash distribution by introducing the slash-elliptical family. An asymmetric version of this family was discussed in the paper by Arslan [3]. Arslan and Genc [4] discussed a symmetric extension of the multivariate slash distribution, and Genc [9] discussed a symmetric generalization of the slash distribution. G´omez et al.[11] used the slash-elliptical family to extend the Birnbaum-Saunders distribution, and Genc [10] introduced a skew extension of the slash distribution based on the beta-normal distribution. Recently, Reyes et al.[24] proposed a modified slash distribution and Rojas et al.[26] extended the slash distribution by considering a beta distribution at the denominator. All these papers illustrate the relevance of slash methodology from a theoretical and practical point of view. In this paper, we focus on univariate symmetrical slash models. The starting point is the 2 stochastic representation of the (standard) slash distribution given by S=Z U1/q ,(1) where q > 0, Z∼N(0,1) and U∼U(0,1) are independent. For q= 1, (1) reduces to the distribution proposed by Rogers and Tukey [25], and for q→ ∞ the N(0,1) distribution is obtained as limit. Following the idea given in (1), Reyes et al.[24] introduced the modified slash (MS) distribution where the uniform distribution in the denominator is replaced by an exponential distribution. In this way, they turned the slash distribution into a more flexible model with the same number of parameters and a possible increase of kurtosis. So the stochastic representation for the MS distribution with location and scale parameters proposed by Reyes et al.[24] is Y=µ+σW V1/q ,(2) where W∼N(0,1) and V∼Exp(2) are independent, µ∈Ris the location parameter, σ > 0 is the scale parameter, and q > 0 is the parameter related to the distribution’s kurtosis. (2) is denoted as Y∼MS(µ, σ, q), and its pdf is fY(y;µ, σ, q) = 2q √2πσ2Z∞ 0 vqe−1 2[(y−µ σ)2v2+4vq]dv, y ∈R.(3) An explicit expression of (3) can be given only for certain values of q. For instance, if q= 1 then the pdf of Y∼MS(µ, σ, 1) is fY(y;µ, σ) =        (8πσ2)−1/2, y =µ, 2σ (y−µ)2 1 √2π−2σe 2σ2 (y−µ)2 |y−µ|Φ−2σ |y−µ|!, y =µ, (4) where Φ(·) denotes the cumulative distribution function (cdf) of the N(0,1) distribution. Another recent extension of the slash model was proposed by Rojas et al.[26]. These authors introduced the extended slash (ES) distribution in which a beta distribution is considered as denominator. So the stochastic representation of the ES distribution with location and scale parameters is Y=µ+σW V,(5) 3 where W∼N(0,1) and V∼Beta(α, β) are independent, µ∈Ris the location parameter, σ > 0 is the scale parameter α > 0 and β > 0 are the parameters related to the distribution’s kurtosis. (5) is denoted as Y∼ES(µ, σ, α, β) and its pdf is fY(y;µ, σ, α, β) = 1 σB(α, β)Z1 0 ϕy−µ σttα(1 −t)β−1dt, (6) with µ∈R,σ, α, β > 0 and B(·,·) denotes the beta function. Rojas et al.[26] showed that ES model exhibits a greater kurtosis than other slash models. The aim of this paper is to improve the modelling of data with heavy tails by increasing the kurtosis of slash models introduced so far. To achieve this end, the Generalized Modified Slash (GMS) distribution is introduced in which V1/q in (2) is replaced by a gamma distribution, Ga(2β, β), with β > 0. This parameterization of the gamma distribution is considered so that the MS model, proposed in Reyes et al.[24], can be considered a particular case of this new extension. The stochastic representation of GMS model is Y=µ+σW V(7) where W∼N(0,1) and V∼Ga(2β, β) are independent, µ∈Ris the location parameter, σ > 0 is the scale parameter, and β > 0 is the parameter related to the distribution’s kurtosis. (7) is denoted as Y∼GMS(µ, σ, β). We highlight that, in this way, the MS model is extended without increasing the number of parameters and achieving a greater degree of kurtosis. The parameters denoted as qand βcontrol the level of kurtosis in MS and GMS models, respectively. Additional benefits of using a GMS model instead of a MS distribution will be given throughout this paper. It is of interest to highlight the advantages of a GMS model versus other slash and t distributions. Throughout this paper, it will be shown that GMS allows us to model more kurtosis than other slash distributions. Moreover, for a GMS model, a hierarchical representation can be given in which, the gamma distribution will be used to describe the stochastic behaviour of the standard deviation of a normal distribution. On the other hand, for slash’s introduced in (1) and (2), they would be an uniform and an exponential distribution, respectively, to the power 1/q,q > 0. The choice of a gamma distribution seems much more 4 natural than the other ones. Details about this point are provided in Property 2 and Remark 1. As for t models, it is obvious that GMS model shares many properties with t-distributions, but they also have relevant differentiating features. We highlight the following ones: - The standard form of a GMS model is defined as X=Z Ga(2β,β),β > 0, and the tdistribution as tν=Z √Ga(2,ν/2)/ν ,ν∈Z+and Z∼N(0,1). In both cases, numerator and denominator are independent. Roughly speaking, it could be said that the main difference is the square root at the denominator of a t-model. However, they also have different uses. Along this paper, we will illustrate this point. - It is also worth mentioning that, in a different context, a variant of our proposal has been simultaneously introduced,by C.H. Nagaraja and H.N. Nagaraja [22]. These authors develop inference methods for central and intermediate quantiles based on the quotient of independent variables, denoted as W(m)=Z Ga(1,m). We cite this fact to highlight the interest of GMS distributions as a model by itself and its different applications in statistical inference. The paper is organized as follows. In Section 2 a closed expression for the pdf of a GMS model is obtained in terms of the confluent hypergeometric function; GMS model is expressed as a scale mixture; the convergence in law to a normal distribution is proven; moments are obtained, with emphasis on the kurtosis coefficient; and comparisons of pdf’s and kurtosis coefficients of slash, MS and GMS distributions are included to illustrate the superiority of our proposal. Section 3 is devoted to moment, iterative an EM maximum likelihood estimation methods. A simulation study has been included to assess the performance of our results. Applications to two real datasets are given in Section 4. There, graphic and descriptive criteria are proposed to show that the new distribution fits better than other slash models. We highlight that simulations and application have been carried out by using R software [23]. Finally, a multivariate version of the distribution GMS is introduced in Section 5 along with its more relevant properties. 5 2 GMS distributions Next the GMS model is introduced by replacing V1/q in (2) by a gamma distribution, Ga(2β, β) with β > 0. First, a closed expression for the pdf of this new model is obtained. This involves the confluent hypergeometric function of the second kind, denoted by U(a, b, z), and whose integral representation is U(a, b, z) = 1 Γ(a)Z∞ 0 ta−1(1 + t)b−a−1e−ztdt (8) for Re[a]>0 and Re[z]>0. Details about this function can be seen in Abramowitz and Stegun [1], p. 505. 2.1 Probability density function Proposition 1 Let Y∼GMS(µ, σ, β)with µ∈R,σ > 0, and β > 0. Then the pdf of Yis fY(y;µ, σ, β) =          1 σ√8π,if y=µ 2β/2 √2π σβ+1ββ+2 |y−µ|β+2 U1 + β 2,3 2,2σ2β2 (y−µ)2,if y=µ. (9) Proof. Given Y∼GMS(µ, σ, β), from (7), we can write Y=µ+σX with X∼ GMS(0,1, β), that is, X=W Vwith W∼N(0,1) and V∼Ga(2β, β) independent. In order to obtain the pdf of X, we consider the transformation X=W Vand Z=V, and by applying the Jacobian technique, we have that the joint pdf of (X, Z) is fX,Z (x, z) = (2β)β √2πΓ(β)zβe−x2z2 2e−2βz ,−∞ <x<∞, z > 0.(10) Marginalizing (10) with respect to Z, the pdf of Xis fX(x) = (2β)β Γ(β)√2πZ∞ 0 zβe−x2z2 2−2βz dz , x ∈R.(11) For x= 0, taking into account that Z∞ 0 zβe−2βzdz =Γ(β+ 1) (2β)β+1 , 6 we obtain that fX(0) = 1 2√2π.(12) For x∈R\{0}, after some algebra, it can be proved that Z∞ 0 zβe−x2z2 2−2βz dz =βΓ(β+ 1) 2β/2|x|β+2 U1 + β 2,3 2,2β2 x2,(13) where U(·) denotes the confluent hypergeometric function whose integral representation was introduced in (8). So, from (11), (12), and (13), we have fX(x) =          1 2√2π,if x= 0 2β/2 √2π ββ+2 |x|β+2 U1 + β 2,3 2,2β2 x2,if x= 0. (14) Since Y=µ+σX, the proposed result follows from (14) and the fact that fY(y) = 1 σfXy−µ σ.□ Figure 1 shows the pdf’s of GMS(0,1, β) distributions for β= 1,3,10 and the pdf of a N(0,2). We observe that for increasing values of βa lower degree of kurtosis is obtained. This appreciation, along with the fact that for β→ ∞ this model approaches to a normal distribution, will be formalized in Corollary 4 and Proposition 3, respectively. 7 −10 −5 0 5 10 0.00 0.05 0.10 0.15 0.20 x Density Beta=1 Beta=3 Beta=10 N(0,2) Figure 1: Comparison of GMS(0,1, β) pdf’s for β= 1 (black line), β= 3 (red line), β= 10 (blue line) and N(0,2) (green line). Next result relates GMS models to MS distributions introduced in Reyes et al.[24]. Specifically, we show that for β=q= 1 both models agree. Corollary 1 For β= 1, the pdf of Y∼GMS(µ, σ, 1) given in (9) reduces to the pdf of a modified slash MS(µ, σ, 1) distribution whose pdf has been given in (4), Reyes et al. [24] Proof. Let Y∼GMS(µ, σ, 1). Then Ycan be written as Y=µ+σX where X∼ GMS(0,1,1). From (9), we have that fX(x) =    1 √8π,if x= 0 1 √π 1 |x|3U3 2,3 2,2 x2,if x= 0. (15) It can be seen in Abramowitz and Stegun [1] that the confluent hypergeometric function can be expanded into simpler functions. Specifically, we have the following expansion of U3 2,3 2,2 x2in terms of the error function U3 2,3 2,2 x2=√2|x|−2√πe2/x2 1−sgn(x)erf √2 x!! (16) where sgn(·) and erf(·) denote the sign and the error function, respectively. 8 Recall that the error function is defined as erf(z) = 2 √πZz 0 e−t2dt, and the relationship with the cdf of a N(0,1) variable, Φ(·), is erf(z) = 2Φ(z√2) −1.(17) So, by applying (16), (17), and the relationship Φ(−z) = 1 −Φ(z) to (15), we obtain fX(x) =    1 √8π,if x= 0 2 x2n1 √2π−2e2/x2 |x|1−Φ2 |x|o,if x= 0.(18) From (18) and since Y=µ+σX, it is straighforward that the pdf of Yagrees with those proposed in Reyes et al.[24] and given in (4). □ Figure 2 compares the shape of the GMS density to the densities of MS and S. distributions for certain values of parameters. −10 −5 0 5 10 0.00 0.05 0.10 0.15 0.20 0.25 0.30 x Density GMS(0,1,2) MS(0,1,2) S(0,1,2) Figure 2: Comparison of GMS(0,1,2) (black line), MS(0,1,2) (blue line), and S(0,1,2) (red line). 9 n= 50 n= 100 n= 200 bµMbσMb βMbµMbσMb βMbµMbσMb βM β=5-0.00260 1.199 48.7 0.002477 1.162 18.9 -0.001318 1.1357 8.94 sd(b θ) 0.41234 0.198 1611.3 0.284744 0.136 233.0 0.203936 0.0963 5.85 qMSE(b θ) 0.41234 0.281 1611.9 0.284754 0.212 233.5 0.203940 0.1664 7.06 β=70.00253 1.110 49.2 -0.000794 1.092 32.1 0.001220 1.0700 14.54 sd(b θ) 0.36814 0.155 371.2 0.258214 0.109 508.3 0.180959 0.0818 51.92 qMSE(b θ) 0.36815 0.190 373.6 0.258215 0.143 508.9 0.180963 0.1077 52.47 β=10 -0.00338 1.059 80.3 0.000514 1.047 54.3 0.000477 1.0377 42.38 sd(b θ) 0.33233 0.137 1367.8 0.236552 0.101 808.7 0.168593 0.0757 851.26 qMSE(b θ) 0.33235 0.149 1369.6 0.236552 0.111 809.9 0.168593 0.0846 851.88 Table 2: Empirical means, sd, and √MSE for the MMEs of µ,σand βin GMS(0,1, β). Table 2 suggests that MM estimators of µand σbehave well when the sample size increases. Although, in general, the standar error and √MSE of b βMM decrease when n increases, this estimator exhibits a large variability. For this reason a second option to initialize MLE procedure is given. These simulations were carried out by using R software. Seed was set equal to 4 for reproducibility and code is available upon request to authors. 3.2 Maximum likelihood estimation Given Y1, . . . , Yna s.r.s from a GMS(µ, σ, β) distribution, the log-likelihood function is l(µ, σ, β) = nβ ln(2) + nβ ln(β)−nln(Γ(β)) −n 2ln(2π)−nln(σ) + n X i=1 ln(G(yi)) (26) with G(yi) = G(yi;µ, σ, β) = Z∞ 0 zβexp −1 2(yi−µ)2z2 σ2+ 4βzdz (27) 16 The maximum likelihood equations are n X i=1 G1(yi) G(yi)= 0 (28) n X i=1 G2(yi) G(yi)=n σ(29) n X i=1 G3(yi) G(yi)=n[Ψ(β)−ln(β)−ln(2) −1] (30) where Ψ(·) denotes the digamma function, that is, the logarithmic derivative of the gamma function, Γ(·), G1(yi) = ∂ ∂µ G(yi), G2(yi) = ∂ ∂σ G(yi),and G3(yi) = ∂ ∂β G(yi). For j > 0, let us denote by ai(j) = ai(µ, σ, j) = Z∞ 0 zjexp −1 2(yi−µ)2z2 σ2+ 4βzdz, bi(j) = bi(µ, σ, j) = Z∞ 0 ln(z)zjexp −1 2(yi−µ)2z2 σ2+ 4βzdz . The maximum likelihood equations can be rewritten in terms of ai(·) and bi(·) as n X i=1 (yi−µ) σ2 ai(β+ 2) ai(β)= 0 (31) n X i=1 (yi−µ)2 σ3 ai(β+ 2) ai(β)=n σ(32) n X i=1 bi(β) ai(β)− n X i=1 ai(β+ 1) ai(β)=n[Ψ(β)−ln(β)−ln(2) −1] (33) Taking xi= (yi−µ)/σ,w1(xi) = ai(β+2) ai(β), and w2(xi) = bi(β)−ai(β+1) ai(β), the previous equations are equivalent to n X i=1 yiw1(xi) = µ n X i=1 w1(xi) (34) n X i=1 x2 iw1(xi) = n(35) 1 n n X i=1 w2(xi) = Ψ(β)−ln(β)−ln(2) −1,(36) 17 which may be solved for µ,σ, and βby using numerical methods such as an iterative Newton-Raphson procedure bµ=Pn i=1 yiw1(xi) Pn i=1 w1(xi)(37) bσ2=σ2 nn X i=1 x2 iw1(xi) (38) b β= exp (Ψ(β)−ln(2) −1−1 n n X i=1 w2(xi)).(39) Since, when sampling from GMS models, there is no guarantee of existence of moment method estimators, following Kafadar [18], to start the iterations we propose, as initial estimators for µthe sample median, bµ0=me, and for σ, the median of the sample absolute deviations, bσ0=medi|yi−bµ0|. As for β, we propose b β0such that it maximizes the profile log-likelihood l(b β0) = maxβ{l(bµ0,bσ0, β)}. 3.3 Observed information matrix Recall that, the Fisher information matrix is given by I(µ, σ, β) = −     jµµ jµσ jµβ jµσ jσσ jσβ jµβ jσβ jββ     , which entries are equal to minus the second partial derivatives of the log-likelihood given in (26) with respect to the model parameters, that is jµµ =∂2 ∂µ2l(µ, σ, β), and so on. Similarly to previous subsection, let us denote by G11(yi) = ∂2 ∂µ2G(yi), G12(yi) = ∂2 ∂µ∂σ G(yi), etc. We have that jµµ = n X i=1 G11(yi)G(yi)−{G1(yi)}2 {G(yi)}2, jµσ = n X i=1 G12(yi)G(yi)−G2(yi)G1(yi) {G(yi)}2 jµβ = n X i=1 G13(yi)G(yi)−G3(yi)G1(yi) {G(yi)}2, jσσ = n X i=1 G22(yi)G(yi)−{G2(yi)}2 {G(yi)}2+n σ2 jσβ = n X i=1 G23(yi)G(yi)−G3(yi)G2(yi) {G(yi)}2, jββ = n X i=1 G33(yi)G(yi)−{G3(yi)}2 {G(yi)}2+n β−nΨ′(β) 18 where G11(yi) = −1 σ2ai(β+ 2) + (yi−µ)2 σ4ai(β+ 4) (40) G12(yi) = −2(yi−µ) σ3ai(β+ 2) + (yi−µ)3 σ5ai(β+ 4) (41) G13(yi) = (yi−µ) σ2[bi(β+ 2) −2ai(β+ 3)] (42) G22(yi) = −3(yi−µ)2 σ4ai(β+ 2) + (yi−µ)4 σ6ai(β+ 4) (43) G23(yi) = (yi−µ)2 σ3[bi(β+ 2) −2ai(β+ 3)] (44) G33(yi) = ci(β)−2ai(β+ 1) −2bi(β+ 1) + 4ai(β+ 2),(45) and ci(j) = ci(µ, σ, j) = Z∞ 0{ln(z)}2zjexp −1 2(yi−µ)2z2 σ2+ 4zβdz . Under regularity conditions, the asymptotic distribution of MLEs, (bµ, bσ, b β), is normal with mean vector (µ, σ, β)′and covariance matrix I(µ, σ, β)−1. Standard large-sample procedures for MLEs can be applied by considering the observed matrix information, I(bµ, bσ, b β), instead of I(µ, σ, β). 3.4 ML Estimation using EM-Algorithm The EM-algorithm (Dempster et al. [7]) is a well-known tool to carry out ML estimation when unobserved (missing) data or latent variables are present. This algorithm enables the computationally efficient determination of the ML estimates when iterative procedures are required. Our proposal is based on the stochastic representation of a GMS distribution given in Proposition 2. Next the implementation of the ML estimation of parameters of GMS distributions by using the EM-algorithm is given. Let Y1, . . . , Ynbe a simple random sample of size nfrom Y∼GMS(µ, σ, β). The parameter vector is θ= (µ, σ, β)⊤, with θ∈Θ⊆R×R2 +. In this setting, next hierarchical 19 representation can be proposed Yi|Vi=vi∼Nµ, σ vi,(46) Vi∼Ga(2β, β), i = 1, . . . , n. (47) Let y= [y1, . . . , yn]⊤and v= [v1, . . . , vn]⊤be observed and unobserved data, respectively. The complete data yc= [y⊤,v⊤]⊤corresponds to the original data yaugmented with v. Let ℓc(θ|yc) and Q(θ|b θ) = E[ℓc(θ|yc)|y,b θ] denote the complete-data log-likelihood function and its expected value, respectively. The hierarchical representation given in (46) and (47) is used to obtain the complete log-likelihood function associated with yc, which can be expressed as ℓc(θ|yc)∝ −nln(σ)− n X i=1 (yi−µ)2 2σ2v2 i+ℓc(β|yc),(48) where ℓc(β|yc)∝ −nln Γ(β) + nβ ln(2β)+(β−1) n X i=1 ln vi−2β n X i=1 vi. Letting b v2 i=E[V2 i|yi,θ=b θ], it follows that the conditional expectation of the complete log-likelihood function has the form Q(θ|b θ)∝ −nln(σ)− n X i=1 (yi−µ)2 2σ2b v2 i+Q(β|b θ),(49) where Q(β|b θ) = E[ℓc(β|yc)|y,b θ] = −nln Γ(β) + nβ ln(2β) + (β−1)S1n−2βS2n,with S1n= Pn i=1 E[ln Vi|yi] and S2n=Pn i=1 E[Vi|yi]. As both quantities S1nand S2nhave no explicit forms in the context of our model, they must be numerically computed. Thus to compute Q(β|b θ) we use a similar approach to that given in Lee and Xu [19] (Section 3.1). Specifically, let {vr;r= 1, ..., R}be a sample randomly drawn from the conditional distribution V|(Y=y, θ=b θ), so the quantity Q(β|b θ) can be approximated as follows: Q(β|b θ)≈1 R R X r=1 ℓc(β|vr). 20 Recall that each iteration of the EM algorithm involves expectation (E) and maximization (M) steps. So the EM-algorithm for the ML estimation of the parameters of the GMS distributions is: E-step. Given θ=b θ(k)= (bµ(k),bσ(k),b β(k))⊤, compute b v2 i (k), for i= 1, . . . , n. By maximizing Q(θ|b θ(k)) over µ,σ, and β, we obtain next results. CM-step I: bµ(k)must be updated as bµ(k+1) = n X i=1 yib v2 i (k) n X i=1 b v2 i (k).(50) CM-step II: bσ(k)must be updated as bσ(k+1) = Pn i=1 yi−bµ(k+1)2b v2 i (k) n  1/2 (51) CM-step III: Obtain b β(k+1) as the solution of 1 + ln(2β(k+1))−Ψ(β(k+1)) = 2S(k) 2n−S(k) 1n n(52) where S(k) 1n=Pn i=1 E[ln Vi|yi,b θ(k)], S2n=Pn i=1 E[Vi|yi,b θ(k)], and Ψ(·) denotes the digamma function, The iterations are repeated until a suitable convergence rule is satisfied, say |ℓ(b θ(k+1))− ℓ(b θ(k))|sufficiently small. As useful starting values required to implement this algorithm can be used the estimates proposed in previous sections. Remark 5 Note that CM-Step III requires an one-dimensional search for the root of β, which can be easily achieved by using the “uniroot” function built in R. As for the standard errors (s.e.) of these estimators, they can be calculated from the observed information matrix. 21 3.5 Simulation study In this subsection a simulation study is carried out in order to illustrate the behaviour of MLEs of the parameters. By using the representation given in (7), it is possible to generate random numbers for the GMS(µ, σ, β) distribution, which leads to the following algorithm. 1. Simulate Wi∼N(0,1), i = 1,2, ..., n. 2. Simulate Vi∼Ga(2β, β), i = 1,2, ..., n. 3. Compute Yi=σWi Vi+µ, i = 1,2, ..., n. It then follows that Y∼GMS(µ, σ, β). Table 3.5 shows results of simulations studies, which illustrate the behaviour of the MLEs for 1000 samples of sizes n= 50, 100, and 200 generated from a population distributed as GMS(µ, σ, β) for different values of µ,σand β, with initial values (0.5,1,1). For each generated sample, MLEs were numerically computed by using a Newton-Raphson procedure. Means, standard deviations (sd), and the squared root of the mean squared error √MSE of simulations are given, along with the empirical coverage (C) of 95% asymptotic confidence intervals for the parameters. 22 nµ σ β bµ sd(bµ)pMSE(bµ)C(bµ)bσ sd(bσ)pMSE(bσ)C(bσ)b β sd(b β)qMSE(b β)C(b β) 50 0,5 1,0 1,0 0,4137 0,3630 0,3567 95,9 1,0574 0,3770 0,3640 96,0 1,1075 0,4719 0,4626 81,9 2,0 0,4998 0,3447 0,3445 92,4 1,0132 0,2162 0,2165 94,0 2,3504 0,8692 0,9368 75,4 3,0 0,5007 0,3379 0,3377 92,4 0,9735 0,1830 0,1848 93,8 3,1130 0,8540 0,8610 74,3 1,0 1,0 1,0 1,0280 0,3445 0,3424 92,2 1,0854 0,3363 0,3439 92,9 1,2331 0,6304 0,6665 77,4 2,0 0,9380 0,6553 0,6545 93,4 2,0368 0,6276 0,6251 97,4 1,0850 0,3156 0,3251 89,6 3,0 1,1525 0,9994 0,9925 95,2 2,8606 0,6667 0,6689 95,4 1,1082 0,1855 0,2118 91,9 1,0 2,0 2,0 1,0382 0,6694 0,6701 92,6 2,0035 0,4256 0,4253 94,9 2,3387 0,8806 0,9430 91,3 2,0 1,9544 0,7498 0,7490 93,3 2,0281 0,4056 0,4054 94,9 2,3668 0,8840 0,9547 92,0 3,0 3,0043 0,7125 0,7122 95,2 2,0066 0,4372 0,4370 96,0 2,2961 0,8530 0,9025 89.8 100 0,5 1,0 1,0 0,4973 0,2306 0,2268 94,4 1,1069 0,3332 0,3448 95,1 1,1083 0,1801 0,2076 88.5 2,0 0,5043 0,2503 0,2501 93,9 1,0129 0,3397 0,3396 95,1 2,2032 0,6380 0,6689 84.0 3,0 0,5017 0,2392 0,2391 91,7 0,9953 0,1273 0,1273 95,2 3,1291 0,6212 0,6458 84.3 1,0 1,0 1,0 0,9197 0,2288 0,2314 94,6 1,0170 0,2850 0,2709 96,8 1,0229 0,2312 0,2205 87.4 2,0 1,0418 0,5025 0,4961 92,9 2,0704 0,3038 0,3070 96,0 1,0588 0,1337 0,1441 84.0 3,0 1,2912 0,8644 0,8891 93,6 3,2206 0,5193 0,5508 95,2 1,0460 0,2505 0,2478 82.2 1,0 2,0 2,0 0,9636 0,4767 0,4774 91,9 2,0279 0,3007 0,3016 96,7 2,2329 0,6219 0,6632 77.2 2,0 1,9693 0,4879 0,4875 95,7 2,0037 0,2921 0,2913 90,5 2,1427 0,6338 0,6480 73,0 3,0 3,0010 0,4694 0,4692 95,0 2,0276 0,7315 0,7317 94,4 2,1592 0,6136 0,6336 81,6 200 0,5 1,0 1,0 0,5737 0,1376 0,1540 94,6 1,0691 0,2103 0,2179 92,9 1,0526 0,1475 0,1542 83,1 2,0 0,4994 0,1685 0,1684 95,8 1,0039 0,1075 0,1075 93,7 2,1048 0,4074 0,4205 86,6 3,0 0,5089 0,1615 0,1614 96,4 1,0068 0,0922 0,0922 88,0 3,1763 0,5964 0,6204 75,3 1,0 1,0 1,0 0,9986 0,1530 0,1515 94,5 1,0928 0,2830 0,2951 91,7 1,1277 0,4434 0,4571 81,4 2,0 0,9527 0,2064 0,1943 94,8 2,4641 0,2167 0,5045 93,1 1,1139 0,1144 0,1545 83,8 3,0 1,0380 0,3945 0,3622 94,7 3,1218 0,4908 0,4643 91,5 1,0972 0,1423 0,1623 84,5 1,0 2,0 2,0 0,9629 0,3344 0,3357 98,3 1,9879 0,2003 0,2002 95,4 2,0517 0,3477 0,3507 92,7 2,0 1,9684 0,3313 0,3314 97,0 2,0266 0,2594 0,2596 96,9 2,1060 0,4979 0,5069 94,8 3,0 3,0153 0,3845 0,3822 95,2 2,0259 0,2276 0,2275 96,5 2,0905 0,4520 0,4580 96,9 Table 3: Mean, standard errors (sd), and √MSE for MLEs, along with empirical coverage (C) for the 95% confidence intervals based on them. The results in Table 3.5 show that the bias becomes smaller as the sample size n increases. In all cases, standard errors and the √MSE’s are quite similar. Moreover, both measures decrease as the sample size increases. Therefore, MLEs are consistent, as one would expect. As for the empirical coverage (C), we highlight that Cimproves if nincreases and Cbehaves better for bµand bσthan for b β. 4 Illustrations Next two real data sets are considered in order to illustrate the applicability of our results. 23 4.1 Illustration I We consider here data coming from an entomology experiment concerning ants. A total of n= 730 ants were placed singly in the center of an arena. The measurements correspond to the initial direction in which they moved in relation to a visual stimulus at an angle of 180 degrees from the zero direction, rounded to the nearest 10 degrees. The data were first presented in Jander [15], later analyzed in Batschelet [5], SenGupta and Pal [27] and Jones and Pewsey [16]. Figure 4 shows the boxplot of these data and Table 4 the statistical summary. ●●●●●●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●●●●●●● ●●●●●●●●●● ●●●●●●●●●● ●●●●●●●●●● ●●●●●●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●● ●●●●●●●●●● 0 50 100 150 200 250 300 350 Direction Figure 4: Boxplot for ants dataset Table 4: Descriptive statistics for ants data set n y S g1g2 730 176.4384 62.6434 −0.2053 4.5879 Figure 4 and Table 4 show the kind of situations in which (symmetric) slash models can be appropriate: symmetric distributions whose tails are much heavier than the normal ones (the typical boxplot would be similar to that given in Figure 4), and a really high sample kurtosis coefficient (g2= 4.5879). Slash models allow to incorporate these heavy tails into 24 the model, that, from the normal point of view would be labelled as outliers. Among the competing slash models previously proposed in literature (Slash (S), Modified Slash (MS), and Extended Slash (ES)), the GMS distributions allow us to model more kurtosis, as it will be next illustrated. We use the AIC (Akaike information criterion), which penalizes the maximized likelihood function by the excess of model parameters (AIC =−2log(lik)+2k, where kis the number of unknown parameters being estimated, see Akaike, 1974). This criterion reveals that the GMS model presents the best fit for the dataset under study. We also computed the Kolmogorov-Smirnov statistic (KSS), whose corresponding values for the fitted models are given in Table 5. KSS also indicates that the best fit is provided by the GMS model. Moreover, Figure 5 depicts the empirical cdf and the estimated cdf’s for S, MS, ES and GMS models on the left hand side and the profile log-likelihood function of βon the right hand side. On the other hand, Figure 6 shows the qqplots for these models. From all these summaries and plots, we can conclude that GMS model provides the best fit to our data. Finally, it is worth pointing out that, by using results in Section 3.1, moment estimators were computed leading to the following values: bµM= 176.438, bσM= 27.368 and b βM= 11.978, which were used as initial estimates for the maximum likelihood approach. Parameter S MS ES GMS µ181.43 (1.269) 181.67 (1.217) 181.321 (0.094) 181.7447 (1.2291) σ16.8 (1.246) 16.7 (0.878) 1.336 (0.108) 14.7227 (0.8849) q 1.17 (0.086) 1.50 (0.034) 1.907 (0.094) β40.084 (4.719) 1.9789 (0.2012) AIC 7950.084 7921.282 7914.642 7911.628 KSS 0.131 0.115 0.105 0.070 Table 5: MLEs for the ants dataset and the corresponding standard errors (in parenthesis), AIC and KSS values 25 [17] Johnson, N.L., Kotz, S., Balakrishnan, N. (1995). Continuous univariate distributions, Vol 1, 2nd edn. New York: Wiley. [18] Kafadar, K. (1982). A biweight approach to the one-sample problem. J. Amer. Statist. Assoc. 77, 416-424. [19] Lee, S., and Xu, L. (2004). Influence analyses of nonlinear mixedeffects models. Comput. Stat. Data Anal. 45, 321-341. [20] Moors, J.J.A. (1988). A quantile alternative for kurtosis, J. R. Stat. Soc. Ser. D, 37, 25–32. [21] Mosteller, F., Tukey,J.W. (1977). Data analysis and regression, Addison-Wesley, Reading, MA. [22] Nagaraja, C.H., Nagaraja, H.N. Distribution-Free Methods for Constructing Confidence Intervals for Quantiles. (submitted for publication). [23] R Core Team (2015). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/. [24] Reyes, J., G´omez, H.W., Bolfarine, H. (2013). Modified slash distribution. Statistics: A Journal of Theoretical and Applied Statistics, 47(5), 929-941. [25] Rogers, W.H., Tukey, J.W. (1072). Understanding some long-tailed symmetrical distributions, Statist. Neerlandica, 26, 211-226. [26] Rojas, M.A., Bolfarine, H., G´omez, H.W. (2014). An extension of the slash-elliptical distribution. Statistics and Operations Research Transactions (SORT), 38(2), 215-230. [27] SenGupta, A., Pal, C. (2001) . On optimal tests for isotropy against the symmetric wrapped stable-circular uniform mixture family. Journal of Applied Statistics 28, 129-143. 32 [28] Wang, J., Genton, M.G. (2006). The multivariate skew-slash distribution. Journal Statistical Planning and Inference, 136, 209-220. 33