A New Weibull-Pareto Distribution: Properties and Applications M. H. TAHIR1 , GAUSS M. CORDEIRO2 , AYMAN ALZAATREH3 , M. MANSOOR4 AND M. ZUBAIR5 1 Department of Statistics, The Islamia University of Bahawalpur, Pakistan Email: mtahir.stat@gmail.com 2 Department of Statistics, Federal University of Pernambuco, Recife, PE, Brazil 3 Department of Mathematics and Statistics, Austin Peay State University, USA 4 Department of Statistics, Punjab College, Model-Town A, Bahawalpur, Pakistan 5 Department of Statistics, Government Degree College Kahrorpacca, Pakistan Abstract. Many distributions have been used as lifetime models. Recently, a generator of distributions called the Weibull-G class was proposed by Bourguignon et al. (2014). We propose a new three-parameter Weibull-Pareto distribution, which can produce the most important hazard rate shapes, namely constant, increasing, decreasing, bathtub and upsidedown-bathtub. Various structural properties of the new distribution are derived including explicit expressions for the moments and incomplete moments, Bonferroni and Lorenz curves, mean deviations, mean residual life, mean waiting time and generating and quantile functions. The Rényi and q entropies are also derived. We obtain the density function of the order statistics and their moments. The model parameters are estimated by maximum likelihood and the observed information matrix is determined. The usefulness of the new model is illustrated by means of two real data sets on Wheaton river flood and bladder cancer. In the two applications, the new model provides better fits than the Kumaraswamy-Pareto, beta-exponentiated Pareto, beta-Pareto, exponentiated-Pareto and Pareto models. Keywords: Hazard function; Likelihood estimation; Moment; Pareto distribution; Weibull-G class Mathematics Subject Classification (2010) 60E05; 62F10; 62N05 1. Introduction The Pareto distribution was pioneered by Vilfraddo Pareto (1896) to explore unequal distribution of wealth and income. It is widely used in modeling actuarial data (e.g. insurance risk) because of its heavy tail properties. It has also wider applications in hydrology, telecommunications, poverty measurement, migration, size of cities and firms, word frequencies, business mortality, service time in queuing theory, etc. A random variable T has the Pareto distribution with two parameters α > 0 and θ > 0, if its 1 cumulative distribution function (cdf) is given by ³ x ´−α Hα,θ (x) = 1 − , θ x > θ. Then, its probability density function (pdf) reduces to ³ α ´ ³ x ´−α hα,θ (x) = , x > θ > 0, x θ (1) α > 0. (2) Keeping in view of its applications, various authors have extended the Pareto distribution to explore properties and for efficient estimation of its parameters. The first parameter induction to the Pareto model called the generalized Pareto (GP) distribution was pioneered by Pickands (1975). Stoppa (1990) proposed another three-parameter generalization of the Pareto distribution known as the exponentiated Pareto (“EP” for short) distribution. It is also named the Stoppa distribution. The three-parameter EP cdf is given by ½ ³ x ´−α ¾λ , x > θ > 0, α, λ > 0. (3) Gα,θ,λ (x) = 1 − θ The density function corresponding to (3) becomes ³ α ´ ³ x ´−α ½ ³ x ´−α ¾λ−1 gα,θ,λ (x) = λ 1− , x θ θ x > θ > 0, α, λ > 0. (4) Gupta et al. (1998) discussed the EP distribution using the Lehmann alternative type I (due to Lehmann, 1953) by taking the λth power of the baseline cdf G(x), say G(x)λ . Akinsete et al. (2008) extended the Pareto distribution by adding two extra shape parameters based on the beta-G class of distributions pioneered by Eugene et al. (2002). The four-parameter beta-Pareto (BP) cdf is given by Z Hα,θ (x) 1 wa−1 (1 − w)b−1 dw = IHα,θ (x) (a, b), x > θ, a, b, α > 0, (5) F (x; a, b, α, θ) = B(a, b) 0 where Hα,θ (x) = 1 − (x/θ)−α is given by (1), Iw (a, b) is the incomplete beta function ratio, and a and b are two additional shape parameters whose role is to govern the skewness and tail weights. The beta-EP (BEP) distribution was defined by Nassar and Nada (2011), Zea et al. (2012) and Mansoor (2013) by using the beta-G class (Eugene et al., 2002) of distributions. So, the five-parameter BEP cdf becomes F (x; a, b, α, θ, λ) = IGα,θ,λ (x) (a, b), x > θ, a, b, α, λ > 0, (6) where Gα,θ,λ (x) is given by equation (3). Further, Mahmoudi (2011) defined the beta generalized Pareto (BGP) distribution and studied some useful properties for modeling extreme value data. Alzaatreh (2012) used the TransformedTransformer (T-X) family of distributions (Alzaatreh et al., 2011, 2013b) to define the gammaPareto (GaP) cdf by Z W (G(x)) π(x) dx, (7) F (x; a, b, θ) = 0 2 where W (G(x)) = − log[1 − G(x)] and π(x) = xa−1 e−x/b /[Γ(a) ba ], a, b > 0. Based on the T-X family, Alzaatreh et al. (2013a) defined the Weibull-Pareto-T-X (WPTX) distribution. If a random variable Z has the Weibull distribution with parameters α and β, then the WPTX cumulative function is given by x α F (x; α, β, θ) = 1 − e−[β log( θ )] , x > θ > 0 α, β > 0. (8) Bourguignon et al. (2013) defined the Kumaraswamy-Pareto (KwP) from the KummaraswamyG (Kw-G) class of distributions proposed by Cordeiro and de Castro (2011). The four-parameter KwP cdf is given by · ½ ³ x ´−α ¾a ¸b F (x; a, b, α, θ) = 1 − 1 − 1 − , θ x > θ > 0, a, b, α > 0. (9) Zagrafos and Balakrishnan (2009) pioneered a versatile and flexible gamma-G class of distributions based on Stacy’s generalized gamma distribution and record value theory. More recently, Bourguignon et al. (2014) proposed the Weibull-G class of distributions influenced by the ZografosBalakrishnan-G class. Let G(x; Θ) and g(x; Θ) denote the cumulative and density functions of a b baseline model with parameter vector Θ and consider the Weibull cdf πW (x) = 1 − e−x (for x > 0) with scale parameter one and shape parameter b > 0. Bourguignon et al. (2014) replaced the argument x by G(x; Θ)/G(x; Θ), where G(x; Θ) = 1 − G(x; Θ), and defined the cdf of their class of distributions, say Weibull-G(b, Θ), by F (x; b, Θ) = b Z 0 h G(x;Θ) G(x;Θ) i b xb−1 e−x dx = 1 − e i h G(x;Θ) b − G(x;Θ) Then, the Weibull-G class pdf is given by ¸ h G(x;Θ) ib · G(x; Θ)b−1 − G(x;Θ) e , f (x; b, Θ) = b g(x; Θ) G(x; Θ)b+1 , x ∈ ℜ, x ∈ ℜ, b > 0. b > 0. (10) (11) It is noteworthy to mention that (10) is a special case of the T-X family proposed by Alzaatreh et al. (2013b) in (7) by taking W (G(x)) = G(x; Θ)/G(x; Θ) and π(x) as the Weibull cdf. In this context, we propose an extension of the Pareto model called the Weibull-Pareto (“WP” for short) distribution based on equations (10) and (11). The proposed distribution is more flexible than the WPTX (Alzaatreh et al., 2013a) model. For example, the hazard function shapes of the WP distribution can be constant, increasing, decreasing, bathtub and upside down bathtub. Alzaatreh et al. (2013a) noted a significant problem to estimate the parameters of the their distribution using maximum likelihood. Further, they showed that the maximum likelihood estimates (MLEs) have considerable bias values. On the other hand, the MLEs of the WP parameters have lower bias values (see Table 1). The paper is outlined as follows. In Section 2, we define the WP distribution and discuss the shapes of the density and hazard rate functions. We provide a mixture representation for its density function in Section 3. Structural properties such as the ordinary and incomplete moments, Bonferroni and Lorenz curves, mean deviations, mean residual life, mean waiting time, generating 3 and quantile functions are derived in Section 4. In Section 5, we determine the Rényi and q entropies. The density of the order statistics is investigated in Section 6. The maximum likelihood estimation of the model parameters is discussed in Section 7. We explore the usefulness of the new model by means of two real data sets in Section 8. Finally, Section 9 offers some concluding remarks. 2. The WP Distribution Inserting (1) in equation (10), the three-parameter WP cdf is defined by x α F (x) = F (x; b, α, θ) = 1 − e−[( θ ) b −1] , x > θ > 0. (12) The pdf corresponding to (12) becomes f (x) = f (x; b, α, θ) = ib−1 b α b α α−1 h³ x ´α −[( xθ ) −1] x − 1 e , θα θ x > θ > 0. (13) Henceforth, a random variable X with density function (13) is denoted by X ∼WP(b, α, θ). The survival function (sf) (S(t)), hazard rate function (hrf) (h(t)), reversed-hazard rate function (rhrf) (r(t)) and cumulative hazard rate function (chrf) (H(t)) of X are given by x α S(x) = e−[( θ ) b −1] , x > θ, ib−1 bα α−1 h³ x ´α x − 1 , x > θ, θα θ ¤b−1 −[( x )α −1]b ¡ ¢ ¡ ¢α £¡ x ¢α −1 e θ b αx xθ θ r(x) = , x > θ, b α −[( xθ ) −1] 1−e h(x) = and H(x) = h³ x ´α θ −1 ib , (14) x > θ, respectively. Lemma 1 provides some relations of the WP distribution with the well-known Weibull and exponential distributions. Lemma 1 (Transformation): (a) If a random variable Y follows the Weibull distribution with shape parameter b and scale parameter one, then the random variable X = θ (Y + 1)1/α follows the W P (b, α, θ) distribution. (b) If a random variable Y follows the exponential distribution, then the random variable X = ¡ ¢1/α θ Y 1/b + 1 follows the W P (b, α, θ) distribution. Plots of the WP density and hazard rate functions for some parameter values are displayed in Figures 1 and 2, respectively. 4 θ = 1.5 3.0 (b) θ=1 2.0 (a) 2.5 b = 0.5 α = 1.2 b = 1.5 α = 2.5 b = 6 α = 1.2 b = 5 α = 1.5 b = 3.5 α = 1.8 1.5 Density 0.0 0.0 0.5 0.5 1.0 1.0 Density 2.0 1.5 b = 0.5 α = 1.5 b = 1.5 α = 2 b = 2 α = 1.5 b = 5 α = 0.9 b = 3.5 α = 0.8 1.0 1.5 2.0 2.5 3.0 1.5 2.0 2.5 x 3.0 x (d) θ = 0.5 3.0 (c) θ = 1.5 1.5 h(x) 1.0 1.5 h(x) 2.0 2.0 2.5 2.5 b=1 α=1 b = 0.22 α = 3.5 b = 1.1 α = 0.4 b = 1.03 α = 0.7 b = 1.7 α = 0.6 2 3 4 5 6 7 0.5 0.0 0.5 1.0 b = 0.35 α = 4.8 b = 0.4 α = 4.5 b = 0.45 α = 4.55 b = 0.48 α = 4.75 b = 0.5 α = 3.9 8 0.5 1.0 1.5 x 2.0 2.5 3.0 x Figure 1: Plots of the density (a)-(b) and hazard rate function (c)-(d) of WP distribution The limiting behaviors of the pdf and hrf of X are given in Theorem 1. Theorem 1: The limit of the pdf as x → ∞ is 0. Also, x → θ+ are given by    0, lim f (x) = lim h(x) = αθ ,  x→θ+ x→θ+  ∞, the limits of the pdf and hrf of X as if b < 1, if b = 1, if b > 1. Further, the limits of the hrf of X as x → ∞ are given by    ∞, if αb < 1, lim h(x) = 1θ , if αb = 1, x→∞   0, if αb > 1. 5 (15) (16) Proof: Most results of this type can be easily shown from equations (13) and (14). We only need to obtain the result in (15). To do this, one can see from the hrf in (14) that ½ ³ x ´−α ¸ ¾b−1 bα ³ x ´αb−1 ³ x ´ αb−1 · bα b−1 lim h(x) = 1− = lim lim x→∞ θ x→∞ θ θ θ x→∞ θ This implies the result in (16). 2.1. Shapes of the density and hazard rate functions We provide two theorems for the shapes of the density and hazard rate functions. Theorem 2: The WP distribution is unimodal and has a unique mode at x = x0 . When b ≤ 1, the mode is x0 = θ and when b > 1, the mode at x0 is the solution of k(x) = 0, where ¾ ³ x ´α ½ h³ x ´α ib+1 k(x) = αb − αb −1 − 1 − α + 1. (17) θ θ Proof: The derivative with respect to x of equation (13) can be reduced to h³ x ´α ib−2 b α d f (x) = α b θα xα−2 −1 e−[(x/θ) −1] k(x). dx θ (18) From (18), the critical point of f (x) are x = θ and x = x0 , where k(x0 ) = 0. Case I: if b ≤ 1, it is easy to show from (13) that f (x) is a decreasing function, and then f (x) has a unique mode at x = θ. Case II: if b > 1, Theorem 1 indicates that limx→θ+ f (x) = 0, which implies that x = θ can not be a mode point. Therefore, the modes of f (x) are the solutions of k(x) = 0. Finally, we need to show that k(x) = 0 has one solution. One can easily see that k ′ (x) < 0 whenever b > 1. Hence, k(x) = 0 has at most one solution. Using the facts from Theorem 1 that limx→θ+ f (x) = 0 and limx→∞ f (x) = 0, we conclude that f (x) must have a unique mode. Theorem 3: The hrf of X possesses the following shapes: (i) Constant failure rate (CFR) whenever α = b = 1, (ii) Increasing failure rate (IFR) whenever αb ≥ 1 and b > 1, (iii) Decreasing failure rate (DFR) whenever αb ≤ 1 and b < 1, (iv) Upsidedown bathtub rate (UBT) whenever αb ≤ 1, b > 1 and α < 1, (v) Bathtub rate (BT) whenever αb ≥ 1, b < 1 and α > 1. Proof: Based on equation (14), the hrf of X can be expressed as ½³ ´ αb−1 · ³ x ´−α ¸¾b−1 x b−1 −1 1− h(x) = α b θ , θ θ x > θ, which implies (i), (ii) and (iii). Now, the derivative of h(x) in (15) is given by ³ x ´α−1 h³ x ´α ib−2 h′ (x) = α b θ−2 w(x), −1 θ θ 6 where w(x) = (αb − 1) ¡ x ¢α θ − (α − 1). ³ ´1/α α−1 The critical value of h(x) is at x0 = θ αb−1 , which is defined only when α(b − 1) > 0. Now, if αb ≤ 1, b > 1 and α < 1, then x0 is defined and w(x) > 0 for all x > x0 and w(x) < 0 for all θ < x < x0 . This proves (iv). Similar argument can be used to prove (v). 3. Mixture Representation In order to obtain a simplified form for the WP pdf, we expand (13) in power series. Then, the WP pdf can be expressed as G(x)b−1 − f (x; b, α, θ) = b g(x) e G(x)b+1 h i G(x) b G(x) . (19) Inserting (1) and (2) in equation (19), we obtain f (x; b, α, θ) = b ³ α ´ ³ x ´−α h x θ h ¡ x ¢−α ib−1 1− θ n ¡ x ¢−α oib+1 × e| 1− 1− θ  − ª b © −α 1−( x ) θ ª © −α 1− 1−( x ) θ {z }. A Consider the quantity A in the last equation. After a power series expansion, it reduces to   h ¡ x ¢−α ibk   ∞   1− θ X (−1)k A= h n ¡ x ¢−α oibk  . k!    k=0 1− 1− θ Combining the last two results, we have ½ · ∞ ³ α ´ ³ x ´−α X ³ x ´−α ¾¸−[b(k+1)+1] ³ x ´−α ¸b(k+1)−1 · (−1)k f (x; b, α, θ) = b . 1− 1− 1− x θ k! θ θ k=0 | {z } B Consider the quantity B in the last equation. After a power series expansion, we obtain · ∞ ³ x ´−α ¸j X Γ(b(k + 1) + j + 1) B= . 1− Γ(b(k + 1) + 1) j! θ j=0 Combining the last two results, we can write f (x; b, α, θ) = ∞ X b Γ(b(k + 1) + j + 1) (−1)k k! j! (bk + b + j)Γ(b(k + 1) + 1) k,j=0 | {z } vk,j × [b(k + 1) + j] | ³ x ´−α ¸bk+b+j−1 ³ α ´ ³ x ´−α · 1− x θ θ {z } gα,θ,b(k+1)+j (x) 7 In a more simplified form, the last equation reduces to f (x; b, α, θ) = ∞ X vk,j gα,θ,b(k+1)+j (x). (20) k,j=0 Equation (20) reveals that the WP density function can be expressed as a mixture of EP densities. So, several of its mathematical properties can be derived from those of the EP distribution. Equation (20) is the main result of this section. 4. Some Structural Properties 4.1Moments The rth moment of X follows from (20) as µ′r r = E(X ) = ∞ X k,j=0 vk,j Z θ ∞ xr gα,θ,b(k+1)+j (x) dx and then (for r ≤ α) µ′r =θ r ∞ X k,j=0 ¡ ¢ vk,j [b(k + 1) + j] B 1 − αr , b[k + 1] + j , (21) R1 where B(a, b) = 0 wa−1 (1 − w)b−1 dw is the beta function. In particular, setting r = 1 in (21), the mean of X reduces to µ′1 = θ ∞ X k,j=0 ¡ ¢ vk,j [b(k + 1) + j] B 1 − α1 , b[k + 1] + j . Further, the central moments (µn ) and cumulants (κn ) of X are easily obtained from (21) as µn = n µ ¶ X n k=0 k (−1) k µ′k 1 µ′n−k and κn = µ′n − n−1 X k=1 µ ¶ n−1 κk µ′n−k , k−1 ′ ′ ′ ′3 respectively, where κ1 = µ′1 . Thus, κ2 = µ′2 − µ′2 1 , κ3 = µ3 − 3µ2 µ1 + 2µ1 , etc. The skewness and kurtosis measures can be calculated from the ordinary moments using well-known relationships. The nth descending factorial moment of X (for n = 1, 2, . . .) is µ′(n) = E[X (n) ] = E[X(X − 1) × · · · × (X − n + 1)] = n X s(n, j) µ′j , j=0 where s(n, j) = (j!)−1 [dj j (n) /dxj ]x=0 is the Stirling number of the first kind. 4.2 Incomplete Moments 8 Let Y be a random variable having the EP distribution with parameters α, θ and λ. Then, the rth incomplete moment of Y (for r < α) is given by Z z ¡ ¢ y r gα,θ,λ (y) dy = λ θr Bz 1 − αr , λ , mr,Y (z) = θ where Bx (a, b) = Rx 0 wa−1 (1 − w)b−1 dw is the incomplete beta function. The rth incomplete moment of X follows from (20) (for r < α) as Z z ∞ X ¡ ¢ mr,X (z) = xr gα,θ,b(k+1)+j (x) dx = θr vk,j [b(k + 1) + j] Bz 1 − αr , b[k + 1] + j . θ (22) k,j=0 The main application of the first incomplete moment is related to the Bonferroni and Lorenz curves. These curves are very useful in several fields such economics, reliability, demography, insurance and medicine. For a given probability π, they are defined by BX (π) = m1,X (q)/(π µ′1 ) and LX (π) = m1,X (q)/µ′1 , respectively, where m1,X (q) comes from (22) with r = 1 and q = QX (π) is determined from (24) (given in Section 4.4). Another important application of the first incomplete moment refers to the R∞ R ∞mean deviations about the mean and the median defined by δ1 = 0 |x − µ| f (x)dx and δ2 = 0 |x − M | f (x)dx, © ª1/α respectively, where µ′1 = E(X) denotes the mean and M = θ 1 + [− log(0.5)]1/b denotes the median, respectively. These measures can be determined from the expressions δ1 = 2µ′1 F (µ′1 ) − 2m1,X (µ′1 ) and δ2 = µ′1 − 2m1,X (M ), where F (µ′1 ) is given by (12) and m1,X (·) is calculated from (22). A third application of the first incomplete moment is to obtain the mean residual life and the mean waiting time. The mean residual life and the mean waiting time are defined by m(t; b, α, θ) = 1/S(t) [1 − m1,X (t)] − t and µ(t; b, α, θ) = t − [m1,X (t)/F (t; b, α, θ)], where S(t) = 1 − F (t) comes from (12). 4.3 Generating Function First, we obtain the moment generating function (mgf) of the EP distribution, say Mα,θ,λ (t). By expanding the binomial term in (4), we can write Mα,θ,λ (t) = α λ ∞ X m (m+1)α (−1) θ m=0 µ λ−1 m ¶ Z ∞ x−(m+1)α−1 et x dx. θ last integral can be computed using Maple. For t ≤ 0 and p > 0 and q > 0, let J(q, p, t) = R ∞ The −p et x dx. We can obtain using this software q x ¸ · etq p Γ(−p, −t q) π csc(p π) p Γ(−p) − + + , J(q, p, t) = (−t) q − t q Γ(p) qt (−t)p+1 q p+1 tq R∞ where Γ(λ, x) = x wλ−1 e−w dw is the complementary incomplete gamma function. Thus, p Mα,θ,λ (t) = α λ ∞ X (−1)m θ(m+1)α m=0 9 µ ¶ λ−1 J(θ, [m + 1]α + 1, t). m Combining (20) and the last result, it follows the mgf of X as M (t) = ∞ X pm J(θ, [m + 1]α + 1, t), (23) m=0 where pm = (−1)m α θ(m+1)α of this section. P∞ k,j=0 [b(k +1)+j] vk,j ¡b(k+1)+j−1¢ m . Equation (23) is the main result 4.4Quantile Function and Simulation By inverting (12), the quantile function (qf) of X follows as QX (u) = F −1 (u) = θ n o1/α 1 + [− log(1 − u)]1/b . (24) Simulating WP random variable is straightforward. Let U be the uniform variate on the unit interval (0, 1]. Thus, using the inverse transformation method, the random variable X=θ n o1/α 1 + [− log(1 − U )]1/b , has pdf given by (13), i.e., X ∼WP(b, α, θ). 5. Rényi and q-Entropies The entropy of a random variable is a measure of the uncertainty variation. The Rényi entropy of X is defined as 1 IR (δ) = log [I(δ)] 1−δ R where I(δ) = ℜ f δ (x) dx, δ > 0 and δ 6= 1. Now, consider h i δ(b−1) G(x) b −δ δ δ δ G(x) G(x) f (x; b, α, θ)] = b g(x) e . G(x)δ(b+1) Using a power series for the exponential function and the generalized binomial expansion in the above result, we obtain f δ (x; b, α, θ) = ∞ X sk,j g δ (x) G(x)kb+δ(b−1)+j , k,j=0 where (−1)k δ k sk,j = bδ k! j! Γ(kb + δ[b + 1] + j) . Γ(kb + δ[b + 1]) Next, inserting (2) and (1) in (25), and then integrating I(δ) = ∞ X k,j=0 δ sk,j (αθ) Z θ ∞ ·³ ´ ¸ δ(α+1) · ³ x ´−α ¸kb+δ(b−1)+j α x −α dx 1− θ θ 10 (25) For δ(α + 1) > 1 and δ(b − 1) > −1, transforming variables and integrating, the last equation becomes ∞ ³ δ(α + 1) − 1 ´ X sk,j (αθ)δ (θ/α) B I(δ) = , kb + δ(b − 1) + j + 1 . (26) α k,j=0 Hence, the Rényi entropy reduces to   ∞ ´ ³ δ(α + 1) − 1 X 1 log  , kb + δ(b − 1) + j + 1  . sk,j (αθ)δ (θ/α) B IR (δ) = 1−δ α k,j=0 The q-entropy, say Hq (f ), is defined by Hq (f ) = 1 log [1 − Iq (f )] , q−1 R where Iq (f ) = ℜ f q (x) dx, q > 0 and q 6= 1. From equation (26), we can easily obtain (for δ(α + 1) > 1 and q(b − 1) > −1)   ∞ ´ ³ q(α + 1) − 1 X 1 log 1 − , kb + q(b − 1) + j + 1  , s⋆k,j (αθ)q (θ/α) B Hq (f ) = q−1 α k,j=0 where s⋆k,j = (−1)k δ k bδ Γ(kb + q[b + 1] + j)/[k! j!Γ(kb + q[b + 1])]. 6. Order Statistics Here, we give the density of the ith order statistic Xi:n , fi:n (x) say, in a random sample of size n from the WP distribution. It is well known that (for i = 1, . . . , n) fi:n (x) = 1 f (x) F i−1 (x) {1 − F (x)}n−i . B(i, n − i + 1) Using the binomial expansion, we can rewrite fi:n (x) as n−i X ¡ ¢ f (x) F (x)i+j−1 . (−1)j n−i fi:n (x) = j B(i, n − i + 1) (27) j=0 Further, we have i+j−1 F (x) = ∞ X k=0 ¡ (−1)k i+j−1 k ¢ · n³ x ´α ob ¸ exp −(k + 1) −1 θ and then by inserting (13) in (27), we obtain fi:n (x) = ∞ X tk+1 f (x; b, α, θ), k=0 11 (28) where tk+1 µ ¶µ ¶ n−i X 1 i+j−1 j+k n − i = (−1) j (k + 1) B(i, n − i + 1) k j=0 and f (x; b, α, θ) = ½ ib−1 ib ¾ h³ x ´α b α (k + 1) α−1 h³ x ´α − 1 − 1 exp −(k + 1) x θα θ θ is the WP density function with parameters b, α and θ. So, the density function of the WP order statistics is a mixture of WP densities. Based on equation (28), we can obtain some structural properties of the WP order statistics from those of the WP properties. 7. Estimation and Information Matrix Here, we discuss maximum likelihood estimation and inference for the WP distribution. Let x1 , . . . , xn be a sample from X ∼WP(b, α, θ) and let Θ = (b, α, θ)⊤ be the vector of the model parameters with known θ (since x > θ). The likelihood function for Θ reduces to n X ℓ = ℓ(Θ) = n log(bα) − n α log(θ) + (α − 1) i=1 log(xi ) + (b − 1) n X i=1 log(zi ) − n X (zib ). (29) i=1 Since θ is assumed known, the score vector can be denoted by U (Θ) = (∂ℓ/∂b, ∂ℓ/∂α)⊤ , where the components corresponding to the model parameters are determined by differentiating (29). By setting zi = [(xi /θ)α − 1], we obtain ∂ℓ ∂b = ∂ℓ ∂α = n n X n X zib log zi log zi − + b i=1 i=1 ¸ n n · n b − 1 X (1 + zi ) log(1 + zi ) b X b−1 n 1X log(1 + zi ) + zi (1 + zi ) log (1 + zi ). − + α α α zi α i=1 i=1 i=1 The maximum likelihood estimates (MLEs) of the model parameters are the solution of the nonlinear equations U (Θ) = 0, which can be solved iteratively. The elements of the observed information matrix Jn (Θ) = {Jrs } (for r, s = b, α) are given by Jbb = − Jαα n h i X n b 2 z (log z ) , − i i b2 i=1 ¸ n · b − 1 X (1 + zi ) [log(1 + zi )]2 n = − 2− α α2 zi2 i=1 − n oi n bXh (1 + zi ) [log(1 + zi )]2 zib−1 + (b − 1) zib−2 (1 + zi ) , α i=1 12 Jb α = ¸ n · n i b X h b−1 1 X (1 + zi ) log(1 + zi ) zi (1 + zi ) (log zi ) {log (1 + zi )} − α zi α i=1 i=1 n i 1 X h b−1 zi (1 + zi ) log (1 + zi ) . α i=1 √ b The distribution of n(Θ − Θ) can be approximated by the multivariate normal N2 (0, J(Θ)−1 ) distribution, where J(Θ) = limn→∞ n−1 Jn (Θ) is the unit information matrix. The estimated bn (Θ)−1 ) distribution of Θ b can be used to construct apasymptotic multivariate normal N2 (0, J proximate confidence intervalsq for the model parameters. The 100(1 − γ)% confidence intervals for p b and α are given by b̂ ± zγ/2 × var(b̂) and α̂ ± zγ/2 × var(α̂), respectively, where the var(·)’s are b −1 corresponding to the model parameters, and zγ/2 is the quantile the diagonal elements of Jn (Θ) − (1 − γ/2) of the standard normal distribution. 7.1. Simulation study We evaluate the performance of the maximum likelihood method for estimating the WP parameters using Monte Carlo simulation for a total of sixteen parameter combinations and the process is repeated 200 times. Two different sample sizes n = 100 and 300 are considered. The MLEs and the standard deviations of the parameter estimates are listed in Table 1. In this simulation study, we assume θ to be unknown and estimate it by the minimum order statistic x(1) . The MLEs of α and b are determined by solving the nonlinear equations U (Θ) = 0, where xi 6= x(1) . From Table 1, we note that the ML method performs well for estimating the model parameters. Also, as the sample size increases, the biases and the standard deviations of the MLEs decrease as expected. Table 1: MLEs and standard deviations for various parameter values Sample size Actual values Estimated values n b α θ b̃ α̃ θ̃ b̃ α̃ θ̃ 100 0.5 0.5 1 0.5045 0.5138 1.0004 0.0560 0.0434 0.0011 0.5 0.5 2 0.5077 0.5092 2.0007 0.0624 0.0422 0.0019 0.5 0.8 1 0.5078 0.7915 1.0066 0.0429 0.0658 0.0098 0.5 0.8 2 0.5139 0.7877 2.0145 0.0456 0.0656 0.0179 0.5 1 1 0.5204 0.9612 1.0203 0.0344 0.0824 0.0200 0.5 1 2 0.5180 0.9570 2.0419 0.0367 0.0803 0.0432 0.5 2 1 0.5877 1.6137 1.1886 0.0492 0.1606 0.1006 0.5 2 2 0.5918 1.6003 2.3880 0.0490 0.1670 0.1992 0.8 0.5 1 0.7794 0.5312 1.0003 0.0884 0.0459 0.0008 0.8 0.5 2 0.7931 0.5216 2.0005 0.0910 0.0449 0.0008 0.8 0.8 1 0.8190 0.8074 1.0049 0.0713 0.0683 0.0060 13 Standard deviations Table 1 (Continued) Sample size Actual values n b α θ b̃ α̃ θ̃ b̃ α̃ θ̃ 100 0.8 0.8 2 0.8313 0.7838 2.0097 0.0759 0.0689 0.0146 0.8 1 1 0.8309 0.9603 1.0129 0.0538 0.0785 0.0136 0.8 1 2 0.8334 0.9591 2.0246 0.0582 0.0747 0.0244 0.8 2 1 0.9472 1.6120 1.1159 0.0843 0.1672 0.0623 0.8 2 2 0.9566 1.5758 2.2380 0.0861 0.1647 0.1273 1 0.5 1 0.9761 0.5295 1.0002 0.1115 0.0421 0.0004 1 0.5 2 0.9870 0.5155 2.0003 0.1168 0.0420 0.0006 1 0.8 1 1.0094 0.8031 1.0042 0.0907 0.0702 0.0055 1 0.8 2 1.0258 0.7868 2.0075 0.0895 0.0603 0.0088 1 1 1 1.0344 0.9834 1.0109 0.0796 0.0805 0.0112 1 1 2 1.0453 0.9587 2.0207 0.0770 0.0870 0.0191 1 2 1 1.1725 1.6353 1.0840 0.0932 0.1658 0.0440 1 2 2 1.1862 1.6028 2.1834 0.1088 0.1694 0.0991 2 0.5 1 1.9345 0.5404 1.0001 0.2125 0.0421 0.0002 2 0.5 2 1.9570 0.5293 2.0002 0.2021 0.0414 0.0004 2 0.8 1 2.0227 0.8127 1.0017 0.1624 0.0636 0.0020 2 0.8 2 2.0287 0.7984 2.0034 0.1741 0.0653 0.0045 2 1 1 2.0497 0.9725 1.0050 0.1395 0.0777 0.0048 2 1 2 2.0667 0.9781 2.0095 0.1511 0.0831 0.0099 2 2 1 2.3585 1.6277 1.0443 0.2108 0.1712 0.0234 2 2 2 2.3629 1.6008 2.0909 0.1908 0.1625 0.0466 0.5 0.5 1 0.4932 0.5140 1.0000 0.0350 0.0262 0.0001 0.5 0.5 2 0.4931 0.5133 2.0001 0.0310 0.0274 0.0002 0.5 0.8 1 0.4991 0.8109 1.002 0.0256 0.0405 0.0023 0.5 0.8 2 0.5042 0.8029 2.0040 0.0263 0.0412 0.0057 0.5 1 1 0.5064 0.9890 1.0059 0.0199 0.0470 0.0059 0.5 1 2 0.5070 0.9806 2.0138 0.0225 0.0462 0.0134 0.5 2 1 0.5490 1.7775 1.1146 0.0269 0.1244 0.0600 0.5 2 2 0.5500 1.7700 2.2181 0.0280 0.1238 0.1172 0.8 0.5 1 0.7794 0.5253 1.0000 0.0510 0.0230 0.0001 0.8 0.5 2 0.7822 0.5173 2.0000 0.0521 0.0263 0.0001 300 Estimated values 14 Standard deviations Table 1 (Continued) Sample size Actual values Estimated values Standard deviations n b α θ b̃ α̃ θ̃ b̃ α̃ θ̃ 300 0.8 0.8 1 0.8028 0.8027 1.0013 0.0412 0.0385 0.0017 0.8 0.8 2 0.8043 0.8038 2.0025 0.0451 0.0431 0.0029 0.8 1 1 0.8087 0.9935 1.0046 0.0315 0.0485 0.0048 0.8 1 2 0.8086 0.9878 2.0082 0.0313 0.0479 0.0092 0.8 2 1 0.8742 1.7872 1.0632 0.0427 0.1150 0.0329 0.8 2 2 0.8722 1.7899 2.1240 0.0441 0.1155 0.0735 1 0.5 1 0.9686 0.5261 1.0000 0.0587 0.0261 0.0001 1 0.5 2 0.9772 0.5172 2.0001 0.0594 0.0245 0.0001 1 0.8 1 0.9974 0.8094 1.0009 0.0469 0.0410 0.0011 1 0.8 2 1.0090 0.7961 2.0022 0.0485 0.0401 0.0028 1 1 1 1.0092 0.9953 1.0036 0.0414 0.0483 0.0035 1 1 2 1.0108 0.9880 2.0063 0.0397 0.0464 0.0061 1 2 1 1.0837 1.8013 1.0485 0.0472 0.1220 0.0235 1 2 2 1.0952 1.7731 2.1040 0.0536 0.1187 0.0529 2 0.5 1 1.9069 0.5385 1.0000 0.1249 0.0245 0.0001 2 0.5 2 1.9389 0.5228 2.0000 0.1159 0.0252 0.0001 2 0.8 1 1.9919 0.8116 1.0004 0.0929 0.0386 0.0006 2 0.8 2 1.9891 0.8078 2.0008 0.0951 0.0406 0.0011 2 1 1 2.0136 0.9989 1.0019 0.0834 0.0471 0.0020 2 1 2 2.0132 0.9897 2.0034 0.0815 0.0442 0.0040 2 2 1 2.1655 1.8139 1.0245 0.1056 0.1251 0.0130 2 2 2 2.1779 1.7910 2.0497 0.1049 0.1095 0.0276 8. Applications In this section, we illustrate the usefulness of the WP distribution. We fit this distribution to two data sets and compare the results with the KwP, BEP, BP, EP and Pareto distributions. 8.1 Wheaton River Data The data consist to the 72 exceedances for the years 1958–1984 (rounded to one decimal place) of flood peaks (in m3 /s) of the Wheaton River near Carcross in Yukon Territory, Canada. The data are listed in Table 1. Choulakian and Stephens (2001), Akinsete et al. (2008), Mahmoudi (2011) and Bourguignon et al. (2013) analyzed the these data using the BP, BGP and KwP distributions. 15 Table 2: Exceedances of Wheaton River flood data 1.7 1.4 0.6 9.0 5.6 1.5 2.2 18.7 2.2 1.7 30.8 2.5 14.4 8.5 39.0 7.0 13.3 27.4 1.1 25.5 0.3 20.1 4.2 1.0 0.4 11.6 15.0 0.4 25.5 27.1 20.6 14.1 11.0 2.8 3.4 20.2 5.3 22.1 7.3 14.1 11.9 16.8 0.7 1.1 22.9 9.9 21.5 5.3 1.9 2.5 1.7 10.4 27.6 9.7 13.0 14.4 0.1 10.7 36.4 27.5 12.0 1.7 1.1 30.0 2.7 2.5 9.3 37.6 0.6 3.6 64.0 27.0 8.2 Bladder Cancer Data We consider an uncensored data set corresponding to the remission times (in months) of a random sample of 128 bladder cancer patients. Bladder cancer is a disease in which abnormal cells multiply without control in the bladder. The most common type of bladder cancer recapitulates the normal histology of the urothelium and is known as transitional cell carcinoma. These data were previously studied by Lemonte (2012), Zea et al. (2012), Lee and Wang (2003) and Lemonte and Cordeiro (2013). Table 3 lists the remission times of the bladder cancer. Table 3: The remission times of the bladder cancer 0.08 1.76 2.75 3.82 5.09 6.76 7.87 10.66 14.24 22.69 0.20 2.02 2.83 3.88 5.17 6.93 7.93 10.75 14.76 23.63 0.40 2.02 2.87 4.18 5.32 6.94 8.26 11.25 14.77 25.74 0.50 2.07 3.02 4.23 5.32 6.97 8.37 11.64 14.83 25.82 0.51 2.09 3.25 4.26 5.34 7.09 8.53 11.79 15.96 26.31 0.81 2.23 3.31 4.33 5.41 7.26 8.65 11.98 16.62 32.15 0.90 2.26 3.36 4.34 5.41 7.28 8.66 12.02 17.12 34.26 1.05 2.46 3.36 4.40 5.49 7.32 9.02 12.03 17.14 36.66 1.19 2.54 3.48 4.50 5.62 7.39 9.22 12.07 17.36 43.01 1.26 2.62 3.52 4.51 5.71 7.59 9.47 12.63 18.10 46.12 1.35 2.64 3.57 4.87 5.85 7.62 9.74 13.11 19.13 79.05 1.40 2.69 3.64 4.98 6.25 7.63 10.06 13.29 20.28 1.46 2.69 3.70 5.06 6.54 7.66 10.34 13.80 21.73 We estimate the model parameters of the distributions by the method of maximum likelihood. There exists many maximization methods in R Packages like NR (Newton-Raphson), BFGS (Broyden-Fletcher-Goldfarb-Shanno), BHHH (Berndt-Hall-Hall-Hausman), SANN (Simulated-Annealing), NM (Nelder-Mead) and L-BFGS-B. Here, the MLEs are computed using LimitedMemory quasi-Newton code for Bound-constrained optimization (L-BFGS-B) and the measures of goodness of fit including the Akaike information criterion (AIC), consistent Akaike information criterion (CAIC), Bayesian information criterion (BIC), Hannan-Quinn information criterion (HQIC), Anderson-Darling (A∗ ) and Cramér–von Mises (W ∗ ) are computed to compare the fitted models. The statistics W ∗ and A∗ are described in details in Chen and Balakrishnan (1995). In general, the smaller the values of these statistics, the better the fit to the data. The required computations are carried out using a script AdequacyModel of the R-package written by Pedro Rafael Diniz Marinho, Cícero Rafael Barros Dias and Marcelo Bourguignon. It is freely available from http://cran.r-project.org/web/packages/AdequacyModel/AdequacyModel.pdf. 16 Tables 4 and 6 list the MLEs and their corresponding standard errors (in parentheses) of the model parameters. The model selection is carried out using the following statistics: 2pn , n−p−1 © £ ¤ª and HQIC = 2 log log(n) k − 2 ℓ̂ , AIC = −2 ℓ̂ + 2p, BIC = −2 ℓ̂ + p log(n), CAIC = −2 ℓ̂ + where ℓ̂ denotes the log-likelihood function evaluated at the MLEs, p is the number of parameters, and n is the sample size. The statistics AIC, CAIC, BIC, HQIC, W ∗ and A∗ are listed in Tables 5 and 7. Table 4: MLEs and their standard errors (in parentheses) for Wheaton river flood data Distribution WP BEP KwP BP EP Pareto a b α λ θ 0.1742 (0.0632) 11.2069 (3.0323) 22.2915 (4.4687) - 4.4363 (0.4409) 83.7393 (251.1542) 31.9168 (31.2434) 12.7688 (12.1744) - 0.0987 (0.0020) 0.3562 (0.0926) 0.1910 (0.06356) 0.1602 (0.1085) 0.7197 (0.0628) 0.1548 (0.0184) 83.8231 (26.3243) 62.4789 (21.4784) - 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 - Table 5: The statistics ℓ̂, AIC, CAIC , BIC , HQIC, A∗ and W ∗ for Wheaton river flood data Distribution WP BEP KwP BP EP Pareto ℓ̂ AIC CAIC BIC HQIC A∗ W∗ -249.3965 -248.0554 -252.4078 -255.8893 -260.9832 -335.2348 502.7930 504.1107 510.8157 517.7787 525.9665 672.4695 502.9695 504.7168 511.1739 518.1369 526.1429 672.5275 507.3184 513.1615 517.6037 524.5667 530.4918 674.7322 504.5926 507.7099 513.515 520.4781 527.7661 673.3693 0.9819 1.0320 1.4604 1.9532 2.5498 1.9146 0.1746 0.1903 0.2619 0.3445 0.4410 0.3381 17 Table 6: MLEs and their standard errors (in parentheses) for bladder cancer data Distribution WP KwP BEP BP EP Pareto a b α λ θ 7.1911 (1.2078) 0.5457 (0.3509) 14.9920 (1.9701) - 3.5707 (0.2477) 66.0736 (67.9676) 30.3072 (16.9047) 23.6299 (44.0328) - 0.1466 (0.0027) 0.1756 (0.0616) 0.2952 (0.1112) 0.1158 (0.1770) 0.8513 (0.0546) 0.2319 (0.0206) 14.3993 (10.6635) 23.9408 (4.7808) - 0.08 0.08 0.08 0.08 0.08 0.08 - Table 7: The statistics ℓ̂, AIC, CAIC, BIC , HQIC, A∗ and W ∗ for bladder cancer data Distribution WP KwP BEP BP EP Pareto ℓ̂ AIC CAIC BIC HQIC A∗ W∗ -407.4665 -409.2100 -409.0659 -418.9310 -432.0897 -539.5911 818.9331 824.4200 826.1318 843.8620 868.1794 1081.1820 819.0298 824.6151 826.4596 844.0571 868.2762 1081.2140 824.6214 832.9525 837.5085 852.3946 873.8678 1084.0260 821.2442 827.8866 830.7540 847.3287 870.4905 1082.3380 0.3650 0.5351 0.5061 1.8946 3.6418 1.8785 0.0563 0.0770 0.0725 0.2930 0.5865 0.2904 In Tables 5 and 7, the WP model is compared with the KwP, BEP, BP and EP models. We note that the WP model has the lowest values of the AIC, BIC, CAIC, HQIC, A∗ and W ∗ statistics among all fitted models. So, the WP distribution could be chosen as the best model for both data sets. The histogram of the data and plots of the estimated pdf and cdf of the WP model are displayed in Figure 2. It is clear that the new distribution provides a better fit to the histogram and therefore could be chosen as the best model for both data sets. 9. Concluding Remarks In this paper, we propose a three-parameter Weibull-Pareto (WP) distribution based on the Weibull-G family of distributions recently introduced by Bourguignon et al. (2014). We derive some of its structural properties and provide explicit expressions for the ordinary and incomplete moments, mean residual life, mean waiting time, generating function, mode and quantile function. We also obtain expressions for the Rényi entropy, q entropy and the density of ith order statistic. The model parameters are estimated by the method of maximum likelihood. The usefulness of the new model is illustrated by means of two real life data sets. The new model provides consistently a better fit than other competitive lifetime models. We hope that the new model will attract wider applications in several areas such as engineering, survival data, economics (income inequality) and others. 18 (b) Estimated cdfs (data set 1) 0.6 0.4 0.03 cdf 0.04 0.05 WP BEP BP KwP EP 0.02 pdf 0.8 0.06 1.0 0.07 (a) Estimated pdfs (data set 1) 0.0 0.00 0.01 0.2 WP BEP BP KwP EP 0 20 40 60 80 0 20 40 x 80 (d) Estimated cdfs (data set 2) 0.04 0.8 0.05 1.0 (c) Estimated pdfs (data set 2) 0.6 cdf 0.4 0.03 WP BEP BP KwP EP WP BEP BP KwP EP 0.0 0.00 0.2 0.01 0.02 pdf 60 x 0 10 20 30 40 50 60 70 0 x 20 40 60 x Figure 2: Plots of the fitted WP, BEP, BP, KwP and EP distributions Acknowledgments The authors would like to thank the Editor-in-Chief, and the referee for constructive comments which greatly improved the paper. References Akinsete, A., Famoye, F., Lee, C. (2008). The beta-Pareto distribution. Statistics 42: 547–563. Alzaatreh, A. (2011). A New Method for Generating Families of Continuous Distributions. Ph.D. thesis, Central Michigan University, Mount Pleasant: Michigan. Alzaatreh, A., Famoye, F., Lee, C. (2012). Gamma-Pareto distribution and its applications. Journal of Modern Applied Statistical Methods 11:78–94. Alzaatreh, A., Famoye, F., Lee, C. (2013a). Weibull-Pareto distribution and its applications. Communications in Statistics–Theory and Methods 42:1673–1691. 19 Alzaatreh, A., Famoye, F., Lee, C. (2013b). A new method for generating families of continuous distributions. Metron 71:63–79. Bourguignon, M., Silva, M.B., Zea, L.M., Cordeiro, G.M. (2013). The Kumaraswamy Pareto distribution. Journal of Statistical Theory and Applications 12:129–144. Bourguignon, M., Silva, R.B., Cordeiro, G.M. (2014). The Weibull–G family of probability distributions. Journal of Data Science 12:53–68. Chen, G., Balakrishnan, N. (1995). A general purpose approximate goodness-of-fit test. Journal of Quality Technology 27:154–161. Choulakian, V., Stephens, M.A. (2001). Goodness-of-fit for the generalized Pareto distribution. Technometrics 43:478–484. Cordeiro, G.M., de Castro, M. (2011). A new family of generalized distributions. Journal of Statistical Computation and Simulation 81:883–893. Eugene, N., Lee, C., Famoye, F. (2002). Beta-normal distribution and its applications. Communications in Statistics–Theory and Methods 31:497–512. Gupta, R.C., Gupta, P.I., Gupta, R.D. (1998). Modeling failure time data by Lehmann alternatives. Communications in Statistics–Theory and Methods 27:887–904. Lee, E.T., Wang, J.W. (2003). Statistical Methods for Survival Data Analysis, 3rd edition, Wiley: New York. Lehmann, E.E. (1953). The power of rank tests. Annals of Mathematical Statistics 24:23–43. Lemonte, A.J. (2012). The beta log-logistic distribution. To appear in Brazilian Journal of Probability and Statistics. Lemonte, A.J., Cordeiro, G.M. (2013). An extended Lomax distribution. Statistics 47:800–816. Mahmoudi, E. (2011). The beta generalized Pareto distribution with application to lifetime data. Mathematics and Computers in Simulation 81:2414–2430. Mansoor, M. (2013). Properties and Extensions of Stoppa Income Size Distribution. Unpublished M.Phil. thesis, Department of Statistics, The Islamia University of Bahawalpur: Pakistan. Nassar, M.M., Nada, N.K. (2011). The beta generalized Pareto distribution. Journal of Statistics: Advances in Theory and Applications 6:1–17. Pareto, V. (1896). La courbe de la rpartition de la richesse. In Recueil publi par la Facult de Droit l.occasion de l.exposition nationale suisse (Eds) C. ViretGenton, Lausanne: Universit de Lausanne, pp. 373–387. Pickands, J. (1975). Statistical inference using extreme order statistics. Annals of Statistics 3:119– 131. 20 Stoppa, G. (1990). A new model for income size distribution. In: C. Dagum and M. Zenga (eds.) Income and Wealth Distribution, Inequality and Poverty, Springer, New York, pp. 33–41. Zea, L.M., Silva, R.B., Bourguignon, M., Santos, A.M., Cordeiro, G.M. (2012). The beta exponentiated Pareto distribution with application to bladder cancer susceptibility. International Journal of Statistics and Probability 1:8–19. Zografos, K., Balakrishnan, N. (2009). On families of beta- and generalized gamma-generated distributions and associated inference. Statistical Methodology 6:344–362. 21