The Power Muth Distribution

Mathematical Modelling and Analysis Publisher: Taylor&Francis and VGTU Volume 22 Number 2, March 2017, 186–201 http://www.tandfonline.com/TMMA https://doi.org/10.3846/13926292.2017.1289481 ISSN: 1392-6292 c Vilnius Gediminas Technical University, 2017 eISSN: 1648-3510 The Power Muth Distribution∗ Pedro Jodr´aa, H´ector Wladimir G´omezb, Mar´ıa Dolores Jim´enez-Gamerocand Mar´ıa Virtudes Alba-Fern´andezd aDpto. de M´etodos Estad´ısticos, Universidad de Zaragoza Mar´ıa de Luna 3, 50018 Zaragoza, Spain bDpto. de Matem´aticas, Universidad de Antofagasta 601 Avda. Angamos, Antofagasta, Chile cDpto. de Estad´ıstica e Investigaci´on Operativa, Universidad de Sevilla Avda. Reina Mercedes s.n., 41012 Sevilla, Spain dDpto. de Estad´ıstica e Investigaci´on Operativa, Universidad de Ja´en Paraje Las Lagunillas s.n., 23071 Ja´en, Spain E-mail(corresp.): [email protected] E-mail: [email protected] E-mail: [email protected] E-mail: [email protected] Received September 29, 2016; revised January 24, 2017; published online March 1, 2017 Abstract. Muth introduced a probability distribution with application in reliability theory. We propose a new model from the Muth law. This paper studies its statistical properties, such as the computation of the moments, computer generation of pseudo-random data and the behavior of the failure rate function, among others. The estimation of parameters is carried out by the method of maximum likelihood and a Monte Carlo simulation study assesses the performance of this method. The practical usefulness of the new model is illustrated by means of two real data sets, showing that it may provide a better fit than other probability distributions. Keywords: Muth distribution, survival analysis, estimation, simulation, data analysis. AMS Subject Classification: 60E05; 62F10; 33B30. 1 Introduction In the last decades, there exists an increasing interest in the development of new parametric distributions with the aim of providing flexible probability models ∗Research in this paper has been partially funded by the research projects: CTM2015- 68276-R MINECO/FEDER,UE (M.V. Alba-Fern´andez), Diputaci´on General de Arag´on –Grupo consolidado PDIE– (P. Jodr´a), FONDECYT 1130495 (H.W. G´omez) and MTM2014-55966-P of the Spanish Ministry of Economy and Competitiveness (M.D. Jim´enez-Gamero). The Power Muth Distribution 187 useful in many different areas. In 1977, Muth [25] introduced a continuous probability distribution in the context of reliability theory. However, this law has been overlooked in the literature until a recent paper by Jodr´a et al. [20], where its main statistical properties were thoroughly studied. To be more precise, a random variable Yis said to have a Muth distribution if its probability density function (pdf) is given by fY(y;α)=(eαy −α) expnαy −1 α(eαy −1) o, y > 0, where α∈(0,1] is a shape parameter. The cumulative distribution function (cdf) of Yis the following FY(y;α)=1−expnαy −1 α(eαy −1) o, y > 0.(1.1) As a natural extension, Jodr´a et al. [20] considered the scaled Muth distribution, which is defined as βY , with β > 0, and they showed the usefulness of the scaled Muth law for modelling rainfall data. In this paper, we introduce a new probability distribution from the Muth law. Specifically, we define the so-called power Muth (PM) distribution by means of the transformation X=βY 1/γ, where β > 0, γ > 0 and Yis the Muth law with parameter α= 1. As it will be seen, the introduction of the new parameter γleads to a rich class of probability distributions for nonnegative random variables with a wide range of values for the asymmetry and kurtosis coefficients, increasing generalized failure rate as well as increasing or bathtub shape failure rate. The choice α= 1 is adopted to avoid unnecessarily increasing the number of parameters since both αand γare shape parameters. From the viewpoint of applications, the new family of distributions may provide a better fit than other probability distributions previously used to this end. Using ordinary results related to the transformation of variables, it is easy to see that the pdf and cdf of Xare given, respectively, by f(x;β, γ) = γ βγxγ−1e(x/β)γ−1exp n(x/β)γ−e(x/β)γ−1o, x > 0, F(x;β, γ)=1−exp n(x/β)γ−e(x/β)γ−1o, x > 0,(1.2) where βis a scale parameter and γis a shape parameter. Figure 1 displays the pdf of the PM distribution for different values of the parameters. Throughout this paper, the PM law with parameters βand γis denoted by PM(β, γ). It is interesting to note that f(x;β, γ) is a well-defined pdf for any β > 0 and γ > 0. However, routine calculations show that lim x→0+f(x;β, γ) =      ∞,if γ∈(0,0.5), 0.5/β, if γ= 0.5, 0,if γ∈(0.5,∞). Commonly used pdf estimators, such as the histogram or kernel based estimators, always yield a finite pdf estimation. Therefore, from a practical point Math. Model. Anal., 22(2):186–201, 2017. 188P. Jodr´a, H.W. G´omez, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez 012345 0.0 0.2 0.4 0.6 0.8 γ =0.5 x 012345 0.0 0.2 0.4 0.6 0.8 γ =0.75 x 012345 0.0 0.2 0.4 0.6 0.8 γ =1 x 0.0 0.5 1.0 1.5 2.0 0 1 2 3 4 γ =2 x 0.0 0.5 1.0 1.5 2.0 0 1 2 3 4 γ =3 x 0.0 0.5 1.0 1.5 2.0 0 1 2 3 4 γ =5 x Figure 1. Pdf of the PM(β, γ) distribution for β= 1 and different values of γ. of view, an unbounded pdf will not provide a reasonable fit for any data set. Because of this reason, we will only consider γ∈[0.5,∞). Nevertheless, most properties studied for the proposed family are valid for any γ > 0. The remainder of this paper is organized as follows. Section 2 studies some statistical properties of the PM distribution. More precisely, the mode, the noncentral moments and moments of the order statistics are expressed in terms of the generalized integro-exponential function. This function has also appeared in the expression of the moments of the order statistics of other distributions such as the Gompertz–Makeham distribution (see Jodr´a [18]). The quantile function is written in terms of the Lambert Wfunction. This function has found nice applications in probability and statistics for calculating certain statistical distances between some discrete distributions (see Adell and Jodr´a [1]) and for the expression of the quantile function of other probability laws (see [16,17,19], among others). The behavior of the (generalized) failure rate function and the mean residual life is also described. The parameter estimation problem is considered in Section 3 and a simulation study is carried out to assess the performance of the maximum likelihood method. In Section 4, an application to two real data sets is presented to illustrate the practical usefulness of the proposed distribution. The Power Muth Distribution 189 2 Statistical properties 2.1 Mode The mode of a continuous probability distribution, denoted by mode(X), is the value at which the pdf has its maximum. The next result shows that the PM(β, γ) distribution is unimodal for γ≥0.5. We first introduce some notation. Let ∆(γ) = γ−1 γ, g(z) = z(ez−r1)(ez−r2) ez−1, where r1= 2 −ϕ,r1=ϕ2and ϕstands for the golden ration. Recall that ϕ can be defined as the positive solution of the equation x2−x−1 = 0, that is, ϕ=1 + √5 2≈1.618034. Therefore, r1≈0.381966 and r2≈2.618034. Clearly, g(0) = g(r2) = 0, g(z)<0, ∀z∈(0, r2), and g(z)>0, ∀z∈(r2,∞). Routine calculations show that gis a strictly convex function for z≥0 with a global minimum at r0∈(0, r2), specifically r0≈0.226874 and g(r0)≈ −1.059945. Proposition 1. Let Xbe a random variable having a PM(β, γ)distribution with β > 0and γ≥0.5. Then, Xhas a unique mode mode(X) = βς1/γ , where: (a) ςis the unique solution in z∈(r0, r2)of the equation g(z) = ∆(γ),(2.1) if γ∈[0.5,1),(b) ς=r2, if γ= 1,(c) ςis the unique solution in z∈(r2,∞) of Eq. (2.1), if γ > 1. Proof. In order to obtain the mode of Xwe have to solve with respect to x the equation (∂/∂x)f(x;β, γ) = 0, which is equivalent to solve the following equation (∂/∂x) log f(x;β, γ) = 0, where log is the natural logarithm. This is tantamount to solve with respect to zEq. (2.1), with z= (x/β)γ. Taking into account that ∆(γ) is a strictly increasing function of γsatisfying ∆(γ)∈[−1,0] if γ∈[0.5,1] and ∆(γ)>0 if γ > 1, as well as the properties of g(z) sketched above, the result follows. ut 2.2 Moments The moments can be written in terms of the generalized integro-exponential function, which is defined by the integral representation (see Milgram [24]) Em s(z) = 1 Γ(m+ 1) Z∞ 1 (log u)me−zuu−sdu, z ∈(−∞,∞),(2.2) where s∈(−∞,∞), m > −1 and Γstands for the gamma function. Math. Model. Anal., 22(2):186–201, 2017. 190P. Jodr´a, H.W. G´omez, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez Proposition 2. Let Xhave the PM(β, γ)distribution. Then, E[Xk] = e βkΓk γ+ 1E k γ−1 0(1) , k = 1,2, . . . . (2.3) Proof. From Eq. (1.2), for any k= 1,2, . . . , we have E[Xk] = Z∞ 0 xkdF(x;β, γ) =Z∞ 0 γ βγxk+γ−1e(x/β)γ−1exp x βγ −e(x/β)γ−1dx =e βkZ∞ 1 (log u)k/γ e−u(u−1)du, where in the last equality we have made the change of variable u= exp {(x/β)γ}. Then, by virtue of Eq. (2.2), we get E[Xk] = e βkΓ(k/γ + 1) nEk/γ −1(1) −Ek/γ 0(1)o. Finally, by taking into account in the above equation the following recurrence formula (see Milgram [24, eqn 2.4]) (1 −s)Em s(z) = zEm s−1(z)−Em−1 s(z), z > 0, s 6= 1, m ≥0,(2.4) where it is assumed E−1 s(z) = e−z, the desired result is obtained. ut The main mathematical properties of the generalized integro-exponential function have been studied by Ozalp and Bairamov [26], who provided an accurate algorithm to evaluate (2.2). As a consequence, from Eq. (2.3) we can compute efficiently the usual statistical measures involving E[Xk]. The behavior of some common measures can be graphically seen in the following figures, where we set β= 1 since βis a scale parameter: Figure 2 displays the mean, standard deviation, asymmetry coefficient (κ=E[(X−E[X])3]/σ3) and kurtosis coefficient (E[(X−E[X])4]/σ4). As can be seen, the PM family has a wide range of values for the coefficients of asymmetry and kurtosis, which provides a great flexibility in modelling data. Moreover, for γ∈[1.12,1.36] the asymmetry coefficient satisfies |κ| ≤ 0.1, leading to quasi-symmetric densities. As suggested in Vargo et al. [30], a diagram of asymmetry coefficient versus kurtosis is also shown in Figure 3. 2.3 Quantile function The PM distribution inherits the variate generation property from the Muth distribution, that is, its quantile function can also be written in closed form. Specifically, it can be expressed explicitly in terms of the Lambert Wfunction (see Corless et al. [8] for a review of the theory and applications of W). In this regard, it is interesting to note that the Lambert Wfunction is implemented in computer algebra systems such as Maple, Mathematica and Matlab and The Power Muth Distribution 191 0 5 10 15 20 1.00 1.10 1.20 mean γ 0 5 10 15 20 0.0 0.2 0.4 0.6 0.8 standard deviation γ 0 10 20 30 40 50 −1.0 0.0 1.0 asymmetry coefficient γ 0 10 20 30 40 50 2.5 3.5 4.5 5.5 kurtosis coefficient γ Figure 2. Mean, standard deviation, asymmetry and kurtosis as functions of γ(β= 1). also in programming languages such as R [29]. Therefore, pseudo-random data from the PM model can be computer-generated in a straightforward manner by virtue of the following result. Proposition 3. The quantile function of the PM(β, γ)distribution is F−1(u;β, γ) = βlog−W−1u−1 e1/γ ,0< u < 1,(2.5) where W−1denotes the negative branch of the Lambert Wfunction. Proof. For any 0 < u < 1, we have to solve with respect to xthe equation F(x;β, γ) = u, which is equivalent to solve FY((x/β)γ; 1) = u, with FYgiven by (1.1). Clearly, the solution of the latter equation is x=βF−1 Y(u; 1)1/γ. Then, taking into account the closed-form expression for F−1 Yprovided in Jodr´a et al. [20, Corollary 2], the result in (2.5) is obtained. ut 2.4 Failure rate function The failure (or hazard) rate function of the PM(β, γ) distribution is h(x;β, γ) = f(x;β, γ) 1−F(x;β, γ)=γ βγe(x/β)γ−1xγ−1, x > 0. Proposition 4. (a) If γ≥1then h(x;β, γ)is increasing in xfor any β > 0. (b) If 0< γ < 1then there exists an x0=x0(β, γ)>0so that h(x;β, γ)is (strictly) decreasing when x<x0and (strictly) increasing when x>x0. Math. Model. Anal., 22(2):186–201, 2017. 192P. Jodr´a, H.W. G´omez, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez −1.5 −1.0 −0.5 0.0 0.5 1.0 3456 asymmetry coefficient kurtosis coefficient Figure 3. Graph of asymmetry coefficient versus kurtosis coefficient. Proof. Let τ= (x/β)γ. The first derivative of hcan be written as follows h0(x;β, γ) = ∂ ∂x h(x;β, γ) = γ βγ(γ−1) (eτ−1) xγ−2+γ2 β2γeτx2(γ−1). Part (a) is evident from the expression of h0. To prove part (b), notice that h0(x;β, γ) = γ βγ(1 −γ) (eτ−1) {c(x)−1}xγ−2 with c(x) = γ βγ(1 −γ) eτ eτ−1xγ. Since for any x > 0, β > 0 and 0 < γ < 1, we have c0(x) = ∂ ∂x c(x) = γ2 βγ(1 −γ) eτ (eτ−1)2(eτ−τ−1)xγ−1>0, limx→0c(x) = 0, limx→∞ c(x) = +∞and γ βγ(1 −γ) (eτ−1) xγ−2>0, then there exists an x0=x0(β, γ)>0 so that h0(x;β, γ)<0 for x < x0and h0(x;β, γ)>0 for x>x0. This proves part (b). ut Therefore, the PM(β, γ) distribution has either increasing failure rate (IFR) when γ≥1 or bathtub-shaped failure rate when 0 < γ < 1. 2.5 Generalized failure rate Lariviere and Porteus [22] introduced the concept of generalized failure rate (GFR) of a continuous non-negative random variable and they showed that the distributions with increasing GFR (IGFR) have useful applications in operations management (see also Lariviere [21]). The GFR of a random variable Xis defined by g(x) = xh(x), where h(x) denotes the failure rate function. It is said that Xis IGFR if g(x) is non-decreasing. The next proposition shows that the PM distribution satisfies this desirable property. Proposition 5. The PM(β, γ)distribution is IGFR. The Power Muth Distribution 193 Proof. The first derivative of the GFR of the PM distribution can be written as follows g0(x;β, γ) = ∂ ∂x g(x;β, γ) = γ2 βγxγ−1 1 βγxγ+ 1e(x/β)γ−1 and clearly g0(x;β, γ)>0 since x > 0, β > 0 and γ > 0. Hence, g(x;β, γ) is a non-decreasing function in x, which implies the result. ut 2.6 Mean residual life Another important reliability measure for non-negative random variables is the mean residual life, which is defined as E(X−x|X > x). It is well-known that if a random variable has IFR, then its mean residual life is decreasing, which occurs when γ≥1 for the PM(β, γ) law as a consequence of Proposition 4. An analytical expression for the mean residual life is given in the following result. Proposition 6. The mean residual life of the PM(β, γ)distribution is r(x;β, γ) = β γ{log u(x)}1 γ−1I(x) (2.6) =β γΓ1 γ{log u(x)}1 γ−1exp {u(x)}E 1 γ−1 0(u(x)) (2.7) with I(x) = Z∞ 0{log(1 + z)}1 γ−1exp {−zu(x)}dz, u(x) = exp {(x/β)γ}. Proof. Let S(x;β, γ)=1−F(x;β, γ). The mean residual life is r(x;β, γ) = E(X−x|X > x) = Z∞ x S(y;β, γ) S(x;β, γ)dy =1 S(x;β, γ)Z∞ x exp y βγ −e(y/β)γ−1dy =β γ u(x)Z∞ u(x) log(t)1 γ−1exp {u(x)−t}dt, where in the last equality we have made the change of variable t= exp {(y/β)γ}. Then, the result in (2.6) follows by making the change of variable v=t−u(x). The result in (2.7) is obtained from (2.6) by taking into account (2.2). ut Figure 4 displays the mean residual life of the PM(β, γ) distribution for β= 1 and several values of γ. As can be seen, for γ∈(0,1) the mean residual life is first increasing up to a point and then decreasing, whereas for γ≥1 it is always decreasing. Math. Model. Anal., 22(2):186–201, 2017. 194P. Jodr´a, H.W. G´omez, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez 0 10 20 30 40 50 0.0 1.0 2.0 γ =0.5 x 0 5 10 15 20 0 2 4 6 8 12 γ =0.75 x 0 1 2 3 4 5 0 20 40 60 80 γ =1 x 0.0 0.2 0.4 0.6 0.8 1.0 0 4000 8000 γ =2 x Figure 4. Mean residual life for β= 1 and γ= 0.5,0.75,1,2. 2.7 Moments of order statistics This section shows that the moments of the order statistics can be expressed in closed form in terms of the generalized integro-exponential function. Let X1, . . . , Xnbe a random sample of size nfrom the PM(β, γ) distribution. Let X1:n≤X2:n≤. . . ≤Xn:nbe the order statistics obtained by arranging Xi,i= 1, . . . , n, in non-decreasing order of magnitude. For any n= 1,2, . . . and k= 1,2, . . . , the kth moment of Xr:n,r= 1, . . . , n, can be computed using the well-known formula (see Arnold et al. [4, p. 108]) E[Xk r:n] = rn rZ∞ 0 xk{F(x;β, γ)}r−1{1−F(x;β, γ)}n−rdF(x;β, γ),(2.8) where Fis defined by (1.2). We state the following. Proposition 7. For any n= 1,2, . . . , the moments of the smallest order statistic X1:nare given by E[Xk 1:n] = enβkΓk γ+ 1E k γ−1 1−n(n), k = 1,2, . . . . Proof. From Eqs. (1.2) and (2.8), together with the change of variable u= exp {(x/β)γ}, we have E[Xk 1:n] = nenβkZ∞ 1 (log u)k/γ une−nudu −Z∞ 1 (log u)k/γun−1e−nudu =nenβkΓk γ+ 1nEk/γ −n(n)−Ek/γ 1−n(n)o, The Power Muth Distribution 201 [17] P. Jodr´a. Computer generation of random variables with Lindley or Poisson– Lindley distribution via the Lambert Wfunction. Math. Comput. Simul., 81(4):851–859, 2010. https://doi.org/10.1016/j.matcom.2010.09.006. [18] P. Jodr´a. On order statistics from the Gompertz–Makeham distribution and the Lambert Wfunction. Math. Model. Anal.,18(3):432–445, 2013. https://doi.org/10.3846/13926292.2013.807316. [19] P. Jodr´a and M.D. Jim´enez-Gamero. A note on the Log- Lindley distribution. Insur. Math. Econ.,71:189–194, 2016. https://doi.org/10.1016/j.insmatheco.2016.09.005. [20] P. Jodr´a, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez. On the Muth distribution. Math. Model. Anal.,20(3):291–310, 2015. https://doi.org/10.3846/13926292.2015.1048540. [21] M.A. Lariviere. A note on probability distributions with increasing generalized failure rates. Oper. Res.,54(3):602–604, 2006. https://doi.org/10.1287/opre.1060.0282. [22] M.A. Lariviere and E.L. Porteus. Selling to the newsvendor: An analysis of priceonly contracts. Manufacturing &Service Operations Management,3(4):293–305, 2001. https://doi.org/10.1287/msom.3.4.293.9971. [23] S.G. Meintanis, M.D. Jim´enez-Gamero and M.V. Alba-Fern´andez. A class of goodness-of-fit tests based on transformation. Comm. Statist. Theory Methods, 43(8):1708–1735, 2014. https://doi.org/10.1080/03610926.2012.673673. [24] M.S. Milgram. The generalized integro-exponential function. Math. Comp., 44(170):443–458, 1985. https://doi.org/10.1090/S0025-5718-1985-0777276-4. [25] E.J. Muth. Reliability models with positive memory derived from the mean residual life function. In C.P. Tsokos and I. Shimi(Eds.), The Theory and Applications of Reliability, volume 2, pp. 401–435, New York, 1977. Academic Press, Inc. [26] N. Ozalp and E. Bairamov. Uniform convergence and computation of the generalized exponential integrals. J. Math. Chem.,49(2):520–530, 2011. https://doi.org/10.1007/s10910-010-9756-5. [27] G. Schwarz. Estimating the dimension of a model. Ann. Statist,6(2):461–464, 1978. https://doi.org/10.1214/aos/1176344136. [28] W. Stute, W. Gonz´alez Manteiga and M. Presedo Quindimil. Bootstrap based goodness-of-fit tests. Metrika,40(1):243–256, 1993. https://doi.org/10.1007/BF02613687. [29] R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org, 2016. [30] E. Vargo, R. Pasupathy and L. Leemis. Moment-ratio diagrams for univariate distributions. J. Qual. Technol.,42(3):276–286, 2010. [31] T. VedoVatto, A.D.C. Nascimento, W.R. Miranda Filho, M.C.S. Lima, L.G.B. Pinho and G.M. Cordeiro. Some computational and theoretical aspects of the exponentiated generalized Nadarajah–Haghighi distribution. Chil. J. Stat.,To appear. Math. Model. Anal., 22(2):186–201, 2017.