A New Bimodal Distribution for Modeling Asymmetric Bimodal Heavy-Tail Real Lifetime Data

We introduced and studied a new generalization of the Burr type X distribution. Some of its properties were derived and numerically analyzed. The new density can be “right-skewed” and symmetric with “unimodal” and many “bimodal” shapes. The new failure rate can be “increasing,” “bathtub,” “J-shape,” “decreasing,” “increasing-constant-increasing,” “reversed J-shape,” and “upside-down (reversed U-shape).” The usefulness and flexibility of the new distribution were illustrated by means of four asymmetric bimodal right- and left-heavy tail real lifetime data.


Introduction
Burr [1] introduced twelve different forms of cumulative distribution functions (CDFs) for modeling real data sets. Among those twelve CDFs, Burr type X (BX) and Burr type XII (BXII) have received special attention. A random variable (RV) W is said to have the BX distribution if its probability density function (PDF), and hence the CDF, are given by: π c 1 ,c 2 (w) = 2c 1 c 2 2 we −(c 2 w) 2 1 − e −(c 2 w) 2 c 1 −1 , and Π c 1 ,c 2 (w) = 1 − e −(c 2 w) 2 c 1 , Respectively, for W > 0 and c 1 , c 2 > 0. For c 1 = 1, we have the standard Rayleigh model. For c 2 = 1, we have the one-parameter BX model (for more details see [2][3][4][5][6][7]). For c 1 = 1, we have the Rayleigh (R) model. Recently, [8] observed that the BX model can be used in modeling strength and general lifetime data. Several aspects of the one-parameter (c 2 = 1) BX distribution were studied by [9], who studied the Bayesian prediction bounds under the BX model, [10] presented the Bayesian approach for prediction with outliers from the BX model, and [11] and [7] investigated the inference for P(Y < X) in the BX model. Many authors have studied the BX distribution and applied it in different areas, e.g., [12] presented an overview about the BX model, [13] proposed a new version of the BX model called the beta Burr type X (BBX) distribution and discussed its relevant mathematical properties along with some real data applications, and [14] presented an estimation of reliability under the one-parameter BX model. [15] Proposed a new generator of distributions based on the BX model. [16] Compounded the Poisson model and the BX model and introduced a new compound generator called the Poisson BX family. Recently, [17] used the Weibull and the BX model to generate a new flexible model and presented a new modified Bagdonavičius-Nikulin goodness-of-fit test for censored validation; [18] presented the quasi Poisson BX-BX distribution along with copula, mathematical properties, and applications; [19] generalized the odd log-logistic BX distribution with an application to the times of failure and running times for a sample of devices from a field-tracking study of a larger system.
By examining the statistical literature, it is noted that all the proposed versions of the BX model have a unimodal density. In this work, we present a new bimodal version of the BX model called the odd Burr BX (OBBX) model based on the family from [20], who merged the generalized odd G (OG) family and the proportional reversed hazard rate family (PRHR) to propose a new wider family called the odd Burr G (OBG) family. The CDF of the OBG family is given by: whereΠ ζ (w) = 1 − Π ζ (w), ν and θ > 0 are two shape parameters, and ζ refers to the parameter vector of any baseline model and the reliability function (RF) of the baseline model. The PDF corresponding to (3) is given by: For θ = 1, the OBG family reduces to the OG family. For ν = 1, the OBG family reduces to the PRHR. Using (3) and (4), the OBBX RF is given by: where Θ = ν, θ, c 1 , c 2 ; O w,c 2 = e −(c 2 w) 2 ; S Θ (w) = 1 − F Θ (w). For θ = 1, the OBBX reduces to the O-BX.
For ν = 1, the OBBX reduces to the PRHR-BX. The PDF corresponding to (5) is given by: The hazard rate function (HRF) for the new model can be obtained from f Θ (w)/S Θ (w).
Many useful distributions are introduced based on the BX model, such as the BX Weibull (BXW) distribution [17], the Burr X Fréchet (BXFr) distribution [21], the Burr X exponentiated Weibull (BXEW) distribution [22], and the Burr X Nadarajah Haghighi (BXNH) distribution [23]. We were motivated to introduce and study the OBBX for the following reasons: 1.
The new density in (6) can be "unimodal and right-skewed," "symmetric and unimodal," and "bimodal density" with many useful shapes (see Figure 1).
In reliability analysis, the OBBX model could be chosen as the best model, especially in modeling asymmetric bimodal failure times data and the asymmetric bimodal right-skewed and heavy-tail survival times data as illustrated in Sections 5.1 and 5.3, respectively. 4.
In medical fields, the OBBX model could be chosen as the best model, especially in modeling the bimodal right-skewed and heavy-tail cancer data, as illustrated in Section 5.2.    The asymptotic behavior of the CDF, PDF, and HRF as W → 0 are respectively given by: The asymptotic behavior of the CDF, PDF, and HRF as W → ∞ are respectively given by: Figure 1 gives some plots of the OBBX PDF for selected parameter values. Figure 2 gives some plots of the OBBX HRF for selected parameter values. Based on Figure 1, the new density can be "right-skewed" and symmetric with "unimodal" and many "bimodal" shapes. Based on Figure 2, the new HRF can take the following forms: "upside-down (reversed U-HRF)" (Figure 2g).
For the simulations of this new model, we obtained the quantile function (QF) of W (by inverting the CDF), say w u = F −1 (u), as: where (7) was used as the random number generator from the OBBX model.

Useful Representations
Based on [20], the PDF in (6) can be expressed as: where: and π c * 1 ,c 2 (w) is the PDF of the BX model. By integrating (8), the CDF of W becomes: where Π c * 1 ,c 2 (w) refers to the CDF of the BX distribution.

Moments and Incomplete Moments
The rth ordinary moment of W is given by: Then, we obtain: where: The variance (V(W)), skewness (S(W)), and kurtosis (K(W)) can be derived easily using the well-known relationships. The rth incomplete moment, say I r (τ), of W can be expressed, from (9), as: Then: where γ(φ, ψ) refers to the incomplete gamma function: The first incomplete moment is given by (11) with r = 1 as: The dispersion index (DisIx), also known as the variance-to-mean ratio, is a measure used to quantify whether a set of observed occurrences are clustered or dispersed compared to a standard statistical model. Therefore, it indicates whether a certain statistical model is suitable for over- (or under-) dispersed data sets and is used widely in ecology as a standard measure for measuring clustering (overdispersion) or repulsion (underdispersion). Thus, the measure can be used to assess whether observed data can be modeled using a Poisson process. For any real data set, when the DisIx is less than 1, the data set is said to be "under-dispersed"; this important condition can relate to occurrence patterns that are more regular than the randomness associated with a Poisson process. Numerical analysis for the DisIx(W) of the new OBBX is presented in Table 1 with useful comments.

Moment Generating Function (MGF)
The MGF M W (τ) = E(exp(τW)) of W can be derived from (8) as: where M c * 1 (τ) is the MGF of the BX model. Then:

Residual Life and Reversed Residual Life Functions
The rth moment of the residual life A r (τ) = E (W − τ) r | w>τ,r=1,2,... . Then, the rth moment of the residual life of W is given by A r (τ) = 1 where a (r,c * 1 ) r). The rth moment of the reversed residual life, say: uniquely determines F(w). Then, we obtain: Then, the rth moment of the reversed residual life of W becomes: 2.5. Numerical Analysis for the Mean, V(W), S(W), K(W), and DisIx(W). Table 1 gives numerical calculations for the mean, V(W), S(W), K(W), and dispersion index (DisIx(W)) for selected parameter values. Based on Table 1, we note that (1) the skewness of the OBBX model can be both positive and negative; (2) the spread for the OBBX kurtosis is much larger, ranging from 1.0005 to 172,196.5; (3) DisIx(W) can be "between 0 and 1" or "more than 1." For more visualization, Figure 3 gives some three-dimensional skewness plots and Figure 4 gives three-dimensional kurtosis plots.

Maximum Likelihood Estimation
The maximum likelihood estimations (MLEs) display desirable properties and can be used for establishing confidence intervals and test statistics. The normal approximation for MLEs in large sample theory is easily handled either numerically or analytically. In this section, we determine the MLEs of the parameters of the OBBX distribution from complete samples only. However, censored samples could be considered in separate works. Let 1 , 2 , …, be an observed random sample

Maximum Likelihood Estimation
The maximum likelihood estimations (MLEs) display desirable properties and can be used for establishing confidence intervals and test statistics. The normal approximation for MLEs in large sample theory is easily handled either numerically or analytically. In this section, we determine the MLEs of the parameters of the OBBX distribution from complete samples only. However, censored samples could be considered in separate works. Let w 1 , w 2 , . . . , w n be an observed random sample from the OBBX model with parameters ν, θ, c 1 , and c 2 . Then, the log-likelihood function for Θ , say L(ν, θ, c 1 , c 2 ), is given by: Symmetry 2020, 12, x FOR PEER REVIEW 11 of 28 from the OBBX model with parameters , , 1 , and 2 . Then, the log-likelihood function for , say ℒ( , , 1 , 2 ), is given by: The equation of ( ) can be maximized either directly by using the R program via the optim function, the SAS program via PROC NLMIXED, or the Ox program using sub-routine MaxBFGS or by solving the nonlinear likelihood equations obtained by differentiating ( , , , ). The score vector components, say: are available from the corresponding author, where: Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods.
where Symmetry 2020, 12, x FOR PEER REVIEW 11 of 28 from the OBBX model with parameters , , 1 , and 2 . Then, the log-likelihood function for , say ℒ( , , 1 , 2 ), is given by: The equation of ( ) can be maximized either directly by using the R program via the optim function, the SAS program via PROC NLMIXED, or the Ox program using sub-routine MaxBFGS or by solving the nonlinear likelihood equations obtained by differentiating ( , , , ). The score vector components, say: are available from the corresponding author, where: Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods. from the OBBX model with parameters , , 1 , and 2 . Then, the log-likelihood function for , say ℒ( , , 1 , 2 ), is given by: The equation of ( ) can be maximized either directly by using the R program via the optim function, the SAS program via PROC NLMIXED, or the Ox program using sub-routine MaxBFGS or by solving the nonlinear likelihood equations obtained by differentiating ( , , , ). The score vector components, say: are available from the corresponding author, where: Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods.
can be maximized either directly by using the R program via the optim function, the SAS program via PROC NLMIXED, or the Ox program using sub-routine MaxBFGS or by solving the nonlinear likelihood equations obtained by differentiating (ν, θ, c 1 , c 2 ). The score vector components, say: Symmetry 2020, 12, x FOR PEER REVIEW 11 of 28 from the OBBX model with parameters , , 1 , and 2 . Then, the log-likelihood function for , say ℒ( , , 1 , 2 ), is given by: The equation of ( ) can be maximized either directly by using the R program via the optim function, the SAS program via PROC NLMIXED, or the Ox program using sub-routine MaxBFGS or by solving the nonlinear likelihood equations obtained by differentiating ( , , , ). The score vector components, say: are available from the corresponding author, where: Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods. are available from the corresponding author, where: Symmetry 2020, 12, x FOR PEER REVIEW 11 of 28 from the OBBX model with parameters , , 1 , and 2 . Then, the log-likelihood function for , say ℒ( , , 1 , 2 ), is given by: Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods. Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods. Setting the nonlinear system of equations to U( ) = U( ) = U( 1 ) = U( 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J( ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods. Setting the nonlinear system of equations to U(ν) = U(θ) = U(c 1 ) = U(c 2 ) = 0 and solving them simultaneously yields the MLEs. These equations cannot be solved analytically, but Newton-Raphson-type algorithms can be used to solve them numerically (see the Appendix A). For the interval estimation of the OBBX model parameters, we require the observed information matrix J(Θ ), which comes as the output using the above maximization procedures. Likelihood ratio tests can be performed for the proposed model in the usual way. Further works could be addressed using different estimation methods to estimate the parameters of the OBBX model, such as bootstrap, least squares, Cramér-von Mises, weighted least squares, Jackknife, Anderson-Darling, Bayesian analysis, and compare the estimators based on these methods.

Graphical Assessment
Graphically and using the biases and mean squared errors (MSEs), we can perform the simulation experiments to assess the finite sample behavior of the MLEs given in Section 4. The assessment was based on the following algorithm (see the Appendix A):

Graphical Assessment
Graphically and using the biases and mean squared errors (MSEs), we can perform the simulation experiments to assess the finite sample behavior of the MLEs given in Section 4. The assessment was based on the following algorithm (see the Appendix A):

Applications
In this section, we provide four applications of the OBBX distribution to empirically show its potentiality. In order to compare the fits of the OBBX distribution with other competing distributions, we consider the Cramér-von Mises (CVM), the Anderson-Darling (AD), and the Kolmogorov-Smirnov (KS) (and its corresponding p-value) estimation methods. These four methods are widely used to determine how closely a specific CDF fits the empirical distribution of a given data set. The smaller the resulting statistics are, the better the fit. The required computations were carried out using

Applications
In this section, we provide four applications of the OBBX distribution to empirically show its potentiality. In order to compare the fits of the OBBX distribution with other competing distributions, we consider the Cramér-von Mises (CVM), the Anderson-Darling (AD), and the Kolmogorov-Smirnov (KS) (and its corresponding p-value) estimation methods. These four methods are widely used to determine how closely a specific CDF fits the empirical distribution of a given data set. The smaller the resulting statistics are, the better the fit. The required computations were carried out using the R software. For data set I: the MLEs and the corresponding standard errors ((SEs) in parentheses) for all the competitive parameters are given in Table 2, the numerical values of the statistics from CVM, AD, and KS (corresponding p-value) are listed in Table 3. For data set II: the MLEs and SEs in for all the competitive parameters are given in Table 4, the numerical values of the statistics from CVM, AD, and KS (corresponding p-value) are listed in Table 5. For data set III: the MLEs and SEs in for all the competitive parameters are given in Table 6, the numerical values of the statistics from CVM, AD, and KS (corresponding p-value) are listed in Table 7. For data set IV: the MLEs and SEs in for all the competitive parameters are given in Table 8, the numerical values of the statistics from CVM, AD, and KS (corresponding p-value) are listed in Table 9.   The total time in test (TTT) plot, box plot, quantile-quantile (Q-Q) plot, and the nonparametric Kernel density estimation (NKDE) plot for data set I are displayed in Figure 9. The estimated PDF (EPDF), estimated CDF (ECDF), probability-probability (P-P), and estimated HRF (EHRF) plots for data set I are displayed in Figure 10. The TTT, box plot, Q-Q and the NKDE plots for data set II are displayed in Figure 11. The EPDF, ECDF, P-P and EHRF plots for data set II are displayed in Figure 12. The TTT, box plot, Q-Q and the NKDE plots for data set II are displayed in Figure 13. The EPDF, ECDF, P-P and EHRF plots for data set III are displayed in Figure 14. The TTT, box plot, Q-Q and the NKDE plots for data set II are displayed in Figure 15. The EPDF, ECDF, P-P and EHRF plots for data set IV are displayed in Figure 16. The dashed line in all the Q-Q plots refers to the safe boundaries for the standard errors. The TTT plots show that the HRFs were "increasing" (data set I), "upside-down" (data set II), "increasing" (data set III), and "increasing" (data set IV). The box plots show that data sets II, III, and IV had some extreme values. The Q-Q plots ensure the results of the box plots. The NKDE plots show that the kernel density estimation was "asymmetric bimodal" (data set I), "asymmetric bimodal with a right-heavy tail" (data sets II and III), and "asymmetric bimodal with a left-heavy tail" (data sets IV). The R codes are given in see Appendix A.     . Estimated probability density function (EPDF), estimated cumulative distribution function (ECDF), probability-probability (P-P), and estimate hazard rate function (EHRF) plots for data set I. Figure 11. The TTT, box, Q-Q, and NKDE plots for data set II. Figure 11. The TTT, box, Q-Q, and NKDE plots for data set II.
Symmetry 2020, 12, x FOR PEER REVIEW 19 of 28 Figure 12. EPDF, ECDF, P-P, and EHRF plots for data set II. Figure 12. EPDF, ECDF, P-P, and EHRF plots for data set II. Figure 12. EPDF, ECDF, P-P, and EHRF plots for data set II. Figure 13. The TTT, box, Q-Q, and NKDE plots for data set III. Figure 13. The TTT, box, Q-Q, and NKDE plots for data set III.

Modeling Cancer Data
This data set represents the remission times (in months) of a random sample of 128 bladder cancer patients as reported in [44]. We compared the fits of the OBBX distribution with other competitive models, namely, the BX, Weibull (W) [45], TrMW, MBW, and transmuted additive W distribution (TrAW) ( [34]) distributions with corresponding densities (for W > 0). Figure 11 presents the TTT, box, Q-Q, and NKDE plots for data set II. Figure 12 shows the EPDF, ECDF, P-P, and EHRF plots for data set II. Table 4 provides the MLEs (SEs) for data set II. Table 5 gives the CVM, AD, and KS (p-value) for data set II. Based on Table 5 and Figure 12, we concluded that the proposed OBBX lifetime model was much better than the W, TrMW, MBW, TrAW, and ETrGR models with small values for Cr [1] and Cr [2] in modeling cancer patient data.

Modeling Survival Times
The second real data set corresponds to the survival times (in days) of 72 guinea pigs infected with virulent tubercle bacilli reported by [46]. We compared the fits of the OBBX distribution with those of other competitive models, namely, the BX, transmuted (TrBX), OLEW, Odd W-W (OWW) [47], gamma exponentiated-exponential (GaEE) [48] distributions with PDFs (for W > 0). Figure 13 presents the TTT, box, Q-Q, and NKDE plots for data set III. Figure 15 shows the EPDF, ECDF, P-P, and EHRF plots for data set III. Table 6 provides the MLEs (SEs) for data set III. Table 7 gives the CVM, AD, and KS (p-value) for data set III. Based on Table 7 and Figure 15, we concluded that the proposed OBBX model was much better than all these models and was a good alternative to these models for modeling the survival times of guinea pigs.

Glass Fibers Data
This data consists of 63 observations of the strengths of 1.5 cm glass fibers, originally obtained by workers at the U.K. National Physical Laboratory. These data are given in [49]. For this data set, we compared the fits of the new distribution with some competitive models, such as the BX, TrBX, OLBX, OLEW, EW, TrW, and the odd log-logistic Weibull (OLLW) models. Figure 15 presents the TTT, box, Q-Q, and NKDE plots for data set IV. Figure 16 shows the EPDF, ECDF, P-P, and EHRF plots for data set IV. Table 8 provides the MLEs (SEs) for data set IV. Table 9 gives the CVM, AD, and KS (p-value) for data set IV. Based on Table 9 and Figure 16, we concluded that the proposed OBBX model was much better than all the competitor models and was a good alternative to these models in modeling glass fibers data.

Concluding Remarks
We introduced and studied a new generalization of the Burr type X distribution called the odd Burr-Burr type X (OBBX) distribution. Some of its properties were derived and numerically analyzed. The new PDF of the OBBX model was re-expressed in terms of the standard BX model. Graphically, we assessed the performance of the estimation method by means of biases and mean squared errors. The usefulness and flexibility of the new distribution was illustrated by means of four real data sets. The superiority of the new model against many other competitor models was proved using the Cramér-von Mises, the Anderson-Darling, and the Kolmogorov-Smirnov (and its corresponding p-value) statistics. The following concluding remarks can be made: The new density can be "unimodal and right-skewed," "symmetric and unimodal," and "bimodal density" with many useful shapes. II.
In the reliability analysis, the OBBX model could be chosen as the best model, especially for modeling the asymmetric bimodal failure times data and the asymmetric bimodal right-skewed and heavy-tail survival times data. IV.
In medical fields, the OBBX model could be chosen as the best model, especially for modeling the bimodal right-skewed and heavy-tail cancer data. V.
In engineering, the OBBX model could be chosen as the best model, especially for modeling the asymmetric bimodal left-skewed and heavy-tail glass fibers data. VI.