A Three-Component Additive Weibull Distribution and Its Reliability Implications

Abstract

We consider a new failure time distribution by combining three Weibull distributions. The proposed distribution exhibits decreasing, increasing, and bathtub-shaped failure rate functions. Some statistical properties of the proposed distribution are presented. The maximum likelihood method is used to estimate the parameters, reliability, and failure rate functions of the proposed distribution, and the uncertainties of the estimates are obtained by using asymptotic confidence intervals. The cross-entropy method, a stochastic optimization procedure, is used to find the maximum of the log-likelihood function. The superiority of the proposed distribution is demonstrated on two benchmark datasets.

1. Introduction

The additive Weibull (AddW) distribution [1] is a distribution in reliability proposed for fitting bathtub-shaped failure rate data. The purpose of the AddW distribution is to provide a single model that can be used to model simultaneously all three phases of a bathtub curve. It combines two Weibull distributions [2], and its cumulative distribution function (CDF) is given by

$F\left(x\right)=1-exp\left[-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }\right],x\ge 0,$

where $\theta$ and $\lambda$ are non-negative and $0<\beta <1<\alpha$. The shape parameters are restricted by $\beta <1<\alpha$ so that its failure rate function is bathtub-shaped only. Later, Lemonte et al. [3] has provided a detailed study on the AddW and has modified the conditions for the parameters as $\theta ,\lambda >0,0<\beta <\alpha$ which allows the AddW to have an increasing failure rate if $1<\beta <\alpha$, a decreasing failure rate if $\beta <\alpha <1$, and a bathtub-shaped failure rate if $\beta <1<\alpha$. Although the AddW distribution has bathtub-shaped failure rate functions, it is not flexible enough for fitting datasets with complex bathtub-shaped failure rates.

Many proposed distributions which result from combining different distributions, rather than purely from two Weibull distributions, have outperformed the AddW distribution. The new modified Weibull (NMW) distribution [4] combines the Weibull distribution and the modified Weibull distribution [5]. It has five parameters and its CDF is given by:

$F\left(x\right)=1-e^{-\left(\alpha x^{\theta }+\beta x^{\gamma }{e}^{\lambda x}\right)},x\ge 0,$

where $\alpha ,\beta ,\gamma ,\theta >0$, and $\lambda \ge 0$. A Bayes study of the improved NMW by using Hamiltonian Monte Carlo simulation is given in [6]. The additive modified Weibull (AMW) distribution [7] combines the modified Weibull distribution and the Gompertz distribution [8]. This distribution has five parameters and its CDF is given by:

$F\left(x\right)=1-e^{-\left(\alpha x^{\theta }{e}^{\gamma x}+e^{\lambda x-\beta }-e^{-\beta }\right)},x\ge 0,$

where $\alpha ,\beta ,\theta >0$, and $\gamma \lambda \ge 0$. The very flexible Weibull distributions [9] results from a linear combination of two logarithms of cumulative hazard functions. Its failure rate function can be monotone, bathtub-shaped, modified bathtub-shaped, or even upside-down bathtub-shaped. The log-normal modified Weibull distribution [10] results from combining the log-normal distribution and the modified Weibull distribution, and this distribution also has five parameters. The additive Chen–Weibull distribution [11] is a four-parameter distribution that combines the Weibull distribution and the Chen distribution. Recently, a generalization of the AddW distribution, called the generalized additive Weibull (GAW) distribution [12], has been proposed. Its CDF is given by:

$F\left(x\right)={\left[1-exp\left(-\alpha x^{\beta }-\gamma x^{\lambda }\right)\right]}^{\theta },x\ge 0,$

where $\alpha ,\theta ,\gamma ,\beta ,\lambda >0$. This distribution comprises a set of distributions, such as modified generalized linear failure rate, exponentiated exponential, exponentiated Weibull, modified Weibull distributions, and additive Weibull, among others.

In addition to the combined distributions that provide a flexible bathtub-shaped failure rate, other distributions also have a bathtub-shaped failure rate function. Some of them can be mentioned as the modified Weibull distribution [5], exponentiated Weibull distribution [13], Hjorth distribution [14], generalized modified Weibull distribution [15], beta modified Weibull distribution [16], upper truncated Weibull distribution [17], new three-parameter exponential-type distribution [18], beta Generalized Weibull distribution [19], Alpha logarithmic transformed Weibull distribution [20], logistic exponential distribution [21], and five-parameter spline distribution [22]. However, these distributions cannot provide bathtub-shaped failure rate functions with a long flat region.

More recently, the logarithmic transformed Weibull distribution was introduced by using a logarithm transformation method [23]. It has constant, increasing, decreasing, unimodal, and unimodal then bathtub-shaped hazard rates. Shakhatreh et al. [24] introduced the generalized extended exponential-Weibull distribution for fitting non-monotonic failure rate data. This model includes the Weibull, generalized Weibull, exponentiated generalized linear exponential, and exponential-Weibull distributions. A new bounded distribution, called bounded weighted exponential distribution is introduced by Mallick et al. [25]. The proposed distribution exhibits increasing and bathtub-shaped failure rate functions.

In this study, we go one step further beyond the AddW distribution by introducing a distribution which combines three Weibull distributions. This distribution will be referred to as the three-component additive Weibull distribution (3CAW). The 3CAW distribution exhibits increasing, decreasing, and bathtub-shaped failure rate functions. It provides a flexible failure rate function, especially a bathtub-shaped failure rate with a long flat region, and it will be demonstrated to provide good fits to given well-known datasets. The 3CAW is useful for modeling either a series system with three independent components where each component follows a Weibull distribution or a component that is affected by three major failure modes, each following a Weibull distribution.

In fact, the combination of three Weibull components has already been proposed by Khalil et al. [26]. The model was named the flexible additive Weibull distribution and its CDF is given by:

$F\left(x\right)=1-e^{-a{x}^b-c{x}^d-g{x}^k},x>0,$

where $a,b,c,d,g,$ and $k>0$. However, it is clear that this model cannot be identifiable, meaning that the same model might be produced by different sets of parameters’ values. The identifiability of a distribution is a very important property of statistical models. It allows us to obtain precise estimates of the parameters’ values. Without identifiability, we cannot estimate the true parameters’ values, even with an infinite number of observations. Although the idea of combining three Weibull distributions was already proposed by Khalil et al. [26], our independent study will differ from [26] in many aspects. Firstly, our proposed model will have a different parameterized form. Secondly, and most importantly, our proposed model will be proved to be identifiable. Thirdly, our study will use a stochastic optimization method that is more likely to find the global maximum since the log-likelihood function of the 6-parameter model is usually highly multimodal. Lastly, in the Applications section, we will carry out the assessments of all the uncertainties of the parameter estimates.

The rest of the paper is organized as follows. Section 2 introduces the 3CAW distribution and its reliability characteristics. Section 3 studies some statistical properties of the 3CAW distribution. Section 4 discusses the methods of maximum likelihood and asymptotic confidence interval estimations for parameters, reliability, and failure rate functions. Section 5 offers the simulation study. Section 6 brings applications of the 3CAW to two well-known data sets. Finally, Section 7 concludes the paper.

2. The 3CAW Distribution

In our previous studies [6,11] we realized that if we combine a model which has a bathtub-shaped failure rate function with a model which has an increasing failure rate function, we might obtain a model which has a very flexible bathtub-shaped failure rate function with a long flat region. In fact, we encountered such distribution in the literature, such as [4,7,10,11]. In this study, the 3CAW is formed by combining three Weibull components, where the first two components form a model which possibly exhibits a bathtub-shaped failure rate function and the last component can provide an increasing failure rate function. The CDF of the 3CAW which combines three Weibull distributions is given by:

(1)$F\left(x\right)=1-exp\left[-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }\right],x>0,$

where $\theta ,\lambda ,\eta >0$, and $0<\alpha <\beta <\gamma$. The restriction of the shape parameters, i.e., $\alpha <\beta <\gamma$, leads to a very important property of the proposed model which will be stated in the following theorem.

Theorem 1

(Identifiability theorem). The 3CAW distribution with the CDF given by Equation (1) is identifiable.

Proof. 

The proof of Theorem 1 is given in Appendix A.    □

The corresponding probability density function (PDF) is given by:

$f\left(x\right)=\left[\alpha \theta {\left(\theta x\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x\right)}^{\beta -1}+\gamma \eta {\left(\eta x\right)}^{\gamma -1}\right]exp\left[-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }\right].$

It has the cumulative failure rate function given in the following form:

$H\left(x\right)={\left(\theta x\right)}^{\alpha }+{\left(\lambda x\right)}^{\beta }+{\left(\eta x\right)}^{\gamma }.$

Its reliability and failure rate functions are, respectively, given by:

$R\left(x\right)=exp\left[-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }\right],$

and

(2)$h\left(x\right)=\alpha \theta {\left(\theta x\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x\right)}^{\beta -1}+\gamma \eta {\left(\eta x\right)}^{\gamma -1}.$

From the CDF of the 3CAW distribution, it is clear that it includes many other distributions as its sub-models. Table 1 shows a list of distributions that can be derived from the 3CAW distribution.

The expression of the failure rate function in Equation (2) is the sum of three Weibull failure rate functions. It is decreasing if $\alpha <\beta <\gamma <1$, increasing if $1<\alpha <\beta <\gamma$, and bathtub-shaped if $\alpha <1<\beta <\gamma$ or $\alpha <\beta <1<\gamma$. Figure 1 shows the plots of the PDF and failure rate function of the 3CAW at selected parameter values. The 3CAW distribution has the PDF and failure rate function with flexible curves, and more importantly it has a long flat region of the bathtub-shaped failure rate function which is very important in reliability and biological studies.

3. Properties of the 3CAW Model

In this section, some statistical properties of the 3CAW distribution including the behavior of the failure rate function, the moments, order statistics, quantiles, and mean residual life are discussed.

3.1. The Failure Rate Function

We study the limiting behavior of $h\left(x\right)$ as x tend to 0 or $+\infty$.

• If $\gamma <1$, then ${lim}_{x\to 0}h\left(x\right)=+\infty$ and ${lim}_{x\to +\infty }h\left(x\right)=0$. In this case, the failure rate function is decreasing;

• If $\alpha >1$, then ${lim}_{x\to 0}h\left(x\right)=0$ and ${lim}_{x\to +\infty }h\left(x\right)=+\infty$. In this case, the failure rate function is increasing;

• If $\alpha <1<\beta$ or $\beta <1<\gamma$ then ${lim}_{x\to 0}h\left(x\right)=+\infty$ and ${lim}_{x\to +\infty }h\left(x\right)=+\infty$. In this case, the failure rate function is bathtub-shaped.

3.2. The Moments

In mathematical analysis we know that:

${\int }_0^{+\infty }{x}^r{e}^{-A{x}^s}dx<+\infty ,\mathrm{for}\mathrm{all}A,s>0\mathrm{and}r\in {\mathbb{Z}}^{+}.$

Therefore, from the expression of the CDF of the 3CAW, we see that the rth non-central moment ${\mu }_r$ (the rth moment about the origin) of the 3CAW always exists for all $r\in {\mathbb{Z}}^{+}$. It can be derived as follows:

(3)$\begin{array}{cc}\hfill {\mu }_r& ={\int }_0^{+\infty }{x}^{r}dF\left(x\right)\hfill \\ \hfill & =-{\int }_0^{+\infty }{x}^{r}d{e}^{-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }}\hfill \\ \hfill & =r{\int }_0^{+\infty }{x}^{r-1}{e}^{-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }}dx,\hfill \end{array}$

for $r=1,2,\dots$. Since (3) can not be obtained in closed form, some suitable numerical method, such as Gauss–Kronrod quadrature method, can be used to obtain this integral.

3.3. Order Statistics

Let a random sample $X_1,X_2,\dots ,X_n$ from the 3CAW distribution and $X_{k:n}$ is the kth order statistic of the sample, then the PDF of $X_{k:n}$ is given by:

$f_{k:n}\left(x\right)=\frac{1}{B(k,n-k+1)}F{\left(x\right)}^{k-1}{(1-F\left(x\right))}^{n-k}f\left(x\right),$

where $B(·,·)$ is the beta function. Since:

${(1-F\left(x\right))}^{n-k}=e^{-(n-k)H\left(x\right)},$

and

$F{\left(x\right)}^{k-1}={(1-e^{-H\left(x\right)})}^{k-1}=\sum _{l=0}^{k-1}\left(\genfrac{}{}{0pt}{}{k-1}{l}\right){(-1)}^l{e}^{-lH\left(x\right)}.$

The PDF of $X_{k:n}$ is derived as:

$\begin{array}{cc}\hfill f_{k:n}\left(x\right)& =\frac{1}{B(k,n-k+1)}\sum _{l=0}^{k-1}\left(\genfrac{}{}{0pt}{}{k-1}{l}\right){(-1)}^{l}h\left(x\right)e^{-(n+l+1-k)H\left(x\right)}\hfill \\ \hfill & =n\left(\genfrac{}{}{0pt}{}{n-1}{k-1}\right)\sum _{l=0}^{k-1}\left(\genfrac{}{}{0pt}{}{k-1}{l}\right)\frac{{(-1)}^l}{n+l+1-k}f(x;\alpha ,\beta ,\gamma ,{\theta }^{\prime },{\lambda }^{\prime },{\eta }^{\prime }),\hfill \end{array}$

where ${\theta }^{\prime }=\theta \sqrt[\alpha ]{n+l+1-k},{\lambda }^{\prime }=\lambda \sqrt[\beta ]{n+l+1-k}$ and ${\eta }^{\prime }=\eta \sqrt[\gamma ]{n+l+1-k}$.

3.4. Quantiles, Median, and Mode

The quantile $x_q$ is the solution of the following non-linear equation

(4)$\begin{array}{c}\hfill {\left(\theta x_q\right)}^{\alpha }+{\left(\lambda x_q\right)}^{\beta }+{\left(\eta x_q\right)}^{\gamma }+log(1-q)=0.\end{array}$

If $q=1/2$, the quantile $x_q$ is the median.

The mode is the value of x for which the PDF $f\left(x\right)$ attains its maximum. Therefore, the mode is the solution of the following non-linear equation:

(5)$\begin{array}{c}\hfill h^{\prime }\left(x\right)R\left(x\right)+h\left(x\right)R^{\prime }\left(x\right)=0\end{array}$

where $R\left(x\right)$ and $h\left(x\right)$ are the reliability failure rate functions, respectively, and $R^{\prime }\left(x\right)$ and $h^{\prime }\left(x\right)$ are their corresponding derivatives with respect to x. Solving (4) and (5) requires some suitable numerical methods.

Data $x_1,\dots ,x_n$ can be simulated from the 3CAW distribution by using a simple algorithm as follows:

• 1:. Generate $u_i\sim U(0,1)$, for $i=1,\dots ,n$;

• 2:. Solve ${\left(\theta x_i\right)}^{\alpha }+{\left(\lambda x_i\right)}^{\beta }+{\left(\eta x_i\right)}^{\gamma }=log\frac{1}{1-u_i}$.

3.5. Mean Residual Life

The mean residual life at time x of a random variable $X\sim$ 3CAW distribution is the expectation of the conditional random variable $X-x|X>x$. It is the mean time to failure of an object that has survived to time x, and it is defined as:

(6)$\begin{array}{cc}\hfill M\left(x\right)& =\mathbb{E}\left[X-x|X>x\right]=\frac{1}{R\left(x\right)}{\int }_x^{+\infty }R\left(u\right)du=\frac{1}{R\left(x\right)}{\int }_0^{+\infty }R(u+x)du\hfill \\ \hfill & =e^{{\left(\theta x\right)}^{\alpha }+{\left(\lambda x\right)}^{\beta }+{\left(\eta x\right)}^{\gamma }}{\int }_0^{+\infty }{e}^{-{\left(\theta (u+x)\right)}^{\alpha }-{\left(\lambda (u+x)\right)}^{\beta }-{\left(\eta (u+x)\right)}^{\gamma }}du.\hfill \end{array}$

In order to obtain the integral in (6), suitable numerical methods need to be used.

4. Parameter Estimation

4.1. Maximum Likelihood Estimation

Let $\mathcal{D}:x_1,\dots ,x_n$ be an observed sample from the 3CAW distribution with unknown parameter vector $\mathsf{\theta }={(\alpha ,\beta ,\gamma ,\theta ,\lambda ,\eta )}^T$. The likelihood function of the 3CAW is defined as:

$\begin{array}{cc}\hfill L\left(\mathsf{\theta }\right)& =\prod _{i=1}^{n}f(x_i;\mathsf{\theta })\hfill \\ \hfill & =\prod _{i=1}^n\left[\left(\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}\right)e^{-{\left(\theta x_i\right)}^{\alpha }-{\left(\lambda x_i\right)}^{\beta }-{\left(\eta x_i\right)}^{\gamma }}\right].\hfill \end{array}$

Then, the log-likelihood function is derived as:

(7)$l=\sum _{i=1}^{n}log\left[\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}\right]-\sum _{i=1}^n\left[{\left(\theta x_i\right)}^{\alpha }+{\left(\lambda x_i\right)}^{\beta }+{\left(\eta x_i\right)}^{\gamma }\right].$

We maximize this function to obtain the maximum likelihood estimates (MLEs) ${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}={(\widehat{\alpha },\widehat{\beta },\widehat{\gamma },\widehat{\theta },\widehat{\lambda },\widehat{\eta })}^T$ of the parameter vector $\mathsf{\theta }={(\alpha ,\beta ,\gamma ,\theta ,\lambda ,\eta )}^T$. This can be written in mathematical form as:

${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}=arg\underset{\theta \in \Theta }{max}l\left(\mathsf{\theta }\right),$

where $\Theta =\{{(\alpha ,\beta ,\gamma ,\theta ,\lambda ,\eta )}^T:\theta ,\lambda ,\eta >0\mathrm{and}0<\alpha <\beta <\gamma \}$.

One way to do so is to first calculate the first-order partial derivatives of the log-likelihood function with respect to the parameters:

$\begin{array}{cc}\hfill \frac{\partial l}{\partial \alpha }& =\sum _{i=1}^n\frac{\theta {\left(\theta x_i\right)}^{\alpha -1}\left(1+\alpha log\left(\theta x_i\right)\right)}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}-\sum _{i=1}^n{\left(\theta x_i\right)}^{\alpha }log\left(\theta x_i\right),\hfill \\ \hfill \frac{\partial l}{\partial \beta }& =\sum _{i=1}^n\frac{\lambda {\left(\lambda x_i\right)}^{\beta -1}\left(1+\beta log\left(\lambda x_i\right)\right)}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}-\sum _{i=1}^n{\left(\lambda x_i\right)}^{\beta }log\left(\lambda x_i\right),\hfill \\ \hfill \frac{\partial l}{\partial \gamma }& =\sum _{i=1}^n\frac{\eta {\left(\eta x_i\right)}^{\gamma -1}\left(1+\gamma log\left(\eta x_i\right)\right)}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}-\sum _{i=1}^n{\left(\eta x_i\right)}^{\gamma }log\left(\eta x_i\right),\hfill \\ \hfill \frac{\partial l}{\partial \theta }& =\sum _{i=1}^n\frac{-{\left(\alpha \theta \right)}^2{\left(\theta x_i\right)}^{\alpha -1}}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}+\alpha \theta \sum _{i=1}^n{\left(\theta x_i\right)}^{\alpha },\hfill \\ \hfill \frac{\partial l}{\partial \lambda }& =\sum _{i=1}^n\frac{-{\left(\beta \lambda \right)}^2{\left(\lambda x_i\right)}^{\beta -1}}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}+\beta \lambda \sum _{i=1}^n{\left(\lambda x_i\right)}^{\beta },\hfill \\ \hfill \frac{\partial l}{\partial \eta }& =\sum _{i=1}^n\frac{-{\left(\gamma \eta \right)}^2{\left(\eta x_i\right)}^{\gamma -1}}{\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1}}+\gamma \eta \sum _{i=1}^n{\left(\eta x_i\right)}^{\gamma },\hfill \end{array}$

and then set these partial derivatives equals zero and solve the system of non-linear equations by using some suitable numerical methods, such as Newton–Raphson or quasi-Newton methods.

Another better way to maximize the log-likelihood function in Equation (7) is to use the stochastic optimization procedure, called cross-entropy (CE) method [29], to maximize the log-likelihood function. The reason why CE works better than gradient-type methods is that the likelihood function is usually highly multimodal (that is, there are many local maxima and minima). For such multi-extremal optimization problems, a gradient-type method would typically find a local maximum rather than a global maximum. CE, on the other hand, is a global optimization technique that does not require local (i.e., gradient) information, and so is more likely to find the global maximum.

We prefer the latter due to its advantages in searching for a global maximum, and it is easy to implement by using the CE package in R or Matlab. All users need to do is to define the log-likelihood function, define the parameter space, specified algorithm’s initial values, and run the algorithm through trial and error until we find the best optimizer which has the largest log-likelihood. The CE algorithm and R code for continuous optimization are, respectively, given in Appendix E and Appendix F.

After obtaining the MLEs of the model parameters, using the invariant property of MLEs, we obtain the MLEs of $R\left(x\right)$ and $h\left(x\right)$ at a certain time x as:

$\begin{array}{cc}\hfill \widehat{R}\left(x\right)& =exp\left[-{\left(\widehat{\theta }x\right)}^{\widehat{\alpha }}-{\left(\widehat{\lambda }x\right)}^{\widehat{\beta }}-{\left(\widehat{\eta }x\right)}^{\widehat{\gamma }}\right],\hfill \\ \hfill \widehat{h}\left(x\right)& =\widehat{\alpha }\widehat{\theta }{\left(\widehat{\theta }x\right)}^{\widehat{\alpha }-1}+\widehat{\beta }\widehat{\lambda }{\left(\widehat{\lambda }x\right)}^{\widehat{\beta }-1}+\widehat{\gamma }\widehat{\eta }{\left(\widehat{\eta }x\right)}^{\widehat{\gamma }-1}.\hfill \end{array}$

4.2. Asymptotic Confidence Interval

To obtain the asymptotic confidence intervals (CIs) for the parameters, and the reliability and failure rate functions, we first calculate the symmetric matrix of second-order partial derivatives.

${\nabla }^{2}l=\left[\begin{array}{cccccc}\frac{{\partial }^{2}l}{\partial {\alpha }^2}& \frac{{\partial }^{2}l}{\partial \alpha \partial \beta }& \frac{{\partial }^{2}l}{\partial \alpha \partial \gamma }& \frac{{\partial }^{2}l}{\partial \alpha \partial \theta }& \frac{{\partial }^{2}l}{\partial \alpha \partial \lambda }& \frac{{\partial }^{2}l}{\partial \alpha \partial \eta }\\ ·& \frac{{\partial }^{2}l}{\partial {\beta }^2}& \frac{{\partial }^{2}l}{\partial \beta \partial \gamma }& \frac{{\partial }^{2}l}{\partial \beta \partial \theta }& \frac{{\partial }^{2}l}{\partial \beta \partial \lambda }& \frac{{\partial }^{2}l}{\partial \beta \partial \eta }\\ ·& ·& \frac{{\partial }^{2}l}{\partial {\gamma }^2}& \frac{{\partial }^{2}l}{\partial \gamma \partial \theta }& \frac{{\partial }^{2}l}{\partial \gamma \partial \lambda }& \frac{{\partial }^{2}l}{\partial \gamma \partial \eta }\\ ·& ·& ·& \frac{{\partial }^{2}l}{\partial {\theta }^2}& \frac{{\partial }^{2}l}{\partial \theta \partial \lambda }& \frac{{\partial }^{2}l}{\partial \theta \partial \eta }\\ ·& ·& ·& ·& \frac{{\partial }^{2}l}{\partial {\lambda }^2}& \frac{{\partial }^{2}l}{\partial \lambda \partial \eta }\\ ·& ·& ·& ·& ·& \frac{{\partial }^{2}l}{\partial {\eta }^2}\end{array}\right],$

where the expressions for second-order partial derivatives are provided in Appendix C. Then, the observed (Fisher) information matrix is defined as

$\mathbf{I}\left(\mathsf{\theta }\right)=-{\nabla }^{2}l.$

The matrix $\mathbf{I}$ can be evaluated at ${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}$, the MLE of $\mathsf{\theta }$, by using the plug-in method

$\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\mathrm{ML}}\right)={\left.-{\nabla }^{2}l\right|}_{\mathsf{\theta }={\widehat{\mathsf{\theta }}}_{\mathrm{ML}}}.$

The approximate variance matrix is obtained by

$\widehat{V}=\left[\begin{array}{cccccc}Var\left(\widehat{\alpha }\right)& Cov(\widehat{\alpha },\widehat{\beta })& Cov(\widehat{\alpha },\widehat{\gamma })& Cov(\widehat{\alpha },\widehat{\theta })& Cov(\widehat{\alpha },\widehat{\lambda })& Cov(\widehat{\alpha },\widehat{\eta })\\ ·& Var\left(\widehat{\beta }\right)& Cov(\widehat{\beta },\widehat{\gamma })& Cov(\widehat{\alpha },\widehat{\theta })& Cov(\widehat{\beta },\widehat{\lambda })& Cov(\beta ,\widehat{\eta })\\ ·& ·& Var\left(\widehat{\gamma }\right)& Cov(\widehat{\gamma },\widehat{\theta })& Cov(\widehat{\gamma },\widehat{\lambda })& Cov(\widehat{\gamma },\widehat{\eta })\\ ·& ·& ·& Var\left(\widehat{\theta }\right)& Cov(\widehat{\theta },\widehat{\lambda })& Cov(\widehat{\theta },\widehat{\eta })\\ ·& ·& ·& ·& Var\left(\widehat{\lambda }\right)& Cov(\widehat{\lambda },\widehat{\eta })\\ ·& ·& ·& ·& ·& Var\left(\widehat{\eta }\right)\end{array}\right]\approx {\mathbf{I}}^{-1}\left({\widehat{\mathsf{\theta }}}_{\mathrm{ML}}\right).$

Then, the CIs for the MLEs can be constructed by using the following approximation [30].

${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}\approx \mathcal{N}\left(\mathsf{\theta },\widehat{V}.\right)$

The $100(1-\xi )\%$ CIs for $\alpha ,\beta ,\gamma ,\theta ,\lambda$, and $\eta$ are, respectively, defined as:

$\begin{array}{cc}\hfill \widehat{\alpha }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\alpha }\right)},& \widehat{\beta }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\beta }\right)},\hfill \\ \hfill \widehat{\gamma }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\gamma }\right)},& \widehat{\theta }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\theta }\right)},\hfill \\ \hfill \widehat{\lambda }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\lambda }\right)},& \widehat{\eta }\pm z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\eta }\right)},\hfill \end{array}$

where $z_{1-\xi /2}$ is the $100(1-\xi /2)$ percentage point of $N(0,1)$.

The delta method which was discussed in [31] and used recently by [32] is used to approximate the variance of $\widehat{R}\left(x\right)$ and $\widehat{h}\left(x\right)$.

$\begin{array}{c}\hfill {\widehat{\sigma }}_{\widehat{R}\left(x\right)}^2={\left[\nabla \widehat{R}\left(x\right)\right]}^T\left[\widehat{V}\right]\left[\nabla \widehat{R}\left(x\right)\right],{\widehat{\sigma }}_{\widehat{h}\left(x\right)}^2={\left[\nabla \widehat{h}\left(x\right)\right]}^T\left[\widehat{V}\right]\left[\nabla \widehat{h}\left(x\right)\right],\end{array}$

where $\nabla \widehat{R}\left(x\right)$ and $\nabla \widehat{h}\left(x\right)$ are, respectively, the gradient of $R\left(x\right)$ and $h\left(x\right)$ with respect to $\alpha ,\beta ,\gamma ,\theta ,\lambda$, and $\eta$ evaluated at $\widehat{\alpha },\widehat{\beta },\widehat{\gamma },\widehat{\theta },\widehat{\lambda }$, and $\widehat{\eta }$. The gradient of $h\left(x\right)$ and $R\left(x\right)$ with respect to $\alpha ,\beta ,\gamma ,\theta ,\lambda$, and $\eta$ are, respectively, given in Appendix B and Appendix D. Then, the $100(1-\xi )\%$ CIs for $R\left(x\right)$ and $h\left(x\right)$ at time x are, respectively, defined as:

$\begin{array}{c}\hfill \widehat{R}\left(x\right)\pm z_{1-\frac{\xi }{2}}\sqrt{{\widehat{\sigma }}_{\widehat{R}\left(x\right)}^2},\widehat{h}\left(x\right)\pm z_{1-\frac{\xi }{2}}\sqrt{{\widehat{\sigma }}_{\widehat{h}\left(x\right)}^2}.\end{array}$

Although the parameters are positive, the CIs may result in a negative lower bound. In this regard, we can replace the negative values by zero. Another way is to use the normal approximation of log-transformed MLE [33]. Then, the $100(1-\xi )\%$ CIs for $\alpha ,\beta ,\gamma ,\theta ,\lambda$, and $\eta$, are, respectively, defined as:

$\begin{array}{cc}\hfill \widehat{\alpha }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\alpha }\right)}}{\widehat{\alpha }}\right],& \widehat{\beta }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\beta }\right)}}{\widehat{\beta }}\right],\hfill \\ \hfill \widehat{\gamma }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\gamma }\right)}}{\widehat{\gamma }}\right],& \widehat{\theta }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\theta }\right)}}{\widehat{\theta }}\right],\hfill \\ \hfill \widehat{\lambda }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\lambda }\right)}}{\widehat{\lambda }}\right],& \widehat{\eta }exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{Var\left(\widehat{\eta }\right)}}{\widehat{\eta }}\right].\hfill \end{array}$

Likewise, the $100(1-\xi )\%$ CIs for $\widehat{R}\left(x\right)$ and $\widehat{h}\left(x\right)$ are, respectively, defined as:

$\begin{array}{c}\hfill \widehat{R}\left(x\right)exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{{\widehat{\sigma }}_{\widehat{R}\left(x\right)}^2}}{\widehat{R}\left(x\right)}\right],\widehat{h}\left(x\right)exp\left[\pm \frac{z_{1-\frac{\xi }{2}}\sqrt{{\widehat{\sigma }}_{\widehat{h}\left(x\right)}^2}}{\widehat{h}\left(x\right)}\right].\end{array}$

5. Simulation Results

We conduct a Monte Carlo simulation study to assess the performances of the proposed estimation methods for the proposed model. We first select five different sets of parameters’ values. For simplicity, we set all scale parameters to unity, i.e., $\theta =\lambda =\eta =1$, in all cases, and vary the values of the shape parameters as given in Table 2. For each set of parameters’ values, we simulate 1000 datasets with sample sizes $n=25,50,100$, and 200, respectively, from the 3CAW distribution by using the algorithm given in Section 3.4. Then, we find the MLEs of the parameters for fitting the 3CAW to each dataset, and calculate the biases and standard deviations (SDs) for the parameter estimates.

The results of the simulation study are listed in Table 2. The outcomes show that the MLEs of the parameters (except $\beta$ and $\gamma$) have low bias and SD in virtually all the cases, demonstrating that the estimated values are close to the parameters’ true values. In the results, we see that $\widehat{\gamma }$ has the largest biases and SDs compared with the other estimators. The biases and SDs of $\widehat{\beta }$ and $\widehat{\gamma }$ increase as we increase the values of $\beta$ and $\gamma$. This can be understandable because, with relatively small sample sizes, it is hard to estimate exactly the value of the Weibull shape parameter when it becomes larger and larger. In that case, the shape of the Weibull distribution does not change much when slightly varying the value of the shape parameter. The biases of $\widehat{\lambda }$ are always negative in all cases which means that the MLE underestimates the parameter $\lambda$. In almost all instances, the biases and SDs decrease as the sample size increases, indicating the consistency of the proposed estimators.

6. Applications

We use two real-life datasets known as the Aarset data [34] and Meeker and Escobar data [33] to investigate the performance of the 3CAW distribution. We also compare this distribution to some distributions, such as the AddW, NMW, AMW, and GAW. The log-likelihood ($\log \left(L\right)$) value, Kolmogorov–Smirnov (K-S) test statistic, Akaike information criterion (AIC), bias-corrected Akaike information criterion (AICc) (the adjusted version of AIC for small dataset), and Bayesian information criterion (BIC) are used for model selection. Among a set of alternative models, the model with the largest $\log \left(L\right)$ value and smallest K-S, AIC, BIC, and AICc values can be considered to be the best approximation model. Notice that K-S and AIC (as well as the AICc and BIC) values do not necessarily follow the same trend (i.e., the model with the smallest K-S may not have the smallest AIC). We calculate the CIs for the MLEs of the parameters, reliability, and failure rate functions of the 3CAW by using the log-transformed MLE method mentioned above.

6.1. Aarset Data

Data in Table 3 represent the times to failure of 50 devices (in weeks) [34]. The empirical scaled TTT-transform plot [35] indicated that this dataset has a bathtub-shaped failure rate because it is first convex and then concave, as shown in Figure 2b, which matches the blue-dashed line in Figure 2a.

The MLEs of the parameters of the 3CAW, GAW, AMW, NMW, and AddW, for fitting to Aarset data, are given in Table 4, and Table 5 provides a numeric comparison between 3CAW, GAW, AMW, NMW, and AddW in terms of $log\left(L\right)$ value, Kolmogorov–Smirnov (K-S) statistic, AIC, BIC, and AICc.

From this table, we can see that the 3CAW has largest $log\left(L\right)$ value and smallest K-S and AIC values. The AMW model performs similarly to the 3CAW and even has smaller BIC and AICc value. This is because it has one less parameter than the 3CAW. However, its model form has a more complicated structure than the 3CAW.

We also provide a visual comparison by plotting the estimated reliability and failure rate functions of the 3CAW, GAW, AMW, NMW, and AddW along with the empirical estimates of reliability and failure rate functions in Figure 3. Once again this figure supports our conclusion that the 3CAW distribution provides the best fit for this dataset because the estimated reliability and failure rate functions of the 3CAW are closer to the corresponding empirical curves. In Figure 3b, the step function is a non-parametric estimation of the failure rate function performed by dividing the time domain into bins of equal width, and then estimating the failure rate in each bin as the number of events in that bin divided by the number of items at risk. This estimate is completed by using the “pehaz” function in the “muhaz” R package [36].

Table 6 provides the MLEs of the parameters and their corresponding CIs for fitting 3CAW to Aarset data. The estimates of the shape parameters show that the first Weibull component has a decreasing failure rate which models the decreasing part of the bathtub curve; the second component has an increasing failure rate with slight change which models the middle part of the bathtub curve; the third component has an increasing failure rate with sharped change which models the sharped wear out of the bathtub curve. The 90% and 95% confidence intervals also support these claims. Table 7 provides the MLEs of the reliability and failure functions and their corresponding $95\%$ confidence intervals at each given failure time. Figure 4 displays the results in Table 7 visually.

6.2. Meeker and Escobar Data

Data in Table 8 represent failure and running times of 30 devices (in thousands of cycles) ([33], p. 383). Two failure modes were observed for this dataset. These data are shown also to have a bathtub-shaped failure rate as indicated by the scaled TTT-transform plot given in Figure 5.

We provide the MLEs of the parameters of the given distributions in Table 9, and the numeric comparison between the given distributions, for fitting to Meeker and Escobar data, in Table 10. Figure 6 displays the visual comparison of the mention distributions by plotting the estimated reliability and failure rate functions. From these results, we see that the 3CAW provides a good fit to Meeker and Escobar data. Again, we can see that AMW performs similar to 3CAW. Even though it has a smaller $log\left(L\right)$ value and a larger K-S value, its other values are still slightly less than that of 3CAW.

Table 11 provides the MLEs of the parameters and their corresponding CIs for fitting 3CAW to Meeker and Escobar data. We also calculate the MLEs of the reliability and failure functions and their corresponding $95\%$ confidence intervals at each given failure time in Table 12. The results in Table 12 are displayed visually in Figure 7.

7. Conclusions

In this study, we have introduced the 3CAW distribution by combining three Weibull distributions. The properties of the model have been studied. The MLE method has been exploited for estimating the parameters and the reliability characteristics of the 3CAW, and the uncertainties of the estimates have been obtained by using asymptotic confidence intervals. The cross-entropy method has been deployed to optimize the log-likelihood function. Two sets of failure data have been used for testing the superiority of the 3CAW. It is interesting to see a Bayesian analysis of this six-parameter distribution by using Hamiltonian Monte Carlo simulation method, and to see more applications of the 3CAW to real datasets arising from practice.

Acknowledgments

The author kindly thanks the anonymous referees and the associate editor for providing constructive comments and suggestions that greatly improved the paper. The author also kindly thanks Dirk P. Kroese, School of Mathematics and Physics, The University of Queensland, for his comments on why CE works better than gradient-type methods. The author also would like to thank Nabendu Pal, Department of Mathematics, University of Louisiana at Lafayette, for his advice on motivation for research, his comments on the paper, and his efforts on teaching and helping researchers in Vietnam from 2013 up till now. Finally, the author thanks TDTU colleagues Cao Xuan Phuong, Vo Xuan Thanh, and Le Ba Khiet for their comments and discussions during the process of writing and revising the paper.

Conflicts of Interest

The author declares no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Plots of (a) PDF and (b) corresponding failure rate functions of the 3CAW.

Figure 2. (a) Types of failure rate identify by shapes of the TTT-plot and (b) the empirical TTT-plot (blue) for Aarset data, where [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.], and its empirical version is [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.] denote the i-th ordered observation.

Figure 3. The estimated (a) reliability and (b) failure rate functions of 3CAW, GAW, AMW, NMW, and AddW for fitting to Aarset data.

Figure 4. The 95% confidence interval for the estimated (a) reliability and (b) failure rate functions (red lines) of 3CAW for fitting to Aarset data. The blue line is the nonparametric estimate of the reliability function.

Figure 5. (a) Types of failure rate identify by shapes of the TTT-plot and (b) the empirical TTT-plot (blue) for Meeker and Escobar data, where u and [Forumla omitted. See PDF.] are explained in Figure 2.

Figure 6. The estimated (a) reliability and (b) failure rate functions of 3CAW, GAW, AMW, NMW, and AddW for fitting to Meeker and Escobar data.

Figure 7. The 95% confidence interval for the estimated (a) reliability and (b) failure rate functions (red lines) of the 3CAW for fitting to Meeker and Escobar data. The blue line is the nonparametric estimate of the reliability function.

Table 1

The sub-models of the 3CAW.

Distribution	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$	$\mathit{\theta }$	$\mathit{\lambda }$	$\mathit{\eta }$	F(x)	Reference
Flexible additive Weibull	b	d	h	$\sqrt[b]{a}$	$\sqrt[d]{c}$	$\sqrt[k]{g}$	$1-e^{-a{x}^b-c{x}^d-g{x}^k}$	Khalil et al. [26]
Additive Weibull	-	-	0	-	-	0	$1-e^{-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }}$	Xie and Lai [1]
Non-linear failure rate	1	-	0	-	-	0	$1-e^{-\theta x-{\left(\lambda x\right)}^{\beta }}$	Salem [27]
Linear failure rate	1	2	0	-	-	0	$1-e^{-\theta x-{\left(\lambda x\right)}^2}$	Bain [28]
Weibull	-	0	0	-	0	0	$1-e^{-{\left(\theta x\right)}^{\alpha }}$	Weibull [2]
Rayleigh	2	0	0	-	0	0	$1-e^{-{\left(\theta x\right)}^2}$	Bain [28]
Exponential	1	0	0	-	0	0	$1-e^{-\theta x}$	Bain [28]

Table 2

Biases and standard deviations of the 3CAW parameter estimates for some selected parameters’ values at different sample sizes.

n	Bias						Standard Deviation
	$\widehat{\alpha }$	$\widehat{\beta }$	$\widehat{\gamma }$	$\widehat{\theta }$	$\widehat{\lambda }$	$\widehat{\eta }$	$\widehat{\alpha }$	$\widehat{\beta }$	$\widehat{\gamma }$	$\widehat{\theta }$	$\widehat{\lambda }$	$\widehat{\eta }$
	$\alpha =0.1,\beta =2,\gamma =10,\theta =1,\lambda =1,\eta =1$
25	0.0042	0.0877	1.3684	0.3379	−0.0195	0.0236	0.0288	1.1173	1.7485	1.3094	0.4039	0.1128
50	−0.0004	0.0645	1.0510	0.2783	−0.0094	0.0058	0.0185	0.9055	1.7630	1.2177	0.2560	0.0588
100	0.0001	0.0832	0.7688	0.3036	−0.0038	−0.0006	0.0126	0.7241	1.5956	1.0784	0.1810	0.0495
200	0.0000	0.0328	0.5359	0.1946	−0.0085	−0.0009	0.0087	0.5193	1.4143	0.9082	0.1399	0.0320
	$\alpha =0.2,\beta =4,\gamma =20,\theta =1,\lambda =1,\eta =1$
25	0.0065	0.1413	1.3990	0.2191	−0.0448	0.0138	0.0548	1.6876	2.1558	1.0436	0.2270	0.0420
50	0.0046	0.0337	0.8171	0.2309	−0.0276	0.0042	0.0376	1.4517	2.2321	0.9281	0.1360	0.0291
100	0.0029	0.0137	0.5407	0.1184	−0.0119	0.0015	0.0251	1.1596	1.9485	0.6962	0.0888	0.0204
200	0.0004	0.0100	0.4517	0.0693	−0.0034	0.0001	0.0181	0.9421	1.7966	0.5369	0.0585	0.0141
	$\alpha =0.4,\beta =6,\gamma =30,\theta =1,\lambda =1,\eta =1$
25	0.0263	0.1721	1.3080	0.1419	−0.0714	0.0075	0.1143	2.0221	2.3351	0.7280	0.2273	0.0274
50	0.0122	0.0848	0.8988	0.0682	−0.0208	0.0030	0.0770	1.6977	2.3738	0.5192	0.1039	0.0174
100	0.0048	0.0987	0.6592	0.0582	−0.0096	0.0002	0.0532	1.4810	2.2169	0.3873	0.0608	0.0126
200	−0.0001	0.0268	0.4164	−0.0049	−0.0028	0.0003	0.0352	1.2605	2.0607	0.2644	0.0395	0.0087
	$\alpha =0.6,\beta =8,\gamma =40,\theta =1,\lambda =1,\eta =1$
25	0.0332	0.2336	1.4777	0.0530	−0.0541	0.0054	0.1727	2.1670	2.4230	0.5059	0.1903	0.0191
50	0.0130	0.1410	0.8352	0.0288	−0.0171	0.0021	0.1190	1.8803	2.5656	0.3610	0.0906	0.0129
100	0.0084	0.1571	0.6749	0.0145	−0.0060	0.0012	0.0817	1.6897	2.4045	0.2641	0.0449	0.0091
200	0.0019	0.0820	0.4934	0.0069	−0.0022	0.0001	0.0546	1.4847	2.1561	0.1843	0.0295	0.0067
	$\alpha =0.8,\beta =10,\gamma =50,\theta =1,\lambda =1,\eta =1$
25	0.0536	0.1830	1.2681	0.0337	−0.0608	0.0049	0.2244	2.2711	2.5538	0.3692	0.1741	0.0167
50	0.0213	0.1323	1.0219	0.0051	−0.0144	0.0017	0.1527	2.0815	2.6315	0.2800	0.0776	0.0107
100	0.0030	0.0653	0.7603	0.0074	−0.0071	0.0008	0.1022	1.8486	2.5559	0.2047	0.0370	0.0077
200	0.0084	0.0476	0.4831	0.0038	−0.0021	0.0004	0.0766	1.6348	2.3628	0.1417	0.0247	0.0053

Table 3

Times to failure of $n=50$ devices put on life test (in weeks) [34].

0.1	0.2	1.0	1.0	1.0	1.0	1.0	2.0	3.0	6.0
7.0	11.0	12.0	18.0	18.0	18.0	18.0	18.0	21.0	32.0
36.0	40.0	45.0	46.0	47.0	50.0	55.0	60.0	63.0	63.0
67.0	67.0	67.0	67.0	72.0	75.0	79.0	82.0	82.0	83.0
84.0	84.0	84.0	85.0	85.0	85.0	85.0	85.0	86.0	86.0

Table 4

The MLEs of parameters for Aarset data.

Model	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\eta }}$
3CAW	0.49874	4.06058	97.88702	0.00748	0.01088	0.01175
GAW	$5.53451\times {10}^{-15}$	7.53436	4.89809	0.05455	324.645	-
AMW	0.0763	90.1357	0.0104	0.4579	1.0604	-
NMW	0.0710	$7.0150\times {10}^{-8}$	0.0160	0.5950	0.1970	-
AddW	0.7025	82.3350	-	0.0162	0.0118	-

Table 5

Log-likelihood, K-S statistic, AIC, BIC, and AICc for Aarset data.

Model	$log\left(\mathit{L}\right)$	K-S (p-Value)	AIC	BIC	AICc
3CAW	−202.53	0.057  (0.997)	417.06	428.53	419.02
GAW	−215.93	0.121  (0.457)	441.86	451.42	443.22
AMW	−203.57	0.068  (0.976)	417.14	426.70	418.50
NMW	−212.90	0.088  (0.838)	435.80	445.36	437.16
AddW	−206.10	0.111  (0.568)	420.19	427.84	421.08

Table 6

MLEs and confidence intervals for the parameters of the 3CAW for fitting to Aarset data.

Parameter	MLE	90% CI	95% CI
$\alpha$ (shape)	0.49874	[0.34723, 0.71637]	[0.32396, 0.76783]
$\beta$ (shape)	4.06058	[1.82322, 9.04351]	[1.56394, 10.5428]
$\gamma$ (shape)	97.8870	[59.1332, 162.039]	[53.6904, 178.465]
$\theta$ (scale)	0.00748	[0.00271, 0.02064]	[0.00223, 0.02507]
$\lambda$ (scale)	0.01088	[0.00911, 0.01300]	[0.00880, 0.01345]
$\eta$ (scale)	0.01175	[0.01168, 0.01182]	[0.01166, 0.01184]

Table 7

MLEs and $95\%$ confidence intervals for the estimated reliability and failure rate functions of the 3CAW for fitting to Aarset data.

x	$\widehat{\mathit{R}}\left(\mathit{x}\right)$	95% CI for $\widehat{\mathit{R}}\left(\mathit{x}\right)$	$\widehat{\mathit{h}}\left(\mathit{x}\right)$	95% CI for $\widehat{\mathit{h}}\left(\mathit{x}\right)$
0.1	0.97278	[0.93988, 1.00683]	0.13765	[0.05822, 0.32544]
0.2	0.96175	[0.92103, 1.00428]	0.09725	[0.04646, 0.20356]
1	0.91666	[0.85448, 0.98335]	0.04340	[0.02593, 0.07264]
2	0.88430	[0.81239, 0.96257]	0.03066	[0.01922, 0.04893]
3	0.86026	[0.78273, 0.94547]	0.02502	[0.01580, 0.03963]
6	0.80840	[0.72113, 0.90624]	0.01769	[0.01092, 0.02866]
7	0.79477	[0.70518, 0.89574]	0.01638	[0.01001, 0.02680]
11	0.74981	[0.65271, 0.86136]	0.01311	[0.00778, 0.02211]
12	0.74025	[0.64155, 0.85414]	0.01258	[0.00743, 0.02130]
18	0.69129	[0.58457, 0.81749]	0.01049	[0.00618, 0.01782]
21	0.67049	[0.56071, 0.80176]	0.00992	[0.00589, 0.01669]
32	0.60416	[0.48815, 0.74775]	0.00939	[0.00549, 0.01606]
36	0.58157	[0.46482, 0.72765]	0.00971	[0.00549, 0.01716]
40	0.55882	[0.44180, 0.70683]	0.01029	[0.00566, 0.01872]
45	0.52941	[0.41247, 0.67949]	0.01140	[0.00621, 0.02096]
46	0.52333	[0.40646, 0.67381]	0.01168	[0.00637, 0.02141]
47	0.51718	[0.40039, 0.66805]	0.01197	[0.00655, 0.02187]
50	0.49822	[0.38178, 0.65017]	0.01296	[0.00722, 0.02327]
55	0.46467	[0.34944, 0.61791]	0.01500	[0.00871, 0.02583]
60	0.42845	[0.31564, 0.58159]	0.01755	[0.01047, 0.02942]
63	0.40541	[0.29477, 0.55757]	0.01935	[0.01146, 0.03264]
67	0.37323	[0.26626, 0.52319]	0.02206	[0.01250, 0.03895]
72	0.33105	[0.22912, 0.47834]	0.02601	[0.01319, 0.05131]
75	0.30499	[0.20566, 0.45230]	0.02870	[0.01333, 0.06177]
79	0.26963	[0.17264, 0.42111]	0.03350	[0.01437, 0.07810]
82	0.23714	[0.14461, 0.38889]	0.06727	[0.02709, 0.16702]
83	0.21538	[0.12866, 0.36056]	0.13852	[0.05370, 0.35734]
84	0.17123	[0.09605, 0.30525]	0.36197	[0.17828, 0.73493]
85	0.08976	[0.04106, 0.19623]	1.05836	[0.56081, 1.99736]
86	0.01296	[0.00176, 0.09529]	3.20497	[1.11010, 9.25304]

Table 8

Failure and running times of a sample of $n=30$ devices (in thousands of cycles) [33].

2	10	13	23	23	28	30	65	80	88
106	143	147	173	181	212	245	247	261	266
275	293	300	300	300	300	300	300	300	300

Table 9

The MLEs of parameters for Meeker and Escobar data.

Models	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	$\widehat{\mathit{\theta }}$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\eta }}$
3CAW	0.55141	4.07734	133.13084	0.00178	0.00287	0.00333
GAW	$4.72286\times {10}^{-8}$	3.06873	3.48132	0.05456	120.502	-
AMW	0.0142	116.9665	0.0019	0.6788	0.3902	-
NMW	0.024	$5.991\times {10}^{-8}$	0.012	0.629	0.056	-
AddW	0.84960	41.25756	-	0.00368	0.00335	-

Table 10

Log-likelihood, K-S statistic, AIC, BIC, and AICc for Meeker and Escobar data.

Model	$log\left(\mathit{L}\right)$	K-S (p-Value)	AIC	BIC	AICc
3CAW	−154.67	0.147  (0.533)	321.33	329.74	324.98
GAW	−175.93	0.146  (0.545)	361.86	368.86	364.36
AMW	−155.58	0.167  (0.374)	321.16	328.17	323.66
NMW	−166.18	0.149  (0.522)	342.36	349.37	344.86
AddW	−162.58	0.168  (0.365)	333.17	340.82	334.05

Table 11

MLEs and confidence intervals for the parameters of the 3CAW for fitting to Meeker and Escobar data.

Parameter	MLE	90% CI	95% CI
$\alpha$ (shape)	0.55141	[0.31596, 0.96231]	[0.28399, 1.07064]
$\beta$ (shape)	4.07734	[0.81092, 20.5010]	[0.59513, 27.9348]
$\gamma$ (shape)	133.130	[71.4311, 248.122]	[63.4005, 279.553]
$\theta$ (scale)	0.00178	[0.00036, 0.00888]	[0.00026, 0.01208]
$\lambda$ (scale)	0.00287	[0.00221, 0.00372]	[0.00211, 0.00391]
$\eta$ (scale)	0.00333	[0.00331, 0.00335]	[0.00331, 0.00335]

Table 12

MLEs and $95\%$ confidence intervals for the estimated reliability and failure rate functions of the 3CAW for fitting to Meeker and Escobar data.

x	$\widehat{\mathit{R}}\left(\mathit{x}\right)$	95% CI for $\widehat{\mathit{R}}\left(\mathit{x}\right)$	$\widehat{\mathit{h}}\left(\mathit{x}\right)$	95% CI for $\widehat{\mathit{h}}\left(\mathit{x}\right)$
2	0.95633	[0.89780, 1.01868]	0.01231	[0.00507, 0.02989]
10	0.89722	[0.81043, 0.99329]	0.00598	[0.00306, 0.01167]
13	0.88219	[0.79072, 0.98425]	0.00532	[0.00271, 0.01045]
23	0.84223	[0.73971, 0.95897]	0.00412	[0.00199, 0.00854]
28	0.82581	[0.71885, 0.94870]	0.00377	[0.00177, 0.00803]
30	0.81970	[0.71105, 0.94495]	0.00366	[0.00170, 0.00786]
65	0.73675	[0.60423, 0.89832]	0.00265	[0.00122, 0.00573]
80	0.70903	[0.56979, 0.88229]	0.00248	[0.00122, 0.00506]
88	0.69527	[0.55344, 0.87345]	0.00242	[0.00122, 0.00481]
106	0.66598	[0.52075, 0.85171]	0.00237	[0.00120, 0.00470]
143	0.60852	[0.46353, 0.79885]	0.00257	[0.00102, 0.00645]
147	0.60224	[0.45748, 0.79281]	0.00261	[0.00102, 0.00672]
173	0.55996	[0.41565, 0.75437]	0.00302	[0.00110, 0.00829]
181	0.54623	[0.40157, 0.74299]	0.00319	[0.00118, 0.00860]
212	0.48863	[0.34233, 0.69744]	0.00406	[0.00181, 0.00910]
245	0.41867	[0.27729, 0.63213]	0.00538	[0.00259, 0.01120]
247	0.41414	[0.27341, 0.62733]	0.00548	[0.00260, 0.01153]
261	0.38169	[0.24618, 0.59180]	0.00619	[0.00257, 0.01494]
266	0.36980	[0.23628, 0.57877]	0.00647	[0.00251, 0.01672]
275	0.34805	[0.21784, 0.55609]	0.00700	[0.00235, 0.02085]
293	0.29242	[0.16971, 0.50387]	0.02534	[0.00796, 0.08069]
300	0.11930	[0.04987, 0.28540]	0.39711	[0.12714, 1.24030]

Appendix A. Proof of Theorem 1

Let ${\mathsf{\theta }}_1=({\alpha }_1,{\beta }_1,{\gamma }_1,{\theta }_1,{\lambda }_1,{\eta }_1)$ where $0<{\alpha }_1<{\beta }_1<{\gamma }_1$ and ${\mathsf{\theta }}_1,{\lambda }_1,{\eta }_1>0$, and ${\theta }_2=({\alpha }_2,{\beta }_2,{\gamma }_2,{\theta }_2,{\lambda }_2,{\eta }_2)$ where $0<{\alpha }_2<{\beta }_2<{\gamma }_2$ and ${\theta }_2,{\lambda }_2,{\eta }_2>0$. To prove the model is identifiable, we show that: If $F_{{\mathsf{\theta }}_1}\left(x\right)=F_{{\mathsf{\theta }}_2}\left(x\right)$ for all $x\in {\mathbb{R}}^{+}$, then ${\mathsf{\theta }}_1={\mathsf{\theta }}_2$. It means the mapping $\mathsf{\theta }\mapsto F_{\mathsf{\theta }}$ is one-to-one.

Indeed, for every $x\in {\mathbb{R}}^{+}$, we have $F_{{\mathsf{\theta }}_1}\left(x\right)=F_{{\mathsf{\theta }}_2}\left(x\right)$, so: (A1)$$e^{-{\left({\theta }_{1}x\right)}^{{\alpha }_1}-{\left({\lambda }_{1}x\right)}^{{\beta }_1}-{\left({\eta }_{1}x\right)}^{{\gamma }_1}}=e^{-{\left({\theta }_{2}x\right)}^{{\alpha }_2}-{\left({\lambda }_{2}x\right)}^{{\beta }_2}-{\left({\eta }_{2}x\right)}^{{\gamma }_2}}.$$

Taking logarithm both side of (A1), we obtain: (A2)$${\left({\theta }_{1}x\right)}^{{\alpha }_1}+{\left({\lambda }_{1}x\right)}^{{\beta }_1}+{\left({\eta }_{1}x\right)}^{{\gamma }_1}={\left({\theta }_{2}x\right)}^{{\alpha }_2}+{\left({\lambda }_{2}x\right)}^{{\beta }_2}+{\left({\eta }_{2}x\right)}^{{\gamma }_2}.$$

If ${\alpha }_1\ne {\alpha }_2$, i.e., ${\alpha }_1<{\alpha }_2$ or ${\alpha }_1>{\alpha }_2$, without lost of generality we assume that ${\alpha }_1<{\alpha }_2$, then dividing both side of (A2) by $x^{{\alpha }_1}$ gives: $${\theta }_1^{{\alpha }_1}+{\lambda }_1^{{\beta }_1}{x}^{{\beta }_1-{\alpha }_1}+{\eta }_1^{{\gamma }_1}{x}^{{\gamma }_1-{\alpha }_1}={\theta }_2^{{\alpha }_2}{x}^{{\alpha }_2-{\alpha }_1}+{\lambda }_2^{{\beta }_2}{x}^{{\beta }_2-{\alpha }_1}+{\eta }_2^{{\gamma }_2}{x}^{{\gamma }_2-{\alpha }_1}.$$

Letting $x\to 0$, we have ${\theta }_1^{{\alpha }_1}=0$, which is a contradiction since ${\theta }_1,{\alpha }_1>0$. Therefore, ${\alpha }_1={\alpha }_2$. Now we have: $${\theta }_1^{{\alpha }_1}+{\lambda }_1^{{\beta }_1}{x}^{{\beta }_1-{\alpha }_1}+{\eta }_1^{{\gamma }_1}{x}^{{\gamma }_1-{\alpha }_1}={\theta }_2^{{\alpha }_1}+{\lambda }_2^{{\beta }_2}{x}^{{\beta }_2-{\alpha }_1}+{\eta }_2^{{\gamma }_2}{x}^{{\gamma }_2-{\alpha }_1}.$$

Letting $x\to 0$, we obtain ${\theta }_1^{{\alpha }_1}={\theta }_2^{{\alpha }_1}$, and hence ${\theta }_1={\theta }_2$. Then (A2) reduces to: $${\left({\lambda }_{1}x\right)}^{{\beta }_1}+{\left({\eta }_{1}x\right)}^{{\gamma }_1}={\left({\lambda }_{2}x\right)}^{{\beta }_2}+{\left({\eta }_{2}x\right)}^{{\gamma }_2}.$$

Using the same reasoning as above, we finally obtain ${\mathsf{\theta }}_1={\mathsf{\theta }}_2$.

Appendix B. The log-Likelihood Function and Its First-Order Partial Derivatives

We rewrite the log-likelihood in a compact form: $$l=\sum _{i=1}^n\left(logh(x_i;\mathsf{\theta })-H(x_i;\mathsf{\theta })\right),$$ where $h(x_i;\theta )$ is the failure rate function and $H(x_i;\mathsf{\theta })$ is the cumulative failure rate function. Then its first-order partial derivatives can be derived as follows: $$\begin{array}{cc}\hfill \frac{\partial l}{\partial {\theta }_k}& =\sum _{i=1}^n\left(\frac{h_{{\theta }_k}(x_i;\mathsf{\theta })}{h(x_i;\mathsf{\theta })}-H_{{\theta }_k}(x_i;\mathsf{\theta })\right),k=\overline{1,6},\hfill \end{array}$$ where: $$\begin{array}{cc}\hfill h(x_i;\mathsf{\theta })& =\alpha \theta {\left(\theta x_i\right)}^{\alpha -1}+\beta \lambda {\left(\lambda x_i\right)}^{\beta -1}+\gamma \eta {\left(\eta x_i\right)}^{\gamma -1},\hfill \\ \hfill h_{\alpha }(x_i;\mathsf{\theta })& =\theta {\left(\theta x_i\right)}^{\alpha -1}(1+\alpha log\left(\theta x_i\right)),\hfill \\ \hfill h_{\beta }(x_i;\mathsf{\theta })& =\lambda {\left(\lambda x_i\right)}^{\beta -1}(1+\beta log\left(\lambda x_i\right)),\hfill \\ \hfill h_{\gamma }(x_i;\mathsf{\theta })& =\eta {\left(\eta x_i\right)}^{\gamma -1}(1+\gamma log\left(\eta x_i\right)),\hfill \\ \hfill h_{\theta }(x_i;\mathsf{\theta })& ={\alpha }^2{\left(\theta x_i\right)}^{\alpha -1},\hfill \\ \hfill h_{\lambda }(x_i;\mathsf{\theta })& ={\beta }^2{\left(\lambda x_i\right)}^{\beta -1},\hfill \\ \hfill h_{\eta }(x_i;\mathsf{\theta })& ={\gamma }^2{\left(\eta x_i\right)}^{\gamma -1},\hfill \\ \hfill H(x_i;\mathsf{\theta })& ={\left(\theta x_i\right)}^{\alpha }+{\left(\lambda x_i\right)}^{\beta }+{\left(\eta x_i\right)}^{\gamma },\hfill \\ \hfill H_{\alpha }(x_i;\mathsf{\theta })& ={\left(\theta x_i\right)}^{\alpha }log\left(\theta x_i\right),\hfill \\ \hfill H_{\beta }(x_i;\mathsf{\theta })& ={\left(\lambda x_i\right)}^{\beta }log\left(\lambda x_i\right),\hfill \\ \hfill H_{\gamma }(x_i;\mathsf{\theta })& ={\left(\eta x_i\right)}^{\gamma }log\left(\eta x_i\right),\hfill \\ \hfill H_{\theta }(x_i;\mathsf{\theta })& =\alpha x_i{\left(\theta x_i\right)}^{\alpha -1},\hfill \\ \hfill H_{\lambda }(x_i;\mathsf{\theta })& =\beta x_i{\left(\lambda x_i\right)}^{\beta -1},\hfill \\ \hfill H_{\eta }(x_i;\mathsf{\theta })& =\gamma x_i{\left(\eta x_i\right)}^{\gamma -1}.\hfill \end{array}$$

Appendix C. The log-Likelihood Function and Its Second-Order Partial Derivatives

The second-order partial derivatives of the log-likelihood function are given as: $$\begin{array}{cc}\hfill \frac{{\partial }^{2}l}{\partial {\theta }_k\partial {\theta }_l}& =\sum _{i=1}^n\left(\frac{h_{{\theta }_k{\theta }_l}(x_i;\mathsf{\theta })h(x_i;\mathsf{\theta })-h_{{\theta }_k}(x_i;\mathsf{\theta })h_{{\theta }_l}(x_i;\mathsf{\theta })}{h^2(x_i;\mathsf{\theta })}-H_{{\theta }_k{\theta }_l}(x_i;\mathsf{\theta })\right),k=\overline{1,6},l=\overline{k,6},\hfill \end{array}$$ where $$\begin{array}{cc}\hfill h_{\alpha \alpha }(x_i;\mathsf{\theta })& =\theta {\left(\theta x_i\right)}^{\alpha -1}log\left(\theta x_i\right)(2+\alpha log\left(\theta x_i\right)),\hfill \\ \hfill h_{\alpha \theta }(x_i;\mathsf{\theta })& =\alpha {\left(\theta x_i\right)}^{\alpha -1}(2+\alpha log\left(\theta x_i\right)),\hfill \\ \hfill h_{\theta \theta }(x_i;\mathsf{\theta })& ={\alpha }^2(\alpha -1)x_i{\left(\theta x_i\right)}^{\alpha -2},\hfill \\ \hfill h_{\beta \beta }(x_i;\mathsf{\theta })& =\lambda {\left(\lambda x_i\right)}^{\beta -1}log\left(\lambda x_i\right)(2+\beta log\left(\lambda x_i\right)),\hfill \\ \hfill h_{\beta \lambda }(x_i;\mathsf{\theta })& =\beta {\left(\lambda x_i\right)}^{\beta -1}(2+\beta log\left(\lambda x_i\right)),\hfill \\ \hfill h_{\lambda \lambda }(x_i;\mathsf{\theta })& ={\beta }^2(\beta -1)x_i{\left(\lambda x_i\right)}^{\beta -2},\hfill \\ \hfill h_{\gamma \gamma }(x_i;\mathsf{\theta })& =\eta {\left(\eta x_i\right)}^{\gamma -1}log\left(\eta x_i\right)(2+\gamma log\left(\eta x_i\right)),\hfill \\ \hfill h_{\gamma \eta }(x_i;\mathsf{\theta })& =\gamma {\left(\eta x_i\right)}^{\gamma -1}(2+\gamma log\left(\eta x_i\right)),\hfill \\ \hfill h_{\eta \eta }(x_i;\mathsf{\theta })& ={\gamma }^2(\gamma -1)x_i{\left(\eta x_i\right)}^{\gamma -2},\hfill \\ \hfill h_{\alpha \beta }(x_i;\mathsf{\theta })& =h_{\alpha \gamma }(x_i;\mathsf{\theta })=h_{\alpha \lambda }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill h_{\alpha \eta }(x_i;\mathsf{\theta })& =h_{\beta \gamma }(x_i;\mathsf{\theta })=h_{\beta \theta }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill h_{\beta \eta }(x_i;\mathsf{\theta })& =h_{\gamma \theta }(x_i;\mathsf{\theta })=h_{\gamma \lambda }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill h_{\theta \lambda }(x_i;\mathsf{\theta })& =h_{\theta \eta }(x_i;\mathsf{\theta })=h_{\lambda \eta }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill H_{\alpha \alpha }(x_i;\mathsf{\theta })& ={\left(\theta x_i\right)}^{\alpha }{log}^2\left(\theta x_i\right),\hfill \\ \hfill H_{\alpha \theta }(x_i;\mathsf{\theta })& =x_i{\left(\theta x_i\right)}^{\alpha -1}(1+\alpha log\left(\theta x_i\right)),\hfill \\ \hfill H_{\theta \theta }(x_i;\mathsf{\theta })& =\alpha (\alpha -1)x_i^2{\left(\theta x_i\right)}^{\alpha -2},\hfill \\ \hfill H_{\beta \beta }(x_i;\mathsf{\theta })& ={\left(\lambda x_i\right)}^{\beta }{log}^2\left(\lambda x_i\right),\hfill \\ \hfill H_{\beta \lambda }(x_i;\mathsf{\theta })& =x_i{\left(\lambda x_i\right)}^{\beta -1}(1+\beta log\left(\lambda x_i\right)),\hfill \\ \hfill H_{\lambda \lambda }(x_i;\mathsf{\theta })& =\beta (\beta -1)x_i^2{\left(\lambda x_i\right)}^{\beta -2},\hfill \\ \hfill H_{\gamma \gamma }(x_i;\mathsf{\theta })& ={\left(\eta x_i\right)}^{\gamma }{log}^2\left(\eta x_i\right),\hfill \\ \hfill H_{\gamma \eta }(x_i;\mathsf{\theta })& =x_i{\left(\eta x_i\right)}^{\gamma -1}(1+\gamma log\left(\eta x_i\right)),\hfill \\ \hfill H_{\eta \eta }(x_i;\mathsf{\theta })& =\gamma (\gamma -1)x_i^2{\left(\eta x_i\right)}^{\gamma -2},\hfill \\ \hfill H_{\alpha \beta }(x_i;\mathsf{\theta })& =H_{\alpha \gamma }(x_i;\mathsf{\theta })=H_{\alpha \lambda }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill H_{\alpha \eta }(x_i;\mathsf{\theta })& =H_{\beta \gamma }(x_i;\mathsf{\theta })=H_{\beta \theta }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill H_{\beta \eta }(x_i;\mathsf{\theta })& =H_{\gamma \theta }(x_i;\mathsf{\theta })=H_{\gamma \lambda }(x_i;\mathsf{\theta })=0,\hfill \\ \hfill H_{\theta \lambda }(x_i;\mathsf{\theta })& =H_{\theta \eta }(x_i;\mathsf{\theta })=H_{\lambda \eta }(x_i;\mathsf{\theta })=0.\hfill \end{array}$$

Appendix D. First-Order Partial Derivatives of the Reliability Function

Given the reliability function, $$R(x;\mathsf{\theta })=exp\left[-{\left(\theta x\right)}^{\alpha }-{\left(\lambda x\right)}^{\beta }-{\left(\eta x\right)}^{\gamma }\right],$$ its first-order partial derivatives with respect to the parameters $\alpha ,\beta ,\gamma ,\theta ,\lambda$, and $\eta$ are, respectively, given by: $$\begin{array}{cc}\hfill R_{\alpha }(x;\mathsf{\theta })& =-{\left(\theta x\right)}^{\alpha }log\left(\theta x\right)R(x;\mathsf{\theta }),\hfill \\ \hfill R_{\beta }(x;\mathsf{\theta })& =-{\left(\lambda x\right)}^{\beta }log\left(\lambda x\right)R(x;\mathsf{\theta }),\hfill \\ \hfill R_{\gamma }(x;\mathsf{\theta })& =-{\left(\eta x\right)}^{\gamma }log\left(\eta x_i\right)R(x;\mathsf{\theta }),\hfill \\ \hfill R_{\theta }(x;\mathsf{\theta })& =-\alpha x{\left(\theta x\right)}^{\alpha -1}R(x;\mathsf{\theta }),\hfill \\ \hfill R_{\lambda }(x;\mathsf{\theta })& =-\beta x{\left(\lambda x\right)}^{\beta -1}R(x;\mathsf{\theta }),\hfill \\ \hfill R_{\eta }(x;\mathsf{\theta })& =-\gamma x{\left(\eta x\right)}^{\gamma -1}R(x;\mathsf{\theta }).\hfill \end{array}$$

Appendix E. Cross-Entropy (CE) Algorithm for Continuous Optimization

In our previous study [37], we have described the CE algorithm in a generic case on how it works. Here, we focus more on the specific case of this study. The objective is to find ${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}$ that maximizes $l\left(\mathsf{\theta }\right)$, (A3)$${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}=arg\underset{\theta \in \Theta }{max}l\left(\mathsf{\theta }\right),$$ where $l\left(\theta \right)$ is the log-likelihood function, and $\Theta =\{{(\alpha ,\beta ,\gamma ,\theta ,\lambda ,\eta )}^T:\theta ,\lambda ,\eta >0\mathrm{and}0<\alpha <\beta <\gamma \}$. For this continuous optimization, let $N_6(\mu ,{\sigma }^2)$ denotes the 6-dimensional normal distribution with independent components, mean vector $\mu =({\mu }_1,\dots ,{\mu }_6)$ and variance vector ${\sigma }^2=({\sigma }_1^2,\dots ,{\sigma }_6^2)$. Then, the CE algorithm is proceeded as follows (Algorithm A1):

Algorithm A1 CE algorithm

1: Choose ${\widehat{\mu }}_0$ and ${\widehat{\sigma }}_0$. Set $t:=0$. 2: while${max}_j\left\{{\tilde{\sigma }}_{tj}\right\}\ge \epsilon$do 3:     $t=t+1$. 4:     Generate ${\mathsf{\theta }}_1,\dots ,{\mathsf{\theta }}_N$ from the $N_6({\widehat{\mu }}_{t-1},{\widehat{\sigma }}_{t-1}^2)$ distribution. 5:     Calculate the performances $l\left({\mathsf{\theta }}_k\right)$ for all k and order them from smallest to largest, $l_{\left(1\right)}\le \dots \le l_{\left(N\right)}$. 6:     Find the sample $(1-\varrho )$-quantile ${\widehat{\gamma }}_t=l_{(\lceil (1-\varrho )N\rceil )}$, where ${10}^{-2}\le \varrho \le {10}^{-1}$. 7:     Let ${\mathcal{E}}_t=\{{\mathsf{\theta }}_k:l\left({\mathsf{\theta }}_k\right)\ge {\widehat{\gamma }}_t\}$ be best performing (=elite) samples. 8:     for $j=1,\dots ,6$ do 9:         ${\tilde{\mu }}_{tj}:=\frac{{\sum }_{{\theta }_k\in {\mathcal{E}}_t}{\theta }_{ij}}{|{\mathcal{E}}_t|}$, 10:         ${\tilde{\sigma }}_{tj}:=\sqrt{\frac{{\sum }_{{\theta }_k\in {\mathcal{E}}_t}{({\theta }_{ij}-{\tilde{\mu }}_{tj})}^2}{|{\mathcal{E}}_t|}}$. 11:     end for 12:     ${\widehat{\mu }}_t:=\alpha {\tilde{\mu }}_t+(1-\alpha ){\widehat{\mu }}_{t-1}$, 13:     ${\widehat{\sigma }}_t:={\beta }_t{\tilde{\sigma }}_t+(1-{\beta }_t){\widehat{\sigma }}_{t-1}$, where $0.5\le \alpha \le 0.9$, ${\beta }_t=\beta -\beta {(1-\frac{1}{t})}^q$, q is an integer (typically between 5 and 10) and $\beta$ is a smoothing constant (typically between 0.8 and 0.99). 14: end while

While applying the CE algorithm, the mean vector ${\widehat{\mu }}_t$ should converge to ${\widehat{\mathsf{\theta }}}_{\mathrm{ML}}$ and the standard deviation vector ${\widehat{\sigma }}_t$ should converge to zero vector. We refer to [38] for more details about the method. To solve (A3) with constraints on the parameters, we can apply the acceptance-rejection (AR) approach. The AR method works as follows: Generate a random vector $\mathsf{\theta }$ from a normal distribution, then accept or reject it depending on whether the sample falls or not in the constrained parameter space [39].

Appendix F. R Code for the Cross-Entropy Method

We provide R code for the cross-entropy method by using the “CEoptim” package [40]. Here, we show the R code for the Aarset data analyzed in Section 6.1. Since the log-likelihood function is highly multi-model and CE method is a stochastic search algorithm, we might need to run the CEoptim function several times in order to obtain the best result.

• Install and invoke the CEoptim package

  install.packages("CEoptim")

  library(CEoptim)

• Input the Aarset data into R

  Aarsetdata <- c(

  136, 246, 255, 376, 421, 565, 616, 617, 652, 655,

  658, 660, 662, 675, 681, 734, 736, 737, 757, 769,

  777, 800, 807, 825, 855, 857, 864, 868, 870, 870,

  873, 882, 895, 910, 934, 943, 1015, 1019

  )

• Define the log-likelihood function

  LogLikelihood <- function(x,t) {

  sum(log(x[1]*x[4]*(t*x[4])^(x[1]-1) + x[2]*x[5]*(t*x[5])^(x[2]-1) +

  x[3]*x[6]*(t*x[6])^(x[3]-1))-((t*x[4])^x[1] + (t*x[5])^x[2]

  + (t*x[6])^x[3]))

  }

• Choose initial values for mean and standard deviation

  mu0 <- c(0,5,90,0,0,0)

  sigma0 <- c(1,10,10,1,1,1)

• Restrict the parameter space

  A=matrix(c(-1,1,0,0,0,0, 0,-1,1,0,0,0, 0,0,-1,0,0,0,

  0,0,0,-1,0,0, 0,0,0,0,-1,0, 0,0,0,0,0,-1), nrow = 6)

  B=c(0,0,0,0,0,0)

• Execute the CEoptim function

  opt <- CEoptim(

  LogLikelihood, f.arg = list(t=Aarsetdata), maximize = T,

  continuous = list(mean=mu0,sd=sigma0, conMat=A, conVec=B,

  smoothMean=0.7, smoothSd=0.7, sdThr=0.00001), rho = 0.01, N=1000

  )

• The optimum, optimizer (MLE), and AIC values are respectively given by

  opt$optimum

  opt$optimizer$continuous

  AIC <- -2*opt$optimum+2*6

References

1. Xie, M.; Lai, C.D. Reliability analysis using an additive Weibull model with bathtub-shaped failure rate function. Reliab. Eng. Syst.; 1995; 52, pp. 87-93. [DOI: https://dx.doi.org/10.1016/0951-8320(95)00149-2]

2. Weibull, W.A. Statistical distribution function of wide applicability. J. Appl. Mech.; 1951; 18, pp. 293-296. [DOI: https://dx.doi.org/10.1115/1.4010337]

3. Lemonte, A.J.; Cordeiroa, G.M.; Ortegab, E.M.M. On the additive Weibull distribution. Commun. Stat. Theory Methods; 2014; 43, pp. 2066-2080. [DOI: https://dx.doi.org/10.1080/03610926.2013.766343]

4. Almalki, S.J.; Yuan, J. A new modified Weibull distribution. Reliab. Eng. Syst.; 2013; 111, pp. 164-170. [DOI: https://dx.doi.org/10.1016/j.ress.2012.10.018]

5. Lai, C.D.; Xie, M.; Murthy, D.N.P. A modified Weibull distribution. IEEE Trans. Reliab.; 2003; 52, pp. 33-37. [DOI: https://dx.doi.org/10.1109/TR.2002.805788]

6. Thach, T.T.; Bris, R. Improved new modified weibull distribution: A bayes study using hamiltonian monte carlo simulation. Proc. Inst. Mech. Eng. O J. Risk Reliab.; 2020; 234, pp. 496-511. [DOI: https://dx.doi.org/10.1177/1748006X19896740]

7. He, B.; Cui, W.; Du, X. An additive modified Weibull distribution. Reliab. Eng. Syst.; 2016; 145, pp. 28-37. [DOI: https://dx.doi.org/10.1016/j.ress.2015.08.010]

8. Gompertz, B. On the nature of the function expressive of the law of human mortality and on the new model of determining the value of life contingencies. Phil. Trans. R. Soc. Lon.; 1825; 115, pp. 513-585.

9. Park, S.; Park, J. A general class of flexible Weibull distributions. Commun. Stat. Theory Methods; 2018; 47, pp. 767-778. [DOI: https://dx.doi.org/10.1080/03610926.2015.1118509]

10. Shakhatreh, M.K.; Lemonte, A.J.; Moreno-Arenas, G. The log-normal modified Weibull distribution and its reliability implications. Reliab. Eng. Syst.; 2019; 188, pp. 6-12. [DOI: https://dx.doi.org/10.1016/j.ress.2019.03.014]

11. Thach, T.T.; Bris, R. An additive Chen-Weibull distribution and its applications in reliability modeling. Qual. Reliab. Eng. Int.; 2021; 37, pp. 352-373. [DOI: https://dx.doi.org/10.1002/qre.2740]

12. Kumar, C.S.; Nair, S.R. On some aspects of a flexible class of additive Weibull distribution. Commun. Stat. Theory Methods; 2018; 47, pp. 1028-1049. [DOI: https://dx.doi.org/10.1080/03610926.2017.1316399]

13. Mudholkar, G.S.; Srivastava, D.K. Exponentiated Weibull family for analysing bathtub failure rate data. IEEE Trans. Reliab.; 1993; 42, pp. 299-302. [DOI: https://dx.doi.org/10.1109/24.229504]

14. Hjorth, U. A reliability distribution with increasing, decreasing, constant and bathtub-shaped failure rates. Technometrics; 1980; 17, pp. 99-107. [DOI: https://dx.doi.org/10.2307/1268388]

15. Carrasco, M.; Ortega, E.M.; Cordeiro, G.M. A generalized modified Weibull distribution for lifetime modeling. Comput. Stat. Data Anal.; 2008; 53, pp. 450-462. [DOI: https://dx.doi.org/10.1016/j.csda.2008.08.023]

16. Silva, G.O.; Ortega, E.M.; Cordeiro, G.M. The beta modified Weibull distribution. Lifetime Data Anal.; 2010; 16, pp. 409-430. [DOI: https://dx.doi.org/10.1007/s10985-010-9161-1]

17. Zhang, T.; Xie, M. On the upper truncated Weibull distribution and its reliability implications. Reliab. Eng. Syst.; 2011; 96, pp. 194-200. [DOI: https://dx.doi.org/10.1016/j.ress.2010.09.004]

18. Lemonte, A.J. A new exponential-type distribution with constant, decreasing, increasing, upside-down bathtub and bathtub-shaped failure rate function. Comput. Stat. Data Anal.; 2013; 62, pp. 149-170. [DOI: https://dx.doi.org/10.1016/j.csda.2013.01.011]

19. Singla, N.; Jain, K.; Sharma, S.K. The beta Generalized Weibull distribution: Properties and applications. Reliab. Eng. Syst.; 2012; 102, pp. 5-15. [DOI: https://dx.doi.org/10.1016/j.ress.2012.02.003]

20. Nassar, M.; Afify, A.Z.; Dey, S.; Kumar, D. A new extension of Weibull distribution: Properties and different methods of estimation. J. Comput. Appl. Math.; 2018; 336, pp. 439-457. [DOI: https://dx.doi.org/10.1016/j.cam.2017.12.001]

21. Ali, S.; Dey, S.; Tahir, M.H.; Mansoor, M. Two-Parameter Logistic-Exponential Distribution: Some New Properties and Estimation Methods. Am. J. Math. Manag. Sci.; 2020; 39, pp. 270-298. [DOI: https://dx.doi.org/10.1080/01966324.2020.1728453]

22. Shih, J.H.; Emura, T. Penalized Cox regression with a five-parameter spline model. Commun. Stat. Theory Methods; 2021; 50, pp. 3749-3768. [DOI: https://dx.doi.org/10.1080/03610926.2020.1772305]

23. Nassar, M.; Afify, A.Z.; Shakhatreh, M.K.; Dey, S. On a new extension of Weibull distribution: Properties, estimation, and applications to one and two causes of failures. Qual. Reliab. Eng. Int.; 2020; 36, pp. 2019-2043. [DOI: https://dx.doi.org/10.1002/qre.2671]

24. Shakhatreh, M.K.; Lemonte, A.J.; Cordeiro, G.M. On the generalized extended exponential-Weibull distribution: Properties and different methods of estimation. Int. J. Comput. Math.; 2020; 97, pp. 1029-1057. [DOI: https://dx.doi.org/10.1080/00207160.2019.1605062]

25. Mallick, A.; Ghosh, I.; Dey, S.; Kumar, D. Bounded Weighted Exponential Distribution with Applications. Am. J. Math. Manag. Sci.; 2021; 40, pp. 68-87. [DOI: https://dx.doi.org/10.1080/01966324.2020.1834893]

26. Khalil, A.; Ijaz, M.; Ali, K.; Mashwani, W.K.; Shafiq, M.; Kumam, P.; Kumam, W. A novel flexible additive Weibull distribution with real-life applications. Commun. Stat. Theory Methods; 2021; 50, pp. 1557-1572. [DOI: https://dx.doi.org/10.1080/03610926.2020.1732658]

27. Salem, S.A. Bayesian estimation of a non-linear failure rate from censored samples type II. Microelectron. Reliab.; 1992; 32, pp. 1385-1388. [DOI: https://dx.doi.org/10.1016/0026-2714(92)90007-8]

28. Bain, L.J. Analysis for the linear failure-rate life-testing distribution. Technometrics; 1974; 16, pp. 551-559. [DOI: https://dx.doi.org/10.1080/00401706.1974.10489237]

29. Rubinstein, R.Y.; Kroese, D.P. The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte Carlo Simulation and Machine Learning; Springer: New York, NY, USA, 2004.

30. Geyer, C.J. Fisher Information and Confidence Intervals Using Maximum Likelihood. Available online: http://www.stat.umn.edu/geyer/old03/5102/notes/fish.pdf (accessed on 10 August 2021).

31. Greene, W. Econometric Analysis; 8th ed. Prentice-Hall: New York, NY, USA, 2020; pp. 118-121.

32. EL-Sagheer, R. Estimation of parameters of Weibull Gamma distribution based on progressively censored data. Stat. Pap.; 2018; 59, pp. 725-757. [DOI: https://dx.doi.org/10.1007/s00362-016-0787-2]

33. Meeker, W.Q.; Escobar, L.A. Statistical Methods for Reliability Data; 2nd ed. John Wiley: New York, NY, USA, 1998.

34. Aarset, M.V. How to identify bathtub hazard rate. IEEE Trans. Reliab.; 1987; 36, pp. 106-108. [DOI: https://dx.doi.org/10.1109/TR.1987.5222310]

35. Bergman, B.; Klefsjou, B. Burn–in models and TTT-transforms. Qual. Reliab. Eng. Int.; 1985; 1, pp. 125-130. [DOI: https://dx.doi.org/10.1002/qre.4680010211]

36. Gentleman, R. Muhaz: A Package for Producing a Smooth Estimate of the Hazard Function for Censored Data. Available online: https://cran.r-project.org/web/packages/muhaz/muhaz.pdf (accessed on 20 August 2021).

37. Thach, T.T.; Bris, R.; Volf, P.; Coolen, F.P.A. Non-linear failure rate: A Bayes study using Hamiltonian Monte Carlo simulation. Int. J. Approx. Reason.; 2020; 123, pp. 55-76. [DOI: https://dx.doi.org/10.1016/j.ijar.2020.04.007]

38. Rubinstein, R.Y.; Kroese, D.P. Simulation and the Monte Carlo Method; 3rd ed. John Wiley & Sons: Hoboken, NJ, USA, 2016; pp. 257-303.

39. Kroese, D.P.; Porotsky, S.; Rubinstein, R.Y. The Cross-Entropy Method for Continuous Multi-Extremal Optimization. Methodol. Comput. Appl. Probab.; 2006; 8, pp. 383-407. [DOI: https://dx.doi.org/10.1007/s11009-006-9753-0]

40. Benham, T.; Duan, Q.; Kroese, D.P.; Liquet, B. CEoptim: Cross-Entropy R Package for Optimization. J. Stat. Softw.; 2017; 76, pp. 1-29. [DOI: https://dx.doi.org/10.18637/jss.v076.i08]