The Kumaraswamy Pareto IV Distribution

We introduce a new model named the Kumaraswamy Pareto IV distribution which extends the Pareto and Pareto IV distributions. The density function is very flexible and can be left-skewed, right-skewed and symmetrical shapes. It has increasing, decreasing, upside-down bathtub, bathtub, J and reversed-J shaped hazard rate shapes. Various structural properties are derived including explicit expressions for the quantile function, ordinary and incomplete moments, Bonferroni and Lorenz curves, mean deviations, mean residual life, mean waiting time, probability weighted moments and generating function. We provide the density function of the order statistics and their moments. The Rényi and q entropies are also obtained. The model parameters are estimated by the method of maximum likelihood and the observed information matrix is determined. The usefulness of the new model is illustrated by means of three real-life data sets. In fact, our proposed model provides a better fit to these data than the gamma-Pareto IV, gamma-Pareto, beta-Pareto, exponentiated Pareto and Pareto IV models.


Introduction
The Pareto distribution and its generalizations are tractable statistical models for scientists, economists and engineers. These distributions cover a wide range of applications especially in economics, finance, actuarial science, risk theory, reliability, telecommunications and medicine. Arnold (1983) established a hierarchy of the Pareto distribution by starting with the classical Pareto (I) to Pareto (IV) distributions by subsequently introducing additional parameters related to location, scale, shape and inequality (Gini index). A very generalized form of the Pareto model is the four-parameter Pareto IV (PIV) distribution defined by the cumulative distribution function (cdf) H µ,α,γ,θ (x) = 1 − 1 + x − µ θ where α, γ, θ > 0, −∞ < µ < ∞ is the location parameter, θ is the scale parameter, γ is the inequality parameter and α is the shape parameter which characterizes the tails of the distribution.
The survival function S(t) and the hazard rate function (hrf) h(t) of Z are given by respectively.
Among these generators, the Kw-G family has received increased attention after the convincing debate on the pitfalls of the beta-G family by Jones (2009).
For a baseline random variable having pdf g(x) and cdf G(x), Cordeiro and de Castro (2011) defined the two-parameter Kw-G cdf by The pdf corresponding to (5) becomes where g(x) = dG(x)/dx and a > 0 and b > 0 are two extra shape parameters whose role are to govern skewness and tail weights.
In this context, we propose an extension of the Pareto-IV model called the Kumaraswamy Pareto-IV ("KwPIV" for short) distribution based on equations (5) and (6). The paper is outlined as follows. In Section 2, we define the KwPIV distribution. The shapes of the density and hazard rate functions are discussed in Section 3. We provide a mixture representation for its density function in Section 4. Mathematical properties such as the quantile function (qf), ordinary and incomplete moments, Bonferroni and Lorenz curves, mean deviations, mean residual life, mean waiting time, probability weighted moments (PWMs) and generating function are derived in Section 5. The density of the order statistics is determined in Section 6. In Section 7, we obtain the Rényi and q entropies. The maximum likelihood estimation of the model parameters is discussed in Section 8. We explore its usefulness by means of three real-life data sets in Section 9. Finally, Section 10 offers some concluding remarks.

The Kumaraswamy Pareto IV distribution
Inserting (2) in (5), the five-parameter KwPIV cdf is defined by where a, b, α > 0 are shape parameters, θ is the scale parameter and γ is the inequality parameter.
The pdf corresponding to (7) becomes Henceforth, we denote by X ∼KwPIV(a, b, α, γ, θ) a random variable having pdf (8). We omit the dependence of the pdf and cdf of X on the parameters.
The Kw-PIV distribution contains as special models the following known and unknown distributions: (i) For a = b = 1, equation (8) gives the PIV density (3); (ii) For a = 1, equation (8) yields the Lehmann type II PIV density, not known in literature yet; (iii) For b = 1, equation (8) gives the exponentiated PIV density, not known in literature yet; (iv) For b = 1, and a is an integer, equation (8) reduces to the distribution of the maximum of a random sample of size "a" from the PIV distribution.
An interpretation of the KwPIV cdf defined by (7) (for a and b positive integers) can be given as follows. Consider that a system is formed by b series independent components and that each component is made up of a parallel independent subcomponents. For j = 1, . . . , b, let X j denote the lifetime of the jth component. Suppose that X j1 , . . . , X ja denote the lifetimes of the subcomponents within the jth component (for j = 1, . . . , b) having a common PIV cdf given by (2). Let X denote the lifetime of the entire system. Then, the cdf of X is given by So, the KwPIV distribution given by (7) represent precisely the time to failure distribution of the entire system.
There may be more than one root to (10). If x = x 0 is a root of (10), then it corresponds to a local maximum, local minimum or a point of inflexion depending on whether The first derivative of log {h(x)} for X takes the form Then, the roots of the following equation are the modes of h(x) There may be more than one root to (11). If x = x 0 is a root of (11), then it corresponds to a local maximum, local minimum or a point of inflexion depending on whether

Mixture representation
Inserting this expansion in equation (8) and, after some algebra, we obtain The quantity B j in the last equation, after a binomial expansion, becomes Combining the last two results, we can write .
The last equation can be rewritten as where Equation (12) is the main result of this section. It gives the KwPIV pdf as a mixture of PIV densities. So, several of its mathematical properties can be derived form those of the PIV distribution. The coefficients v +,k depend only on the generator parameters.

Some mathematical properties
Established algebraic expansions to determine some mathematical properties of the KwPIV distribution can be more efficient than computing those directly by numerical integration of (8), which can be prone to rounding off errors among others.

Quantile function
The qf of X is obtained by inverting (7) as Simulating the KwPIV random variable is straightforward. If U is a uniform variate on the unit interval (0, 1), then the random variable X = Q(U ) follows (8).

Moments
The rth moment of X can be written from (12) as Using (4), we obtain (for r < α/γ) Setting r = 1 in (14), we have the mean µ 1 of X. The central moments (µ n ) and cumulants (κ n ) of X follow from (14) as The skewness and kurtosis measures can be calculated from the ordinary moments using wellknown relationships.

Incomplete moments
The answers to many important questions in economics require more than just knowing the mean of the distribution, but its shape as well. The rth incomplete moment of X follows from (12) (for r < α/γ) as where B z (a, b) = z 0 w a−1 (1 − w) b−1 is the incomplete beta function. The main application of the first incomplete moment refers to the Bonferroni and Lorenz curves. These curves are very useful in several fields. For a given probability π, they are defined by B(π) = m 1 (q)/(π µ 1 ) and L(π) = m 1 (q)/µ 1 , respectively, where m 1 (q) follows from (15) with r = 1 and q = Q(π) is calculated from (13). The plots of Bonferroni curve with fixed parameters b = 1, α = 2, γ = 1.5 and θ = 2, and Lorenz curve with fixed parameters α = 2, γ = 1.5 and θ = 2 are given in Figure 3.
The amount of scatter in a population is measured to some extent by the totality of deviations from the mean and median defined by is the mean and M = Q(0.5) is the median. These measures can be determined from δ 1 = 2µ 1 F (µ 1 ) − 2m 1 (µ 1 ) and δ 2 = µ 1 − 2m 1 (M ), where F (µ 1 ) comes from (7).
A further application of the first incomplete moment is related to the mean residual life and the mean waiting time given by , respectively, where F (·; ·) and S(·; ·) = 1 − F (·; ·) are obtained from (7).

Probability weighted moments
The (s, r)th PWM of X (for r ≥ 1, s ≥ 0) is formally defined by For some specific distributions, the relations between the PWMs and the parameters are of a simpler analytical structure than those between the conventional moments and the parameters. The simpler analytical structure suggests that it may be possible to derive the relations between the parameters and the PWMs even though it may not be possible to derive the relations between the parameters and the conventional moments. Consider Using the binomial expansion twice in the above equation, we obtain Inserting (12) and (17) in (16), and after some algebra, ρ s,r can be expressed as By using (4), we obtain (for r < α/γ) ρ r,s = α θ r ∞ n,k=0 (n + k + 1) η (s) n,k B (rγ + 1, (n + k + 1)α − rγ) .

Generating function
We obtain the moment generating function (mgf) M X (t) of X from (12) as By expanding the binomial terms, we can write For t < 0, M (t) can be expressed as which is the main result of this section.

Reliability
Let X and Y be two continuous and independent random variables with cdfs F 1 (x) and F 2 (y) and pdfs f 1 (x) and f 2 (y), respectively. The reliability parameter R = P (X > Y ) is defined by Suppose that X ∼KwPIV(a 1 , b 1 , α 1 , γ, θ) and Y ∼KwPIV(a 2 , b 2 , α 2 , γ, θ). We can use the mixture representation (12) to obtain the power series expansion for the cdf of X as Let w +,m = ∞ l=0 (−1) l+m b 2 l la 2 m . By using again (12), the reliability of the KwPIV model (for fixed values of γ and θ) can be expressed as

The Kumaraswamy Pareto IV Distribution
Here v +,k is obtained from (12) in terms of the parameters a 1 and b 1 .

Order statistics
Here, we provide the density of the ith order statistic X i:n , f i:n (x) say, in a random sample of size n from the KwPIV distribution. By suppressing the parameters, we have (for i = 1, . . . , n) Following (17), we can write Then, by inserting (12) in (18), we obtain where where v +,k is given in (12) and g(x; α(p + k + 1), γ, θ) denotes the PIV density function with parameters α(p + k + 1), γ and θ. So, the density function of the KwPIV order statistics is a mixture of PIV densities. Based on equation (19), we can obtain some structural properties of X i:n from those of the PIV properties.
The L-moments are analogous to the ordinary moments but can be estimated by linear combinations of order statistics. They exist whenever the mean of the distribution exists, even though some higher moments may not exist, and are relatively robust to the effects of outliers. Recently, L-moments have been noticed as appealing alternatives to the conventional moments. Based upon the moments (20), we can determine closed-form expressions for the L-moments of X as finite weighted linear combinations of KwPIV expected order statistics given by (for s ≥ 1) The first four L-moments are: λ 1 = E(X 1:1 ), λ 2 = 1 2 E(X 2:2 − X 1:2 ), λ 3 = 1 3 E(X 3:3 − 2X 2:3 + X 1:3 ) and λ 4 = 1 4 E(X 4:4 − 3X 3:4 + 3X 2:4 − X 1:4 ). We can easily obtain the λ's for X from (20) with r = 1.

Rényi and q-entropies
The entropy of a random variable X is a measure of the uncertain variations. The Rényi entropy is defined by where I(δ) = f δ (x) dx, δ > 0 and δ = 1.
Applying the power series used before twice to the last result, we obtain where Inserting (2) and (3) in equation (21) and integrating, we obtain Computing the last integral for δ >min{γ/(γ − 1), γ/(α + γ)}, Hence, the Rényi entropy reduces to The q-entropy, say H q (f ), is defined by where I q (f ) = f q (x) dx, q > 0 and q = 1.
From equation (22), we can easily obtain

Estimation
Here, we consider the estimation of the unknown parameters of the new distribution by the method of maximum likelihood. Let x 1 , . . . , x n be a random sample of size n from the KwPIV distribution given by (8). The log-likelihood function for the vector of parameters Θ = (a, b, α, γ, θ) can be expressed as Then, we can write as The components of the score vector U (Θ) are given by Setting these equations to zero and solving them simultaneously yields the MLEs of the five parameters. For interval estimation of the model parameters, we require the 5 × 5 observed information matrix J(Θ) = {U rs } (for r, s = a, b, α, γ, θ) given in Appendix A. Under standard regularity conditions, the multivariate normal N 5 (0, J( Θ) −1 ) distribution can be used to construct approximate confidence intervals for the model parameters. Here, J( Θ) is the total observed information matrix evaluated at Θ. Then, the 100(1 − γ)% confidence intervals for a, b, α, γ and θ are given byâ ± z α * /2 × var(â),b ± z α * /2 × var(b),α ± z α * /2 × var(α) γ ± z α * /2 × var(γ) andθ ± z α * /2 × var(θ), respectively, where the var(·)'s denote the diagonal elements of J( Θ) −1 corresponding to the model parameters, and z α * /2 is the quantile (1 − α * /2) of the standard normal distribution.
The likelihood ratio (LR) statistic can be used to check if the KwPIV distribution is strictly "superior" to the PIV distribution for a given data set. Then, the test of H 0 : a = b = 1 versus H 1 : H 0 is not true is equivalent to compare the KwPIV and PIV distributions and the LR statistic becomes w = 2{ ( a, b, α, γ, θ) − (1, 1,α,γ,θ)}, where a, b, α, γ and θ are the MLEs under H 1 andα,γ andθ are the estimates under H 0 .

Applications
Here, we illustrate the usefulness of the KwPIV model by means of three applications to real data sets. The MLEs and goodness-of-fit measures for some fitted models are computed to compare them.

Application 1: Actuarial science data
This data set has been studied by Balakrishnan, Leiva, Sanhueza, and Cabrera (2009) to show the flexibility of the mixture inverse Gaussian model, which describes the distributional behavior of the mortality of retired people on disability of the Mexican Institute of Social Security. The data set corresponding to lifetimes (in years) of retired women with temporary disabilities, which are incorporated in the Mexican insurance public system and who died during 2004, are: 22, 24, 25(2), 27, 28, 29(4), 30, 31(6), 32 (7)

Application 3: Sum of skin folds data
The maximum log-likelihood (ˆ ) and measures of goodness of fit including Akaike information criterion (AIC), consistent Akaike information criterion (CAIC), Anderson-Darling (A * ) and Cramér-von Mises (W * ) are computed to compare the fitted models. The statistics A * and W * 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 in the R-language introduced by R (2009).
The numerical values of statistics AIC, CAIC, A * and W * are listed in Tables 1, 2 and 3.  Tables 4, 5 and 6 list the MLEs and their corresponding standard errors (in parentheses) of the model parameters.
In Tables 1, 2 and 3, we compare the KwPIV model with the GPIV, GP, BP, PIV and EP models. We note that the KwPIV model gives the lowest values for the AIC, CAIC, A * and W * statistics among all fitted models. So, the KwPIV model could be chosen as the best model for these three data sets. The histogram of the three data sets and their estimated pdfs and cdfs for the fitted models are displayed in Figures 3 and 4. These plots indicate that the KwPIV distribution gives a better fit to the histograms and therefore could be chosen to explain these data sets.

Conclusions
In this paper, we propose a new five-parameter Kumaraswamy Pareto IV distribution (Kw-PIV) distribution. We study some of its structural properties including an expansion for the density function and explicit expressions for the quantile function, ordinary and incomplete moments, mean residual life, mean waiting time, probability weighted moments and generating functions. Explicit expressions for the order statistics, Rényi entropy and q entropy are also derived. The maximum likelihood method is employed for estimating the model parameters.
We also obtain the observed information matrix. We fit the KwPIV model to three real life data sets to show the usefulness of the proposed distribution. The new model provides consistently better fit than other lifetime models such as the gamma-Pareto IV, gamma-Pareto, beta-Pareto, exponentiated Pareto and Pareto IV models. We hope that the proposed model may attract wider application in areas such as economics, finance, actuarial science, engineering, survival and lifetime data, among others.