Modified Chen distribution: Properties, estimation, and applications in reliability analysis

Abbreviations

1.   Introduction

In survival analysis, the hazard rate (HR) function plays a critical role in studying a product's lifecycle. A bathtub-shaped HR function illustrates a pattern over time that begins with a high initial HR, followed by a period of low, steady HR, and then an increase later in the lifecycle. This pattern is relevant in fields such as product reliability, warranty analysis, infrastructure planning, healthcare systems, manufacturing, aerospace, aviation, and telecommunications. Understanding the HR pattern is essential for determining optimal maintenance strategies. For instance, electronic components or machinery may exhibit a high initial HR due to manufacturing defects or early-life failures, which then decreases to a stable rate during their useful life. Toward the end of the lifecycle, the HR may rise again due to wear and aging. Similarly, manufacturing equipment may experience high early HRs from issues like malfunctioning components or calibration errors, stabilize, and then increase again due to wear and tear or obsolescence. Recognizing this pattern helps in scheduling preventive maintenance, optimizing production efficiency, and reducing downtime ^[1,2,3,4].

Practitioners analyze statistical probability distributions to determine whether a device exhibits bathtub curve behavior, with a focus on warranty periods and optimal maintenance schedules. Classical probability distributions, like the Weibull, lognormal, and exponential, are often employed for their mathematical flexibility. However, these distributions may not accurately fit data that exhibit non-monotonic behavior, such as the bathtub curve. This limitation can introduce uncertainty in data interpretation and may lead to inconclusive results in the analysis due to the specific characteristics of these distributions.

The Chen distribution (CD) ^[5], introduced by Chen in 2000, is a life distribution widely used in survival analysis and modeling. It is defined by two parameters and is adaptable for various applications. The model also provides closed-form confidence intervals for the shape parameter and joint confidence regions for both parameters. However, the CD has limitations in accurately modeling certain survival datasets, particularly due to its asymmetric HR, which does not conform to a bathtub shape, and its lack of a scale parameter.

To address these limitations, researchers have developed modified versions of the CD to improve its flexibility in describing diverse data types. For instance, Xie et al. ^[6] introduced the modified Weibull extension distribution (MWED) by incorporating a scale parameter. Pappas et al. ^[7] further expanded on this by adding a shape parameter to the MWED, resulting in the Marshall-Olkin extended CD. Sarhan and Apaloo ^[8] introduced the exponentiated MWED, a four-parameter distribution that extends the MWED. Chaubey and Zhang ^[9] proposed the exponentiated CD (ExpCD), a three-parameter extension of the CD. Recently, Reis et al. ^[10] presented the gamma-CD (GCD), and Bhatti et al. ^[11] developed another flexible model, also named the extended CD (ECD), derived from the generalized Burr-Hatke differential equation. Additionally, Anafo et al. ^[12] proposed the modified ECD (MECD), and Acquah et al. ^[13] introduced the new ECD (NECD), a three-parameter distribution.

While these modifications have improved data-fitting capabilities for specific applications, they have failed to accurately represent complex hazard behaviors, particularly non-monotonic patterns such as bathtub-shaped HRs. One significant issue is their inability to model the bathtub curve accurately, often producing hazard functions that resemble 'V' or 'J' shapes, which incorrectly suggest a short operational life for the device. This misrepresentation biases the data, leading to inadequate estimates of device reliability. The research gap lies in the absence of a three-parameter distribution capable of effectively modeling all three phases of a product's lifecycle-early failure, stable operation, and wear-out-without introducing excessive mathematical complexity or practical limitations. Therefore, it is essential to explore alternative methodologies or develop hybrid distributions that better align with the characteristics of the bathtub curve.

In recent years, several models have been developed to represent reliability by combining two distinct classical distributions, even when they belong to the same distribution family. Examples of these models include the additive Chen-Weibull distribution ^[14], Weibull-CD ^[15], additive Chen-Gompertz distribution ^[16], additive CD ^[17], and additive Chen-Perks distribution ^[18]. These models are capable of providing a bathtub-shaped hazard function. However, by merging two classical distributions, each with at least two parameters, these models typically involve four or five parameters, increasing their complexity and making parameter estimation more challenging. This highlights the need for a distribution that can effectively handle system reliability with fewer parameters, balancing flexibility and simplicity.

In this article, we present a novel, flexible three-parameter distribution known as the modified CD (MCD). The MCD addresses critical shortcomings of existing models by incorporating an additional parameter that significantly enhances flexibility while maintaining interpretability. This novel three-parameter distribution overcomes the absence of a scale parameter in the original CD and enables the modeling of both monotonic and non-monotonic hazard functions, including the three phases of the bathtub curve-early failure, stable operation, and wear-out. The MCD's highly adaptable HR function exhibits a bathtub shape, offering a versatile and practical solution readily applicable to real-world scenarios. Its simplicity facilitates ease of parameter estimation, making it an attractive option for reliability engineers, and it consistently provides a better fit compared to the other three-parameter models. This unique combination of flexibility, simplicity, and practical applicability distinguishes the MCD from previous modifications, making it a valuable contribution to the field of reliability analysis.

The article is organized as follows: Section 2 introduces the new distribution. Section 3 delves into the statistical properties of the distribution. Section 4 details the parameter estimation methods, including maximum likelihood estimators, two bootstrap confidence intervals, and Bayesian estimates via the Markov chain Monte Carlo (MCMC) method under balanced loss functions. Section 5 presents the results of a simulation study. In Section 6, the applications of the MCD to two real datasets are discussed, highlighting its importance and flexibility. Finally, Section 7 provides the conclusions and future work regarding the MCD.

2.   The MCD

In this section, we introduce a new three-parameter MCD, obtained by modifying the cumulative distribution function (CDF) of the CD ^[5]. The CDF of the MCD, parameterized by the vector $\underline{\mathbf{\psi}} = (\alpha, \, \gamma, \, \lambda)$, is defined as follows:

where $\alpha$ is the location parameter, allowing the distribution to adjust its position; $\gamma$ is the scale parameter, which adjusts the scale of the distribution, effectively stretching or compressing the survival function along the time axis; and $\lambda$ is the shape parameter, which governs the overall shape of the distribution, influencing how rapidly the hazard rate changes.

The probability density function (PDF) and the survival function $S(x)$ are defined as follows:

and

The HR and cumulative hazard functions are expressed as follows:

and

Different shapes of the PDF of the MCD are shown in Figure 1. The PDF of the MCD can be unimodal, decreasing, or increasing. Equation (2.4) demonstrates that the HR function increases when $\lambda \geq 1$ and displays a bathtub shape when $0 < \lambda < 1$. For $0 < \lambda < 1$, the HR function presents a bathtub shape if $h^{'}(x^{*}; \underline{\mathbf{\psi}}) = 0$, with $x = x^{*}$ being the root of the equation:

The HR function decreases when $h^{'}(x^{*}; \underline{\mathbf{\psi}}) < 0$ for $x < x^{*}$ and increases when $h^{'}(x^{*}; \underline{\mathbf{\psi}}) > 0$ for $x > x^{*}$. Therefore, the HR function exhibits different behaviors depending on the value of $\lambda$, summarized as follows:

● For $\lambda\geq1$, the HR function increases steadily throughout its entire range (refer to Figure 2(a)).

● For $0 < \lambda < 1$, the HR function exhibits a bathtub-like shape (refer to Figure 2(b)–(d)).

3.   Statistical properties of the MCD

3.1.   Quantiles and mode

The $p$th quantile ($x_{p}$) of the MCD$(\underline{\mathbf{\psi}})$ is calculated as the real solution to the following nonlinear equation:

By setting $p = 0.25, \, 0.5, \, 0.75$, we can determine the first, second, and third quartiles of the MCD.

The mode of the MCD is obtained by solving the following nonlinear equation:

Equations (3.1) and (3.2) lack simple closed form solutions, necessitating the use of numerical methods.

3.2.   Moments

Theorem 3.1. Given that $X$ follows the MCD$(\underline{\mathbf{\psi}})$, the $r$th moment of $X$, for $r = 1, \, 2, \, \cdots$ is defined as

where $\mu'_{m+r, \textbf{CD}}$ is the $(m+r)$th noncentral moment of the CD.

Proof. See Appendix A.

Using the first four ordinary moments of the MCD, the measures of skewness ($SK$) and kurtosis ($KU$) can be determined as follows:

and

Table 1 presents the values of the first four moments, variance, $SK$, and $KU$ for the MCD, taking into account different combinations of $\alpha, \, \gamma$, and $\lambda$.

3.3.   Incomplete moments

Theorem 3.2. Assume $X$ is a continuous random variable following the MCD$(\underline{\mathbf{\psi}})$. The incomplete moment of $X$ is then given by:

where

and

Proof. See Appendix A.

3.4.   Mean time to failure

Theorem 3.3. Given that $X$ follows the MCD$(\underline{\mathbf{\psi}})$, the mean time to failure (MTTF) of the MCD is expressed as:

where $\Gamma(.)$ denotes the gamma function.

Proof. See Appendix A.

3.5.   Mean residual life

Theorem 3.4. The mean residual life (MRL) of the MCD$(\underline{\mathbf{\psi}})$ is expressed as

Proof. See Appendix A.

Table 2 presents the outcomes of a Monte Carlo simulation for the MRL $(M_{X}(t))$ and MTTF. We generated $N = 2000$ random samples with sizes of $50,150,250,350,$ and $450$ from the MCD, using various parameter values. The $M_{X}(t)$ was evaluated at time points $t = 0.1, 0.25, 0.4, 0.9,$ and $1.5$. It was observed that the expected remaining lifetime (MRL) of an individual or system approaches the MTTF as $t$ approaches zero, i.e., $M_{X}(t)\rightarrow$ MTTF as $t\rightarrow 0$, confirming that $\text{MTTF} = M_{X}(t)|t = 0$. Additionally, we noticed that $M_{X}(t)$ decreases as $t$ increases for a fixed $n$, and the mean squared errors (MSEs) (in parentheses) decreases as $n$ increases (see Figure 3). Moreover, the MTTF remained relatively stable across different sample sizes. Figure 4 also presents the trace plots and density estimates for MRL and MTTF produced by the Monte Carlo simulation.

3.6.   Rényi entropy

Theorem 3.5. Given that $X$ follows the MCD$(\underline{\mathbf{\psi}})$, the Rényi entropy of $X$ is given by

Proof. See Appendix A.

3.7.   Order statistics

Theorem 3.6. Consider an ordered sample $\{X_{i}\}_{i = 1}^{n}, \, n\geq1$ from the MCD$(\underline{\mathbf{\psi}})$. The $r$th moments of the $l$th order statistic are given as

where $\mu'_{r}(x; \alpha', \gamma, \lambda)$ is provided in Eq (3.3), with the parameters $\alpha' = (n+j+1-l)\, \alpha, \, \gamma, \, \lambda$.

Proof. See Appendix A.

4.   Parameter estimation

4.1.   Maximum likelihood estimation

In this subsection, we estimate the parameters $\alpha, \, \gamma,$ and $\lambda$ of the MCD using maximum likelihood estimators (MLEs), along with the functions $S(x)$ and $h(x)$. Additionally, we calculate asymptotic confidence intervals (ACIs) for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$ by leveraging the normality property of the corresponding MLEs.

Consider $\underline{x} = (x_{1}, \ldots, x_{n})$ as an observed sample drawn from the MCD with an unknown parameter vector $\underline{\mathbf{\psi}} = (\alpha, \, \gamma, \, \lambda)^{T}$. Based on Eq (2.2), the likelihood function of MCD can be expressed as follows:

The log-likelihood function can then be derived as:

To obtain the MLEs of $\underline{\mathbf{\psi}}$, we begin by computing the first partial derivatives of $\ell(D|\underline{\mathbf{\psi}})$ with respect to $\alpha, \, \gamma$, and $\lambda$. The resulting likelihood equations are as follows:

and

By setting these expressions to zero, the MLEs for the parameters $\hat{\underline{\mathbf{\psi}}} = (\hat{\alpha}, \, \hat{\gamma}, \, \hat{\lambda})^T$ can be determined by solving the system of nonlinear Eqs (4.3)–(4.5) for $\alpha, \, \gamma$, and $\lambda$. Since these equations do not have an analytical solution, numerical methods, such as the Newton–Raphson method, are necessary to find the estimates. Mathematica software provides subroutines for nonlinear optimization problems, and we utilized the $FindRoot[\{{eqn_{1}}, {eqn_{2}}, \ldots\}, \{\{x, x_{0}\}, \{y, y_{0}\}, \ldots\}]$ package for this purpose.

The MLEs for $S(x)$ and $h(x)$ are given by:

and

4.1.1.   ACI

To obtain the CIs for the parameters $\alpha, \, \gamma$, and $\lambda$, we require the distributions of the MLEs $\hat{\alpha}, \, \hat{\gamma}$, and $\hat{\lambda}$. Since these MLEs do not have closed-form solutions, determining their exact distributions is not feasible. Therefore, we derive ACIs based on the asymptotic normality of these parameters. Under standard regularity conditions, $\sqrt{n}\, (\underline{\mathbf{\psi}}-\hat{\underline{\mathbf{\psi}}})$ is asymptotically distributed as a multivariate normal $N_{3} (0, \, J^{-1}(\underline{\mathbf{\psi}}))$, where

where the formulas for the second partial derivatives are provided in Appendix B. The approximate variance-covariance matrix can be evaluated at $\hat{\underline{\mathbf{\psi}}} = (\hat{\alpha}, \, \hat{\gamma}, \, \hat{\lambda})^T$, the MLE of $(\hat{\alpha}, \, \hat{\gamma}, \, \hat{\lambda})^T$:

Thus, the $(1-\tau)100\%$ symmetric approximate normal CIs for $\underline{\mathbf{\psi}} = (\alpha, \, \gamma, \, \lambda)$ are given by

where $i = 1, \, 2, \, 3$ and $z_{\frac{\tau}{2}}$ is the upper $\frac{\tau}{2}$ point of the standard normal distribution.

To derive the ACIs for $S(x)$ and $h(x)$, we apply the delta method to compute their variances, as described by Greene ^[19]. Let $A_{1} = \Big(\begin{array}{ccc} \frac{\partial\, S(x)}{\partial\, \alpha} & \frac{\partial\, S(x)}{\partial\, \gamma} & \frac{\partial\, S(x)}{\partial\, \lambda} \end{array}\Big)$ and $A_{2} = \Big(\begin{array}{ccc} \frac{\partial\, h(x)}{\partial\, \alpha} & \frac{\partial\, h(x)}{\partial\, \gamma} & \frac{\partial\, h(x)}{\partial\, \lambda} \end{array}\Big)$, where $\frac{\partial\, S(x)}{\partial\, \alpha}, \, \frac{\partial\, S(x)}{\partial\, \gamma}, \, \frac{\partial\, S(x)}{\partial\, \lambda}, \, \frac{\partial\, h(x)}{\partial\, \alpha}, \, \frac{\partial\, h(x)}{\partial\, \gamma},$ and $\frac{\partial\, h(x)}{\partial\, \lambda}$ are the first derivatives of $S(x)$ and $h(x)$ with respect to the parameters $\alpha, \, \gamma,$ and $\lambda$, respectively. The approximate asymptotic variances of $\hat{S}(x)$ and $\hat{h}(x)$ are given by

where $A^{T}_{k}$ is the transpose of $A_{k},$ for $k = 1, \, 2.$ These results yield the ACIs for $S(x)$ and $h(x)$ as

4.2.   Bootstrap CIs

In this subsection, we present two parametric bootstrap methods for constructing CIs for $\alpha, \, \gamma, \, \lambda, \, S(x)$, and $h(x)$. These methods are the percentile bootstrap-p (Boot-P) CI, introduced by Efron ^[20], and the bootstrap-t (Boot-T) CI, proposed by Hall ^[21]. The following algorithms outline the process of estimating CIs using both approaches.

4.2.1.   Parametric Boot-P

Algorithm 1: Boot-P method

1. Begin by using the original dataset $\underline{x} = X_{1:n}, \, X_{2:n}, \, \ldots, \, X_{n:n}$ to compute $\hat{\alpha}, \, \hat{\gamma},$ and $\hat{\lambda}$ by maximizing Eqs (4.3)–(4.5).

2. Generate a bootstrap sample $\underline{x}^* = X_{1:n}^*, \, X_{2:n}^*, \, \ldots, \, X_{n:n}^*$ by performing resampling with replacement.

3. Calculate the MLEs from the bootstrap sample and denote this estimate as $\hat{\varphi}^*$, where $\varphi = (\alpha, \, \gamma, \, \lambda, \, S(x), \, h(x)).$

4. Repeat Steps 2 and 3 NBoot times, producing a series of estimates $\hat{\varphi}^*_{1}, \hat{\varphi}^*_{2}, \, \ldots, \, \hat{\varphi}^*_{ \text{NBoot}}$, where $\hat{\varphi}^*_{l} = \big(\hat{\alpha}^{*}_{l}, \, \hat{\gamma}^{*}_{l}, \, \hat{\lambda}^{*}_{l}, \, \hat{S}^{*}_{l}(x), \, \hat{h}^{*}_{l}(x)\big)$ for $l = 1, \, 2, \, \ldots, \, \text{NBoot}$.

5. Arrange the series of estimates $\hat{\varphi}^*_{l}, \, l = 1, \, 2, \, \ldots, \, \text{NBoot}$ in ascending order, and then compute $\hat{\varphi}^*_{(1)}, \hat{\varphi}^*_{(2)}, \, \ldots, \, \hat{\varphi}^*_{(\text{NBoot})}$.

6. Define $\widehat{G}_{1}(u) = P(\hat{\varphi}^*\le u)$ as the CDF of $\hat{\varphi}^*$. Idintify $\hat{\varphi}_{ \text{NBootP}} = \widehat{G}^{-1}_{1}(u)$ for a given value of $u$. The approximate bootstrap-P $100(1-\tau)\%$ CIs for $\hat{\varphi}$ are then given by:

4.2.2.   Parametric Boot-T

Algorithm 2: Boot-T method

1. Follow Steps 1 through 3 as described in the Boot-P method.

2. Compute the $T^{*\varphi}$ statistic defined as: $T^{*\varphi} = \dfrac{(\hat{\varphi}^{*}-\hat{\varphi})}{\sqrt{\widehat{ \text{Var}}(\hat{\varphi}^{*})}}$, where $\widehat{ \text{Var}}(\hat{\varphi}^{*})$ is determined using Eq (4.6).

3. Repeat Steps 1 and 2 for NBoot iterations, calculating $T^{*\varphi}_{1}, \, T^{*\varphi}_{2}, \ldots, \, T^{*\varphi}_{ \text{NBoot}}.$

4. Arrange the sequence $T^{*\varphi}_{1}, \, T^{*\varphi}_{2}, \ldots, \, T^{*\varphi}_{ \text{NBoot}}$ in ascending order to obtain $T^{*\varphi}_{(1)}, \, T^{*\varphi}_{(2)}, \ldots, \, T^{*\varphi}_{ \text{(NBoot)}}.$

5. Define $\widehat{G}_{2}(u) = P(T^*\le u)$ as the CDF of $T^*$. For a given $u$, set $\hat{\varphi}_{ \text{NBootT}} = \hat{\varphi}+\widehat{G}^{-1}_{2}(u)\, \sqrt{\widehat{ \text{Var}}(\hat{\varphi}^{*})}.$ The approximate bootstrap-T $100(1-\tau)\%$ CIs for $\hat{\varphi}$ are then:

4.3.   Bayesian estimation

In this subsection, we derive Bayesian estimates and the corresponding CIs for the unknown parameters $\alpha, \gamma$, and $\lambda$, as well as for the functions $S(x)$ and $h(x)$. Additionally, we apply Bayes estimation using the MCMC method under balanced loss functions, including balanced squared error loss (BSEL), balanced linear-exponential loss (BLINEXL), and balanced general entropy loss (BGEL).

4.3.1.   Bayesian estimation using the MCMC method

In Bayesian analysis, specifying a loss function is crucial for determining the optimal estimate of an unknown parameter. To comprehensively compare Bayes estimates, we use three types of loss functions: BSEL, BLINEXL, and BGEL. For further research on the balanced loss function, refer to references ^{[22,23,24,25,26]}.

Assume that the parameters $\alpha, \gamma$, and $\lambda$ are independent random variables, and according to the literature ^[27,28], they follow gamma prior distributions:

As a result, the joint prior distribution $\pi(\underline{\mathbf{\psi}})$ for $\alpha, \, \gamma$ and $\lambda$ is

The posterior distribution of $\alpha, \, \gamma,$ and $\lambda,$ denoted as $\pi^*(\underline{\mathbf{\psi}}|D),$ can be derived by integrating the likelihood function from Eq (4.1) with the joint prior distribution from Eq (4.7). By applying Bayes' theorem, the posterior distribution $\pi^*(\underline{\mathbf{\psi}}|D)$ for $\underline{\mathbf{\psi}}|D$ is given by:

where $\int_{\underline{\mathbf{\psi}}}\, L(D|\underline{\mathbf{\psi}})\times\pi(\underline{\mathbf{\psi}})\, d\underline{\mathbf{\psi}}$ represents the normalizing constant of the posterior distribution of $\underline{\mathbf{\psi}}$, also known as the marginal distribution of $D$. Since this constant is not required for Bayesian inference using MCMC methods, the posterior distribution is typically expressed as:

By substituting Eqs (4.1), (4.7), and (4.8), the joint posterior distribution of $\underline{\mathbf{\psi}}|D$ is given by:

The MCMC technique is utilized to generate samples from Eq (4.9). These samples are then used to compute the Bayes estimates of $\alpha, \, \gamma$ and $\lambda$, as well as related functions like $S(x)$ and $h(x)$, and to construct the CIs. The Gibbs within Metropolis sampler is applied to perform the MCMC technique, which requires deriving the full set of conditional posterior distributions. The marginal posterior density for $\alpha, \, \gamma,$ and $\lambda$ is derived from Eq (4.9) as follows:

and

The conditional posterior density of $\alpha$ follows a gamma distribution with a shape parameter of $(n+c_{1})$ and a scale parameter of $\Big(d_{1}-\sum _{i = 1}^n \left(2-{\rm e}^{\gamma\, x_i}-{\rm e}^{x_i^{\lambda }}\right)\Big)$. Consequently, samples of $\alpha$ can be generated using any gamma distribution generation method. Although the conditional posteriors of $\gamma$ and $\lambda$ do not follow standard forms, they resemble a normal distribution, as shown in Figure 5, making Gibbs sampling impractical. Therefore, the Metropolis-Hastings (M-H) sampler is necessary for implementing the MCMC methodology. We will outline the steps involved in the M-H within the Gibbs sampling method as presented in Algorithm 3.

Algorithm 3: M-H within Gibbs sampling

1. Begin by selecting initial values for the chain $(\alpha^{0}, \, \gamma^{0}, \, \lambda^{0})$ and define $M$ as the burn-in period.

2. Set the iteration counter to $j = 1$.

3. Generate $\alpha^{j}$ from a Gamma distribution with parameters $\Big(n+c_{1}, \, d_{1}-\sum _{i = 1}^n \left(2-{\rm e}^{\gamma\, x_i}-{\rm e}^{x_i^{\lambda }}\right)\Big)$.

4. Using the M-H algorithm, generate $\gamma^{(j)}$ and $\lambda^{(j)}$ from Eqs (4.11) and (4.12) utilizing normal proposal distributions $N\big(\gamma^{(j-1)}, \, \mathrm{Var}({\gamma})\big)$ and $N\big(\lambda^{(j-1)}, \, \mathrm{Var}({\lambda})\big),$ where $\mathrm{Var}({\gamma})$ and $\mathrm{Var}({\lambda})$ are obtained using Eq (4.6).

● Generate proposals $\gamma^{*}$ from $N\big(\gamma^{(j-1)}, \, \mathrm{Var}({\gamma})\big)$ and $\lambda^{*}$ from $N\big(\lambda^{(j-1)}, \, \mathrm{Var}({\lambda})\big)$.

● Calculate the acceptance probabilities:

● Generate $u$ from a Uniform $(0, 1)$ distribution.

● If $u\leq\rho_{1}$, accept the proposal and set $\gamma^{(j)} = \gamma^{*}$; otherwise, set $\gamma^{(j)} = \gamma^{(j-1)}$.

● If $u\leq\rho_{2}$, accept the proposal and set $\lambda^{(j)} = \lambda^{*}$; otherwise, set $\lambda^{(j)} = \lambda^{(j-1)}$.

5. Compute $S(x)$ and $h(x)$ as follows:

6. Set $j = j+1$.

7. Repeat Steps 3–6 for $N$ iterations to generate

8. Calculate the Bayes estimates of $\varphi = \big(\alpha, , \gamma, , \lambda, , S(x), , h(x)\big)$ after the burn-in period $M$ as:

9. Determine the Bayes estimates of $\varphi$ under the BSEL, as introduced by Jozani ^[23], as follows:

10. Obtain the Bayes estimates of $\varphi$ under the BLINEXL, as presented by Jozani ^[22], as follows:

11. Compute the Bayes estimates of $\varphi$ using the BGEL, as outlined by Jozani ^[25], as follows:

12. To calculate the CIs for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$, first sort the generated values $\alpha^{(l)}, \, \gamma^{(l)}, \, \lambda^{(l)}, \, S^{(l)}(x),$ and $h^{(l)}(x)$ for $l = 1, \, \ldots, N$ in ascending order as follows: $\Big\{\alpha^{(1)} < \ldots < \alpha^{(N)}\Big\}$, $\Big\{\gamma^{(1)} < \ldots < \gamma^{(N)}\Big\}$, $\Big\{\lambda^{(1)} < \ldots < \lambda^{(N)}\Big\}$, $\Big\{S^{(1)}(x) < \ldots < S^{(N)}(x)\Big\}$, and $\Big\{h^{(1)}(x) < \ldots < h^{(N)}(x)\Big\}$. Then, the $(1-\tau)100\%$ CIs for $\varphi = (\alpha, \, \gamma, \, \lambda, \, S(x), \, h(x))$ are given by $\Big(\varphi_{\big(N\, \dfrac{\tau}{2}\big)}, \, \, \varphi_{\big(N\, (1-\dfrac{\tau}{2})\big)}\Big).$

5.   Simulation study

This section outlines a simulation study evaluating the performance of the MLEs, along with two bootstrap CIs (Boot-P and Boot-T) and Bayesian estimators using MCMC under BSEL, BLINEXL, and BGEL, to estimate the parameters of the MCD, as well as $S(x)$ and $h(x)$. The following steps are employed for the simulation study:

1. Set initial values for the parameters $\alpha, \, \gamma$, and $\lambda$, and assume the samples sizes $n = 50, \, 100, \, 150,$ and $200$.

2. Generate a random sample $x_{1}, \ldots, x_{n}$ of size $n$ from Eq (3.1).

3. Replicate each sample $N = 1000$ times.

4. Calculate the bias and MSE for the parameters $\alpha, \, \gamma$, and $\lambda$, as well as for $S(x)$ and $h(x)$, according to the procedures detailed in Subsection 4.1.

5. Apply Algorithm 1 of the Boot-B to calculate the bias and MSE for the parameters $\alpha, \, \gamma$, and $\lambda$, as well as for $S(x)$ and $h(x)$.

6. Use Algorithm 2 of the Boot-T to compute the bias and MSE for the parameters $\alpha$, $\gamma$, and $\lambda$, as well as for $S(x)$ and $h(x)$.

7. Use Algorithm 3 to calculate the bias and MSE for the parameters $\alpha$, $\gamma$, and $\lambda$, as well as $S(x)$ and $h(x)$, utilizing chosen hyperparameter values to generate posterior samples from each marginal posterior distribution.

8. Compute the bias and MSEs using the results from Step 7 based on BSEL, BLINEXL, and BGEL.

All simulations were performed using Wolfram Mathematica, generating 11000 MCMC samples and discarding the first 1000 samples as the 'burn-in' period.

The bias and MSE of the MCD parameter estimates across various sample sizes using MLE, Bayes MCMC, Boot-P, and Boot-T are presented in Table 3. The table indicates that the MLEs and Bayes MCMC demonstrate similar efficiency. Additionally, Boot-P outperforms Boot-T by exhibiting smaller MSE for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$. Tables 4–8 provide numerical results on bias and MSE for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$ based on Bayes estimates under BSEL, BLINEXL, and BGEL. The results indicate that Bayes estimates under these loss functions maintain low bias and MSE across all scenarios. Furthermore, both bias and MSE tend to decrease as the sample size increases, indicating the consistency of the estimates. When $\omega = 0$, Bayes estimates yield better results for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$ by achieving smaller MSEs. Bayes estimates under BLINEXL with $c = 0.7$ consistently produce better estimates with smaller MSE when $\omega = 0$, while BGEL with $q = 0.7$ offers improved estimates with smaller MSE when $\omega = 0.9$. Finally, Figure 6 presents trace plots and density estimates derived from the Monte Carlo simulation conducted for the MCD, utilizing $n = 450$ and $t = 0.9$, with parameters $\alpha = 0.1$, $\gamma = 0.3$, and $\lambda = 0.7$.

6.   Applications

This section illustrates the application of the MCD using two real datasets from reliability engineering. The first dataset is the Aarset data ^[36], which contains the lifetimes of fifty devices. The second dataset is the Meeker-Escobar data ^[37], which represents the failure and operating times of thirty devices. These datasets were chosen because their underlying distributions exhibit a characteristic bathtub shape (see Figures 7 and 10). This makes them widely recognized as benchmark datasets in the literature for evaluating the fit of distributions with a bathtub-shaped HR function. Additionally, they enable practitioners to effectively utilize the HR function for predictive maintenance and reliability analysis in engineering applications. We compared the MCD with competing models listed in Table 9, using various metrics such as Log-likelihood ($\ell$), Kolmogorov-Smirnov ($\text{K-S}$) statistics with their corresponding P-values, Anderson-Darling ($A^*$), Cramér-von Mises ($W^*$), and several information criteria including the Akaike information criterion (AIC), Bayesian information criterion (BIC), and Hannan-Quinn information criterion (HQIC). All computations were carried out using Wolfram Mathematica 12.3 software.

6.1.   The lifetime of fifty electronic devices

The Aarset dataset ^[36], which represents the failure times of fifty electronic devices, has been extensively analyzed in the literature, with the latest studies referenced in ^[38,39,40]. As demonstrated by the scaled total time on test transform (TTT-transform) plot in Figure 7, the dataset exhibits a bathtub-shaped HR. Table 10 provides estimates for MLEs and Bayes MCMC, along with their $95\%$ CIs, for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$ applied to the Aarset data. Additionally, Table 11 presents estimates for Boot-P and Boot-T methods, including their corresponding $95\%$ CIs, for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$.

Bayes estimates derived from BSEL, BLINEXL, and BGEL functions, with various values of $c, \, q,$ and $\omega$ for the parameters $\alpha, \, \gamma, \, \lambda$, as well as $S(x)$ and $h(x),$ are summarized in Table 12. Further, Table 13 describes the MRL at specific time points and the Rényi entropy at different $\rho$ values for the fitted MCD. Observations from Table 13 indicate that the MRL tends to increase and then decrease as time progresses, reflecting its inverse relationship with the bathtub-shaped HR of the MCD. The calculated MTTF values are 46.4929, 46.9939, 44.9965, and 41.0331 for the MLE, Bayes MCMC, Boot-P, and Boot-T methods, respectively.

Table 14 compares the MLEs, $\ell$, K-S statistics, P-values, $A^*$, $W^*$, $AIC$, $BIC$, and $HQIC$ for the MCD and other competitive models. The MCD shows the smallest values for $\text{K-S}$, $A^*$, $W^*$, $AIC$, $BIC$, and $HQIC$, and the highest $\ell$ and P-value, highlighting its superior fit to the Aarset data compared to other models. Figure 7 shows the empirical and fitted scaled TTT-transform plot for the MCD, while Figure 8 depicts the survival, hazard, and cumulative hazard functions of the MCD in comparison to competitive models. Figure 8 offers a graphical depiction of the devices' behavior over time. Figure 8(a), (b) illustrates the lifetimes of the Aarset data in relation to the MCD's survival equation and competing models, which are crucial in reliability engineering for estimating lifetimes and maintenance times to enhance product life. Notably, the MCD shows a better fit with Kaplan-Meier empirical reliability compared to other distributions, suggesting it more accurately describes lifetimes, maintenance times, and MTTF. Figure 8(c), (d) depicts the failure patterns, showing that the MCD aligns more closely with Kaplan-Meier empirical behavior compared to other distributions, suggesting it provides a more accurate representation of device failures. This information is valuable for product engineers aiming to improve designs, reduce costs, and estimate maintenance expenses. Finally, Figure 9 presents boxplots for the Aarset data and samples generated from competing distributions, showing that the MCD's boxplot more accurately reflects the range and variability of the Aarset data compared to other models. Consequently, the graphical results (Figures 8 and 9) support the numerical findings, endorsing the MCD as a suitable model for analyzing and predicting device failure times.

6.2.   The lifetime of thirty electronic devices

The second application, referred to as the Meeker-Escobar data ^[37], includes the failure and running times of thirty electronic devices. This dataset has been extensively analyzed by numerous researchers ^[39,41]. As indicated by the scaled TTT-transform plot in Figure 10, the dataset exhibits a bathtub-shaped HR. Table 17 presents MLEs, Bayes MCMC estimates, and their $95\%$ CIs for the parameters $\alpha, \, \gamma, \, \lambda$, as well as $S(x)$ and $h(x),$ based on the Meeker-Escobar data. Additionally, Table 15 provides estimates for the Boot-P and Boot-T methods, along with their corresponding $95\%$ CIs for $\alpha, \, \gamma, \, \lambda, \, S(x),$ and $h(x)$.

Table 16 summarizes the results of Bayes estimates using the BSEL, BLINEXL, and BGEL functions, with varying values of $c$ for BLINEXL, $q$ for BGEL, and different values of $\omega$ for $\alpha, \, \gamma, \, \lambda$, as well as $S(x)$ and $h(x)$ applied to the Meeker-Escobar data. Furthermore, Table 18 presents the MRL at specific time points and the Rényi entropy for different values of $\rho$ for the fitted MCD. The table shows that the MRL for the MLE, Bayes MCMC, Boot-P, and Boot-T methods decreases over time, consistent with the empirical MRL curve. The calculated MTTF values are 182.437,180.078,112.624, and 112.547 for the MLE, Bayes MCMC, Boot-P, and Boot-T methods, respectively.

Table 19 compares the MLEs, $\ell$, $\text{K-S}$ statistics, P-values, $A^*$, $W^*$, $AIC$, $BIC$, and $HQIC$ for the MCD and other competitive models. The MCD shows the lowest values for $\text{K-S}$, $A^*$, $W^*$, $AIC$, $BIC$, and $HQIC$, along with the highest $\ell$ and P-value, indicating a better fit to the Meeker-Escobar data compared to the other models. Figure 10 shows the empirical and fitted scaled TTT-transform plot for the MCD, while Figure 11 illustrates the survival, hazard, and cumulative hazard functions of the MCD in comparison to competitive models fitted to the Meeker-Escobar data. Figure 11(a), (b) highlights the lifetimes of the Meeker-Escobar data relative to the MCD's survival equation and competing models. The Kaplan-Meier empirical reliability indicates that the MCD provides a superior fit compared to other analyzed distributions, offering a more accurate representation of lifetimes, maintenance times, and MTTF. Figure 11(c), (d) illustrates the failure patterns, showing that the MCD aligns more closely with empirical behavior compared to other distributions, providing a more accurate representation of device failures. Finally, Figure 12 shows boxplots for the Meeker-Escobar data and samples from competing distributions, with the MCD's boxplot closely resembling the range and variability of the Meeker-Escobar data, indicating a more accurate representation of device failure times compared to other fitted distributions. Thus, Figures 11 and 12 reinforce the numerical results, supporting the MCD as a suitable model for analyzing and predicting device failure times. This suggests that the MCD is more appropriate for modeling device failures over its lifespan, making it a valuable tool for reliability engineering practitioners.

7.   Conclusions

In this paper, we introduced and examined a novel lifetime distribution called the MCD for use as a reliability and survival model. The HR function of this distribution is straightforward and capable of encompassing both increasing and bathtub-shaped HRs. We studied several statistical properties of this distribution, such as moments, MTTF, MRL, Rényi entropy, and order statistics. These properties establish a robust mathematical framework for understanding the MCD's behavior and its real-world applications. The unknown model parameters, along with the survival and hazard functions, were estimated using maximum likelihood, two parametric bootstrap methods (Boot-P and Boot-T), and Bayesian methods via MCMC with SEL, BSEL, BLINEXL, and BGEL loss functions. The MCD's flexibility with various estimation methods and loss functions makes it a versatile and powerful tool for reliability practitioners. Additionally, ACIs for the parameters, as well as survival and hazard functions, were obtained using various methods. A simulation study was conducted to evaluate the performance of the proposed methods. It found that MLE and Bayesian MCMC exhibited similar efficiency, whereas the Boot-P method consistently outperformed Boot-T by demonstrating lower MSE across diverse parameter combinations and sample sizes. To demonstrate the flexibility of the MCD, we analyzed two real-world reliability datasets (Aarset and Meeker-Escobar). In both cases, the MCD demonstrably outperforms competing distributions, evidenced by superior fit statistics (AIC, BIC, HQIC, K-S) and graphical analysis. This superior performance underscores the MCD's ability to provide more accurate and reliable estimates of device lifetimes, MTTF, and failure patterns, ultimately leading to better decision-making in maintenance planning, warranty analysis, and system design.

In future research, it may be beneficial to apply the MCD to censored datasets, such as those involving type-Ⅱ progressive censoring, joint progressive type-Ⅱ censoring schemes, and generalized hybrid censoring schemes. Additionally, the MCD could be adapted for use in accelerated life tests, a method commonly employed in reliability engineering to gather information quickly by accelerating a variable that significantly influences the lifespan of the devices under study.

Appendix

A.   Proofs

Proof of Theorem 3.1. The $r$th noncentral moment of the MCD$(\underline{\mathbf{\psi}})$ can be expressed using Eq (2.3) as follows:

By applying the Taylor expansion for $e^{x}$, the $r$th noncentral moment of the MCD$(\underline{\mathbf{\psi}})$ can be expressed as a linear combination of the moments of the CD, as shown below.

where $\mu'_{m+r, \textbf{CD}}$ represents the $(m+r)$th noncentral moment of the CD.

Proof of Theorem 3.2. Incomplete moments are determined using the following equation:

By applying the series expansions for ${\rm e}^{-\alpha\, {\rm e}^{\gamma\, x}}$ and ${\rm e}^{-\alpha\, {\rm e}^{x^{\lambda}}}$, we derive the following results:

Additionally, by using the series expansions for ${\rm e}^{j\, x^{\lambda}}$ and ${\rm e}^{x^{\lambda}}$, we obtain

After solving the integrals mentioned above, we obtain the result stated in Theorem 3.2.

Proof of Theorem 3.3. The MTTF of the MCD is obtained as

Proof of Theorem 3.4. The MRL of the MCD is expressed as

Proof of Theorem 3.5. The Rényi entropy of $X$ for the MCD$(\underline{\mathbf{\psi}})$ is defined as follows:

By substituting Eq (2.2) into Eq (A.1), we get

Applying the binomial expansion to the function $\Big(\gamma\, {\rm e}^{\gamma\, x}+\lambda\, x^{\lambda-1}\, {\rm e}^{x^{\lambda}}\Big)^{\rho}$ yields:

Utilizing the Taylor series expansion, we obtain the following result:

Solving the integrals above leads to the result presented in Theorem 3.5.

Proof of Theorem 3.6. Consider an ordered sample $\{X_{i}\}_{i = 1}^{n}, \, n\geq1$ from the MCD, with its PDF given by (2.2) and its CDF by (2.1). The PDF of the $l$th order statistic, denoted as $f_{l:n}(x)$, is defined as:

which can be expressed as

Substituting Eqs (2.1) and (2.2) into (A.2) leads to

Consequently,

Here, $f(x; \alpha', \gamma, \lambda)$ denotes the PDF of the MCD with parameters $\alpha' = (n+j+1-l)\, \alpha, \, \gamma, \, \lambda$. By applying Eqs (3.3) and (A.3), the $r$th moment of the $l$th order statistics is derived as shown in (3.6).

B.   Elements of the Fisher information matrix for the MCD

The second partial derivatives of the log-likelihood function for the MCD are as follows:

Use of Generative-AI tools declaration

The author declare that he has not used artificial intelligence tools in the creation of this article.

Acknowledgments

The author thanks the editor and the referees for their valuable comments and suggestions which improved greatly the quality of this paper.

Conflict of interest

There is no conflict of interest declared by the author.