A New Extended Model with Bathtub-Shaped Failure Rate: Properties, Inference, Simulation, and Applications

: Theoretical and applied researchers have been frequently interested in proposing alternative skewed and symmetric lifetime parametric models that provide greater ﬂexibility in modeling real-life data in several applied sciences. To ﬁll this gap, we introduce a three-parameter bounded lifetime model called the exponentiated new power function (E-NPF) distribution. Some of its mathematical and reliability features are discussed. Furthermore, many possible shapes over certain choices of the model parameters are presented to understand the behavior of the density and hazard rate functions. For the estimation of the model parameters, we utilize eight classical approaches of estimation and provide a simulation study to assess and explore the asymptotic behaviors of these estimators. The maximum likelihood approach is used to estimate the E-NPF parameters under the type II censored samples. The efﬁciency of the E-NPF distribution is evaluated by modeling three lifetime datasets, showing that the E-NPF distribution gives a better ﬁt over its competing models such as the Kumaraswamy-PF, Weibull-PF, generalized-PF, Kumaraswamy, and beta distributions. θ = 0.5), S-II( η = 1.1, ζ = 1.5, θ = 1.5), S-III( η = 2.1, ζ = 3.5, θ = 0.5), S-IV( η = 1.1, ζ = 1.5, θ = 2.5), and S-V( η = 4.1, ζ = 1.5, θ = 3.9). The results indicate that the moments and variance decrease gradually, though skewness falls between 0 and 1, and kurtosis might be negative as per the different combinations of the parameters.


Introduction
Over the past three decades, the promising attention of researchers towards the development of new generalized models has increased to explore the hidden characters of baseline models. New generalized models open new horizons to address real-world problems and to provide an adequate fit to the complex and asymmetric random phenomena. Hence, various models have been constructed and studied in the literature. One of the simplest and most handy lifetime models induced in the statistical literature is the Lehmann Type I (L-I) and Type II models [1]. In the literature, the L-I model is most often discussed in favor of the power function (PF) distribution. The L-I model is simply the exponentiation of any baseline model, and it can be specified by the following cumulative distribution function (CDF): F(x; α, ξ) = G α (x; ξ), α > 0, x ∈ R.
Gupta et al. [2] utilized the L-I approach to define a generalized form of the exponential distribution. On the other hand, Cordeiro et al. [3] proposed a dual transformation to the , and its associated PDF reduces to f (x; α, β) = β(1+α) (1−x) , where x ∈ (0, 1), −1 < α< ∞, β >0 are the scale and shape parameters, respectively. Mathematics 2021, 9, 2024 3 of 32 To this end, we define the E-NPF distribution that can model asymmetric and bathtubshaped failure rate phenomena. For this, following Gupta et al. [2], we add up a shape parameter to the baseline NPF model. The newly added parameter advances the tail weight of the E-NPF density and starts providing a consistently better fit over its competitors.
A random variable X is said to follow the E-NPF distribution, say X ∼E-NPF(η, ζ, θ), if its CDF takes the form and its corresponding PDF reduces to where x ∈ (0, 1), −1 < η < ∞, is a scale and ζ, θ > 0 are the two shape parameters, respectively. The two special models of the E-NPF distribution are listed in Table 1.

Shape
Several shapes for the E-NPF density and failure rate functions are displayed in Figures 1 and 2 for various choices of the parameters. Note that the possible shapes of the PDF corresponding to the parameter η, which controls the scale of the distribution, along with the two shape parameters ζ and θ, which control the shapes of the distribution, including increasing, symmetric, upside-down bathtub, decreasing, and J shapes. These shapes are presented in Figure 1a,b. Furthermore, Figure 2a,b presents the hazard rate function (HRF) shapes, including increasing, U, and bathtub shapes. These flexible HRF shapes are suitable for both the monotonic and non-monotonic hazard rate behaviors, which are most likely to appear in real-time situations (see Pu et al. [27] and Oluyede et al. [28]). Such kinds of shapes are often observed in non-stationary lifetime phenomena.

Linear Representation
Linear combination provides a simple way to explore the mathematical properties of the model. For this reason, binomial expansion is utilized. It is given as (1) and (2) are given, respectively, by

Linear Representation
Linear combination provides a simple way to explore the mathematical properties of the model. For this reason, binomial expansion is utilized. It is given as Infinite linear combinations of CDF and PDF in Equations (1) and (2) are given, respectively, by and f E−NPF (x; η, ζ, θ) = ζθ(1 + η) where The expression in Equation (4) will be adopted in the forthcoming calculations of several mathematical properties of the E-NPF distribution.

Reliability Characteristics of the E-NPF Model
Analyzing and predicting the lifetime of a component are important roles in reliability engineering. Hence, some useful and well-established reliability measures are accessible in the literature. The survival function of X, which represents the probability that a component will survive till time x, takes the form In reliability theory, the function that measures the failure rate of a component in a particular period t is also referred to as the force of mortality or HRF. The HRF of the E-NPF distribution has the form It is well known that most mechanical parts/components of some systems follow a bathtub-shaped hazard rate phenomenon. The cumulative hazard rate (CHR) function is expressed by H(x) = − log(S(x)). The CHR function of X is given by The reverse HRF (RHRF) is expressed by W(x) = f (x)/F(x). The RHRF of the E-NPF model takes the form The Mills ratio (MR) is expressed by M(x) = S(x)/ f (x). The MR of X reduces to The odd function is expressed by O(x) = F(x)/S(x) and it is defined, for the E-NPF model, by

Quantile Function, Skewness, and Kurtosis
Hyndman and Fan [29] introduced the concept of quantile function (QF). The q-th QF of the E-NPF distribution can be adapted to generate random numbers and is obtained by inverting its CDF (1). The QF is defined by q = F x q = P X ≤ x q , q ∈ (0, 1). Then, the QF of X follows as To derive the 1st quartile, median, and 3rd quartile of X, one may place q = 0.25, 0.5, and 0.75, respectively, in Equation (14).
The Bowley [30] skewness, say B, and Moors [31] kurtosis, say M, of the E-NPF distribution can be calculated using Equation (14) with the following two formulas: These descriptive measures depend on quartiles and octiles and can provide more robust estimates than the traditional measures of skewness and kurtosis. Additionally, B and M are less sensitive to the outliers and work more effectively for the deficient moment distributions. Some possible shapes of the skewness and kurtosis for various choices of η, ζ and θ are presented in Figure 3.

Moments and Associated Measures
Moments have a useful role in distribution theory, to address the significant characteristics of a probability distribution such as mean, variance, skewness, and kurtosis. Theorem 1. If ~E-NPF( , , ), with , , > 0, then the r-th ordinary moment ( ) of is given by

Moments and Associated Measures
Moments have a useful role in distribution theory, to address the significant characteristics of a probability distribution such as mean, variance, skewness, and kurtosis. Theorem 1. If X ∼E-NPF(η, ζ, θ), with η, ζ, θ > 0, then the r-th ordinary moment (µ r ) of X is given by Proof. Using Equation (2), µ r can be written as After some algebra, we get Simple computations on the last expression lead to the final form of µ r as follows: The moment formula (15) is supportive in the development of some useful statistical measures. For instance, the mean of X follows with r = 1 in Equation (15). The 2nd, 3rd, 4th, and higher-order moments of X are formulated by replacing r = 2, 3, and 4 in Equation (15), respectively. Additionally, the Fisher index (F.I. = (Var(X)/E(X))) plays a significant role in discussing the variability in X. The negative moment of X is simply derived through substituting r by -w in (15). Moreover, one can calculate the skewness (τ 1 = µ 2 3 /µ 3 2 ) and kurtosis (τ 2 = µ 4 /µ 2 2 ) of X by integrating Equation (15). Well-established relationships between the central moments (µ s ) and cumulants (K s ) of X may easily be defined through µ r . The moment generating (MG) function of X, M X (t), is defined by The MG function of X follows as Table 2 displays some numerical results of the first four moments µ 1 , µ 2 , µ 3 , µ 4 , variance ε 2 , τ 1 , and τ 2 for some choices of η, ζ and θ. The results are presented for S-I(η = 0.1, ζ = 1.2, θ = 0.5), S-II(η = 1.1, ζ = 1.5, θ = 1.5), S-III(η = 2.1, ζ = 3.5, θ = 0.5), S-IV(η = 1.1, ζ = 1.5, θ = 2.5), and S-V(η = 4.1, ζ = 1.5, θ = 3.9). The results indicate that the moments and variance decrease gradually, though skewness falls between 0 and 1, and kurtosis might be negative as per the different combinations of the parameters.
The conditional survivor/residual life function is the probability that the life of a component, say x, will survive in an additional interval at t.
It is given by The residual life function of X has the form Furthermore, the reverse residual life is obtained by S R(t) (t|x) = S(x−t) S(t) . It is derived for X by the formula

Entropy
In general, entropy is defined as the system's disorderedness, uncertainty, or a measure of entanglement.
Hence, the simplified form of H φ (X) follows as where

Stress-Strength Reliability
Let X 1 and X 2 be defined to discuss the strength and stress of a component, respectively. The reliability, say R, of X is defined (for X 2 < X 1 ) by R = P(X 2 < X 1 ). Theorem 2. Let X 1 ∼E-NPF(η, ζ, θ 1 ) and X 2 ∼E-NPF(η, ζ, θ 2 ) be independent random variables, then the reliability R reduces to Proof. The reliability R is given by Hence, R reduces to On simplification, Equation (21) takes the form After simple computations, R can be expressed in terms of θ 1 and θ 2 as The last expression illustrates that the reliability R of the E-NPF distribution is an increasing function of θ 1 and a decreasing function of θ 2 .

Stochastic Ordering
Over the past 40 years, the concept of stochastic ordering has engaged scientists and practitioners and is quite useful in economics, reliability theory, queuing theory, insurance, and ecology, among other fields of science.
Let X 1 and X 2 be the two continuous, nonnegative, and univariate random variables, with Q 1 (X) and Q 2 (X) being the CDFs with corresponding PDFs q 1 (x) and q 2 (x), respectively. The random variable X 1 is smaller than X 2 under the following constraints: The four stochastic orders are connected to each other based on the following implications [33].
Proof. Using Equation (2), we have and Hence, we can write By taking derivative w.r.t x on both sides, we get d dx Then, d dx is a decreasing function for all θ 1 < θ 2 .

Order Statistics
Order statistics (OS) and their moments are considered noteworthy measures in quality control, reliability analysis, and life testing. Let X 1 , X 2 , . . . , X n be a random sample of size n that follows the E-NPF distribution and X (1:n) < X (2:n) < X (3:n) < . . . < X (n:n) be the corresponding OS.
The PDF of X (i:n) is By replacing (1) and (2) in the above equation, the PDF of X (i) takes the form The above equation is adopted to compute the w-th moment OS of the E-NPF distribution. Furthermore, the minimum and maximum OS of X follow directly from the above equation with i = 1, 2, respectively. The w-th moment OS, E X w OS , of X follows as The CDF of X (i:n) is defined (for i = 1, 2, 3, . . . , n) by For the E-NPF model, the CDF of X (i:n) reduces to

Estimation under Complete Samples
In this section, we estimate the E-NPF parameters η, ζ, and θ using eight frequentist approaches under complete samples.

Maximum Likelihood Estimation under Complete Samples
In this section, we estimate the parameters of the E-NPF distribution using the maximum likelihood (ML) method. Let X 1 , X 2 , . . . , X n be a random sample of size n from the E-NPF distribution, then the likelihood function of X, L (Θ), takes the form The log-likelihood function, (Θ) = log L(Θ), of X is By taking partial derivatives of Equation (24), we get Mathematics 2021, 9,2024 12 of 32 The ML estimators (MLE)η MLE ,ζ MLE andθ MLE of the E-NPF distribution can be obtained by maximizing Equation (24) or by solving the above equations simultaneously. These equations cannot provide analytical solutions for the MLEs and the optimum values of η, ζ, and θ. Consequently, the Newton-Raphson-type algorithm is an appropriate way in the support of the MLE. This can be done by using different programs, namely R (optim function), Mathematica, or SAS (PROC NLMIXED), or by solving the nonlinear likelihood equations determined by differentiating (Θ).

Maximum Product of Spacing
The maximum product of spacing (MPS) method, as an approximation to the Kullback-Leibler information measure, is a good alternative to the MLE method.
η, ζ, θ , for i = 1, 2, . . . , n + 1, be the uniform spacing of a random sample from the E-NPF distribution, where with respect to η, ζ, and θ, or, equivalently, by maximizing the logarithm of the geometric mean of sample spacings.
The MPSE of the E-NPF parameters can be obtained by solving the nonlinear equations defined by (25).

The Cramér-von Mises Estimators
The Cramér-von Mises estimators (CVME) as a type of minimum distance (MD) estimator have less bias than the other MD estimators. The CVME are obtained based on the difference between the estimates of the CDF and the empirical distribution function. The CVME of the E-NPF parameters,η CV ME ,ζ CV ME , andθ CV ME , can be obtained by minimizing with respect to η, ζ, and θ. Furthermore, the CVME follow by solving the nonlinear equations where ∆ 1 (·|η, ζ, θ), ∆ 2 (·|η, ζ, θ), and ∆ 3 (·|η, ζ, θ) are provided in (25).

Method of Percentile Estimation
Then, the percentile estimators (PCE) of the parameters of the E-NPF distribution are obtained by minimizing the following function: with respect to η, ζ, and θ, where x q is the QF of the E-NPF distribution defined in (14).

Estimation under Type II Censored Samples
In this section, we estimate the E-NPF parameters η, ζ, and θ under complete samples using the MLE method.
Let x 1 , x 2 , . . . , x n be a random sample of size n from the E-NPF distribution; if in the type II censoring scheme we observe only the first r order statistics, then the likelihood function of X, L * (Θ) takes the form f (x i:r:n )[1 − F(x r:r:n )] n−r , x 1:r:n ≤ x 2:r:n ≤ · · · ≤ x r:r:n , where C is a constant that does not depend on the parameters and x 1:r:n , x 2:r:n , . . ., x r:r:n are the censored data. The log-likelihood function without the constant term can be written as follows: By taking partial derivatives of Equation (26), we get The type II censored ML estimators (CMLE)η CMLE , , and of the E-NPF distribution can be obtained by maximizing Equation (26) or by solving the above equations simultaneously by any numerical method as in Section 4.1.1.

Simulation Study
In this section, we perform the simulation study of the E-NPF distribution for complete samples using eight estimation methods and under the type II censored samples for the MLE method.

Simulation Results under Complete Samples
In this section, the performance of the estimates is discussed by the following algorithm.
Step 4: Results of the average of absolute value of biases ( Bias Θ ), and the average of mean relative errors (MRE),      We used R software (version 4.1.0) [34]. Furthermore, Tables 3-6 show the rank of each of the estimators among all the estimators in each row, which is the superscript indicators, and the ∑ Ranks, which is the partial sum of the ranks for each column in a certain sample size. From the results of the 36 combinations, we observe that all the estimation methods show the property of consistency for all parameter combinations, except the Cramér-von Mises method at the combinations: Θ = (η = 1.5, ζ = 0.5, θ = 2.75) T and Θ = (η = 4.00, ζ = 0.50, θ = 2.75) T for the parameter η.
Table S1 (in Supplementary Materials) shows the partial and overall rank of the estimators. From Table S1, and for the parametric values, we can conclude that the AD method outperforms all the other methods (its overall score of 236). Therefore, we confirm the superiority of ADE for the E-NPF distribution.
For visual illustration, we display the MSE of η, ζ, and θ graphically to show that the MSE decrease with an increase in n as expected. Figures 4 and 5 show the MSE of the three parameters based on the values in Tables 5 and 6, respectively. The plots illustrate that the MSE of the parameters decrease gradually with the increase in n. methods show the property of consistency for all parameter combinations, except the Cramér-von Mises method at the combinations: = ( = 1.5, = 0.5, = 2.75) and = ( = 4.00, = 0.50, = 2.75) for the parameter .
Table S1 (in Supplementary Materials) shows the partial and overall rank of the estimators. From Table S1, and for the parametric values, we can conclude that the AD method outperforms all the other methods (its overall score of 236). Therefore, we confirm the superiority of ADE for the E-NPF distribution.
For visual illustration, we display the MSE of , , and graphically to show that the MSE decrease with an increase in as expected. Figures 4 and 5 show the MSE of the three parameters based on the values in Tables 5 and 6, respectively. The plots illustrate that the MSE of the parameters decrease gradually with the increase in .  Table 5 (a-c).
(a) (b) (c) Figure 5. Plots of the mean square errors (MSE) of the E-NPF parameters for the eight methods of estimation using the values in Table 6 (a-c).

Simulation Results under Type II Censored Samples
In this section, the performance of the estimates is explored. We consider that the random sample , , … of sizes = 20, 40, 60, 80, 100, and 150 were generated from the QF in Equation (14), and the values of are chosen to be 80% of . We choose different values for the parameters = −0.50, 0.75, 4.00 , = 0.50, 3.00 and = 0.40, 1.60 and replicate the process = 5000 times.  Table 5 (a-c).   Supplementary Materials) shows the partial and overall rank of the estimators. From Table S1, and for the parametric values, we can conclude that the AD method outperforms all the other methods (its overall score of 236). Therefore, we confirm the superiority of ADE for the E-NPF distribution.
For visual illustration, we display the MSE of , , and graphically to show that the MSE decrease with an increase in as expected. Figures 4 and 5 show the MSE of the three parameters based on the values in Tables 5 and 6, respectively. The plots illustrate that the MSE of the parameters decrease gradually with the increase in .  Table 5 (a-c).
(a) (b) (c) Figure 5. Plots of the mean square errors (MSE) of the E-NPF parameters for the eight methods of estimation using the values in Table 6 (a-c).

Simulation Results under Type II Censored Samples
In this section, the performance of the estimates is explored. We consider that the random sample , , … of sizes = 20, 40, 60, 80, 100, and 150 were generated from the QF in Equation (14), and the values of are chosen to be 80% of . We choose different values for the parameters = −0.50, 0.75, 4.00 , = 0.50, 3.00 and = 0.40, 1.60 and replicate the process = 5000 times.  Table 6 (a-c).

Simulation Results under Type II Censored Samples
In this section, the performance of the estimates is explored. We consider that the random sample x 1 , x 2 , . . . x n of sizes n = 20, 40, 60, 80, 100, and 150 were generated from the QF in Equation (14), and the values of r are chosen to be 80% of n. We choose different values for the parameters η = {−0.50, 0.75, 4.00}, ζ = {0.50, 3.00} and θ = {0.40, 1.60} and replicate the process N = 5000 times.
In each setting, we obtain the mean of the estimates (ME), ME = 1 N ∑ N i=1 Θ and the corresponding |Bias(Θ)|, MSE and MRE. These results are displayed in Table 7.       Table 7, it is observed that the MSEs decrease as the sample size increases in all the cases under the complete sample, as we discussed in Section 5.1. In the case of type II censored sample, as the number of failures r increases, the MSE decrease in all the cases. Furthermore, the MEs tend to the true parameter values as the sample size increases. Thus, we can say that the MLEs of the parameters η, ζ, and θ under the two schemes are asymptotically unbiased and consistent.

Results and Discussion
In this section, we report the application of the E-NPF distribution in applied sciences by analyzing three suitable lifetime datasets. The first dataset relates to oceanography, and it represents the synthetic aperture radar (SAR) image for modeling oil slick visibility in the ocean [35]. The second dataset from the reliability engineering field consists of 20 observations about failure times of mechanical components [36]. The third dataset is discussed by Ahsanullah et al. [8], and it contains measurements on petroleum rock samples. It is worth mentioning that the observations of the three datasets are bounded in the interval (0, 1). The E-NPF distribution is compared with its competing models, which are present in Table 8, based on some criteria such as the Akaike information criterion (AIC), Cramér-von Mises (CM), Anderson-Darling (AD), and Kolmogorov-Smirnov (KS) test with its p-value (KS p-value) statistics. However, using the statistical software R under the package Adequacy Model [37] is considered. In addition, Some choices of descriptive statistics are presented in Table 9. Table 8. List of some competitive models along with their CDFs.

Competitive Models
The CDF of Competitive Models Author(s) Mood et al. [38] Topp-Leone (TL) Weibull-PF (WPF) Tahir et al. [11] Kumaraswamy-PF (Kum-PF) Mustapha Type II (MT-II) Iqbal et al. [26] Exponentiated-PF (EPF) Tables 10-12 illustrate the estimates of the parameters, standard errors (S.E.), and goodness-of-fit statistics. For the three datasets, the E-NPF distribution has the lowest values of AIC, CM, DA, and KS measures and the largest KS p-value among all models studied. That is, the E-NPF distribution provides a superior fit than other models for the three datasets.  Furthermore, the empirically fitted plots of the PDF, CDF, Kaplan-Meier survival, and probability-probability for the three datasets are presented in Figure 6, Figure 7 and Figure 8, respectively. These plots confirm the close fit of the E-NPF distribution to the three datasets. The numerical results in this section are calculated using the statistical software R under the package Adequacy Model [37]. Furthermore, the empirically fitted plots of the PDF, CDF, Kaplan-Meier surv and probability-probability for the three datasets are presented in Figures 6-8, res tively. These plots confirm the close fit of the E-NPF distribution to the three datasets. numerical results in this section are calculated using the statistical software R under package Adequacy Model [37].

Conclusions
In this article, we develop a bounded lifetime model that exhibits a bathtub-shap failure rate. The proposed distribution is called the exponentiated new power function NPF) distribution. Numerous mathematical and reliability measures are derived in plicit expressions. Some classical methods of estimation are adopted to estimate the NPF parameters. Moreover, the maximum likelihood is utilized to estimate the E-N parameters under the type II censored samples. Two Monte Carlo simulation studies

Conclusions
In this article, we develop a bounded lifetime model that exhibits a bathtub-shaped failure rate. The proposed distribution is called the exponentiated new power function (E-NPF) distribution. Numerous mathematical and reliability measures are derived in explicit expressions. Some classical methods of estimation are adopted to estimate the E-NPF parameters. Moreover, the maximum likelihood is utilized to estimate the E-NPF parameters under the type II censored samples. Two Monte Carlo simulation studies are carried out to investigate the asymptotic performances of the estimates under complete and type II censored samples. The most efficient and consistent results of the E-NPF distribution are explored by modeling three real-life datasets related to the fields of oceanography, reliability engineering, and petroleum engineering. It is hoped that in the future the E-NPF distribution will be quite helpful for researchers and will be considered a better choice against the baseline model.

Supplementary Materials:
The following is available online at https://www.mdpi.com/article/ 10.3390/math9172024/s1: Table S1: Partial and overall ranks of all estimation methods for various combinations of Θ.