A robust bathtub-shaped failure time model for a two-component system with applications to complete and censored reliability data

ABSTRACT This article proposes a flexible additive model that adequately describes complex reliability and survival data. The proposed methodology, referred to as the flexible exponential power-Gompertz (FEPG4) distribution, is able to characterize the behavior of a complex system whose failure times have bathtub-shaped, with clear burn-in and wear-out change points and a low, yet lengthy, flat middle segment as its underlying failure rate distribution. We discuss some properties of the model. Parameter inferences are proposed under maximum likelihood and Bayesian techniques. We determine the Bayes estimators of the FEPG4 parameters and used Hamiltonian Monte Carlo for posterior simulations. Extensive simulation experiments are performed to validate the proposed estimators. For assessing the potential of the FEPG4, the model is compared with other recent bathtub distributions constructed via the same approach on devices’ failure and running times (censored and non-censored case) and failure times of some devices, each with the bathtub failure rate. Seven parametric and non-parametric selection criteria and other supporting plots are utilized for comparison purposes. Findings indicate that the FEPG4 model might be the best alternative for representing device failure times, particularly when the bathtub-shaped failure rate of the available data clearly illustrates its three phases.


Introduction
In reliability and survival research, the failure rate (FR) function fundamentally characterizes the failure time data (Shakhatreh et al., 2019;Kosky et al., 2021).It may be applied to describe the lifetime distribution for life occurrence, which can be utilized to explain the quality and reliability of systems.As the number of complex systems increases in reliability engineering, several failure time data display non-monotone FR patterns, most notably the bathtub-shaped FR, which corresponds to the three separate phases of the system.The FR drops during the early life (burn-in) phase, remains virtually constant throughout the mid-life (useful life) phase, and rises during the wear-out period.Collins and Warr, (2019) characterized the bathtub-shaped curve as the reliability engineering rule of thumb for complex components and systems.There are many interesting examples of failure time data with bathtub-shaped FR patterns and other non-monotone FR forms.A few examples are the times to failure of electronic devices (EDs) (Aarset, 1987, Kosznik-biernacka, 2007;Méndez-gonzález et al., 2022;Méndez-González et al., 2019), death times of psychiatric patients (Klein and Moeschberger, 2003), early cable-joint failure times (Tang et al., 2015), and time to failure of sugarcane harvesting machines (Prataviera et al., 2018).In previous studies, considerable attempts have been made to model failure time data using various non-monotone FRs, and many lifetime models are developed for investigating the reliability of complex systems (see, for example Fan andHsu, 2015 andShakhatreh et al., 2019).
In particular, numerous Weibull and non-Weibull generalized models have been proposed for modeling lifetime data with non-monotone FRs from various viewpoints, notably reliability engineering.Few of these are Zhang and Xie, 2011, Chaubey and Zhang, 2015, Saboor et al., 2016, Bakouch and Abd El-Bar, 2017, Shakil et al., 2018, Lu and Chiang, 2018, Ahmad and Ghazal, 2020), Muhammad and Liu, 2021, Méndez-gonzález et al., 2022, Wang et al., 2022, and Muhammad et al., 2022).However, it was learned that a significant proportion of these distributions have bathtub FR shapes but lack the FR's relatively constant phase (Shakhatreh et al., 2019).This flat zone is crucial, or may be the most critical phase for reliability modeling (Kosky et al., 2021), since it describes the useful life span of the component or system.Therefore, it is vital to have a model that can accurately represent this flat region.
A considerable focus has been placed on modeling series systems reliability among the other multi-component systems constructed with various configurations.Each component of the series system has an underlying FR that is seldom specified directly.Still, system operators can make implicit assessments of the FR using the lifetime data of the components, recorded during system operation, or or by performing accelerated life test (ALT) on components separately (Ikonen et al., 2020).Consequently, there is a need for a model that can adequately account for the individual component FR in determining the system FR.
Various models are presented for analyzing the reliability of a series system with two independent components whose failure times reflect monotone or non-monotonic FR.In this context, each component follows a separate lifetime model.The additive modified Weibull (AMW) (He et al., 2016), log-normal modified Weibull (LNMW) (Shakhatreh et al., 2019), improved new modified Weibull (INMW) (Thach and Briš, 2020), exponentiated additive Weibull (EAddW) (Ahmad and Ghazal, 2020), additive Chen Weibull (ACW) (Thach and Briš, 2021), flexible additive Chen-Gompertz (FACG) (Abba et al., 2022), and additive Perks (Méndez-González et al., 2022) distributions are some recent related models in this regard.These models are created by combining two classical models under the assumption that the failure times of the first component match one of the models and the failure times of the second component suit the other model.The technique yields adaptable distributions for complex series systems' reliability, especially for data sets with bathtubshaped FRs (Shakhatreh et al., 2019;Abba et al., 2022).The AMW, LNMW, INMW, EAddW, ACW, and FACG models were shown to provide good fits to the well-known failure times of EDs by Aarset (Aarset, 1987) or the device failure and running times reported by Meeker and Escobar (Meeker and Escobar, 1998), each having a bathtub-shaped FR.However, the majority of these studies largely concentrate on non-censored scenarios, with the FACG being the only model recently demonstrated for fitting right-censored data.Furthermore, AMW, EAddW, and ACW distributions were fitted to early cable-joint failures (Tang et al., 2015) and death times of psychiatric patients (Klein and Moeschberger, 2003), each with bathtub FR.The findings favored the FACG model and suggested that the fitted FR curves of the AMW, EAddW, and ACW models were unable to identify the bathtub-shaped FR's long-flat phase and the wear-out change point for the early cable failure data (Abba et al., 2022).
Hence, the possibility of current models to identify the low and yet long-flat segment of the bathtub-shaped FRs and its two change points for various series systems failure times remains a discoverable task.Thus, the need to search for more suitable lifetime models capable of effectively describing the system's reliability.To our understanding, the ACW (Tarvirdizade andAhmadpour, 2021, Thach andBriš, 2021) and FACG (Abba et al., 2022) models represent the most recent Chen distribution (Chen, 2000) extensions defined based on complex system interpretation.However, because the Chen model component in these models lacks a scale parameter, it may be problematic to monitor the spread of the distribution.Besides that, there has been no study combining exponential power and Gompertz models to develop additive distribution.
This work intends to introduce an alternate model to the current literature by establishing a new lifetime model for modeling the reliability of a two-component series system or devices with one or two failure modes whose failure times conform to bathtub-shaped FR.The proposed model termed the flexible exponential power-Gompertz (FEPG4) distribution, is developed by combining the exponential power (EP) distribution (Smith and Bain, 1975) with a special variant of the Gompertz (G) distribution.The FEPG4 is defined based on a series system consisting of two non-dependent components or failure times having two independent failure modes.Each failure is attributed to one of the two failure types, and the failure time is supposed to follow either the EP or G distribution.The FEPG4 distribution could be used instead of ACW and FACG models, particularly when the role of the scale parameter from each model component is required.The study used Gompertz to strengthen the FEPG4's ability to adequately describe the FR's low and long-flat useful period segment.It also helps the proposed model to provide a clear change point at the wear-out phase of the bathtub-shaped FR function (FRF) when modeling failure time data with such feature.To accommodate the lack of scale parameter in the Chen model component of the ACW and FACG models, we used the exponential power model as the other component of the proposed distribution.
The main motivations of this article include: (i) The new FEPG4 distribution, which is a combination of the EP and G models, is substantially more flexible than several current hybrid bathtub distributions.
(ii) The FEPG4 model includes a scale parameter in each of the model's components, making it a viable alternative to the ACW and FACG models, where the role of scale parameters is required.
(iii) For the FEPG4 parameters, we propose the Hamiltonian Monte Carlo (HMC) technique and the maximum likelihood method, with HMC being insensitive to correlated parameters and allowing for more robust parameter estimation.
In addition to the newly presented model, we used the HMC to simulate data from the non-closed form posterior distribution derived from the Bayesian estimation technique.The HMC algorithm is a potential tool for sophisticated models with several parameters, as they give more precise findings and adjustable settings (Thach and Briš, 2021).We also proposed maximum likelihood estimators for the model's parameters.The appropriateness of the proposed FEPG4 model is demonstrated by fitting two complete and one right-censored data with bathtub-shaped failure rates.The data sets used include the failure and running times of 30 devices (Meeker and Escobar, 1998) (censored and uncensored) and failure times of 30 EDs (Kosznik-biernacka, 2007).Several goodness-of-fit statistics, including the Akaike information criterion (AIC), bias-corrected Akaike information criterion (AICc), Bayesian information criterion (BIC), Kolmogorov Smirnov (KS), Anderson-DarlingAnderson-Darling (AD), and Cramér-von Mises (CVM), acknowledged that the FEPG4 model offers a better fit for each of the three illustrations among other recent additive models.The examples also highlight how the FRF of the FEPG4 model uncovers the low and yet long-flat segment of the bathtub-shaped curves, particularly in the first case.
Beyond the direct fitting of failure time distributions to various lifetime data, there are several facets of reliability in which lifetime distributions are used to build diverse reliability models.Some highlights of these areas and potential uses of the proposed FEPG4 model in this context are provided below.Silva et al., (2010) and Méndez-González et al., (2016) determined that when two or more types of failures are observed during the ALT of EDs or when the EDs are subjected to time-varying stress, the FR of the EDs lifetimes is non-monotone.This development has resulted in the growing application of bathtub and other non-monotone FR distributions in constructing several reliability models.For example, the exponentiated Weibull (Mudholkar and Srivastava, 1993), Beta-Weibull (Famoye et al., 2005), Beta generalized Weibull (Singla et al., 2012), alpha power transformed Weibull (Dey et al., 2017), and alpha power transformation exponentiated Perks (Méndez-González et al., 2022) distributions were utilized in conjunction with inverse power law (IPL) model for developing a life-stress model (LSM) (Ali et al., 2019;Méndez-González et al., 2022;Méndez-gonzález et al., 2022;Méndez-González et al., 2017;Méndez-González et al., 2019).Even though the LSM created using these distributions provided a good fit for the case study data.Nevertheless, with the three-phase bathtub-shaped FRF of the proposed FEPG4 model and its good direct fitting performance, the new methodology, if implemented, could provide a more accurate description of the behavior of EDs in the implementations of LSM.A location-scale regression model based on the FEPG4 distribution may also be suited to analyze and predict early cable-joint failures linked with power cables, given that the failure times of power system assets are believed to follow the bathtub-shaped FR (Klutke et al., 2003).This results from fitting similar bathtub-shaped models and Weibull distribution by Abba et al., (2022) and Wang et al., (2022) on early cable-joint failure times (Tang et al., 2015).The study demonstrates that the widely used Weibull distribution in this reliability aspect is insufficient to describe the bathtub shape FR of the cable-joint failure data by a wide margin.Thus, the suggested FEPG4 and similar models may provide more accurate cable-joint failure predictions than the standard Weibull model.Similar to the recent work of Wang et al (Wang et al., (2022), the FEPG4 or its sub-model can also be a contender for establishing a multi-component stress--strength model (SSM) when the stress and strength variables adhere to bathtub-shaped FRs.This is crucial, as the majority of previous studies on SSM employed a lifetime model with less capacity to accommodate data with a distinct three-phase bathtub-shaped FR.
We introduce the proposed model and some of its properties in Section 2. Parameter inferences are provided in Section 3. In Section 4, we validate the proposed estimation approaches using simulation.In Section 5, we demonstrate the appropriateness of the FEPG4 model for modeling reliability data.Section 6 serves as the conclusion of the study.

The FEPG4 model and some of its properties
This section begins with a brief overview of the exponential power (EP) distribution and the special case of the Gompertz (G) distribution, followed by the definition of the proposed model and some of its basic functions.We offer a comprehensive assessment of the model's FRF.At the end of the section, we present some properties of the proposed FEPG4 model, including the quantile function, incomplete moments, mean residual life (MRL), mean time to failure (MTTF), and stress--strength reliability.Proofs of the incomplete moments, MRL, MTTF, and stress--strength reliability, along with other properties of the model, are detailed in the supplementary materials.Smith and Bain, (1975) established the EP distribution as a life-testing model.The model is applicable in situations when failure time data exhibit either a U-shaped bathtub or an exponentially increasing FR.In contrast to the Chen distribution (Chen, 2000), which is a modified version of the model, the EP distribution contains a scale parameter, making it more applicable in reliability context.Numerous extensions for modeling various problems have been proposed (see, for example Shakil et al., (2018) and Korkmaz et al., (2021).
The G distribution is a well-known model with an increasing FRF that has been widely investigated for modeling a variety of data (see, for example El-Gohary et al., (2013) and Bakouch and Abd El-Bar, (2017).It has been identified as a potential candidate for constructing hybrid bathtub distributions (Abba et al., 2022;He et al., 2016).On the other hand, numerous studies have modified the distribution to accommodate non-monotone FR data (e.g Bakouch and Abd El-Bar, (2017) and Shama et al., (2022).
Assume that the EP random variable Z 1 has the reliability function (REF) given as The FEPG4 random variable T is constructed as the minimum between the two independent random variables Z 1 and Z 2 , that is, and therefore, T,FEPG4ðμÞ, where μ ¼ ðα; τ; θ; υÞ 0 is the FEPG4 parameter vector.Then, we define the REF of T as: where τ > 0 is the FEPG4 shape parameter and α � 0; θ > 0, and υ � 0 are the scale parameters of the model.From Eq. (1), the FEPG4 cumulative distribution function (CDF) is given as Whereas the associated probability density function (PDF) of the proposed model is obtained as We define the FRF of the FEPG4 distribution from Eq.s (1)-( 2) as The FEPG4 model is applied to two independent components in a series system, with one component's lifetime following the EP distribution and the other component's lifetime following the G model.It can also be used for a component that is mainly affected by dual failure mechanisms, one of which follows the EP distribution and the other the G distribution.In this instance, the component might fail due to either the first or second failure mode.
Figure 1 shows many potential shapes of FEPG4's density function.Depending on the parameter values, we can see, among other PDF patterns, a decreasing, unimodal, and decreasing-increasingdecreasing (modified bathtub) shapes.

Three-parameter flexible exponential power-Gompertz (FEPG3) distribution
When the scale parameter α ¼ 1, we obtained an interesting special case of Eq. ( 2); that is, the three-parameter flexible exponential power-Gompertz (FEPG3) distribution, represented by T,FEPG3ðτ; θ; υÞ.In terms of FRF and real-life applications, the FEPG3 model is equally versatile as the FEPG4 distribution when α is very close to unity (α ! 1).Its FRF is expressed as

FRF of the FEPG4 model
The FRF is a crucial feature of a probability distribution that illustrates how the likelihood of failure increases with time or age and is relevant to the vast majority of applications, especially when dealing with failure time data.Utilizing prior knowledge of the shape of the FRF facilitates the selection of a suitable model.When the attributes that may affect the lifetime of a system or individual vary over time, it is essential to analyze the problem using the FRF.Similar to other generalized distributions, it is easy to establish the pattern of the FEPG4's FRF, which is dependent on τ, the only shape parameter of the model.It is straightforward to deduce from Eq. ( 4) that the FRF increases when τ � 1 and bathtub when τ < 1.When h 0 ðt � Þ ¼ 0, the FRF appears bathtub-shaped for τ < 1, if and only if t ¼ t � is the root of the following equation: It can be verified that the FRF is decreasing if h 0 ðtÞ < 0 for t < t � and is increasing if h 0 ðtÞ > 0 for t > t � .Various shapes of the FRF plots are presented in Figure 2(a-d) for distinct values of parameters.The proposed distribution has increasing FR shape and various bathtub forms for the FRF and is versatile enough to represent complex lifetime data with clear bathtub-shaped FR. Figure 2(c,d) depicts several forms of bathtub FRs with or without a low and yet long-flat phase (useful time) of the FEPG4 distribution.Remarkably, the proposed FEPG4 model exhibits the bathtub-shaped characteristics mentioned in Section 1, and may thus be used to explain burn-in procedures for enhancing system reliability, along with other reliability applications.Further analyses on the bathtub-shaped FR reflect that the shapes of the FRF are only defined by τ.Moreover, α and θ determine the length of the flat region (see, Figure 2(d)), and υ stabilizes the low and yet constant zone as it approaches zero (see, Figure 2(c)).
Typically, a total time on test (TTT) plot is used to determine a data set's FR distribution, thereby guiding researchers in model selection.The reader is referred to (Abba et al., 2022) for a recent reexplanation of the strategy.

Quantile function
The quantile function of the FEPG4 model is determined as a real solution to the following nonlinear equation.Equation ( 6) has no closed-form solution in t q ; hence, we employ a numerical approach to calculate the quantile.At q ¼ 1=2, the solution to Equation ( 6) simplifies to the median.A straightforward approach for simulating data from the FEPG4ðα; τ; θ; υÞ model is shown below.
Step 1 Simulate q i ,Uð0; 1Þ; i ¼ 1; 2; . . .; n, and then Step 2 Take some initial values for μ and then apply Step 1 to determine for t i in

Incomplete moments
Theorem 2.1 The r th incomplete moment of the FEPG4 random variable T is given as , and y dy is the incomplete generalized integro-exponential function (Astorga et al., 2020).Proof.See the supplementary material.It is observed from Theorem 2.1 that, E n s z 1 ; z ¼ z 2 ð Þ converges to a generalized integro-exponential function (Milgram, 1985),

MRL of the FEPG4 model
Theorem 2.2 The MRL for the FEPG4 random variable T is given as 1 ðln yÞ jτ y À 1 e À e À θ y dy.
Proof.See the supplementary material.

Stress--strength reliability
Theorem 2.4 Suppose T and Y are two independent strength and stress random variables, which, respectively, follow FEPG4 and Weibull distributions, with parameter vectors ðα; τ; θ; υÞ 0 and ðλ; τÞ 0 .Then, for the strength T of a component with CDF F 1 ðtÞ in Eq. ( 2), that is subjected to the stress Y with CDF F 2 y ð Þ ¼ 1 À e À λy τ ; y > 0; λ; τ > 0, the stress--strength reliability of the component is given as Proof.See the supplementary material.

Parameter estimation
In this part, we offer maximum likelihood (ML) and Bayesian estimation techniques for the FEPG4 parameters.

Maximum likelihood estimation for uncensored data
For an n-size random sample t 1 ; t 2 ; . . .; t n drawn from the FEPG4 distribution having the parameter vector μ ¼ α; τ; θ; υ ð Þ 0 , the uncensored log-likelihood function of μjt i ; i ¼ 12; . . .; n is constructed from the PDF f ðtÞ (3) as Thus, the associated score function for each FEPG4 parameter from (8) are where and hðt i Þ is the FRF given in Eq. ( 4).Due to the nature of the score functions U α ðμÞ; U τ ðμÞ; U θ ðμÞ; and U υ ðμÞ, we suggest computing the parameter numerically using statistical software.For the optimization, we used the optim and maxLik package (Henningsen and Toomet, 2011) in R statistical software.
Additionally, we can determine the ML estimates (MLEs) of RðtÞ and hðtÞ by applying the invariant property of the MLE (Casella and Berger, 2001).Thus, the MLEs of RðtÞ and hðtÞ are, respectively, given by

Asymptotic confidence intervals
It is conceivable but difficult to derive the closed-form expressions for the exact distributions of the MLEs.Due to this, we may obtain approximate normal confidence intervals for α; τ; θ, and υ when n becomes sufficiently large (n !1), that is, Þa four-dimensional normal distribution with zero mean and I À 1 variance-covariance matrix, where The entries of IðμÞ represent the second partial derivatives of ,ðμÞ, with regard to α; τ; θ; and υ, and are, respectively, given as follows: where and g α ðt i Þ; g τ ðt i Þ; g θ ðt i Þ and g υ ðt i Þ are provided in Section 3.1.Consequently, the estimated variancecovariance matrix may then be evaluated at the MLE of α; τ; Therefore, the 100ð1 À �Þ% asymptotic confidence intervals for α; τ; θ and υ are expressed as where Z� 2 represents the standard normal distribution � th upper percentile.

Maximum likelihood estimation for censored data
In a similar manner, we construct the log-likelihood function of the FEPG4 model for rightcensored data, as follows: Suppose (t i ; γ i ), i ¼ 1; 2; . . .; n is a censored random sample, where t i is a failure or survival time for γ i ¼ 1 and censored time for γ i ¼ 0. The log-likelihood function of the FEPG4 model is constructed from the reliability function RðtÞ (1) and FRF hðtÞ (4) as Thus, the associated score function for each FEPG4 parameter from ( 9) are where g α ðt i Þ; g τ ðt i Þ; g θ ðt i Þ and g υ ðt i Þ are given in Section 3.1.We adopt similar optimization techniques discussed in Section 3.1 for estimating the FEPG4 parameters.

Bayesian estimation
Bayesian inference is a robust method of estimation that has been used to estimate parameters in a variety of domains, particularly when dealing with complex situations.In this study, we developed the Bayesian framework for estimating the parameters of the FEPG4 model.The posterior distribution of μ is given by πðμjDÞ, where the prior πðμÞ reflects the distribution of μ before the data D : t 1 ; t 2 ; ; t n are obtained.The posterior distribution of πðμjDÞ is defined using Bayes theorem as with ð μ LðDjμÞπðμÞdμ representing the marginal distribution of D. The likelihood function of FEPG4 for censored/uncensored data is expressed as , γ i is the censoring indicator as specified in Section 3.1 for censored data, and γ i ¼ 1 for all i's for uncensored data.Let the FEPG4 parameters α; τ; θ, and υ be non-dependent random variables.The choice of the prior is very sensitive, which lies in the availability of prior information on the values of the parameters.However, with no prior belief for the values of these parameters, we assume the random variables α; τ; θ, and υ to each follow nondefinite diffuse prior, that is, the gamma distribution.Additionally, the choices of gamma priors are encouraged following its promising Bayesian results in similar recent studies (Abba and Wang, 2023, Abba et al., 2022, Thach and Briš, 2021).Let πðμ k Þ denote the gamma distribution with parameters ðλ k ; γ k Þ -also called hyper-parameters, with density given by Thus, α,πðαjλ 1 ; γ 1 Þ; τ,πðτjλ 2 ; γ 2 Þ; θ,πðθjα 3 ; γ 3 Þ and υ,πðυjλ 4 ; γ 4 Þ, and then the typical joint prior is the product of the four non-dependent priors The joint posterior density of μjD is defined as 10).The marginal posterior density of α; τ; θ, and υ can be obtained from (11) as From the marginal posterior densities in Eq. ( 12), the Bayesian estimators of α; τ; θ, υ, RðtÞ, and hðtÞ under the square error loss function (SELF) are, respectively, given as With no closed-form of the marginal posterior densities πðαjDÞ; πðτjDÞ; πðθjDÞ; and πðυjDÞ available, we could not analytically compute the Bayesian estimates (BEs) through the posterior means.Thus, we use Markov chain Monte Carlo (MCMC) methods.Here, we adopt the Hamiltonian Monte Carlo (HMC) approach, an MCMC method alternative to the frequently used Metropolis algorithm and Gibbs sampler.Unlike the Metropolis algorithm in highly correlated or high-dimensional parameters, the HMC suppresses the local random walk behavior and is not sensitive to correlated parameters, which makes it much more efficient at exploring posteriors in models with correlated parameters (Stan Development Team, 2018).It also enables faster mixing and is rapidly moving through the target distribution.A recent excellent step-bystep re-explanation of the HMC algorithm can be found in (Abba andWang, 2023, Thach andBriš, 2020).In this paper, we used a user-friendly R interface to Stan (Rstan (Stan Development Team, 2020)) to sample from the posterior distribution using the No-U-Turn sampler (NUTS) (Hoffman and Gelman, 2014).NUTS is a variant of HMC that automatically update the HMC's parameters.Therefore, for a sample μ s ; s ¼ 1; 2; . . .; N f g obtained from πðDjμÞ, the approximate BEs of FEPG4 parameters α; τ; θ, and υ, as well as the reliability function RðtÞ and FRF hðtÞ are computed by where } is the warm-up/burn-in period, representing the number of iterations before the stationary samples are determined.It is recommended in practice to run for b parallel chains (b ¼ 3; 4 or 5) to enable the assessment of the sampler convergence (Stan Development Team, 2018).Thus, the posterior means for b parallel chains are computed as follows: For each chain, we monitored the convergence of the HMC samples using the Gelman and Rubin's (Gelman and Rubin, 1992)'s potential scale reduction factor (PSRF, b R), which is the weighted average of the between-chain and within-chain variances.The HMC sample converges if b R < 1:1.To compute the (1 À δ)100% Bayesian credible intervals (BCIs) of the parameters, expressed as ½μ k; δ 2 ðNÀ sÞ ; μ k;ð1À δ 2 ÞðNÀ sÞ �, we sort the HMC samples in ascending order as μ k;1 � μ k;2 � . . .� μ k;NÀ s ; k ¼ 1; . . .; 4. We can also compute the approximate (1 À δ)100% highest posterior density (HPD) interval of μ k , estimated as the interval that contains the highest posterior density.

Simulation results
In this part, the performance of the estimators derived using maximum likelihood and Bayesian estimation methodologies for estimating the parameters of the proposed distribution are evaluated.For simplicity, we choose the three-parameter FEPG3 model, which is as flexible as the FEPG4 model when α approaches unity.We utilize the following steps for the simulation studies: (i) set initial values for the parameterτ; θ; and υ, and assume the sample size n; (ii) apply Equation ( 6) to simulate a random sample of size n from the FEPG3 distribution; (iii) use the samples drawn from (ii) to optimize the log-likelihood function given in (9) to determine the estimates of τ; θ; and υ; (iv) repeat steps (i), (ii) and (iii) N ¼ 1000 times, and calculate the average estimates for all the parameters, biases and mean square errors (MSEs) of the estimates; and (v) for each case, we also used the N ¼ 1000 samples of size n to compute the Bayes estimates with 1000 iterations under gamma priors and counted out the first 20% as warm-up samples.We carried out all the simulations in R with 10 different sample sizes ranging from 30 to 300.
Performance comparisons of the estimators for each sample size were made based on the bias and MSE obtained from (iv). Figure 3(a-h) provides the pictorial views of the simulation experiments for various parameter settings and different sample sizes.Whereas the numerical results are provided in Tables 1-2.Obviously, the biases and MSEs approach zero as the n increases, which explains the proposed methods' consistency.Similar findings can be made from the numerical results.In addition, the estimates for some parameter settings are shown to have converged more quickly to the true values than others (see, for instance, (Figure 3  (a-d)).

Applications to reliability data
In this section, we demonstrate the breadth of applications of the proposed FEPG4 model and its sub-model, the FEPG3, by comparing its performance on three instances of reliability data to five related models, which are constructed using the same methodology.They include the AMW (He et al., 2016), INMW (Thach and Briš, 2020), EAddW (Ahmad and Ghazal, 2020), ACW (Thach and Briš, 2021), and FACG (Abba et al., 2022) (see, Table 3).The three sets of data include the failure and running times of 30 devices (Meeker and Escobar, 1998) (censored and uncensored) and the failure times of 30 devices reported by Kosznik-biernacka, (2007).Estimation techniques introduced in Section 3 are utilized for all the computations.All the chosen data sets are characterized by bathtub-shaped FR, with the former data serving as a benchmark for evaluating the adequacy of fit for models with bathtub-shaped FRF.We use the well-known measures, including the negative log-likelihood value ( À ,), AIC, AICc, BIC, KS, AD, and CVM to select the appropriate model for fitting the non-censored data.Whereas for the censored data, we only used À ,, AIC, AICc, and BIC criteria, since KS, AD and CVM do not account for censoring information.In addition, the reliability, FRF, and MRL curves, along with box plots for the fitted models, are employed to help in the comparisons.For the Bayesian analysis, we report the BEs and their related 95% BCIs and HPDs, along with other supporting plots for the FEPG4 parameter estimates.

Failure times of thirty devices
The data contains the failure and running times (FRTs) for n ¼ 30 device samples from a large system field tracking study (Meeker and Escobar, 1998).During the process, two failure modes were discovered: electric surge and wear-out, which can be modeled with the proposed FEPG4 and FEPG3 models.We shall call it the Meeker--Escobar data.The data has been studied by numerous authors, including Shakhatreh et al., (2019), Thach and Briš, (2020), and Thach and Briš, (2021), among others.Nevertheless, most authors treated the right-censored values as actual failure times and modeled the data as uncensored.In 2020, Ikonen et al., (2020) used the data in its full scale (as censored data) to demonstrate the use of statistical bathtub failure rate models in selective maintenance optimization.The supplementary material provides the data, and Figure 6(a) shows a scaled TTT-transform plot of the data's bathtub-shaped FR.The scaled TTT-transform plot was obtained using TTT function available in AdequacyModel package (Diniz-Marinho et al., 2016) in R. We analyze both the uncensored and censored cases of the data in Sections 5.1.1 -5.1.2to demonstrate the applicability of the suggested FEPG4 and FEPG3 models over the aforementioned distributions.
Table 3.Some recent bathtub-shaped failure rate models used for comparison.

Model
The displays of the fitted models' reliability functions, FRFs, and MRL on the data are shown in Figure 6(b-d).The reliability curves of the trained models and the KM estimate can be found in Figure 6(b).From this plot, we may deduce that the INMW, ACW, FACG FEPG3, FEPG4-ML, and FEPG4-BE reliability functions all have curves closely resembling the KM step function and compete with one another.The FEPG3, FEPG4-ML, FEPG4-BE, AMW, INMW, FACG, and ACW models' estimated FR curves closely mirrored the step function of the observed data (see, Figure 6(c)).The results are manifested in their FRs' early and wearout stages, which have clear and sharp change points.Although the AMW, INMW, ACW, and FACG have somewhat lower FRs in the early stage, the fitted FEPG3, FEPG4-ML, and FEPG4-BE have the lowest rates of failure in the wear-out phase and reasonably have longer flat FRs over the mid-time period than the other fitted distributions.Besides that, the proposed FEPG4 and FEPG3 models' FR curves are quite comparable to the step function and present the lowest FR values at the minimum failure time ðt 1 ¼ 0:1Þ.Thus, it can lead to a better burn-in strategy.It is worth noting that the FEPG3 and FEPG4 distributions can effectively describe the data.As a result, we believe that all of the graphical results (Figures 4-6) support the numerical results and that either the FEPG4 or FEPG3 models may be utilized to investigate and predict device reliability.

Censored case
Table 6 contains the estimated parameter values (MLEs), À ,, AIC, BIC, and AICc for the FEPG3, FEPG4 and other fitted additive distributions.Compared to the proposed FEPG3 and FEPG4 distributions, the five-parameter INMW model has the smallest values of À , with a nonsignificant margin of less than 0.3, as seen in Table 6.According to the AIC, AICc, and BIC results, the FEPG3 and FEPG4 were ranked first and second, respectively, indicating that they fit the devices censored failure and running time data better than the other trained additive models.
We use the same procedure as in Section 5.1 with hyper-parameter values of ðλ 1 ¼ 0:22; γ 1 ¼ 0:5Þ; ðλ 2 ¼ 0:395; γ 2 ¼ 0:5Þ; ðλ 3 ¼ 200; γ 3 ¼ 0:6Þ; and ðλ 4 ¼ 2:6e À 4 ; γ 4 ¼ 0:1Þ, to simulate the posterior samples of α; τ; θ, and υ, respectively, for the censored Meeker--Escobar data.Figure 7 shows the trace plots and marginal density estimates of the parameters generated by the HMC.The trace plots visualize that the three parallel HMC chains converge for each parameter to the same target distribution.In addition, the PSRF for each BE is b R < 1:1, assuring chain convergence (see, Table 7).The posterior samples produce good BEs under SELF with their respective densities symmetrically spread around the central values.Table 7 gives the BEs, median, 95% BCI and HPD intervals for the α; τ; θ and υ.We can see that the values of MLEs and BEs of the proposed FEPG4 in Tables 6 and 7 are quite identical in υ but slightly different in α; τ and θ.
Figure 8(b-d) displays the non-parametric estimates and curves of the fitted models' reliability functions, FRFs, and MRL on the censored data.The reliability curves of the trained   8(b).When compared to the KM estimate, it is obvious that each of the fitted reliability curves presents, to some extent, a decent match to the KM.The fitted FEPG4 distribution for the BEs (FEPG4-BE) (Figure 8(c)) has produced a more realistic failure rate curve than the other competing models, with both early and wear-out phase change points plainly visible.The non-parametric failure rate of the censored Meeker--Escobar data does not appear to have a long flat zone as the non-censored case, and and hence is similar to the fitted FEPG4-BE failure rate curves.Similarly, the FEPG4-BE works well in describing the MRL and resembles the MRL of the observed sample better than the other models (see, Figure 8(d)).Therefore, we conclude that the numerical findings are in line with Figures 7-8.Thus, the FAEG4 model can be used to investigate and predict the reliability of the censored failure and running time of devices under the Bayesian method.

Thirty FTs of some certain devices
These data represent the failure times (FTs) for n ¼ 30 devices, as reported by Kosznik-biernacka, (2007).For simplicity, the data was named Sylwia data and is presented in the supplementary material.Kosznik-biernacka, (2007) and Tarvirdizade and Ahmadpour, (2021), respectively, employed the data to demonstrate the robustness of Makeham's generalized and Weibull-Chen distributions for fitting FTs data with bathtub FRs.The scaled TTT-transform plot in Figure 11(a) displays a small convex curve followed by a large concave curve, which likely indicates a bathtubshaped FR with a brief flat phase.Consequently, the applicability of the proposed FEPG4 model to the Sylwia data.
Table 8 provides the estimated parameter values (MLEs), À ,, AIC, BIC, AICc, and KS, AD, and CVM for the five fitted models.It can be observed from Table 8 that the introduced FEPG4 model has the least values for the four parametric measures, À ,, AIC, AICc, and BIC, among the competing models.In terms of the similarity between the empirical and fitted distributions, the INMW model has the greatest p-values of KS and AD statistics, followed by the FEPG4 and FACG models.Contrary to KS and AD statistics, the proposed FEPG4 model has the highest p-value of CVM statistic.For an overall picture of the models' performance, we ranked each of the seven criteria, and Figure 9   average, indicating that it provides a better fit for the Sylwia data than the FACG (second), INMW (third), and other trained models.
The trace and marginal density plots of the posterior samples from each of the three chains for the FEPG4 parameters could be observed in Figure 10(a,b).The trace plots reveal that the three  parallel chains of the HMC algorithm converge to the same target distribution for each parameter.Further to that, the PSRF, b R is smaller than 1:1 for each as stated in Table 9, and therefore, ensuring convergence of the chain (see, Figure 10(b)).The densities present relatively acceptable BEs under the SELF as estimates of α; τ; θ; and υ are distributed fairly symmetrically around the midpoints.
The posterior means (BEs), median, and associated 95% credible and HPD intervals for the α; τ; θ and υ are presented in Table 9.Both intervals accommodate their respective estimates adequately.The posterior means and medians of the parameters are comparable and look similar to their respective MLEs.
Figure 11(b-d) displays the reliability functions, FRFs, and MRLs plots for the fitted models on Sylwia data.provide an excellent estimation of reliability curves.It is more evident from the FR curves that the three best models (FEPG4, FACG, and INMW) from numerical findings provide the best description of the FR curves for the data (see, Figure 11(c)).Figure 11(d) Table 8.MLEs, À ,, AIC, BIC, AICc, and KS, AD, and CVM along with their p-values between parenthesis for the considered models on Sylwia data; where the superscript numbers are the ranks for each of the seven criteria.

Model
Optimized  demonstrates the MRLs plots of the fitted distributions and the non-parametric MRL estimate of the data.The MRL curve of the proposed FEPG4 (FEPG4-MLE and FEPG4-BE) model best matches the non-parametric estimate from early through phases.Consequently, the numerical results are similar to the findings presented in Figures 9-11, and we are confident that the FEPG4 model is competent for the analysis and prediction of the device's reliability.

Compatibility of the FEPG3 and FEPG4 Models
If a generated sample from a trained model and the observed data compare favorably, the model is regarded acceptable for the studied data (Gupta et al., 2008;Thach and Briš, 2020).Thus, we use a box plot to compare the samples produced from the trained models with each associated observed data.In each instance, the generated samples from the fitted models have the same sample size as the observed data.
Figure 12 presents box plots of the observed data and predictive samples from the learned distributions for the three illustrations.The Figure 12(a) shows that the FEPG4 Therefore, by these findings and the previous numerical and graphical results, we sum up that our proposed FEPG4 and its submodel, the FEPG3, work better for analyzing failure time data sets whose underlying failure rate distributions are bathtub-shaped with three clear segments as described in Section 1.The FEPG4 distribution's exceptional potential suggests that the model is enough to enable reliability analysis practitioners to determine more precise optimal burn-in procedures or maintenance times that are closely aligned with the reality of the system under consideration.The proposed model may also be appropriate for practitioners interested in developing parametric reliability models such as life-stress and stress--strength models for devices or systems whose failure time distributions are defined by bathtub-shaped FR.Other areas of possible implementation of FEPG4 distribution are discussed in Section 6.

Conclusions and possible future research
This paper introduces a model based on exponential power and Gompertz distributions, referred to as flexible exponential power-Gompertz (FEPG4) distribution.The model is constructed to provide better descriptions of failure time data resulting from a series system with two independent components or a system with one or two failure modes, which due to its complex nature, exhibits bathtub-shaped FR.Classical maximum likelihood and Bayesian are used for the model parameter estimation.The FEPG4 Bayes estimators of the parameters are exploited with Hamiltonian Monte Carlo for posterior simulations, and an extensive simulation experiment shows that the proposed estimators work well.We verified the performance of FEPG4 and its submodel, the FEPG3 model, on failure and running time data (censored and non-censored cases) and EDs failure times, each having bathtub-shaped FR.Based on the description of the FRs, a comparison study was conducted with five more bathtub-shaped methodologies: the AMW, INMW, ACW, EAddW, and FACG models.To understand how well each model describes the data under each case study, the optimized parameters, selection criteria, reliability plots, FR and MRL curves, and box plots are provided for each model.The findings revealed that the FEPG4 (and FEPG3) model might be the best option to represent device failure times, most importantly, when the bathtub-shaped FR of the data at hand clearly shows the burn-in and wear-out change points with a low, yet lengthy, flat middle segment.
The future study proposed for the FEPG4 or its submodel, the FEPG3 distribution, is to construct an LSM for analyzing non-monotone FTs of EDs recorded from an ALT under a voltage stress profile.This proposal is to be achieved by combining FEPG4 and IPL models.Given the versatility of the FEPG4 model, it may provide a strong life-stress connection that may more accurately assess the impact of the parameters on the device's lifetime.We are besides interested in establishing a location-scale regression model based on the FEPG4 distribution.The resulting model may be suited for analysis and prediction of analyzing and predicting early cablejoint failures linked with power cables, given that failure times of power system assets are thought to follow the bathtub-shaped FR.Another possible future investigation derived from this study is to propose several Bayesian estimators of the FEPG4 parameters under different loss functions with various informative and non-informative priors.
It is important to note the proposed methodology's shortcomings.The FEPG4 model is inapplicable for studying the reliability of devices whose failure times follow an upside-down bathtub FR.When modeling data with U-shaped bathtub FR, where the two change points are rarely detected, the proposed model may not be the best solution.That is, it may not be appropriate, for example, for modeling the failure times of 60 electrical appliances placed in ALT as described by Lawless, (2003), in which the FR's change points could not possibly be identified.Hence, other bathtub distributions with fewer parameters could be utilized in this case.

Figure 1 .
Figure 1.Some possible shapes of the FEPG4 density function.

Figure 2 .
Figure 2. Increasing and bathtub FR shapes of the FEPG4 distribution at distinct values of parameters.

Figure 3 .
Figure 3. Performance comparison based on biases and MSEs of the FEPG3 parameters for several chosen values of the parameters at 10 distinct sample sizes.
Figure 6(d) shows the empirical MRL functions for all models, which further support the results.It is clear that the proposed FEPG4 and FEPG3 distributions' fitted MRL and FR functions are close to the empirical MRL and FR values and present the reciprocal relationship between them.

Figure 4 .
Figure 4. Line plots for the ranks of performance measures for the fitted models on non-censored Meeker-Escobar data.

Figure 5 .
Figure 5. (a) Trace plots and (b) posterior density estimates of the parameters obtained by training FEPG4 to non-censored Meeker-Escobar data using HMC sampling for three parallel chains.
depicts the associated line plots for each model.The figure indicates that the FEPG4 model occupies the first place on

Figure 7 .
Figure 7. (a) Trace plots and (b) posterior density estimates of the parameters obtained by training FEPG4 to censored Meeker-Escobar data using HMC sampling for three parallel chains.
Figure 11(b) displays the reliability curves of the trained models as well as the non-parametric Kaplan-MeierKaplan-Meier (KM) estimate of the observed data.Except for the AMW reliability curve, the curves of the other fitted models have closely followed the KM step-function, as seen in Figure 11(b).The FEPG4-MLE and FEPG4-BE also

Figure 9 .
Figure 9. Line plots for the ranks of performance measures for the fitted models on Sylwia data.

Figure 10 .
Figure 10.(a) Trace plots and (b) marginal posterior density plots for three parallel chains of the FEPG4 parameters on Sylwia data.

Figure 12 .
Figure 12.Box plots for the (a) non-censored Meeker-Escobar, (b) censored Meeker-Escobar, and (c) Sylwia data (raw data) with separate generated samples from the trained FEPG3, FEPG4, and other distributions under the three illustrations.

Table 1 .
Estimators' performance comparison based on biases and MSEs of the FEPG3 parameters for several chosen values of the parameters at five distinct sample sizes.

Table 2 .
Estimators' performance comparison based on biases and MSEs of the FEPG3 parameters for several chosen values of the parameters at five distinct sample sizes.

Table 5 .
Posterior means, median with their 95% BCI and HPD, and scale reduction factor ( b R) for the FEPG4 parameters on noncensored Meeker--Escobar data.

Table 7 .
Posterior means, median with their 95% BCI and HPD, and scale reduction factor ( b R) for the FEPG4 parameters on censored Meeker--Escobar data.

Table 9 .
Posterior means, median with their 95% BCI and HPD, and scale reduction factor ( b R) for the FEPG4 parameters on Sylwia data.