A New Probability Distribution for Modeling Failure and Service Times: Properties, Copulas and Various Estimation Methods

In this paper, a new generalization of the Pareto type II model is introduced and studied. The new density can be “right skewed” with heavy tail shape and its corresponding failure rate can be “ J -shape”, “decreasing” and “upside down (or increasing-constant-decreasing)”. The new model may be used as an “under-dispersed” and “over-dispersed” model. Bayesian and non-Bayesian estimation methods are considered. We assessed the performance of all methods via simulation study. Bayesian and non-Bayesian estimation methods are compared in modeling real data via two applications. In modeling real data, the maximum likelihood method is the best estimation method. So, we used it in comparing competitive models. Before using the the maximum likelihood method, we performed simulation experiments to assess the ﬁnite sample behavior of it using the biases and mean squared errors.


Introduction
The Pareto type II (also called as Lomax) distribution is a heavy-tail probability distribution used in business, actuarial science, biological sciences, engineering, economics, income and wealth inequality, queueing theory, size of cities, and Internet traffic modeling (see Lomax [41]). It has been applied to model data obtained from income and wealth (Harris [34] and Atkinson and Harrison [10]), firm size (Corbellini et al. [31]) and reliability and life testing (Hassan and Al-Ghamdi [35]). The Pareto type II model is known as a special model form of Pearson type VI distribution and has also considered as a mixture of exponential (Exp) and gamma (Gam) distributions. The Pareto type II model belongs to the family of "decreasing" hazard rate function (HRF) and considered as a limiting model of residual lifetimes at great age (Balkema and de Hann [19] and Chahkandi and Ganjali [21]). The Pareto type II distribution has been suggested as heavy tailed alternative to the Exp, Weibull (W) and Gam distributions (Bryson [20]). For details about relation between the Pareto type II model and the Burr family and Compound Gamma (CGam) model see Tadikamalla [62] and Durbey [24].
The main aim of this work is to provide a flexible extension of the Pareto type II distribution using the odd Burr-G (OB-G) family defined by Alizadeh et al. [6]. Bayesian and Non-Bayesian estimation methods such as uch as the maximum likelihood estimation (MLE), ordinary least square estimation (OLSE), weighted least square estimation Due to Alizadeh et al. [6], the CDF of the OB-G family is given by where G ζ (z) = 1 − G ζ (z) . The PDF corresponding to (3) is given by For θ = 1, the OB-G family reduces to the Odd G (O-G) family (see Gleaton and Lynch [42]). For v = 1, the OB-G family reduces to the proportional reversed hazard rate (PRHR) family. The odd Burr Pareto type II (OBP) survival function (SF) is given by where ∇ β (z) = 1 + z β and S ξ (z) = 1 − F ξ (z) | (ξ=v,θ,α,β) . For θ = 1, the OBP reduces to the odd Pareto type II. For v = 1, the OBP reduces to the proportional reversed hazard rate Pareto type II. The PDF corresponding to (5) is given by The HRF for the new model can be derived from f ξ (z) /S ξ (z). Many useful Pareto type II extensions can be found in Tahir et al. [64] (Weibull Pareto type II distribution), Cordeiro et al. [23] (the one parameter Pareto type II system of densities), Altun et al. [12] (Odd log-logistic Pareto type II), Altun et al. [12] (Zografos-Balakrishnan Pareto type II distribution), Elbiely and Yousof [33] (Weibull generalized Pareto type II, Rayleigh generalized Pareto type II and Exponential generalized Pareto type II distributions), Yousof et al. [66] (Topp Leone Poisson Pareto type II distribution), Goual and Yousof [43] (Pareto type II inverse Rayleigh), Gad et al. [27] (Burr type XII Pareto type II, Pareto type II Burr type XII and Pareto type II Pareto type II distributions), Yadav et al. [65] (Topp Leone Pareto type II distribution) and Ibrahim and Yousof [36] (Poisson Burr X generalized Pareto type II and Poisson Rayleigh generalized Pareto type II distributions). To illustrate the flexibility of the new PDF and its corresponding HRF we present Figure 1. Figure 1(a) gives some PDF plots for some selected parameters value. Figure 1(b) gives some HRF plots for some selected parameters value. Based on Figure 1(a) the OBP density can be right skewed with heavy tail shape. Based on Figure 1(b) the OBP HRF can be "Jshape" (v = 100, θ = 1, α = 5, β = 0.55), "decreasing" (v = 0.75, θ = 7.5, α = 0.1, β = 0.1) and "upside down (or increasing-constant-decreasing)" (v = 1.5 = 1.5, α = β = 1).

BOBP type via Ali-Mikhail-Haq copula
Under the stronger Lipschitz condition, the joint CDF of the archimedean Ali-Mikhail-Haq copula can expressed as the corresponding joint PDF of the archimedean Ali-Mikhail-Haq copula can expressed as we can derive the joint CDF and the joint PDF of the BOBP type via Ali-Mikhail-Haq copula.

Asymptotics and quantile function
In mathematical analysis, the asymptotic analysis is used for describing the limiting behavior of some functions. Asymptotics derivations for the CDF, PDF and HRF can be obtained for the new model. The asymptotics of the CDF, PDF and HRF as z → 0 are given by The asymptotics of CDF, PDF and HRF as z → ∞ are given by For simulation of this new model, we obtain the quantile function (QF) of Z (by inverting (5)), say z u = F −1 (u), as Equation (7) is used for simulating the new model.

Useful representations
Due to Alizadeh et al. (2017a), the PDF in (6) can be expressed as where where G s • ,α,β (z) is the CDF of the Pareto type II distribution with power parameter s • .

Moments and incomplete moments
The r th ordinary moment of Z is given by then we obtain w2−1 dm. Setting r = 1, 2, 3 and 4 in (10), we have where E(Z) = µ ′ 1 is the mean of Z. The r th incomplete moment, say I r (t), of Z can be expressed, from (9), as The first incomplete moment given by (11) with r = 1 as The index of dispersion IxDis or the variance to mean ratio can derived as IxDis(z) = µ 2 /µ ′ 1 . It is a measure used to quantify whether a set of observed occurrences are clustered or dispersed compared to a standard statistical model. Table 1 we noted that, the Skew(Z) of the OBP distribution can range in the interval ( −69216.87, 98983.16). The spread for the Kur(Z) of the OBP model is ranging from 74.17124 to ∞. The IxDis(Z) for the OBP model can be in (0, 1) an dalso > 1 so it may be used as an "under-dispersed" and "over-dispersed" model. Figure 2 gives some three dimensional skewness plots for parameter α. The flexibility of the skewness of the new model are proved in Figures 2 and 3 using some some three dimensional plots. The flexibility of the kurtosis of the new model are proved in Figures 4 and 5 using some some three dimensional plots. Figure 3 gives some three dimensional skewness plots for parameter θ. Figure 4 gives some three dimensional kurtosis plots for parameter α. Finally, Figure 5 gives three dimensional kurtosis plots for parameter θ.

Some generating functions (GF)
The moment generating function (MGF) can be derived using (8) as  where M s • (t; α, β) is the MGF of the ExpP model, then  The first r derivatives of M Z (t), with respect to t at t = 0, yield the first r moments about the origin, i.e.,  The generating function GF (CGF) is the logarithm of the MGF. Thus, r th cumulant, say κ r , can be obtained from ) ] | (t=0, and r=1,2,3,...) .
The 1 st cumulant is the mean (κ 1,Z = µ ′ 1 ), the 2 nd cumulant is the variance, and the 3 rd cumulant is the same as the 3 rd central moment κ 3,Z = µ 3,Z . But 4 th and higher-order cumulants are not equal to central moments, that being said In some cases theoretical treatments of problems in terms of cumulants are simpler than those using moments. In particular, when two or more RVs are statistically independent, the r th order cumulant of their sum is equal to the sum of their r th order cumulants. Moreover, the cumulants can be also obtained from κ r,

Reversed residual life function
The n th moment of the reversed residual life, say V n, Then, the n th moment of the reversed residual life of Z becomes

Estimation
In this Section, non-Bayesian and Bayesian estimation methods are considered. In first subsection, we will consider four non-Bayesian estimation methods such as the MLE, OLSE, WLSE and the KE methods. In the second subsection, the Bayesian estimation method under the squared error loss function (SELF).

Non-Bayesian estimation methods
. . , z n be a random sample from size n from the OBP distribution with parameters v, θ, α and β. Let ξ be the 4 × 1 parameter vector. For determining the MLE of ξ, we have the log-likelihood function The components of the score vector, , are available if needed. Setting U (v) = U (θ) = U (α) = U (β) = 0 and solving them simultaneously yields the MLE ξ = ( v, θ, α, β) . To solve these equations, it is usually more convenient to use nonlinear optimization methods such as the quasi-Newton algorithm to numerically maximize ℓ. For interval estimation of the parameters, we obtain the 4 × 4 observed information matrix J(ξ) = {∂ 2 ℓ(ξ)/∂m ∂w}symbol| (m,w=v,θ,α,β) .

OLS Let
denotes the CDF of OBP model and let z 1 < z 2 < · · · < z n be the n ordered RS. The OLSEs are obtained upon minimizing then, we have The LSEs are obtained via solving the following non-linear equations WRT v, θ, α and β The WLSEs are obtained by solving , ξ) and ς β (z [i:n] , ξ) defined above.

KE method
The Kolmogorov estimates (KEs) v, θ, α and β and β of v, θ, α and β are obtained by minimizing the function

Bayesian estimation
Assume the gamma priors of the parameters v, θ, α and β of the following forms Assume that the parameters are independently distributed. The joint prior distribution can be written as The posterior distribution π (v, θ, α, β|Z) of the parameters is defined as Under squared error loss function, the Bayesian estimators of v, θ, α and β are the means of their marginal posteriors and defined as Repeat steps 2 − 5, M = 100000 times to get the samples of size M from the corresponding posteriors of interest. Obtain the Bayesian estimates of v, θ, α and β using the following formulae respectively, where M 0 (≈ 50000) is the burn-in period of the generated MCMC.

Simulation studies for comparing estimation methods
A numerical simulation is performed in to compare the classical estimation methods. The simulation study is based on N=1000 generated data sets from the OBP version where n = 50, 100, 150 and 300 and The estimates are compared in terms of their 1-Bias BIAS(ξ); 2-Root mean-standard error RMSE(ξ); 3-the mean of the absolute difference between the theoretical and the estimates "D-abs" and 4-the maximum absolute difference between the true parameters and estimates "D-max". Where  Tables 2, 3 and 4 we note that: 1-The BIAS(ξ) tend to zero when n increases which means that all estimators are non-biased. 2-The RMSE(ξ) tend to zero when n increases which means incidence of consistency property.

Comparing classical methods under failure times
The first real data set (data set I) represents the data on failure times of 84 aircraft windshield given in Murthy et al. [58].  Figure 6 gives probability-probability (P-P) plots for comparing all methods under the failure times data set. From Table 5, the MLE method is the best method with W ⋆ =0.05840 and A ⋆ =0.56749 then Bayesian method with W ⋆ =0.06309 and A ⋆ =0.63239.

Comparing classical methods under service times
The second real data set (data set II) represents the data on service times of 63 aircraft windshield given in Murthy et al. [58].  [18] and Ibrahim and Yousof [36]. Figure 7 gives P-P plots for comparing all methods under the failure times data set. From Table 6, the MLE method is the best method with W ⋆ =0.09443 and A ⋆ =0.57283 then the OLS method with W ⋆ =0.09443 and A ⋆ =0.71474.

Applications
In this section, we provide two real life applications to two real data sets to illustrate the importance and flexibility of the OBP model. We compare the fit of the OBP with some well-known competitive models (see Table 7)  [22] For checking the normality, the Quantile-Quantile (Q-Q) plot is sketched. For exploring the HRF for real data, the total time test (TTT) plot is provided. For exploring the initial shape of real data nonparametrically, kernel density estimation (KDE) is provided. Figures 12 and 13 give the box plot, normal Q-Q plot, TTT plot and nonparametric KDE for data set I and data set II respectively. In general, the smaller the values of these statistics, the better the fit to the data. The required computations are obtained by using the "maxLik" and "goftest" sub-routines in R-software. For failure times data: the analysis results of are listed in Tables 8 and 9. Table 8 gives the MLEs and standard errors (SEs) for failure times data. Table  9 gives the −l and goodness-of-fits statistics for failure times data. For service times data: the analysis results of are listed in Tables 10 and 11. Table 10 gives the MLEs and SEs for service times data. Table 11 give the −l and goodness-of-fits statistics for the service times data. Based on Tables 9 and 11, we note that the OBP model gives the lowest values for the AIC, CAIC, BIC, HQIC, A * and W * among all fitted models. Hence, it could be chosen as the best model under these criteria. Figure 14 gives the estimated PDF (EPDF) plot, estimated CDF (ECDF) plot, probability-probability (P-P) plot and estimated HRF (EHRF) plot for data set I. Figure 15 gives the EPDF plot, ECDF plot, P-P plot and EHRF plot for data set II. From Figures 14 and 15, it is noted that the OBP provides a very adequate fits to the empirical functions.     Figure 16. EPDF, ECDF, P-P plot and EHRF for data set I.

Conclusions
A new lifetime model called the odd Burr Pareto type II (OBP) model is introduced and studied. The major motivation for the practicality of the OBP model is based on the wider importance of the standard Pareto type II model. TheOBP density can be "right skewed" with heavy tail shape. The OBP failure rate can be "J-shape", "decreasing" and "upside down (or increasing-constant-decreasing)". The new OBP density can be expressed as as a mixture of the exponentiated Pareto type II model. The skewness of the OBP distribution can range in the interval (−69216.87, 98983.16). The spread for the kurtosis of the OBP model is ranging from 74.17124 to ∞. The index of dispersion for the OBP model can be in (0, 1) an dalso > 1 so it may be used as an "under-dispersed" and "over-dispersed" model. Bayesian and non-Bayesian estimation methods are considered. We assessed the performance of all methods via simulation study. Bayesian and many non-Bayesian estimation methods (such as uch as the maximum likelihood estimation method, ordinary least square estimation method, weighted least square estimationmethod and the Kolmogorov estimation method) are compared in modeling real data via two applications. In modeling real data, the maximum likelihood method is the best estimation method. So, we used it in comparing competitive models. Before using the the maximum likelihood method, we performed simulation experiments to assess the finite sample behavior of it using the biases and mean squared errors. The OBP model could be chosen as the best model among Pareto type II, exponentiated Pareto type II, Kumaraswamy Pareto type II, Macdonald Pareto type II, beta Pareto type II, gamma Pareto type II, odd log-logistic Pareto type II, reduced odd log-logistic Pareto type II, reduced Burr-Hatke Pareto type II, reduced OBP and special generalized mixture Pareto type II distribution in modeling the "failure times" and the "service times" data sets.