A Novel Lomax Extension with Statistical Properties, Copulas, Different Estimation Methods and Applications

A new compound Lomax model is proposed and analyzed. The novel distribution is derived based on compounding the zero truncated Poisson distribution and the exponentiated exponential Lomax distribution. The new density can be “monotonically left skewed,” “monotonically right skewed” and “symmetric” with various useful shapes. The new hazard rate can be “upside down bathtub-increasing,” “bathtub (U-shape),” “monotonically decreasing,” “increasing-constant” and “monotonically increasing.” Relevant statistical properties are derived. We briefly describe different estimation methods, namely the maximum likelihood, Cramér-von-Mises, ordinary least squares, weighted least square, Anderson–Darling, right tail Anderson–Darling and left tail Anderson–Darling. Monte Carlo simulation experiments are performed for comparing the performances of the proposed methods of estimation for both small and large samples. For facilitating the mathematical modeling of the bivariate real data sets, we derive some new corresponding bivariate distributions. Graphical simulation study is performed for assessing the finite sample behavior of the estimators using the maximum likelihood method. Two applications are provided for illustrating the applicability of the new model.

A random variable (RV) is said to have the Lomax distribution if its cumulative distribution function (CDF) is given by where θ > 0 refers to the shape parameter. The above CDF is a special case from the Burr type XII (BXII) model. So, many useful details about the Lomax model along with its relationship with other related models can be found in Burr [11][12][13], Lomax [33], Burr and Cislak [14], Harris [25], Rodriguez [34], Tadikamalla [35] and Yadav et al. [44]. In this paper, we propose and study a new compound version Lomax (L) distribution using zero truncated Poisson (ZTP) distribution. Suppose that a system has N (a discrete random variable) subsystems functioning independently at a given time where N has ZTP distribution with parameter a and the failure time of ith component Y i |i 1, 2, . . . (say), independent of N . It is the conditional probability distribution of a Poisson-distributed random variable (RV), given that the value of the RV is not zero. The probability mass function (PMF) of N is given by P a (N n) a n exp(−a) (1 + n) 1 − exp(−a) | a∈R + .
Note that for ZTP RV, the expected value E(N |a) and variance V(N |a) are, respectively, given by and Suppose that for each sub-device, the failure time (i.e., ith sub-device) has the exponentiated exponential Lomax (EEL) and having the following CDF where b > 0 is the shape parameter. For b 1, the exponentiated exponential Lomax model reduces to exponential Lomax model. For θ 1, the exponentiated exponential Lomax model reduces to exponentiated exponential model. For b θ 1, the exponentiated exponential Lomax model reduces to exponential model. Let Y i denotes the failure time of the ith subsystem and let Z min{Y 1 , Y 2 , . . . , Y N }. Then, the conditional CDF of Z given N is Therefore, the unconditional CDF of X, as described in Aryal and Yousof [5], Korkmaz et al. [31], Alizadeh et al. [6], can be expressed as The CDF in (2) is called the Poisson exponentiated exponential Lomax (PEEL) model. The corresponding probability density function (PDF) can be derived as A RV Z having PDF (3) will be denoted by Z~PEEL (a, b, θ). The PDF in (3) is said to be "concave PDF" if for any Z 1 ∼ PEEL (a 1 , b 1 , θ 1 ) and Z 2 ∼ PEEL (a 2 , b 2 , θ 2 ) the PDF satisfies If the function f (εz 1 If the function f (εz 1 is convex for every ψ > 0. If f (εz 1 + εz 2 ) and g(εz 1 + εz 2 ) are "convex PDF," then [ f (εz 1 + εz 2 ) + g(εz 1 + εz 2 )] is also "convex PDF." If f (εz 1 + εz 2 ) and g(εz 1 + εz 2 ) are "convex PDF," then [ f (εz 1 + εz 2 ) and g(εz 1 + εz 2 )] is also "convex PDF." If the function − f (εz 1 + εz 2 ) is "convex PDF," then the function  Figure 1 (left plot) gives some plots of the PEEL PDF. Figure 1 (right plot) gives some plots of the PEEL HRF for some selected values of the parameters. Based on Fig. 1 (left plot), it is noted that the PDF of the PEEL can be "monotonically left skewed," "monotonically right skewed" and "symmetric" with various useful shapes. Based on Fig. 1 (right plot), it is noted that the HRF of the PEEL can be "upside down bathtubincreasing," "bathtub (U-shape)," "monotonically decreasing," "increasing-constant" and "monotonically increasing." The PEEL model could be useful in modeling the "asymmetric monotonically increasing HRF" real data sets as illustrated in Sect. 6 First, for facilitating the mathematical modeling of the bivariate RVs, we derive some new bivariate PEEL (BPEEL) type distributions using "Farlie-Gumbel-Morgenstern copula" (FGMCp) Morgenstern [37], Farlie [20], Gumbel [22], Gumbel [23], Johnson [29] and Johnson [30], modified FGMCp which contains four internal types, " Clayton copula (CCp)" (see Nelsen [39] for details), "Renyi's entropy copula (RECp)" (Pougaza and Djafari [40] and "Ali-Mikhail-Haq copula (AMHCp)" [4]. The multivariate PEEL (MPEEL) type can be easily derived based on the Clayton copula. However, future works may be allocated to study these new models. Additionally, we briefly describe different estimation methods, namely the maximum likelihood, Cramér-von-Mises, ordinary least square, weighted least square, Anderson-Darling, right tail Anderson-Darling and left tail Anderson-Darling. Monte Carlo simulation experiments are performed for comparing the performances of the proposed methods of estimation for both small and large samples. The above-mentioned estimation methods are compared also using two real data sets.
Second, we briefly considered and then described different estimation methods, namely the maximum likelihood estimation (MLE) method, Cramér-von-Mises estimation (CVM) method, ordinary least square estimation (OLS) method, weighted least square estimation (WLSE) method, Anderson-Darling estimation (ADE) method, right tail Anderson-Darling estimation (RTADE) method, left tail Anderson-Darling estimation (LTADE) method. These methods are used in estimation process of the unknown parameters.
Generally, in statistical modeling of the failure times of the aircraft windshield data, the PEEL model is compared with many common Lomax extensions such as the special generalized mixture Lomax, the odd log-logistic Lomax, the Burr-Hatke Lomax, the transmuted Topp-Leone Lomax, the Gamma Lomax, the Kumaraswamy Lomax, the McDonald Lomax, the exponentiated Lomax and the proportional reversed hazard rate Lomax under the Akaike information criteria, consistent-information criteria, Bayesian information criteria and Hannan-Quinn information criteria. The PEEL model provided better fits in modeling failure times of the aircraft windshield data. In statistical modeling of the service times of the aircraft windshield, the PEEL model is compared with many common Lomax extensions such as the special generalized mixture Lomax, the odd log-logistic Lomax, the Burr-Hatke Lomax, the transmuted Topp-Leone Lomax, the Gamma Lomax, the Kumaraswamy Lomax, the McDonald Lomax, the exponentiated Lomax and the proportional reversed hazard rate Lomax under the Akaike information criteria, consistent-information criteria, Bayesian information criteria and Hannan-Quinn information criteria. The PEEL model provided better fits in modeling service times of the aircraft windshield data.
In this paper, after studying the main statistical properties and presenting some bivariate type extensions, we briefly considered and then described different estimations. Monte Carlo simulation experiments are performed for comparing the performances of the proposed methods of estimation for both small and large samples.

Expanding the New Density
We present a simple useful expansion for the new PDF given on (3) in terms of the exponentiated Lx (exp-L) model. Using the obtained expansion, we derive the main mathematical properties of the new PDF of the PEE model. Note thatThen, the PDF in (4) can be expressed as Considering the power series and applying (5) to the quantity B b(h+1), (z) in (4), we get Inserting the expansion of the quantity B b(h+1), (z) into (4), then, the PDF of the PEEL can be expressed as Expanding the quantity C (l+1) (z), we can write Inserting the result of (7) into (6), the PEEL density can be reduced to Expanding (1 + z) −τ −2 via generalized binomial expansion, we get Inserting (9) in (8), the PDF of the PEEL can be summarized as where is the PDF of the exp-L model with power τ + d + 1 and ε τ ,d is a constant where Similarly, the CDF of the PEEL model can also be expressed as where H τ +d+1 (z; θ) 1 − (1 + z) −θ τ +d+1 is the CDF of the exp-L family with power τ + d + 1.

Moments
The calculations of this subsection involve several special functions, including the complete beta function B(l 1 , l 2 ) 3 ) , and the upper incomplete gamma function. Noting that Let Z be a RV having the PEEL (a, b, θ) model. Then, the pth moment of the RV Z is

Moment Generating Function (MGF)
Clearly, the MGF can be derived from Eq. (10) as

Incomplete Moments
The pth incomplete moments, say I p,Z (t), of the RV Z can be obtained from (10) as Then, the pth incomplete moments can be written as the 1st incomplete moments can be written as is the median of the RV Z , and I 1,Z (t) is the first incomplete moment is given by These results for I 1,Z (t) can be directly applied for calculating the Bonferroni (Bon( )) and Lorenz (Lor ( )) curves defined, for a certain given probability, say , by Bon( ) I 1,Z (Q( ))/ μ 1,Z and Lor ( ) I 1,Z (Q( ))/μ 1,Z , respectively.

Residual Life (RL) and Reversed Residual Life (RRL)
The jth moment of the RL of the RV Z can be obtained from w j, which can also be written as Then, which can also be expressed as Then, For j 1, we obtain the mean waiting time (MWT) which also called the mean inactivity time (MIT) which can be derived where ε τ ,d,l (W, 1) ε τ ,d  Table 1 gives some numerical calculations for the mean (E(Z )), variance (V (Z )), skewness (S(Z )) and kurtosis (K(Z )) for PEEL distribution. Based on results listed in Table 1, it is noted that S(Z ) ∈ (0.9, 1730.1) and K(Z ) ranging from 3.525878 to 3,332,135.

BPEEL Type via CCp
Let us assume that X 1 ∼ PEEL(a 1 , b 1 , θ 1 ) and X 2 ∼ PEEL(a 2 , b 2 , θ 2 ). The CCp depending on the continuous marginal functions U 1 − U and V 1 − V can be considered as Then, the BPEEL type distribution can be obtained from (5). A straightforward multivariate PEEL (mdimensional extension) via CCp can be easily derived analogously. The m-dimensional extension via CCp which is function operating in [0, 1] m , and in that case, z i is not a value in [0, 1] necessarily.

BPEEL type via RECp
Following Pougaza and Djafari [40], the RECp can be derived as where the values z 1 and z 2 are in order to guarantee that C(U, V) is a copula. Then, the associated CDF of the BPEEL will be It is worth mentioning that in [18], the authors emphasize that this copula does not show a closed form and numerical approaches become necessary.

BPEEL type via FGMCp
Considering the FGMCp (see ( [10][11][12][13][14][15]), the joint CDF can be written as where the continuous marginal function U ∈ (0, 1), V ∈ (0, 1) and where which is "grounded minimum condition" and and which is "grounded maximum condition ." The grounded minimum/maximum conditions are valid for any copula. Setting U The joint PDF can be derived from or from where the two function and are densities corresponding to the joint CDFs and .

BPEEL type via modified FGMCp
The modified formula of the modified FGMCp can be written as The following four types can be derived and considered:

BPEEL type via AMHCp
Under the "stronger Lipschitz condition," the joint CDF of the Archimedean AMHCp can be written as the corresponding joint PDF of the Archimedean AMHCp can be expressed as ,we have and

Estimation Methods
In this Section, we briefly describe and consider different classical estimation methods, namely the MLE method, CVM method, OLS method, WLSE method, ADE method, RTADE method, left tail LTADE. All these methods are discussed in the statistical literature with more details. In this work, we may ignore some of its derivation details for avoiding repetition.

The ML Method
Let This system can only be solved numerically for the complicated models using some common iterative algorithms such as the "Newton-Raphson" algorithm.

The CVM Method
The CVM estimates (CVMEs) of the parameters a, b and θ are obtained via minimizing the following expression with respect to the parameters a, b and θ , respectively, where The CVME of the parameters a, b and θ is obtained by solving the following nonlinear equations

The OLS Method
Let F a,b,θ z [i,m] denote the CDF of PEEL model and let z 1 < z 2 < · · · < z m be the m ordered random sample. The OLS estimates (OLSEs) are obtained upon minimizing and equivalently The OLSEs are obtained via solving the following nonlinear equations

The WLS Method
The WLS estimates (WLSEs) are obtained by minimizing the function WLSE (a, b, θ) with respect to a, b and θ where and The WLSEs are obtained by solving

The AD Method
The AD estimates (ADEs) of a, b and θ are obtained by minimizing the function The parameter estimates of a, b and θ follow by solving the nonlinear equations

The RTAD Method
The RTAD estimates (RTADEs) of a, b and θ are obtained by minimizing

The LTAD Method
The LTAD estimates (LTADEs) of a, b and θ are obtained by minimizing The parameter estimates of δ, θ and β are obtained by solving the nonlinear equations

Simulations for Comparing 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 PEEL version where m 50, 100, 150 and 300 and The estimates are compared in terms of their bias and the root mean-standard error (RMSE). The mean of the absolute difference between the theoretical and the estimates (D-abs) and the maximum absolute difference between the true parameters and estimates (D-max) are also reported. Tables 2, 3 and 4 give the simulation results. From Tables 2, 3 and 4, we note that the RMSE tends to zero when m increases which means incidence of consistency property.

Comparing Methods
Two applications to real data sets are considered for comparing the estimation methods. The data set I represents the data on failure times of 84 aircraft windshields. The data set II represents the data on service times of 63 aircraft windshields. The two real data sets were reported by Murthy et al. [38]. The required computations are carried out using the MATHCAD software. In order to compare the estimation methods, we consider the Cramér-von Mises (CVM) and the Anderson-Darling (AD) statistics. These two statistics are widely used to determine how closely a specific CDF fits the empirical distribution of a given data set. The results are given in Tables 5 and 6. From Table 5, we conclude that the ML method is the best method with CVM* 0.06444 and AD* 0.64651. From Table 6, we conclude that the RTAD method is the best method with CVM* 0.10075 and AD* 0.61025. However, all other methods performed well.

Applications
In this Section, we consider the same two real data sets of Murthy et al. [38] for applications to show the flexibility and the importance of the family presented under the L case. The fits of the PEEL are compared with many common Lomax Table 3 Simulation results for parameters a  Table 4 Simulation results for parameters a  XI. Lomax (L) Lomax [33]. XII. Proportional reversed hazard rate Lomax (PRHRL).
These two data sets are considered by matching their properties and the shapes of the PDF of the new model (see Fig. 1 (right plot)). By examining Fig. 1 (the right panel), it is noted that the new PDF can be "symmetric" and also "asymmetric right skewed function" with variable shapes. Additionally, by examining the initial density shapes of the two real data sets, it is seen that the initial densities are "semi symmetric" PDFs. Furthermore, the HRF of the new family can be "upside down bathtub-increasing," "bathtub (U-shape)," monotonically decreasing, "increasing-constant" and "monotonically increasing" (see Fig. 1 (left plot)). Many other symmetric and asymmetric useful data sets can be found in Yousof et al. [41], Aryal et al. [9], Altun et al. [7]. For model comparison, some competitive models using a certain real data set (sets), we first need to explore the data. Exploring real data set can be used either numerically or graphically or with both techniques. In this section, we will consider many graphical techniques such as the skewness-kurtosis plot (or the Cullen and Frey plot) for exploring initial fit to the theoretical distributions such as normal, uniform, exponential, logistic, beta, lognormal and Weibull distributions (see Fig. 2 top left (1st data) and Fig. 2 top right (2nd data)). Bootstrapping is applied and plotted for more accuracy. Cullen and Frey plot just compares distributions in terms of squared skewness and kurtosis. This is a good summary but still only a summary of the properties of a distribution. The scattergram plots are also given in Fig. 2 middle left and bottom left for the 1st data and Fig. 2 middle right and bottom right for the 2nd data.
The "normality" of the two real data sets is checked using the "Quantile-Quantile" (Q-Q) plot (Figs. 3, 4 (top right plots)). The initial HRFs shape is explored by using the "total time in test (TTT)" tool (Figs. 3, 4 (bottom left plots)). The "nonparametric Kernel density estimation (NKDE)" tool is used for exploring the initial PDF shape (Figs. 3, 4 (top left plots)). The outliers are checked by the "box plot" (Figs. 3, 4 (bottom right plots)). Based on Figs. 3 and 4 (top left plots), it is seen that the NKDE is bimodal and semi-symmetric functions. Based on Figs. 3, 4 (top right plots), it is seen that the "normality" nearly exists for the two data sets (bottom left plots). It is  shown that the HRF is " monotonically increasing HRF" for the two data sets. From Figs. 3 and 4 (bottom right plots), it is observed that no extreme values were spotted.
The following goodness-of-fit (GOF) statistics are used for comparing competitive models: I. The "Akaike information" (AICr). II. The "consistent-AIC" (CAICr). III. The "Bayesian-IC" (BICr). IV. The "Hannan-Quinn-IC" (HQICr). Tables 7 and 9 give the MLEs and the corresponding SEs for the two data sets, respectively. Tables 8 and 10 give the four GOF tests for the two data sets, respectively. Figures 5 and 6 give fitted PDFs, the Probability-Probability (P-P) plots, Kaplan-Meier Survival (KMS) plot and estimated HRF (E-HRF) plot for the two data sets, respectively. Based on Tables 3 and 5

Conclusions
In this work, a new compound Lomax model called the Poisson exponentiated exponential Lomax distribution is proposed and analyzed. The Poisson exponentiated exponential Lomax distribution is derived based on compounding the zero truncated Poisson distribution and the exponentiated exponential Lomax distribution. The  new density can be "monotonically left skewed," "monotonically right skewed" and "symmetric" with various useful shapes. The new hazard rate can be "upside down bathtub-increasing," "bathtub (U-shape)," "monotonically decreasing," "increasingconstant" and "monotonically increasing." Relevant statistical properties such as ordinary moments, incomplete moments, moments of residual life, moments of the reversed residual life and mean deviation are derived. For facilitating the mathematical modeling of the bivariate real data sets, we derive some new bivariate Poisson exponentiated exponential Lomax distributions using "Farlie-Gumbel-Morgenstern copula, modified Farlie-Gumbel-Morgenstern copula which contains four internal types," Clayton copula," "Renyi's entropy copula (RECp)" and "Ali-Mikhail-Haq copula." However, future works may be allocated to study these new models. After studying the main statistical properties and presenting some bivariate type extensions, we briefly considered and then described different estimation methods, namely the maximum likelihood, Cramér-von-Mises estimation, ordinary least square, weighted least square, Anderson-Darling, right tail Anderson-Darling and left tail Anderson-Darling. These methods are used in estimation process of the unknown parameters. Monte Carlo simulation experiments are performed for comparing the performances of the proposed methods of estimation for both small and large samples. Furthermore, two applications are provided for illustrating the applicability of the Poisson exponentiated exponential Lomax model. The Kernel density plot, the "Quantile-Quantile plot, the total time in test plot and box plot are provided and analyzed. Based on two real data sets, the Poisson exponentiated exponential Lomax    Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.