The Chen–Perks Distribution: Properties and Reliability Applications

Abstract

In this paper, a statistical distribution is presented that possesses the ability to describe failure rates exhibiting both monotonic and non-monotonic behaviors, and the bathtub curve, which represents the performance of a device in reliability engineering. The proposed distribution is based on the sum of the hazard functions of the Chen distribution and the Perks distribution, thus presenting the Chen–Perks distribution (CPD). Statistical properties of the CPD focused on reliability engineering are presented to make the model attractive to practitioners of the discipline. The parameters of the CPD were calculated via the maximum likelihood estimator. On the other hand, a comparative analysis was conducted in three study cases to determine the behavior of the CPD relative to other distributions that can describe failure times with the shape of a bathtub curve. The results show that the CPD can offer competitive results, which practitioners can consider when conducting reliability analysis.

1. Introduction

Today, reliability engineering recognizes one of its most well-known concepts in the field: the representation of device behavior over time, thus expressing the stages of the product’s useful life. This concept is related to the bathtub curve, which consists of three parts: infant mortality (decreasing failure rate), normal life (constant failure rate), and wear-out (increasing failure rate). In practice, reliability engineering uses statistical distributions with which it is possible to represent how a device behaves under the time variable. However, the distributions most used in the analysis, such as the Weibull (WD), exponential, and log-normal, do not have the flexibility to represent non-monotonic behavior within their properties, as expressed by the bathtub curve. Those mentioned above can represent a problem in product design and production stages since the warranty and maintenance times are limited or do not approach the real behavior given the limitations exhibited by “classic” reliability distributions.

To solve the problem of accurately representing the lifetimes of devices, alternative statistical distributions and methodologies have been presented, with which it is possible to represent the lifetimes of a device with non-monotonic behavior. Aarset [1] proposed that it is possible to empirically determine the behaviors of the failure times of devices through the total test on time (TTT). With the TTT, it is possible to determine whether the experimental data exhibit increasing, decreasing, or bathtub curve behaviors. With the above, some authors initiated modifications of the most popular distributions in the reliability analysis. For example, in [2,3,4,5], the authors took the properties of mathematical flexibility offered by the WD to establish variants with extra parameters that allow the WD to represent non-monotonic data.

Other methodologies, such as those presented by Mahdavi and Kundu [6], proposed the alpha power transformation (APT), which allows distributions to add one more scale parameter to the base form of the distribution by exponentiating the cumulative distribution function (CDF). Kaushik and Nigam [7] proposed using the DUS transformation; the base distribution is exponentiated by the probability density function (PDF) and the CDF, which enables savings in the calculation time, as it does not introduce any new parameters other than the parameters involved in the baseline distribution. Related works using the two transformations above can be consulted at [8,9,10,11,12,13]. Lee et al. [14] presented another alternative that proposed a generalized use of the Beta distribution to modify the representations of the hazard functions of other popular distributions and obtain behaviors close to the bathtub curve. Specifically, the work focused on the introduction of the Beta–Weibull distribution.

Xie and Lai [15] established that it is possible to obtain statistical distributions with properties of representing data—with bathtub curve behavior—through the additive methodology. This methodology adds the hazard rate functions (HRFs) of two equal or different distributions. Xie introduced the additive Weibull distribution (AWD) with four parameters to evaluate this methodology, obtaining competitive results concerning other distributions. Thach [16] extended Xie’s methodology by proposing a triple sum of hazard functions. For this, Thach established the triple additive WD, which is an excellent fit to the bathtub curve; however, makes it complex as it has too many parameters to estimate. Other distributions based on the additive methodology can be seen in [17,18,19,20,21,22].

One of the most significant problems presented by the previous works is that they do not manage to represent the curve of the bathtub closely since the shapes resemble a “V,” “J,” or “U,” which means that the operational life of the device is short. Therefore, the information is biased, meaning the device estimates are inadequate. Thus, there is a need to explore alternatives in methodologies or hybrid distributions that establish behaviors closer to that marked by the assumptions of the bathtub curve.

Chen [23] proposed a distribution of only two parameters, with which it is possible to represent non-monotonic data, consequently forming the Chen distribution (ChD). One of its characteristics is given its mathematics; the shape parameter has confidence intervals with a closed form. The weakness of ChD lies in the poor flexibility of the model, as it lacks a scale parameter, which has been one of the impediments when considering the ChD as an alternative in reliability analyses. Thanh Thach and Briš [24] and Khan et al. [25] modified the ChD to analyze reliability and survival. The proposals were based on WD due to the popularity and mathematical flexibility the WD can offer to the ChD by taking a scale parameter from the WD. In both cases, the Chen–Weibull distribution (CWD) was proposed. At the same time, the use of distributions in the actuarial area has become popular because the behavior exhibited by a living being is similar to a bathtub shape. The Perks distribution (PD), proposed by Perks [26], has been used with some classical reliability distributions with satisfactory results used to describe lifetimes. Zeng et al. [27] proposed Perks4 and Perks5, with which it is possible to determine non-monotonic behaviors close to the bathtub curve by combining the HRF. Singh [28] proposed the additive Perks–Weibull (APW) distribution with applications in reliability. The APW takes advantage of the PD representation before non-monotonic data, making them more flexible with the WD. Méndez-González et al. [29] proposed the additive Perks distribution (APD), which involves four parameters. For this, Méndez modified the failure rate function of one of the PDs to provide flexibility in the wear-out stage of the device under analysis. Méndez-González et al. [30] extended the PD through the APT, thus proposing the alpha-exponentiated Perks distribution (AEXP). In studies with real data from the latest generation medical devices, the AEXP is shown to be a valid alternative to the reliability analysis within electronic devices, which are more associated with the behavior of the bathtub curve.

Therefore, this paper proposes a new statistical distribution based on the ChD and PD, i.e., the Chen–Perks distribution (CPD). The motivation for this new distribution lies in the following.

• We present an alternative statistical distribution with reliability applications that possess the ability to describe non-monotonic behaviors, such as those exhibited by the bathtub curve.
• The new CPD distribution, which combines the ChD and PD, is significantly more flexible than the current hybrid bathtub distributions presented in the literature review.
• We provide the ChD with a scale parameter through the PD, with which it is possible to establish life-stress models. By adding this scale parameter, CPD can be used within accelerated life tests (ALTs), offering practical modeling for reliability engineers.
• We establish an attractive distribution for reliability engineers to conduct analyses, considering the benefits of modeling formed from an actuarial point of view.

The CPD is based on the sum of the HRFs of the ChD and PD. By employing additive methodology, the individual statistical properties of distributions in question are brought together to represent diverse types of behavior, such as the bathtub curve. To complement the modeling of the CPD, reliability-oriented statistical and application properties are analyzed and presented. The parameters of the CPD were calculated via the maximum likelihood estimator (MLE) implemented in RStudio. To evaluate the results from the proposed model, the CPD was compared with other distributions used in reliability and with comparable properties to the CPD in studies with data that have non-monotonic behaviors. To issue a recommendation to readers about the models compared, we considered the parameter estimation, the Akaike information criteria (AIC), the Bayesian information criteria (BIC), the Kolmogorov–Smirnov test (K-S), and the p-value.

Finally, this paper is organized as follows. Section 2 presents the general equations of the CPD for the reliability analysis. Section 3 shows the measures of the central tendency of APD. Section 4 presents the moments and incomplete moments. Section 5 presents the order statistics. In Section 6, the mean residual lifetime function is estimated. Section 7 presents the Rényi entropy. Section 8 presents the likelihood function to calculate the parameters proposed in Section 2. Section 9 presents the case studies of the paper. The last section provides the concluding remarks and future work about the proposed model.

2. The CPD Model Construction

Let $h_{ChD}\left(x\right)=\beta \theta x^{\theta -1}{e}^{x^{\theta }}$ and $h_{PD}\left(x\right)=\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}$ be the HRF of ChD and PD, respectively. Hence, to form the new HRF of CPD based on the additive methodology, it is defined as follows:

$$h_{CPD}\left(x\right)=h_{ChD}\left(x\right)+h_{PD}\left(x\right)$$

(1)

where $x>0$, $\beta >0,\theta >0$ and $\alpha >0$ are the shape parameters and $\lambda >0$ is the scale parameter. In Figure 1, different representations of the HRF of the CPD are shown, considering what was obtained in Equation (1).

The PDF of the CPD can be defined from Equation (1) and is described as follows:

$$f_{CPD}\left(x\right)=h_{CPD}\left(x\right)·exp\left({\int }_0^x{h}_{CPD}\left(v\right)dv\right)$$

$$=\left(\beta \theta x^{\theta -1}{e}^{x^{\theta }}+\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}\right)·\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}$$

(2)

Figure 2 shows some shapes the PDF takes for some values.

Other Reliability Equation

Based on Equation (2), it is possible to determine other essential equations to determine the device’s reliability in normal operating conditions.

The survival function $S\left(x\right)$ of CPD from Equation (2) can be written as follows:

$$S_{CPD}\left(x\right)=\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}$$

(3)

The CDF $F\left(x\right)$ of CDP can be calculated as follows:

$$F_{CPD}\left(x\right)=1-\left[\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}\right]$$

(4)

Finally, the cumulative hazard function $H\left(x\right)$, from Equation (3) can be written as follows:

$$H_{CPD}\left(x\right)=\beta \left(-1+e^{x^{\theta }}\right)-ln\left(1+\alpha \right)+ln\left(1+\alpha e^{\lambda x}\right)$$

(5)

3. CPD Measures of Central Tendency

This section establishes some statistical elements that can complement the reliability analysis.

3.1. Quantile

The p-$th$ quantile $q_p$ of the random variable x for the CPD based on Equation (4) can be written as follows:

$$p=1-\left[\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{q^{\theta }}-1\right)}}{1+\alpha e^{\lambda q}}\right]$$

(6)

As Equation (6) can be appreciated, it does not have a closed-form solution. This non-solution makes it necessary to make approximations through numerical methods programmed in specialized software.

3.2. Mode

The mode of the CPD can be obtained from the first derivative of the PDF calculated in Equation (2), which must be equal to zero. By performing the operation marked in the previous statement, we obtain the following:

$$\frac{\left(1+\alpha \right)·\theta \left\{\alpha \left[a·e^{\lambda x}+b\right]·\beta \left(1+\alpha e^{\lambda x}\right)e^{x^{\theta }}{x}^{\theta -1}+e^{\lambda x}\left[c\right]\lambda x\right\}{e}^{-\beta \left(e^{x^{\theta }}-1\right)}}{x{\left[1+\alpha ·e^{\left(\lambda x\right)}\right]}^3}=0$$

(7)

where

$$\begin{array}{c}\hfill a=\theta x^{\theta }\left(\beta e^{x^{\theta }}-1\right)+\lambda x-\theta +1\\ \hfill b=\theta \left(\beta x^{\theta }{e}^{x^{\theta }}-x-1\right)+1\\ \hfill c=\alpha \left\{\beta \theta x^{\theta }{e}^{x^{\theta }}+\lambda x\right\}+\beta \theta x^{\theta }{e}^{x^{\theta }}-\lambda x\end{array}$$

As can be seen, Equation (7) does not have a closed solution, so it is necessary to use numerical methods programmed in specialized software to approximate a solution of the mode.

4. Moments, Moment-Generating Function, and Incomplete Moments

This section calculates the moments and the incomplete moments. These functions can determine some aspects of the device, such as the mean time between failure (MTTF), the variance, skewness, kurtosis, and the Bonferroni and Lorenz curves.

4.1. Moments

Theorem 1. 

Let X be a continuous random variable with the CPD with $\alpha ,\beta ,\lambda ,\theta \in \mathbb{R}>0$. The $rth$ moment of X is then given by:

$$E\left(X^r\right)=r\left(1+\alpha \right)\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}\mathsf{\Gamma }\left(\frac{r+𝚥+k}{\theta k}\right)$$

(8)

where $\mathsf{\Gamma }\left(.\right)$ represents the gamma function.

Proof. 

The reliability function S(x) can calculate the moment function. Mathematically, moments are defined as follows:

$$E\left(X^r\right)={\int }_0^{\infty }r·x^{r-1}s\left(x\right)dx$$

By substituting Equation (3), in the above equation, the following is obtained:

$$E\left(X^r\right)=r\left(1+\alpha \right)·{\int }_0^{\infty }{x}^{r-1}{e}^{-\beta \left(e^{x^{\theta }}-1\right)}·{\left(1+\alpha e^{\lambda x}\right)}^{-1}dx$$

(9)

Equation (9) can be solved by series expansion. For this, we define the following series:

$$\begin{array}{c}\hfill {\left(1+\alpha e^{\lambda x}\right)}^{-1}=\sum _{𝚤=0}^{\infty }{\left(-\alpha \right)}^{𝚤}{e}^{\lambda 𝚤x}.\\ \hfill e^{\lambda 𝚤x}=\sum _{𝚥=0}^{\infty }\frac{{\left(\lambda 𝚤\right)}^{𝚥}}{\mathsf{\Gamma }\left(𝚥+1\right)}\end{array}$$

By developing the integral by the series defined above, we obtain the following:

$$E\left(X^r\right)=r\left(1+\alpha \right)\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}{\int }_0^{\infty }{x}^{r+𝚥+k-1}{e}^{x^{\theta k}}dx$$

□

Finally, the remaining integral can be solved by defining the gamma function. Therefore, when applying said definition, it is possible to arrive at what is obtained in Equation (8).

MTTF, Variance, Skewness, and Kurtosis of CPD

Based on the obtained in Theorem 1, it is possible to obtain other equations that are very useful within the reliability analysis.

The MTTF or the mean can be obtained from the first moment of Equation (9). When performing the relevant calculation, the MTTF can be expressed as follows:

$$E\left(X\right)=MTTF=\left(1+\alpha \right)\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)$$

(10)

As is well known, the variance of a PDF can be defined as follows:

$$Var\left(X\right)=E\left(X^2\right)-{\left(E\left(X\right)\right)}^2$$

Therefore, by taking Equation (10), the variance of the CPD can be estimated as follows:

$$Var\left(X\right)=(1-{\alpha }^2)·\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\left[\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{2+𝚥+k}{\theta k}\right)\right)-{\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 1\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)\right)}^2\right]$$

(11)

The skewness of a PDF can be estimated as follows:

$${\gamma }_1=\frac{E\left(X^3\right)-3\mu {\sigma }^2-{\mu }^3}{{\sigma }^3}$$

If we consider the following:

$$E{\left(x\right)}^3=3(1+\alpha )·\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\left[\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{3+𝚥+k}{\theta k}\right)\right]$$

and taking into account the $E\left(x\right)$ and $Var\left(x\right)$ obtained in Equations (10) and (11), the skewness of CPD is given by:

$${\gamma }_1=\sum _{\mathit{ı}=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\left[\frac{\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{3+𝚥+k}{\theta k}\right)-3\left[\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)·\mathsf{\Gamma }\left(\frac{2+𝚥+k}{\theta k}\right)\right]+2{\left[\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)\right]}^3}{{\left[\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{2+𝚥+k}{\theta k}\right)\right)-{\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 1\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)\right)}^2\right]}^{3/2}}\right]$$

(12)

Finally, the kurtosis of a PDF is defined mathematically as follows:

$${\gamma }_2=\frac{E\left(X\right){}^4}{\sigma {}^4}$$

If we consider the following:

$$E{\left(X\right)}^4=4(1+\alpha )·\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\left[\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{4+𝚥+k}{\theta k}\right)\right]$$

and take the obtained in Equation (11), the kurtosis of CDP can be estimated as follows:

$${\gamma }_2=4(1+\alpha )·\sum _{𝚤=0}^{\infty }\sum _{𝚥=0}^{\infty }\sum _{k=0}^{\infty }\left[\frac{\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }(𝚥+1)\mathsf{\Gamma }(k+1)}·\mathsf{\Gamma }\left(\frac{4+𝚥+k}{\theta k}\right)}{{\left[\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{2+𝚥+k}{\theta k}\right)\right)-{\left(\frac{-{\alpha }^{𝚤}{\left(\lambda 1\right)}^{𝚥}{e}^{\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{𝚥+k}{\theta k}\right)\right)}^2\right]}^2}\right]$$

(13)

4.2. Incomplete Moments

Theorem 2. 

Let X be a continuous random variable with the CPD with $\alpha ,\beta ,\lambda ,\theta \in \mathbb{R}>0$. The incomplete moment of X is then given by:

$$H_s\left(y\right)=\left(1+\alpha \right)·\sum _{𝚤=𝚥=0}^{\infty }\sum _{k=𝓁=0}^{\infty }\sum _{m=n=0}^{\infty }\frac{2^{𝚤}{\left(-\alpha \right)}^{𝚤}{\left(-1\right)}^{2k}{\alpha }^m{\beta }^k{\left(𝚤𝓁\right)}^{𝚥}{k}^{𝓁}{\left(𝓁m\right)}^n}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)\mathsf{\Gamma }\left(𝓁+1\right)}\mathsf{\Gamma }\left[\frac{𝚥-\beta 𝓁+n+\theta +s}{\theta },y^{\theta }\right]$$

(14)

Proof. 

The incomplete moments are calculated from the following equation:

$$\begin{array}{c}\hfill H_s\left(y\right)={\int }_0^y{x}^{s}f\left(x\right)dx\\ \hfill ={\int }_0^y{x}^s·\left(\beta \theta x^{\theta -1}{e}^{x^{\theta }}+\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}\right)·\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}dx\\ \hfill =\left(1+\alpha \right){\int }_0^y{x}^s·\left(\frac{e^{-\beta \left(e^{x^{\theta }}-1\right)}·\left[\beta \theta e^{x^{\theta }}\left(1+\alpha e^{\lambda x}\right)x^{\theta -1}+\alpha e^{\lambda x}\right]}{{\left(1+\alpha e^{\lambda x}\right)}^2}\right)dx\end{array}$$

By solving the previous equation through series, we can obtain the following:

$$\begin{array}{c}\hfill H_s\left(y\right)=\left(1+\alpha \right)·\sum _{𝚤=𝚥=0}^{\infty }\frac{{\left(-2\alpha \right)}^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}}{\mathsf{\Gamma }\left(𝚥+1\right)}{\int }_0^y{e}^{-\beta \left(e^{x^{\theta }}-1\right)}\left[\beta \theta e^{x^{\theta }}\left(1+\alpha e^{\lambda x}\right)x^{\theta -1}+\alpha e^{\lambda x}\right]·x^{s+𝚥}dx\\ \hfill =\left(1+\alpha \right)·\sum _{𝚤=𝚥=0}^{\infty }\sum _{k=𝓁=0}^{\infty }\frac{{\left(-2\alpha \right)}^{𝚤}{\left(\lambda 1\right)}^{𝚥}{\left(-1\right)}^{2k}{\beta }^k{k}^{𝓁}{e}^{-\beta }}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)\mathsf{\Gamma }\left(𝓁+1\right)}{\int }_0^y\left[\beta \theta e^{x^{\theta }}\left(1+\alpha e^{\lambda x}\right)x^{\theta -1}+\alpha e^{\lambda x}\right]x^{s+𝚥-𝓁\beta }dx\\ \hfill =\left(1+\alpha \right)·\sum _{𝚤=𝚥=0}^{\infty }\sum _{k=𝓁=0}^{\infty }\sum _{m=n=0}^{\infty }\frac{2^{𝚤}{\left(-\alpha \right)}^{𝚤}{\left(-1\right)}^{2k}{\alpha }^m{\beta }^k{\left(𝚤𝓁\right)}^{𝚥}{k}^{𝓁}{\left(𝓁m\right)}^n}{\mathsf{\Gamma }\left(𝚥+1\right)\mathsf{\Gamma }\left(k+1\right)\mathsf{\Gamma }\left(𝓁+1\right)}{\int }_0^y\left[\beta \theta e^{x^{\theta }}{x}^{\theta -1}+\alpha e^{\lambda x}\right]x^{s+𝚥-𝓁\beta +n}dx\end{array}$$

By solving the last integral through the gamma function, we can arrive at the result obtained in Theorem 2. □

5. Order Statistics

Order statistics are usually used in reliability analyses for devices with redundant systems. This class of systems requires ensuring full operation in the event of different failures, which is why order statistics is commonly used in critical systems. Order statistics are defined mathematically as follows:

$$f_{k:n}\left(x\right)=n\left(\genfrac{}{}{0pt}{}{n-1}{k-1}\right)F{\left(x\right)}^{k-1}{\left[1-F\left(x\right)\right]}^{n-k}f\left(x\right)$$

(15)

By taking the results obtained in Section 2 and substituting them in Equation (13), the order statistics of the CPD can be written as follows:

$$\begin{array}{c}\hfill f_{k:n}\left(x\right)=n\left(\genfrac{}{}{0pt}{}{n-1}{k-1}\right){\left\{1-\left[\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}\right]\right\}}^{k-1}{\left\{\left[\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}\right]\right\}}^{n-k}\\ \hfill \left(\beta \theta x^{\theta -1}{e}^{x^{\theta }}+\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}\right)·\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}\end{array}$$

By developing the series and applying algebraic reductions to the previous equation, the order statistics of the CPD are described as follows:

$$f_{k:n}\left(x\right)=n\left(\genfrac{}{}{0pt}{}{n-1}{k-1}\right)h\left(x\right)·\sum _{𝚤=0}^{n-k}\sum _{𝚥=0}^{k-1}{\left[{\left(k-𝚤\right)}^{𝚥}·S\left({x_{𝚤}}_{𝚥}\right)\right]}^{𝚤+𝚥}$$

(16)

6. Mean Residual Lifetime

In scenarios where the expected lifetime for a given device is computed, a time t as a function of t is denoted as MRL. This information makes it possible to characterize the operation of the devices under normal operating conditions. The MRL is written as follows:

$$M\left(t\right)=E\left(x-t|x+t\right)=\frac{1}{s\left(x\right)}·{\int }_0^{\infty }s\left(x+t\right)dx.$$

(17)

By taking Equation (3) and substituting it into Equation (17), the following is obtained:

$$M\left(t\right)=\frac{1+\alpha e^{\lambda x}}{e^{-\beta \left(e^{x^{\theta }}-1\right)}}·{\int }_0^{\infty }{e}^{-\beta \left(e^{{\left(x+t\right)}^{\theta }}-1\right)}{\left(1+\alpha e^{\lambda \left(x+t\right)}\right)}^{-1}dx.$$

By developing the series and solving the integral, the MRL of the CPD can be approximated at the rate of:

$$M\left(t\right)=\frac{1+\alpha e^{\lambda t}}{e^{-\beta \left(e^{t^{\theta }}-1\right)}}·\sum _{𝚤=𝚥=0}^{\infty }\frac{-{\alpha }^{𝚤}{\left(\lambda 𝚤\right)}^{𝚥}}{\mathsf{\Gamma }\left(𝚥+1\right)}·E(-\beta e^t+\theta ),$$

(18)

Mills Ratio

The Mills ratio (ML) is a unique technique for describing reliability because of its connection to the failure rate. The ML for the CPD can be expressed as follows:

$$M\left(x|\beta ,\theta ,\lambda ,\alpha \right)=\frac{s\left(x|\beta ,\theta ,\lambda ,\alpha \right)}{f\left(x|\beta ,\theta ,\lambda ,\alpha \right)}$$

$$M\left(x|\beta ,\theta ,\lambda ,\alpha \right)=\frac{x\left(1+\alpha e^{\lambda x}\right)}{x^{\theta }{e}^{x^{\theta +\lambda x}}\alpha \beta \theta +\beta \theta x^{\theta }{e}^{x^{\theta }}+\alpha \lambda e^{\lambda x}}$$

(19)

The odd function (OF) is defined by:

$$\widehat{O}\left(x|\beta ,\theta ,\lambda ,\alpha \right)=\frac{F\left(x|\beta ,\theta ,\lambda ,\alpha \right)}{s\left(x|\beta ,\theta ,\lambda ,\alpha \right)}$$

The OF of x for CPD is given by:

$$\widehat{O}\left(x|\beta ,\theta ,\lambda ,\alpha \right)=\frac{\left(-\alpha -1\right)e^{-\beta \left(e^{x^{\theta }}-1\right)}+\alpha e^{\lambda x}+1}{\left(1+\alpha \right)e^{-\beta \left(e^{x^{\theta }}-1\right)}}$$

(20)

7. Rényi Entropy

Entropy within the reliability analysis quantifies a system’s diversity, uncertainty, or randomness. One of the applications is associated with life-stress models that are obtained from accelerated life (ALT) analyses.

In this manuscript, we will discuss the Rényi entropy, which is defined mathematically as follows:

$$R\left(x\right)=\frac{1}{1-p}log{\int }_0^{\infty }f{\left(x\right)}^{p}dx<0<p<1$$

(21)

By taking Equation (2) and substituting it in Equation (21), the following is obtained:

$$R\left(x\right)=\frac{1}{1-p}log{\int }_0^{\infty }{\left[\left(\beta \theta x^{\theta -1}{e}^{x^{\theta }}+\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}\right)·\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x^{\theta }}-1\right)}}{1+\alpha e^{\lambda x}}\right]}^{p}dx$$

In the first instance, we develop $f{\left(x\right)}^p$. For this, we will use the following series ${\left(a+b\right)}^n={\sum }_{𝚤=0}^p{a}^{p-𝚤}{b}^{𝚤}$, thus

$$\begin{array}{c}\hfill f{\left(x\right)}^p=\sum _{𝚤=0}^p\left(\genfrac{}{}{0pt}{}{p}{𝚤}\right){\left[\beta \theta x^{\theta -1}{e}^{x^{\theta }}\right]}^{p-𝚤}·{\left(\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}\right)}^{𝚤}·{\left(1+\alpha \right)}^p·e^{-\beta p\left(e^{x^{\theta }}-1\right)}·{\left(1+\alpha e^{\lambda x}\right)}^{-p}\\ \hfill =\sum _{𝚤=0}^p\sum _{𝚥=k=0}^{\infty }\sum _{𝓁=m=0}^{\infty }\left(\genfrac{}{}{0pt}{}{p}{𝚤}\right)\frac{{\left[\beta \left(p-𝚤\right)\theta \right]}^{𝚥}{𝚤}^{-p-𝚤+𝓁}{\left(\lambda -p-𝚤+𝓁\right)}^m{\left(1+\alpha \right)}^p}{\mathsf{\Gamma }\left(k+1\right)}·{\int }_0^{\infty }{\left(\alpha \lambda e^{\lambda x}\right)}^{𝚤}·e^{-\beta p\left(e^{x^{\theta }}-1\right)}·x^{𝚥+\theta +k\theta +m-1}dx\end{array}$$

Solving the last integral, we obtain the Rényi entropy for the CPD:

$$\begin{array}{c}\hfill R\left(x\right)=\frac{1}{1-p}log\left\{\sum _{𝚤=0}^p\sum _{𝚥=k=𝓁=0}^{\infty }\sum _{m=n=s=0}^{\infty }\left(\genfrac{}{}{0pt}{}{p}{𝚤}\right)\frac{{\left[\beta \left(p-𝚤\right)\theta \right]}^{𝚥}{𝚤}^{-p-𝚤+𝓁}{\left(\lambda -p-𝚤+𝓁\right)}^m{\left(1+\alpha \right)}^p{\alpha }^{𝚤+n}{\left(\lambda 1\right)}^n{\left(\beta p\right)}^s}{\mathsf{\Gamma }\left(k+1\right)}·\mathsf{\Gamma }\left(\frac{𝚥+\theta \left(k+1\right)+m+n+s}{s\theta }\right)\right\}\end{array}$$

(22)

8. Parameter Estimation for CPD

This section establishes equations with which the parameters of the CPD can be estimated. For each parameter estimation of the CPD, the MLE technique is used. Let $x_1,x_2\dots x_m$ be a random sample from the lifetime distribution with PDF $f\left(x\right)$ based on a sample of size m. The likelihood function of Equation (2) is written as follows:

$$L=\prod _{𝚤=0}^m\left[\left(\beta \theta x_{𝚤}^{\theta -1}{e}^{x_{𝚤}^{\theta }}+\frac{\alpha \lambda e^{\lambda x_{𝚤}}}{1+\alpha e^{\lambda x_{𝚤}}}\right)·\frac{\left(1+\alpha \right)·e^{-\beta \left(e^{x_{𝚤}^{\theta }}-1\right)}}{1+\alpha e^{\lambda x_{𝚤}}}\right]$$

By obtaining the natural algorithm from the equation described above, the log-likelihood $\mathsf{\Lambda }$ function becomes:

$$\mathsf{\Lambda }=m·ln\left(1+\alpha \right)·\sum _{𝚤=0}^{m}ln\left[\beta \theta x_{𝚤}^{\theta -1}{e}^{x_{𝚤}^{\theta }}+\frac{\alpha \lambda e^{\lambda x_{𝚤}}}{1+\alpha e^{\lambda x_{𝚤}}}\right]-\beta ·\sum _{𝚤=0}^m\left[e^{x_{𝚤}^{\theta }}-1\right]-\sum _{𝚤=0}^{m}ln\left[1+\alpha e^{\lambda x_{𝚤}}\right]$$

(23)

As is known, the MLE needs the partial derivatives of Equation (23) concerning each parameter defined in the CPD. Given the above, the estimation of the parameters $\beta ,\theta ,\alpha ,\lambda$ are defined as follows:

$$\frac{\partial \mathsf{\Lambda }}{\partial \beta }=m+\sum _{𝚤=1}^m\left[\frac{\theta x_{𝚤}^{\theta -1}{e}^{x_{𝚤}^{\theta }}\left(1+\alpha e^{\lambda x_{𝚤}}\right)}{\beta \theta x_{𝚤}^{\theta -1}{e}^{x_{𝚤}^{\theta }}\left(1+\alpha e^{\lambda x_{𝚤}}\right)+\alpha \lambda e^{\lambda x_{𝚤}}}-e^{x_{𝚤}^{\theta }}\right]$$

(24)

$$\frac{\partial \mathsf{\Lambda }}{\partial \theta }=\sum _{𝚤=1}^m\left[\frac{\beta \theta x_{𝚤}^{\theta -1}ln\left(x_{𝚤}\right)e^{x_{𝚤}^{\theta }}+\theta x_{𝚤}^{2\theta -1}\beta ln\left(x_{𝚤}\right)e^{x_{𝚤}^{\theta }}+\beta x_{𝚤}^{\theta -1}{e}^{x_{𝚤}^{\theta }}}{h\left(x_{𝚤}\right)}-\beta x_{𝚤}^{\theta }ln\left(x_{𝚤}\right)e^{x_{𝚤}^{\theta }}\right]$$

(25)

$$\frac{\partial \mathsf{\Lambda }}{\partial \alpha }=\frac{m}{1+\alpha }-\sum _{𝚤=1}^m\left[\frac{e^{\lambda x_{𝚤}}\left(e^{\lambda x_{𝚤}}h\left(x_{𝚤}\right)\alpha +h\left(x_{𝚤}\right)-\lambda \right)}{{\left(1+\alpha e^{\lambda x_{𝚤}}\right)}^{2}h\left(x_{𝚤}\right)}\right]$$

(26)

$$\frac{\partial \mathsf{\Lambda }}{\partial \lambda }=\sum _{𝚤=1}^m\left[\frac{\left(-e^{2\lambda x_{𝚤}}{x}_{𝚤}\alpha \lambda +e^{\lambda x_{𝚤}}\left(\lambda x_{𝚤}-x_{𝚤}h\left(x_{𝚤}\right)+1\right)\left(1+\alpha e^{\lambda x_{𝚤}}\right)\right)\alpha }{{\left(1+\alpha e^{\lambda x_{𝚤}}\right)}^{2}h\left(x_{𝚤}\right)}\right]$$

(27)

In turn, the elements of the Fisher matrix can be estimated by obtaining the partial derivatives of each parameter concerning the results obtained in Equation (24) to Equation (27). These elements are presented in Appendix A.

9. Case Study

In this section, three case studies are shown, where the CPD is put to the test and compared with other hybrid distributions. For this, the following are considered:

• Distributions with different construction methodologies are contemplated to establish a representative analysis. The distributions considered are AWD [15], the modified Weibull extension (MWE) [4], Perks4 [27], APW [28], APD [29], AEXP [30], and CWD [24,31]
• The case studies do not consider censored data of any kind.
• The parameter estimation of the distributions presented in

  Table 1. HRF of some recent Bathtub-shaped models used for all case studies.

  Statistical Distribution	h(x)
  AWD	$\alpha \lambda x^{\lambda -1}+\beta \theta x^{\theta -1}$
  MWE	$\alpha (\beta +\lambda x)x^{\beta -1}{e}^{\lambda x}$
  Perks4	$\frac{\theta +e^{\left(\beta x+\alpha \right)}}{1-e^{\left(-\beta x+\lambda \right)}}$
  APW	$\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}+\theta \beta x^{\beta -1}$
  APD	$\frac{\alpha \lambda e^{\lambda x}}{1+\alpha e^{\lambda x}}+\frac{\beta \theta e^{-\theta x}}{1+\beta e^{-\theta x}}$
  AEXP	$\frac{ln\left(\alpha \right)\beta \lambda \theta \left(1+\beta \right)e^{\lambda x}{\left(1-\frac{1+\beta }{1+\beta e^{\lambda x}}\right)}^{\theta -1}{\alpha }^{{\left(1-\frac{1+\beta }{1+\beta e^{\lambda x}}\right)}^{\theta }}}{{\left(1+\beta e^{\lambda x}\right)}^2\alpha \left(1-{\alpha }^{{\left(1-\frac{1+\beta }{1+\beta e^{\lambda x}}\right)}^{\theta }-1}\right)}$
  CWD	$\alpha \beta x^{\beta -1}{e}^{x^{\beta }}+\lambda \theta x^{\theta -1}$

  was performed in RStudio through the MaxLik library. For this, the gradient and Hessian of each of the distributions are considered.

9.1. Case Study 1. Reliability Analysis of Eighteen Electronic Devices

This case study focuses on the reliability analysis of eighteen electronic devices. Wang et al. [32] presented the lifetimes; the above-mentioned data can be consulted in

Table 2. The lifetime of eighteen devices by Wang et al. [32].

Lifetimes
5	11	21	31	46	75	98	122	145
165	196	224	245	293	321	330	350	420

.

The analysis of this case study begins by empirically verifying the behaviors of the data presented in Table 1. For this, we will use the TTT plot; the results of this graph can be seen in Figure 3. The results from the TTT plot of case study 1 suggest that the data have non-monotone behavior.

Table 3. Parameter values and statistics estimated for the lifetime data presented byWang et al. [32].

Model	Parameters				Statistics
	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\lambda }$	$\mathit{\theta }$	Loglik	AIC	BIC	K-S	AD	CVM
CPD	$2.606\times {10}^{-5}$ (1.874 $\times {10}^{-4})$	$8.692\times {10}^{-3}$ (1.08 $\times {10}^{-2})$	$3.068\times {10}^{-2}$ (1.956 $\times {10}^{-2})$	0.874 (0.237)	$-107.441$	$222.882$	$226.443$	0.044 (0.975)	0.068 (0.983)	0.024 (0.965)
AWD	$4.366\times {10}^{-7}$ (4.147 $\times {10}^{-6})$	$1.521\times {10}^{-2}$ (1.745 $\times {10}^{-2})$	2.578 (1.457)	0.714 (0.301)	$-109.451$	$226.902$	$230.463$	0.101 (0.801)	0.257 (0.811)	0.0336 (0.800)
MWE	13.472 (9.124)	0.745 (0.145)	$2.541\times {10}^{-3}$ (7.114 $\times {10}^{-4})$	$---$	$-111.478$	$228.956$	$231.627$	0.161 (0.672)	0.411 (0.655)	0.062 (0.682)
PERKS4	−9.282 (2.225)	$1.441\times {10}^{-2}$ (2.225 $\times {10}^{-3})$	−0.422 (1.022 $\times {10}^{-2})$	$2.821\times {10}^{-3}$ (3.221 $\times {10}^{-5})$	$-108.125$	$224.250$	$227.811$	0.095 (0.945)	0.117 (0.951)	0.011 (0.942)
APW	$1.012\times {10}^{-4}$ (2.014 $\times {10}^{-3})$	0.799 (0.325)	$2.722\times {10}^{-2}$ (4.321 $\times {10}^{-2})$	$1.322\times {10}^{-2}(1.722$)	$-108.201$	$224.402$	$227.963$	0.059 (0.950)	0.075 (0.960)	0.035 (0.946)
APD	$0.170$ (3.868 $\times {10}^{-2})$	$0.221$ (6.144 $\times {10}^{-2})$	$7.362\times {10}^{-3}$ (5.172 $\times {10}^{-3})$	$5.532\times {10}^{-2}$ (2.962 $\times {10}^{-2})$	$-108.373$	$224.746$	$228.308$	0.081 (0.948)	0.119 (0.953)	0.0113 (0.939)
AEXP	2.627 (0.614)	0.990 (0.114)	$8.623\times {10}^{-3}$ (1.772 $\times {10}^{-3})$	0.807 (0.409)	$-109.963$	$227.925$	$231.488$	0.401 (0.324)	0.511 (0.344)	0.381 (0.301)
CWD	$6.884\times {10}^{-2}$ (4.25 $\times {10}^{-3})$	$0.228$ (1.242 $\times {10}^{-3})$	$0.224$ (2.335 $\times {10}^{-2})$	$4.638\times {10}^{-2}$ (4.221 $\times {10}^{-5})$	$-110.544$	$229.088$	$232.650$	0.131 (0.825)	0.354 (0.876)	0.052 (0.866)

shows the results for each distribution in the comparative analysis. In the parameter estimation part, the values shown in parentheses reflect the standard errors of the parameters. The values shown in parentheses show the p-values for the criteria in the statistical part.

Figure 4 shows the reliability graphs relevant to the analysis of the times presented in Table 2. The PDF, reliability, hazard, and cumulative hazard graphs are computed based on the estimates obtained from each parameter in Table 3. The results presented in Figure 4 show the behaviors of the data in Table 2 under the different distributions analyzed. In the first instance, Figure 4a shows the PDF of each statistical distribution and how it fits the histogram of the lifetimes. As can be seen, the CPD offers an extremely competitive fit relating to the histogram, which means that the CPD expresses the real behavior of the device in the field more closely. This aspect is essential for practitioners since the risk of information loss or bias is minimal concerning other distributions with properties similar to the CPD. Figure 4b shows how lifetimes behave when analyzed from the perspective of the reliability graph. At this point, we can observe that the CPD concerning the empirical reliability of Kaplan–Meier is associated with more points relating to the other distributions under analysis. From a practical point of view, the CPD will be able to more closely determine the MTTF that the product will exhibit under the designed operating conditions. This last affirmation can be supported by the estimations of Table 2 and the first moment of each distribution. Figure 4c shows how device failures behave under the different distributions. Figure 4c shows the lifetimes of case study one. The results show that the CPD offers a representation closer to the bathtub curve than the other distributions with comparable properties, showing that non-monotonic behaviors exist. However, the shapes they take are irregular. The conclusion that a practitioner can obtain at this point is that the CPD can be a good option to determine maintenance strategies closer to the real behavior of the piece under analysis. Figure 4d shows the cumulative hazard plot, confirming consistency between the model (the CPD) and the lifetimes established in this case. Further, Figure 4e presents the MRL plots of the fitted distributions and the non-parametric MRL. We learn that the MRL curve of the proposed CPD has the best non-parametric estimate match.

Based on the mathematical and graphical evidence presented in Table 3 and Figure 4, it can be determined that the CPD more closely describes the device’s behavior in operating environments. The competitive advantage that CPD offers for this specific device is that failures are described in the form of a bathtub curve, which allows practitioners to draw better conclusions to improve the performance of the device under analysis.

9.2. Case Study 2. Reliability Analysis for Aarset Data

This case study used the data obtained by Aarset [1]. The data consisted of analyzing the lifespan behaviors of fifty devices.

Table 4. The lifespan behaviors of eighteen devices by Aarset [1].

Lifetimes
0.1	0.2	1	1	1	1	1	2	3	6
7	11	12	18	18	18	18	18	21	32
36	40	45	46	47	50	55	60	63	63
67	67	67	67	72	75	79	82	82	83
84	84	84	85	85	85	85	85	86	86

shows the lifespan behaviors of the devices under analysis. As in the previous case study, Figure 5 shows the TTT plot of the data presented by Aarset. The behaviors observed in the Aarset data suggest that the devices under analysis exhibit a behavior similar to that of a bathtub curve.

Table 5. Parameter values and statistics estimated for the lifetime data presented by Aarset [1].

	Parameters				Statistics
Model	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\lambda }$	$\mathit{\theta }$	Loglik	AIC	BIC	K-S	AD	CVM
CPD	$1.214\times {10}^{-16}$ (1.874 $\times {10}^{-4})$	$4.895\times {10}^{-2}$ (1.08 $\times {10}^{-2})$	$0.436$ (1.956 $\times {10}^{-2})$	0.258 (0.101)	$-201.222$	$410.444$	$418.092$	0.055 (0.978)	0.065 (0.988)	0.030 (0.960)
AWD	$8.451\times {10}^{-9}$ (1.498 $\times {10}^{-9})$	$9.121\times {10}^{-2}$ (3.821 $\times {10}^{-2})$	$4.279$ (4.912 $\times {10}^{-2})$	$0.456$ (9.831 $\times {10}^{-2})$	$-221.441$	$450.874$	$458.974$	0.112 (0.451)	0.571 (0.466)	0.271 (0.449)
MWE	13.736 (5.088)	$0.588$ (7.814 $\times {10}^{-2})$	$9.014\times {10}^{-3}$ (3.141 $\times {10}^{-2})$	$---$	$-231.657$	$469.296$	$475.032$	0.155 (0.175)	1.424 (0.165)	0.491 (0.170)
PERKS4	−71.432 (9.124)	0.839 (0.111)	$-6.211\times {10}^{-2}$ (1.321 $\times {10}^{-4})$	$1.534\times {10}^{-2}$ (1.194 $\times {10}^{-3})$	$-205.610$	$419.112$	$424.824$	0.079 (0.901)	0.083 (0.912)	0.075 (0.881)
APW	$7.158\times {10}^{-17}$ (1.112 $\times {10}^{-19})$	$0.688$ (1.025 $\times {10}^{-2})$	$0.443$ (1.921 $\times {10}^{-2})$	$5.320\times {10}^{-2}$ (4.221 $\times {10}^{-4})$	$-212.870$	$433.750$	$441.440$	0.091 (0.804)	0.095 (0.824)	0.086 (0.782)
APD	$1.748\times {10}^{-3}$ (6.538 $\times 10-4)$	0.695 (0.242)	$9.049\times {10}^{-2}$ (4.294 $\times {10}^{-3})$	$9.683\times {10}^{-2}$ (3.609 $\times {10}^{-2})$	$-219.224$	$446.449$	$454.096$	0.080 (0.832)	0.221 (0.855)	0.088 (0.824)
AEXP	2.652 (0.335)	$1.162\times {10}^{-3}$ (4.057 $\times {10}^{-3})$	$8.751\times {10}^{-2}$ (3.645 $\times {10}^{-2})$	0.221 (0.107)	$-227.040$	$462.081$	$469.728$	0.089 (0.753)	0.095 (0.791)	0.081 (0.725)
CWD	$4.215\times {10}^{-2}$ (1.021 $\times {10}^{-5})$	$0.278$ (2.224 $\times {10}^{-2})$	$1.1183\times {10}^{-2}$ (3.995 $\times {10}^{-4})$	86.231 (2.414)	$-205.350$	$418.710$	$426.360$	0.070 (0.955)	0.077 (0.971)	0.067(0.940)

shows the result of the estimates for the parameters of the models under analysis. With the information obtained in Table 5, practitioners can derive essential conclusions on the behavior of the device under analysis.

For this case study, we note that the CPD more closely describes the behaviors of the data established in Table 4. This statement is based on the behaviors of the statistics shown in Table 5, where the CPD shows lower AIC and BIC. The K-S, AD, and CVM statistics show competitive p-values concerning the other distributions in the comparative study. Those values are reflected directly in the behavior graphs shown in Figure 6.

Figure 6a shows the behaviors of the PDFs for each distribution under analysis. As can be seen concerning the histogram of the data presented in Table 4, the CPD offers an excellent adjustment to the device’s behavior when it is below the designed operating levels. Figure 6b shows the reliability behavior of the device throughout its useful life. For this specific case, it is observed that the distributions exhibit different behaviors with respect to the empirical reliability of Kaplan–Meier. For the CPD, it is observed that it fits with a form comparable to that presented by the empirical reliability. Figure 6c represents the useful life of the device operating under normal operating conditions. In this graph, practitioners can see that the CPD describes the data behaviors presented in Table 4 as a bathtub curve. This statement can be validated by the shape that the failures take in a non-parametric way. We should note that the CPD is not the only one that demonstrates the data as a bathtub curve. However, the CPD is the only distribution that adjusts to the non-parametric fault line. Figure 6d shows how cumulative failures behave over the life of the device. The distributions under analysis show competitive results, but the CPD fits more points in the non-parametric curve. Finally, Figure 6e presents the empirical form that the MRL takes and how the MRL adjusts to the different distributions under analysis. In the same way, as in the previous reliability graphs, the CPD exhibits competitive results.

For this case study, the graphic and statistical evidence shows that the CPD can be considered a viable option for reliability analysis. As established in this case study, the CPD adjusted the non-parametric curve of the hazard function to a great extent, which suggested that there will not be any bias in the information obtained by the CPD.

9.3. Case Study 3. Reliability Analysis for Sylwia Data

This case study focused on analyzing the data provided by Sylwia [33], which shows the thirty failure times obtained during the reliability tests.

Table 6. Parameter values and statistics estimated for the lifetime data presented by Sylwia [33].

Lifetimes
0.0094	0.0500	0.4064	4.6307	5.1741	5.8808	6.3348	7.1645	7.2316	8.2604
9.2662	9.3812	9.5223	9.8783	9.9346	10.019	10.408	10.479	11.076	11.325
11.528	11.923	12.029	12.074	12.184	12.355	12.538	12.805	13.462	13.853

presents the results obtained by Sylwia. In the same way, as in the previous cases, Figure 7 shows the empirical behaviors of the data through the TTT plot of the data presented in Table 6.

Table 7. Parameter values and statistics estimated for the lifetime data presented by Sylwia [33].

	Parameters				Statistics
Model	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\lambda }$	$\mathit{\theta }$	Loglik	AIC	BIC	K-S	AD	CVM
CPD	$2.701\times {10}^{-5}$ (1.601 $\times {10}^{-4})$	$4.780\times {10}^{-2}$ (3.057 $\times {10}^{-2})$	0.949 (0.460)	0.311 (0.217)	$-68.801$	$145.602$	$151.206$	0.069 (0.989)	0.024 (0.938)	$0.160\left(0.942\right)$
AWD	$0.345$ (1.313 $\times {10}^{-2})$	$3.932\times {10}^{-2}$ (2.336 $\times {10}^{-2})$	0.519 (0.223)	0.717 (0.200)	$-75.335$	$158.670$	$164.274$	0.125 (0.725)	0.092 (0.425)	0.442 (0.436)
MWE	$1.947\times {10}^{-2}$ (1.401 $\times {10}^{-2})$	0.247 (0.164)	$0.317$ (6.035 $\times {10}^{-2})$	$---$	$-73.221$	$152.442$	$156.646$	0.109 (0.865)	0.040 (0.785)	0.301 (0.701)
PERKS4	−4.994 (0.776)	$0.357$ (7.838 $\times {10}^{-2})$	$5.787$ (1.015 $\times {10}^{-2})$	$-0.193$ (8.645 $\times {10}^{-2})$	$-205.610$	$419.112$	$424.824$	0.079 (0.901)	0.083 (0.912)	0.075 (0.881)
APW	$3.535\times {10}^{-4}$ (7.264 $\times {10}^{-4})$	0.311 (0.188)	0.753 (0.171)	$8.760\times {10}^{-2}$ (5.089 $\times {10}^{-2})$	$-71.715$	$151.430$	$157.034$	0.106 (0.872)	0.045 (0.655)	0.232 (0.725)
APD	$9.704\times {10}^{-4}$ (3.278 $\times {10}^{-4})$	$0.110$ (6.757 $\times {10}^{-2})$	$0.677$ (3.076 $\times {10}^{-2})$	6.494 (3.791)	$-70.994$	$149.989$	$155.592$	0.090 (0.893)	0.030 (0.711)	0.221 (0.745)
AEXP	1.901 (0.223)	$2.259\times {10}^{-9}$ (4.236 $\times {10}^{-10})$	1.572 (0.654)	$0.137$ (1.063 $\times {10}^{-3})$	$-70.609$	$149.218$	$154.882$	0.092 (0.914)	0.033 (0.731)	0.225 (0.756)
CWD	$4.861\times {10}^{-2}$ (1.020 $\times {10}^{-2})$	0.197 (0.100)	$8.890\times {10}^{-2}$ (3.142 $\times {10}^{-3})$	5.6345 (2.336)	$-70.561$	$149.122$	$154.726$	0.095 (0.925)	0.035 (0.755)	0.228 (0.794)

shows the estimates obtained for each parameter of the distributions under study. The results of this specific case study show that the distributions under analysis exhibit results that are close to each other. Despite this, the CPD has proven to be a good option since the statistics in Table 7 show a better fit and lower AIC and BIC relating to the other distributions. The above can be established by the reliability graphs shown in Figure 8.

Figure 8a shows the behaviors of the different PDFs for the distributions under analysis. In this case study, it is observed that some of the distributions under analysis offer remarkably close results. However, the CPD shows a slightly closer fit to the histogram, so practitioners could consider the CPD as an alternative analysis. In turn, Figure 8b shows the shape of the reliability curve relating to the empirical Kaplan–Meier curve. In this case, the reliability behavior of the equipment throughout its useful life is slightly better described with the CPD. Regarding the shape of the failure rate presented in Figure 8c, it can be observed that in all cases, the distributions under examination exhibit non-monotonic behavior close to the non-parametric curve established in the graph cited before. Nevertheless, based on the results obtained in Table 7, the CPD is our best option, as it is the closest to describing the failure times of the device under analysis. Regarding the graphs that show the cumulative failure and the MRL, which can be seen in Figure 8c,d, respectively, the results between the distributions are close for this case study. In both cases, the distribution behaviors are similar; thus, it is necessary to refer to the statistics obtained in Table 7, so that practitioners can choose the appropriate distribution based on the application and device under analysis.

In conclusion, for this case study, the following can be established. The graphic results show that the distributions under analysis offer comparable results, even though all distributions can model the behavior of failures in a non-monotonic way. Nevertheless, the statistical decision criteria obtained in Table 7 show that the CPD obtains a slight advantage regarding the competing distributions.

It should be noted that the efficiency of the distributions established in the comparative analysis could vary significantly based on the nature of the data. For example, if the data were increasing or decreasing monotones, in this manuscript, they were not considered within the study cases since the purpose was to assess one of the most deeply rooted concepts of reliability, the description of the failures in the devices under the bathtub shape.

10. Conclusions and Future Work

In this paper, a hybrid statistical distribution with applications in reliability engineering was presented. The proposed distribution was based on the additive methodology, which added the hazard function of the ChD and the PD, thus forming the CPD. The motivation for using PD is due to the similarity in the behavior of certain actuary systems with reliability engineering. The proposed distribution consists of four parameters, distributed in three shape parameters and one scale parameter. One of the most important characteristics of CPD is the ability to characterize failure times in a non-monotonic way. Non-monotonic failure times usually occur in devices whose physical constitution is based on semiconductors. To make the CPD attractive to reliability engineering practitioners, statistical properties such as measures of central tendency, order statistics, moments, MRL, and entropy were developed. In turn, the base functions of the CPD were demonstrated with data gathered from the ALT, with voltage and temperature, as well as the combination of temperature with a non-thermal variable, established as stress signals. The parameters of CPD were estimated via MLE; for this, an RStudio code was developed.

The CPD was tested in three case studies, focusing on devices whose lifetimes had the ability to be non-monotonic. In each case study, the CPD was compared with statistical distributions with reliability applications and properties such as the CPD. The distributions considered in the comparative analysis can also model non-monotonic failure behavior. The results obtained in the case studies showed that the CPD offers competitive results and that the practitioners can consider this distribution within the reliability analysis.

In future work, it may be beneficial to conduct a Bayesian analysis of the CPD and verify its behavior under the same analysis technique involving other distributions with equivalent properties. On the other hand, the PERKS5 distribution can be considered, which is derived from the same family as the PD considered in this manuscript. Finally, given the nature of PD used in actuarial applications, CPD can be conceptualized in other analyses not solely focused on reliability engineering.