Statistical Analysis of the Lifetime Distribution with Bathtub-Shaped Hazard Function under Lagged-Effect Step-Stress Model

Abstract

In survival analysis, applying stress is often used to accelerate an experiment. Stress can be discontinuous, and the step-stress model is applied widely due to its flexibility. However, in reality, when new stress is applied, it often does not take effect immediately, but there will be a lagged effect. Under the lagged-effect step-stress model, the statistical inference of the Chen distribution is discussed. The Chen distribution is an important life distribution as its risk function is bathtub-shaped with certain parameters. In this paper, the maximum likelihood estimators are presented and the Newton–Raphson algorithm is used. According to the form of risk function under this model, the explicit expressions of least squares estimators are obtained. The calculation methods of asymptotic confidence intervals and coverage probabilities are proposed by using the observed Fisher matrix. Finally, to evaluate the performance of the above estimation methods, a Monte Carlo simulation study is provided.

1. Introduction

1.1. Chen Distribution

In survival analysis, hazard function plays an important role in studying the life phenomenon of a product. For many products, their failure rates decrease first, then keep at a constant level, and increase finally. Such failure rate is like a bathtub, and this life distribution is widely used in electronic, machinery, and medical fields. For example, some drugs do not work well for children and the elderly, but they work well for middle-aged people. In other words, the failure rate of drugs is relatively high in childhood but gradually decreases with age. Then, the failure rate remains low in middle age for some time and eventually increases with age. One of the life distributions with such hazard function is the Chen distribution, which was first proposed by ref. [1]. It is a two-parameter lifetime distribution with the bathtub-shaped or increasing hazard function and can model the real data well.

Ref. [1] proposed confidence intervals and joint confidence regions for the Chen distribution’s parameters under Type-II censoring. Ref. [2] investigated a simple method to conduct the statistical test and obtain the exact confidence interval of the Chen distribution’s shape parameter, which can also be applied to models under Type-II right censoring. Based on Type-II right-censored samples of the Chen distribution, ref. [3] later discussed several test statistics for an exact hypothesis test concerning the shape parameter. Ref. [4] obtained the point estimations and interval estimations for the parameters under a Type-II censored model. Ref. [5] proposed an extended maximum spacing method to estimate parameters of the Chen distribution. Under hybrid censoring, ref. [6] discussed the maximum likelihood estimations and several asymptotic confidence intervals. They also used the Lindley method, and the Tierney and Kadane method, to calculate Bayes estimates. Based on the Chen distribution, ref. [7] analyzed the stress–strength reliability under progressive Type-II censoring and generalized it to the proportional hazard family. Under progressively censored samples of the Chen distribution, ref. [8] discussed maximum likelihood estimates, different Bayes estimates, asymptotic confidence intervals, and prediction intervals. Based on data from the Chen distribution, ref. [9] developed simplified forms of the single moments and covariances. The estimates of the shape parameters as well as the prediction of the records are also proposed.

A Chen (Chen($\beta$,$\lambda$)) random variable T with two positive parameters $\beta$ (≥0) and $\lambda$ (≥0) has the following probability density function (pdf):

$$\begin{array}{c}\hfill f(t;\beta ,\lambda )=\lambda \beta t^{\beta -1}{e}^{t^{\beta }}exp\left\{\lambda \left(1-e^{t^{\beta }}\right)\right\},t>0.\end{array}$$

(1)

The cumulative distribution function and the survival function are, respectively, given by:

$$\begin{array}{c}\hfill F(t;\beta ,\lambda )=1-exp\left\{\lambda \left(1-e^{t^{\beta }}\right)\right\},t>0.\end{array}$$

(2)

$$\begin{array}{c}\hfill S(t;\beta ,\lambda )=1-F(t;\beta ,\lambda )=exp\left\{\lambda \left(1-e^{t^{\beta }}\right)\right\},t>0.\end{array}$$

(3)

Accordingly, the hazard function is:

$$\begin{array}{c}\hfill h(t;\beta ,\lambda )=\frac{f(t;\beta ,\lambda )}{S(t;\beta ,\lambda )}=\lambda \beta t^{\beta -1}{e}^{t^{\beta }},t>0.\end{array}$$

(4)

The shape of the pdf varies with the parameters and the characteristics are summarized as follows: (1) If $0<\beta <1$: the pdf will decrease throughout or decrease first and then increase when $0<\lambda <1$; the pdf will decrease throughout when $\lambda \ge 1$. (2) If $\beta =1$: the pdf will be unimodal when $0<\lambda <1$; the pdf will decrease throughout when $\lambda \ge 1$. (3) If $\beta >1$: the pdf will always be unimodal no matter which value $\lambda$ takes. Different plots of pdf are shown in Figure 1, Figure 2, Figure 3 and Figure 4, respectively.

Take the derivative of $h(t;\beta ,\lambda )$ with respect to t, then $h^{\prime }(t;\beta ,\lambda )=\lambda \beta t^{\beta -2}{e}^{t^{\beta }}[(\beta -1)+\beta t^{\beta }]$. Thus, the hazard function shows different shapes when $\beta$ differs and the properties are as follows: (1) The hazard function is bathtub-shaped when $0<\beta <1$. (2) The hazard function increases throughout when $\beta \ge 1$. The corresponding plots are shown in Figure 5 and Figure 6, respectively.

1.2. Step-Stress Model with Lagged Effect

Nowadays, due to the development of science and technology, the life of a product is getting longer and longer, and waiting for the product to fail will cause a great waste of time and manpower. Therefore, some measures need to be taken to accelerate product failure. Applying stress is a common means to accelerate the experiment in life test and reliability analysis, which can reduce time waste and other related costs. Stress can be voltage, temperature, oxygen, etc. There are three common stress-application schemes: constant-stress model, step-stress model, and progressive-stress model. Under the constant-stress model, the stress remains unchanged until the products fail. The increase in stress in the progressive-stress model is linear and continuous. In the step-stress model, the stress can be changed, but it does not have to be changed continuously, and sudden change is allowed. In this paper, a simple step-stress model is considered: at first, the initial stress level lasts for some time, and then at a given time, the stress level increases and remains unchanged until all products fail.

The cumulative exposure model (CEM) is a commonly used step-stress model, which assumes that the remaining life of the product is only associated with the cumulative exposure experienced previously and current stress. Ref. [10] first proposed the concept of the cumulative exposure model. Ref. [11] used the CEM to analyze step-stress data of the Weibull distribution and presented the maximum likelihood estimation and interval estimation under this method. Ref. [12] then presented the optimum scheme of the model, including the optimum duration of the first stress, the optimum proportion failing, and the asymptotic variance. Ref. [13] took into account the multiplier effect of stress, calculating the maximum likelihood estimation of the Weibull family of functions and the Fisher information matrix. Under CEM, ref. [14] discussed the maximum likelihood estimation and interval estimations of the exponential distribution with Type-I hybrid censoring. Ref. [15] later considered the Type-II censoring and independent competing risks in the model. Ref. [16] proposed the optimal life tests of the Weibull distribution using the Bayesian method under the model and used two algorithms to optimize it. Under CEM, Ref. [17] discussed the maximum likelihood estimation of the Weibull distribution with Type-I progressive hybrid censoring. Based on Type-II progressive hybrid censoring, Ref. [18] discussed statistical inference and optimal design on a step-stress partially accelerated life test for a hybrid system in the presence of masked data.

Although the CEM is widely used, ref. [19] pointed out that the hazard function is discontinuous when the stress level changes. That is, the impact of stress change is instantaneous. In reality, when the stress level changes, it often does not take effect immediately, but there exists a lag period. The CEM is unreasonable and inappropriate in this case. To solve this problem, the cumulative risk model (CRM) is proposed, which takes into account the lagged effect. In this model, the risk function is continuous, and it is supposed that the lagged effect causes a linear risk function in the intermediate period, which is more consistent with reality. Ref. [19] first proposed the concept of the cumulative risk model, and discussed the maximum likelihood estimation and least squares estimation of the model under exponential distribution. Ref. [20] combined the CRM with the degradation test model for data analysis. Ref. [21] took into account competing risks under the exponential distribution. In addition to calculating the maximum likelihood estimation, it also used three methods to calculate the confidence interval and coverage probabilities. Ref. [22] later extended this model to the Weibull distribution, and took the competing risks into account. Under masked data, ref. [23] also introduced competing risks based on the CRM. Ref. [24] applied the CRM to fuzzy lifetime data. Ref. [25] calculated the maximum likelihood estimation, the least square estimation, and Bayesian estimation under a Weibull cumulative risk model.

Many studies on the step-stress model consider the CEM, but the CRM is more in line with reality. In addition, most of the existing research on the CRM only involves the exponential distribution or Weibull distribution. From the point of view of the hazard function, although the Weibull distribution is widely used in survival and reliability analysis, its hazard function can only be monotonic or constant. Compared with the Weibull distribution, the Chen distribution has a hazard function that can not only be monotonous but also show the shape of the bathtub, which is important in practical fields. Statistical analysis based on the Chen distribution can make applications of the CRM deeper and wider. In this article, the Chen distribution and step-stress with lagged effect model are both considered, which is of great significance in theory and practice.

It is assumed that lifetime under the initial stress level obeys the Chen distribution. The stress level changes at ${\tau }_1$, which starts to take effect at ${\tau }_2$ due to the lagged effect, and the parameters of the Chen distribution change at ${\tau }_2$ as well. From ${\tau }_1$ to ${\tau }_2$, the hazard functions under these two stress levels are connected by a linear function.

The rest of the paper is arranged as follows. Some basic calculations and derivations of the model are shown in Section 2. In Section 3, the maximum likelihood estimation and least square estimation under the CRM are given. In Section 4, the asymptotic confidence intervals and coverage probabilities are discussed by using the large sample theory. To evaluate the performance of the estimators, the simulation results are presented in Section 5. Section 6 considers a special case where only one parameter changes when stress level changes. Section 7 is the summary of the article.

2. Model Description

Assume that the lifetime under the initial stress obeys Chen$({\beta }_1,{\lambda }_1)$. The new stress is applied at ${\tau }_1$, and it starts to take effect at ${\tau }_2$ (${\tau }_1$ and ${\tau }_2$ are known). The lifetime under the new stress obeys Chen$({\beta }_2,{\lambda }_2)$. From ${\tau }_1$ to ${\tau }_2$, the hazard function is linear and denoted as $a+bt$ (here, a and b are parameters).

The Chen hazard functions under the initial stress level and the second level are:

$$h_1\left(t\right)=\begin{array}{c}\hfill {\lambda }_1{\beta }_1{t}^{{\beta }_1-1}{e}^{t^{{\beta }_1}},t>0,\end{array}$$

(5)

$$h_2\left(t\right)=\begin{array}{c}\hfill {\lambda }_2{\beta }_2{t}^{{\beta }_2-1}{e}^{t^{{\beta }_2}},t>0.\end{array}$$

(6)

Under the CRM, the hazard function is given by:

$$h\left(t\right)=\left\{\begin{array}{c}\hfill {\lambda }_1{\beta }_1{t}^{{\beta }_1-1}{e}^{t^{{\beta }_1}},0<t<{\tau }_1\\ \hfill a+bt,{\tau }_1\le t<{\tau }_2\\ \hfill {\lambda }_2{\beta }_2{t}^{{\beta }_2-1}{e}^{t^{{\beta }_2}},t\ge {\tau }_2\end{array}\right..$$

(7)

Note that when ${\tau }_1$ = ${\tau }_2$, the hazard function $h_0\left(t\right)$ can be written as follows, which is the hazard function of the CEM as well:

$$h_0\left(t\right)=\left\{\begin{array}{c}\hfill {\lambda }_1{\beta }_1{t}^{{\beta }_1-1}{e}^{t^{{\beta }_1}},0<t<{\tau }_1\\ \hfill {\lambda }_2{\beta }_2{t}^{{\beta }_2-1}{e}^{t^{{\beta }_2}},t\ge {\tau }_1\end{array}\right.$$

(8)

In the CRM, we assume that ${\tau }_1\ne {\tau }_2$.

To make sure that the hazard function is continuous at ${\tau }_1$ and ${\tau }_2$, the following equations must be satisfied:

$$\left\{\begin{array}{c}\hfill {\lambda }_1{\beta }_1{{\tau }_1}^{{\beta }_1-1}{e}^{{\tau }_1^{{\beta }_1}}=a+b{\tau }_1\\ \hfill {\lambda }_2{\beta }_2{{\tau }_2}^{{\beta }_2-1}{e}^{\tau 2_2^{{\beta }_2}}=a+b{\tau }_2\end{array}\right.$$

(9)

According to (9), ${\lambda }_1$ and ${\lambda }_2$ can be solved as:

$$\left\{\begin{array}{c}\hfill {\lambda }_1=\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}\\ \hfill {\lambda }_2=\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}\end{array}\right..$$

(10)

The cumulative hazard function $H\left(t\right)$ under the model can be obtained by using the formula $H\left(t\right)={\int }_0^{t}h\left(x\right)dx$ and replacing the parameters ${\lambda }_1$ and ${\lambda }_2$ according to (10).

$$H\left(t\right)=\left\{\begin{array}{c}\hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{t^{{\beta }_1}}-1),0<t<{\tau }_1\\ \\ \hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)+a(t-{\tau }_1)+\frac{b}{2}(t^2-{\tau }_1^2),{\tau }_1\le t<{\tau }_2\\ \\ \hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)+a({\tau }_2-{\tau }_1)+\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\\ \hfill +\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}(e^{t^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}}),t\ge {\tau }_2\end{array}\right.$$

(11)

The survival function $S\left(t\right)$ under the model can be given as follows by the formula $S\left(t\right)=e^{-H\left(t\right)}$:

$$S\left(t\right)=\left\{\begin{array}{c}\hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{t^{{\beta }_1}}-1)\right\},0<t<{\tau }_1\\ \\ \hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\\ \\ \hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\\ \hfill \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}(e^{t^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\right\},t\ge {\tau }_2\end{array}\right.$$

(12)

According to the formula $f\left(t\right)=h\left(t\right)S\left(t\right)$, the probability density function $f\left(t\right)$ of the lifetime under the CRM is as follows:

$$f\left(t\right)=\left\{\begin{array}{c}\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{{\tau }_1}^{{\beta }_1-1}}{t}^{{\beta }_1-1}{e}^{t^{{\beta }_1}}exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{t^{{\beta }_1}}-1)\right\},0<t<{\tau }_1\hfill \\ \\ (a+bt)exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\hfill \\ \\ \frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{{\tau }_2}^{{\beta }_2-1}}{t}^{{\beta }_2-1}{e}^{t^{{\beta }_2}}exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\hfill \\ \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}(e^{t^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\right\},t\ge {\tau }_2\hfill \end{array}\right.$$

(13)

Thus, the corresponding cumulative distribution function $F\left(t\right)$ under the CRM is obtained by:

$$F\left(t\right)=\left\{\begin{array}{c}\hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{t^{{\beta }_1}}-1)\right\},0<t<{\tau }_1\\ \\ \hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\\ \\ \hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\\ \hfill \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}(e^{t^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\right\},t\ge {\tau }_2\end{array}\right..$$

(14)

Based on the above analysis, the estimations of the parameters can be given in the following section.

3. Point Estimation

3.1. Maximum Likelihood Estimation

Assume that $t_1,t_2,\cdots ,t_n$ are the failure times under the model. Among them, $n_1$ products fail during the first stress application (before ${\tau }_1$), $n_2$ products fail in the lag period (from ${\tau }_1$ to ${\tau }_2$), $n_3$ products fail during the second stress application (after ${\tau }_2$), and $n_1+n_2+n_3=n$.

The maximum likelihood estimation method (MLE) is a classical point estimation method and is widely used in estimating parameters. According to the theory of maximum likelihood estimation, the likelihood function can be written as follows. Denote it as $L({\beta }_1,{\beta }_2,a,b)$.

$$\begin{array}{c}\hfill L({\beta }_1,{\beta }_2,a,b)=\prod _{i=1}^{n}f\left(t_i\right)\end{array}$$

(15)

Plug (13) into (15), and the likelihood function can be expressed as:

$$\begin{array}{cc}\hfill L({\beta }_1,{\beta }_2,a,b)=& \prod _{i=1}^{n_1}\left[\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{{\tau }_1}^{{\beta }_1-1}}{t_i}^{{\beta }_1-1}{e}^{t_i^{{\beta }_1}}exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{t_i^{{\beta }_1}}-1)\right\}\right]\hfill \\ \hfill & \times \prod _{i=n_1+1}^{n_1+n_2}\left[(a+b{t}_i)exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-a(t_i-{\tau }_1)-\frac{b}{2}(t_i^2-{\tau }_1^2)\right\}\right]\hfill \\ \hfill & \times \prod _{i=n_1+n_2+1}^n[\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{{\tau }_2}^{{\beta }_2-1}}{t_i}^{{\beta }_2-1}{e}^{t^{{\beta }_2}}exp\left\{-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\hfill \\ \hfill & \times exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{{\beta }_1}}-1)-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}(e^{t_i^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\right\}].\hfill \end{array}$$

(16)

Based on the form of $L({\beta }_1,{\beta }_2,a,b)$, it can be seen that when $n_2=0,n_1=0$ or $n_2=0,n_3=0$, the maximum likelihood estimates (MLEs) do not exist. In the following, it is assumed that $n_i>0$.

The log-likelihood function $l({\beta }_1,{\beta }_2,a,b)$ is given by:

$$\begin{array}{cc}\hfill l({\beta }_1,{\beta }_2,a,b)=& lnL({\beta }_1,{\beta }_2,a,b)\hfill \\ \hfill =& n_1\left[ln(a+b{\tau }_1)-{{\tau }_1}^{{\beta }_1}\right]+({\beta }_1-1)\sum _{i=1}^{n_1}ln\left(\frac{t_i}{{\tau }_1}\right)+\sum _{i=1}^{n_1}{t_i}^{{\beta }_1}-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}\sum _{i=1}^{n_1}(e^{t_i^{{\beta }_1}}-1)\hfill \\ \hfill & +\sum _{i=n_1+1}^{n_1+n_2}ln(a+b{t}_i)-(n_2+n_3)\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{{\tau }_1}^{{\beta }_1}}-1)-a\sum _{i=n_1+1}^{n_1+n_2}(t_i-{\tau }_1)\hfill \\ \hfill & -\frac{b}{2}\sum _{i=n_1+1}^{n_1+n_2}(t_i^2-{\tau }_1^2)+n_3\left[ln(a+b{\tau }_2)-{{\tau }_2}^{{\beta }_2}\right]+({\beta }_2-1)\sum _{i=n_1+n_2+1}^{n}ln\left(\frac{t_i}{{\tau }_1}\right)+\sum _{i=n_1+n_2+1}^n{t_i}^{{\beta }_2}\hfill \\ \hfill & -a({\tau }_2-{\tau }_1)n_3-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)n_3-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}\sum _{i=n_1+n_2+1}^n(e^{t_i^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}}).\hfill \end{array}$$

(17)

In order to maximize the $l({\beta }_1,{\beta }_2,a,b)$, take partial derivatives in (17) with respect to ${\beta }_1$, ${\beta }_2$, a, and b. The results are as follows:

$$\begin{array}{cc}\hfill \frac{\partial l({\beta }_1,{\beta }_2,a,b)}{\partial {\beta }_1}=& -\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}\left\{-(\frac{1}{{\beta }_1}+ln{\tau }_1+{\tau }_1^{{\beta }_1}ln{\tau }_1)\left[\sum _{i=1}^{n_1}(e^{t_i^{{\beta }_1}}-1)+(n_2+n_3)(e^{{{\tau }_1}^{{\beta }_1}}-1)\right]\right.\hfill \\ \hfill & \left.+(n_2+n_3)(e^{{{\tau }_1}^{{\beta }_1}}{{\tau }_1}^{{\beta }_1}ln{\tau }_1)+\sum _{i=1}^{n_1}{e}^{t_i^{{\beta }_1}}{t}_i^{{\beta }_1}ln{t}_i\right\}-n_1{\tau }_1^{{\beta }_1}ln{\tau }_1+\sum _{i=1}^{n_1}ln\left(\frac{t_i}{{\tau }_1}\right)\hfill \\ \hfill & +\sum _{i=1}^{n_1}{t_i}^{{\beta }_1}ln{t}_i,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \frac{\partial l({\beta }_1,{\beta }_2,a,b)}{\partial {\beta }_2}=& -\frac{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}\left\{-(\frac{1}{{\beta }_2}+ln{\tau }_2+{\tau }_2^{{\beta }_2}ln{\tau }_2)\left[\sum _{i=n_1+n_2+1}^n(e^{t_i^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\right]\right.\hfill \\ \hfill & \left.+\sum _{i=n_1+n_2+1}^n(e^{t_i^{{\beta }_2}}{t}_i^{{\beta }_2}ln{t}_i-e^{{\tau }_2^{{\beta }_2}}{\tau }_2^{{\beta }_2}ln{\tau }_2)\right\}-n_3{\tau }_2^{{\beta }_2}ln{\tau }_2+\sum _{i=n_1+n_2+1}^{n}ln\left(\frac{t_i}{{\tau }_1}\right)\hfill \\ \hfill & +\sum _{i=n_1+n_2+1}^n{t_i}^{{\beta }_2}ln{t}_i,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{\partial l({\beta }_1,{\beta }_2,a,b)}{\partial a}=& \frac{n_1}{a+b{\tau }_1}-\frac{e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}\left[\sum _{i=1}^{n_1}(e^{t_i^{{\beta }_1}}-1)+(n_2+n_3)(e^{{{\tau }_1}^{{\beta }_1}}-1)\right]+\sum _{i=n_1+1}^{n_1+n_2}\frac{1}{a+b{t}_i}\hfill \\ \hfill & -\sum _{i=n_1+1}^{n_1+n_2}(t_i-t_1)+\frac{n_3}{a+b{\tau }_2}-({\tau }_2-{\tau }_1)n_3-\frac{e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-1}}\sum _{i=n_1+n_2+1}^n(e^{t_i^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}}),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{\partial l({\beta }_1,{\beta }_2,a,b)}{\partial b}=& \frac{n_1{\tau }_1}{a+b{\tau }_1}-\frac{e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-2}}\left[\sum _{i=1}^{n_1}(e^{t_i^{{\beta }_1}}-1)+(n_2+n_3)(e^{{{\tau }_1}^{{\beta }_1}}-1)\right]+\sum _{i=n_1+1}^{n_1+n_2}\frac{t_i}{a+b{t}_i}\hfill \\ \hfill & -\frac{1}{2}\sum _{i=n_1+1}^{n_1+n_2}(t_i^2-t_1^2)+\frac{n_3{\tau }_2}{a+b{\tau }_2}-\frac{1}{2}({\tau }_2^2-{\tau }_1^2)n_3-\frac{e^{-{\tau }_2^{{\beta }_2}}}{{\beta }_2{{\tau }_2}^{{\beta }_2-2}}\sum _{i=n_1+n_2+1}^n(e^{t_i^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}}).\hfill \end{array}$$

(21)

By making the above functions equal to 0 simultaneously, the maximum likelihood estimates of ${\beta }_1$, ${\beta }_2$, a, and b can be solved. However, explicit solutions cannot be given because the forms of the equations are complex and nonlinear. Therefore, some numerical techniques, such as the Newton–Raphson algorithm, can be used to calculate approximate estimates of parameters. This can be realized by using the $optim$ function in R software.

3.2. Least Squares Estimation

Observing the form of cumulative hazard function (11), we notice that it is a linear function of a and b when assuming other parameters are known. As a result, least squares estimation (LSE) can be used to estimate a and b.

For a dataset size of n, if we estimate the probability of the i-th failure time by its relative frequency, using the non-parametric estimation, the fitted cumulative density function $\widehat{F}\left(t_i\right)$ can be obtained by:

$$\begin{array}{c}\hfill \widehat{F}\left(t_i\right)=\widehat{P}(t\le t_i)=\frac{i-1}{n}.\end{array}$$

(22)

According to the formula $H\left(t\right)=-ln(1-F(t\left)\right)$, the fitted cumulative hazard function $\widehat{H}\left(t_i\right)$ is:

$$\begin{array}{c}\hfill \widehat{H}\left(t_i\right)=ln\left(\frac{n}{n-i+1}\right).\end{array}$$

(23)

Based on the above analysis, when the parameters ${\beta }_1$ and ${\beta }_2$ are known, the least squares estimates of a, b can be obtained by minimizing the least squares distance between $H\left(t\right)$ and $\widehat{H}\left(t\right)$. Denote the least squares distance function as $Q(a,b)$, and it is given by:

$$\begin{array}{cc}\hfill Q(a,b)=& \sum _{i=1}^n(H\left(t_i\right)-\widehat{H}\left(t_i\right))\hfill \\ \hfill =& \sum _{i=1}^{n_1}{[(k_{1}a+k_{2}b)(e^{t_i^{{\beta }_1}}-1)-ln\left(\frac{n}{n-i+1}\right)]}^2\hfill \\ \hfill & +\sum _{i=n_1+1}^{n_1+n_2}{[(k_{1}a+k_{2}b)(e^{{\tau }_1^{{\beta }_1}}-1)+a(t_i-{\tau }_1)+\frac{b}{2}(t_i^2-{\tau }_1^2)-ln\left(\frac{n}{n-i+1}\right)]}^2\hfill \\ \hfill & +\sum _{i=n_1+n_2+1}^n[(k_{1}a+k_{2}b)(e^{{\tau }_1^{{\beta }_1}}-1)+a({\tau }_2-{\tau }_1)+\frac{b}{2}({\tau }_2^2-{\tau }_1^2)+(k_{3}a+k_{4}b)(e^{t^{{\beta }_2}}-e^{{\tau }_2^{{\beta }_2}})\hfill \\ \hfill & -ln\left(\frac{n}{n-i+1}\right){]}^2\hfill \end{array}$$

(24)

where:

$$\begin{array}{c}\hfill k_1=\frac{1}{e^{{\tau }_1^{{\beta }_1}}{\beta }_1{{\tau }_1}^{{\beta }_1-1}},k_2=\frac{1}{e^{{\tau }_1^{{\beta }_1}}{\beta }_1{{\tau }_1}^{{\beta }_1-2}},k_3=\frac{1}{e^{{\tau }_2^{{\beta }_2}}{\beta }_2{{\tau }_2}^{{\beta }_2-1}},k_4=\frac{1}{e^{{\tau }_2^{{\beta }_2}}{\beta }_2{{\tau }_2}^{{\beta }_2-2}}.\end{array}$$

(25)

For the given ${\beta }_1$ and ${\beta }_2$, the analytic expression of least square estimates $\widehat{a}({\beta }_1,{\beta }_2)$ and $\widehat{b}({\beta }_1,{\beta }_2)$ can be obtained by taking the derivative of $Q(a,b)$. The results are as follows:

$$\begin{array}{c}\hfill \widehat{a}({\beta }_1,{\beta }_2)=\frac{B_1{C}_2-B_2{C}_1}{A_1{B}_2-B_1^2}\\ \hfill \widehat{b}({\beta }_1,{\beta }_2)=\frac{B_1{C}_1-A_1{C}_2}{A_1{B}_2-B_1^2}\end{array}$$

(26)

where $A_1,B_1,C_1,B_2$, and $C_2$ are concerned with ${\beta }_1,{\beta }_2,a,b,{\tau }_1,{\tau }_2,t_i$ and are shown specifically in the Appendix A.

Note that if ${\beta }_1$ and ${\beta }_2$ are assumed to be unknown, we can plug $\widehat{a}({\beta }_1,{\beta }_2)$ and $\widehat{b}({\beta }_1,{\beta }_2)$ into the log-likelihood function $l({\beta }_1,{\beta }_2,a,b)$. Thus, the log-likelihood function is only concerned with ${\beta }_1$ and ${\beta }_2$ (denote it as $l({\beta }_1,{\beta }_2)$), which makes it more conducive to calculate the maximum likelihood estimates. By maximizing $l({\beta }_1,{\beta }_2)$, the estimates of ${\beta }_1$ and ${\beta }_2$ can be obtained. Using (10) and (26), the estimates of ${\lambda }_1,{\lambda }_2,a$, and b can be calculated as well.

The least squares estimates (LSEs) of the parameters calculated in this section can also be used as the initial iterative values when calculating the maximum likelihood estimates in the previous section.

4. Interval Estimation

4.1. Observed Fisher Information Matrix

Based on the large-sample theory, when the sample size n is large enough, the inverse of the Fisher information matrix can be used as the approximation of the variance–covariance matrix. Denote the Fisher information matrix as I.

$$I=E\left(\begin{array}{cccc}-\frac{{\partial }^{2}l}{\partial a^2}& -\frac{{\partial }^{2}l}{\partial a\partial b}& -\frac{{\partial }^{2}l}{\partial a\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial a\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial b\partial a}& -\frac{{\partial }^{2}l}{\partial b^2}& -\frac{{\partial }^{2}l}{\partial b\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial b\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial a}& -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial b}& -\frac{{\partial }^{2}l}{\partial {\beta }_1^2}& -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial a}& -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial b}& -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial {\beta }_2^2}\end{array}\right)$$

(27)

The specific elements of I are provided in the Appendix A.

Since it is difficult to calculate the above expectations, the observed Fisher information matrix is often used as a substitute for the Fisher matrix, which does not take expectations but takes the parameter values as the maximum likelihood estimates. Denote it as O.

$$O={\left(\begin{array}{cccc}-\frac{{\partial }^{2}l}{\partial a^2}& -\frac{{\partial }^{2}l}{\partial a\partial b}& -\frac{{\partial }^{2}l}{\partial a\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial a\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial b\partial a}& -\frac{{\partial }^{2}l}{\partial b^2}& -\frac{{\partial }^{2}l}{\partial b\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial b\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial a}& -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial b}& -\frac{{\partial }^{2}l}{\partial {\beta }_1^2}& -\frac{{\partial }^{2}l}{\partial {\beta }_1\partial {\beta }_2}\\ -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial a}& -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial b}& -\frac{{\partial }^{2}l}{\partial {\beta }_2\partial {\beta }_1}& -\frac{{\partial }^{2}l}{\partial {\beta }_2^2}\end{array}\right)}_{|a=\widehat{a},b=\widehat{b},{\beta }_1={\widehat{\beta }}_1,{\beta }_2={\widehat{\beta }}_2}$$

(28)

Therefore, the approximated variance–covariance matrix of $\widehat{a},\widehat{b},\widehat{{\beta }_1}$, and $\widehat{{\beta }_2}$ is given by:

$$\left(\begin{array}{cccc}\widehat{Var}\left(\widehat{a}\right)& \widehat{Cov}(\widehat{a},\widehat{b})& \widehat{Cov}(\widehat{a},{\widehat{\beta }}_1)& \widehat{Cov}(\widehat{a},{\widehat{\beta }}_2)\\ \widehat{Cov}(\widehat{b},\widehat{a})& \widehat{Var}\left(\widehat{b}\right)& \widehat{Cov}(\widehat{b},{\widehat{\beta }}_1)& \widehat{Cov}(\widehat{b},{\widehat{\beta }}_2)\\ \widehat{Cov}({\widehat{\beta }}_1,\widehat{a})& \widehat{Cov}({\widehat{\beta }}_1,\widehat{b})& \widehat{Var}\left({\widehat{\beta }}_1\right)& \widehat{Cov}({\widehat{\beta }}_1,{\widehat{\beta }}_2)\\ \widehat{Cov}({\widehat{\beta }}_2,\widehat{a})& \widehat{Cov}({\widehat{\beta }}_2,\widehat{b})& \widehat{Cov}({\widehat{\beta }}_2,{\widehat{\beta }}_1)& \widehat{Var}\left({\widehat{\beta }}_2\right)\end{array}\right)=O^{-1}.$$

(29)

As the maximum likelihood estimators have asymptotic normality under regularity condition, it can be known that $(\widehat{a},\widehat{b},{\widehat{\beta }}_1,{\widehat{\beta }}_2)$ obeys the quaternion normal distribution approximately. Its mean vector is $(a,b,{\beta }_2,{\beta }_2)$ and the variance–covariance matrix is $O^{-1}$. Based on the above analysis, the asymptotic confidence intervals of $\widehat{a}$, $\widehat{b}$, ${\widehat{\beta }}_1$, and ${\widehat{\beta }}_2$ can also be calculated. In the next section, the specific implementation steps and the calculation method of coverage probabilities are given.

4.2. Asymptotic Confidence Interval

When given a set of initial parameters ${\beta }_1,{\beta }_2,{\lambda }_1$, and ${\lambda }_2$, the following steps can generate sample data and compute the confidence intervals and coverage probabilities.

Step 1: Generate n random numbers that are independent and identically distributed in a Uniform distribution $U(0,1)$. Then, invert $F\left(t\right)$ in (14) to generate the survival time $t_i$. The corresponding function is as follows:

$$t_i=\left\{\begin{array}{c}\hfill {\left[ln(1-\frac{{\beta }_1{{\tau }_1}^{{\beta }_1-1}ln(1-u_i)}{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}})\right]}^{\frac{1}{{\beta }_1}},0<u_i<F\left({\tau }_1\right)\\ \\ \hfill \frac{-a+\sqrt{a^2-2b[ln(1-u_i)+\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{{\tau }_1}^{{\beta }_1}}-1)-a{\tau }_1-\frac{b}{2}{\tau }_1^2]}}{b}F\left({\tau }_1\right)\le u_i<F\left({\tau }_2\right)\\ \\ \hfill {\left[ln(e^{-{\tau }_2^{{\beta }_2}}-\frac{ln(1-u_i)+\frac{(a+b{\tau }_1)e^{-{\tau }_1^{{\beta }_1}}}{{\beta }_1{{\tau }_1}^{{\beta }_1-1}}(e^{{{\tau }_1}^{{\beta }_1}}-1)+a({\tau }_2-{\tau }_1)+\frac{b}{2}({\tau }_2^2-{\tau }_1^2)-ln(1-u_i)}{(a+b{\tau }_2)e^{-{\tau }_2^{{\beta }_2}}}{\beta }_2{{\tau }_2}^{{\beta }_2-1})\right]}^{\frac{1}{{\beta }_2}}\\ \hfill u_i\ge F\left({\tau }_2\right)\\ \end{array}\right.$$

(30)

Step 2: Use the data $t_i$ generated from $Step1$ and the log-likelihood function (17) to calculate the MLEs of a, b, ${\beta }_1$, and ${\beta }_2$. Denote them as $\widehat{a}$, $\widehat{b}$, $\widehat{{\beta }_1}$, and $\widehat{{\beta }_2}$. Calculate the MLEs of ${\lambda }_1$ and ${\lambda }_2$ via the equation (10) and denote them as ${\widehat{\lambda }}_1$ and ${\widehat{\lambda }}_2$.

Step 3: Use the data $t_i$ generated from $Step1$ and the MLEs from $Step2$ to calculate the observed Fisher information matrix O.

Step 4: Invert O matrix to calculate the asymptotic variance–covariance matrix, and denote it as A. Obtain the asymptotic variance of ${\beta }_1$ and ${\beta }_2$ as $\widehat{var}$(${\widehat{\beta }}_1$) and $\widehat{var}$(${\widehat{\beta }}_2$).

Step 5: Based on the theory of the Delta method ([26]), the asymptotic variance of ${\widehat{\lambda }}_1$ and ${\widehat{\lambda }}_2$ can be calculated using the following equations:

$$\begin{array}{c}\hfill \widehat{var}\left({\widehat{\lambda }}_1\right)=C_{1}A{C}_1^T\\ \hfill \widehat{var}\left({\widehat{\lambda }}_2\right)=C_{2}A{C}_2^T\end{array}$$

(31)

where:

$$\begin{array}{c}\hfill C_1=(\frac{\partial {\lambda }_1(a,b,{\beta }_1)}{\partial a},\frac{\partial {\lambda }_1(a,b,{\beta }_1)}{\partial b},\frac{\partial {\lambda }_1(a,b,{\beta }_1)}{\partial {\beta }_1},0)\\ \hfill C_2=(\frac{\partial {\lambda }_2(a,b,{\beta }_2)}{\partial a},\frac{\partial {\lambda }_2(a,b,{\beta }_2)}{\partial b},0,\frac{\partial {\lambda }_2(a,b,{\beta }_2)}{\partial {\beta }_2})\end{array}$$

(32)

$C_1^T$ and $C_2^T$ are the transpose matrices of $C_1$ and $C_2$, respectively. Further, the specific expressions of $C_1$ and $C_2$ are shown in the Appendix A.

$Step6$: The lower and upper bounds of the $100(1-\alpha )$% confidence intervals for ${\beta }_1,{\beta }_2,{\lambda }_1,{\lambda }_2$ are given by:

$$\begin{array}{c}\hfill {{\widehat{\beta }}_i}^L=min\{{\widehat{\beta }}_i-u_{\frac{\alpha }{2}}\sqrt{\widehat{var}\left({\widehat{\beta }}_i\right)},0\}{{\widehat{\beta }}_i}^U={\widehat{\beta }}_i+u_{\frac{\alpha }{2}}\sqrt{\widehat{var}\left({\widehat{\beta }}_i\right)},i=1,2\\ \hfill {{\widehat{\lambda }}_i}^L=min\{{\widehat{\lambda }}_i-u_{\frac{\alpha }{2}}\sqrt{\widehat{var}\left({\widehat{\lambda }}_i\right)},0\}{{\widehat{\lambda }}_i}^U={\widehat{\lambda }}_i+u_{\frac{\alpha }{2}}\sqrt{\widehat{var}\left({\widehat{\lambda }}_i\right)},i=1,2\end{array}$$

(33)

where $u_q$ is the q-quantile of a standardized normal distribution.

Step 7: Repeat the foregoing steps 999 times to obtain the coverage probabilities as CPrs.

$$\begin{array}{c}\hfill CPr\left({\beta }_1\right)=\sum _{j=1}^{999}\frac{I({{\widehat{\beta }}_{1j}}^L<{\beta }_1<{{\widehat{\beta }}_{1j}}^U)}{999}\end{array}$$

(34)

where $I({{\widehat{\beta }}_{1j}}^L<{\beta }_1<{{\widehat{\beta }}_{1j}}^U)$ is the indicator function. ${{\widehat{\beta }}_{1j}}^L$ and ${{\widehat{\beta }}_{1j}}^U$ are the j-th results of the ${\beta }_1$’ lower and upper bounds of the $100(1-\alpha )$% confidence intervals.

In the same way, we can obtain the CPrs of ${\beta }_2,{\lambda }_1$, and ${\lambda }_2$.

5. Simulation Results and Analysis

In this section, the simulation results under different sample sizes (n) and different parameters are presented using the method given in the previous section with the R program.

The simulation results of the MLEs, LSEs, 95% and 99% confidence intervals, and the corresponding coverage probabilities are given by Monte Carlo simulations, which evaluate the performance of the estimation methods. By comparing the mean, bias, and mean square error of MLEs and LSEs, the advantages and disadvantages of the two methods are compared.

Based on the characteristics of the Chen distribution’s hazard function, different parameters are chosen to generate random numbers, and the results are listed in

Table 1. The results of MLEs, LSEs, interval estimates, and CPrs when ${\beta }_1=0.7,{\beta }_2=0.9$, ${\lambda }_1=0.5016$, ${\lambda }_2=1.0015,{\tau }_1=0.5,{\tau }_2=1,a=-0.85,b=3.3$.

n	$\mathbf{Par}$	$\mathbf{Method}$	$\mathbf{Mean}$	$\mathbf{Bias}$	$\mathbf{MSE}$	${\mathbf{LB}}_{95\%}$	${\mathbf{UB}}_{95\%}$	${\mathbf{CPr}}_{95\%}$	${\mathbf{LB}}_{99\%}$	${\mathbf{UB}}_{99\%}$	${\mathbf{CPr}}_{99\%}$
50	${\lambda }_1$	MLE	0.5469	0.0454	0.0247	0.2580	0.8359	0.944	0.1666	0.9273	0.984
		LSE	0.5165	0.0149	0.0173
	${\lambda }_2$	MLE	1.0326	0.0311	0.1454	0	2.5273	0.964	0	3.0001	0.981
		LSE	0.6894	−0.3121	0.0137
	${\beta }_1$	MLE	0.7597	0.0597	0.0473	0.3832	1.1362	0.947	0.2641	1.2553	0.988
		LSE	0.7294	0.0294	0.0191
	${\beta }_2$	MLE	1.1487	0.2487	0.2595	0.2605	2.0368	0.901	-0.0204	2.3178	0.964
		LSE	1.1840	0.2840	0.0512
100	${\lambda }_1$	MLE	0.5181	0.0165	0.0102	0.3220	0.7142	0.948	0.2599	0.7762	0.993
		LSE	0.5286	0.0270	0.0089
	${\lambda }_2$	MLE	1.1001	0.0986	1.6755	0.1456	2.0545	0.953	0	2.3564	0.984
		LSE	0.7872	−0.2143	0.0114
	${\beta }_1$	MLE	0.7238	0.0238	0.0178	0.4708	0.9768	0.951	0.3907	1.0568	0.988
		LSE	0.7327	0.0327	0.0095
	${\beta }_2$	MLE	1.0009	0.1009	0.1027	0.4122	1.5896	0.917	0.2259	1.7758	0.970
		LSE	1.0530	0.1530	0.0117
200	${\lambda }_1$	MLE	0.5108	0.0092	0.0047	0.3736	0.6480	0.960	0.3302	0.6914	0.992
		LSE	0.5171	0.0155	0.0049
	${\lambda }_2$	MLE	0.9841	−0.0173	0.1017	0.3514	1.6168	0.969	0.1512	1.8170	0.995
		LSE	0.8632	−0.1382	0.0072
	${\beta }_1$	MLE	0.7143	0.0143	0.0089	0.5382	0.8904	0.953	0.4825	0.9462	0.988
		LSE	0.7210	0.0210	0.0041
	${\beta }_2$	MLE	0.9707	0.0707	0.0439	0.5673	1.3740	0.930	0.4397	1.5016	0.977
		LSE	0.9881	0.0881	0.0046

,

Table 2. The results of MLEs, LSEs, interval estimates, and CPrs when ${\beta }_1=0.7,{\beta }_2=0.9$, ${\lambda }_1=0.5016$, ${\lambda }_2=1.0015,{\tau }_1=0.5,{\tau }_2=0.7,a=-2.0235,b=5.6470$.

n	$\mathbf{Par}$	$\mathbf{Method}$	$\mathbf{Mean}$	$\mathbf{Bias}$	$\mathbf{MSE}$	${\mathbf{LB}}_{95\%}$	${\mathbf{UB}}_{95\%}$	${\mathbf{CPr}}_{95\%}$	${\mathbf{LB}}_{99\%}$	${\mathbf{UB}}_{99\%}$	${\mathbf{CPr}}_{99\%}$
50	${\lambda }_1$	MLE	0.5406	0.0390	0.0278	0.2367	0.8444	0.947	0.1406	0.9406	0.981
		LSE	0.5392	0.0376	0.0221
	${\lambda }_2$	MLE	1.0617	0.0602	0.3453	0.2177	1.9056	0.955	0	2.1726	0.988
		LSE	0.8517	−0.1498	0.0387
	${\beta }_1$	MLE	0.7529	0.0529	0.0489	0.3675	1.1382	0.944	0.2456	1.2602	0.979
		LSE	0.7409	0.0409	0.0194
	${\beta }_2$	MLE	1.0346	0.1346	0.1296	0.3951	1.6742	0.930	0.1928	1.8765	0.983
		LSE	1.0550	0.1550	0.0041
100	${\lambda }_1$	MLE	0.5205	0.0189	0.0130	0.3122	0.7288	0.931	0.2464	0.7947	0.981
		LSE	0.5247	0.0231	0.0124
	${\lambda }_2$	MLE	1.0101	0.0086	0.1397	0.4483	1.5719	0.953	0.2706	1.7496	0.984
		LSE	0.8970	−0.1044	0.0212
	${\beta }_1$	MLE	0.7259	0.0259	0.0189	0.4647	0.9872	0.952	0.3820	1.0698	0.990
		LSE	0.7302	0.0302	0.0107
	${\beta }_2$	MLE	0.9733	0.0733	0.0555	0.5470	1.3996	0.927	0.4122	1.5345	0.978
		LSE	0.9928	0.0928	0.0012
200	${\lambda }_1$	MLE	0.5099	0.0083	0.0058	0.3649	0.6549	0.949	0.3191	0.7007	0.987
		LSE	0.5235	0.0219	0.0063
	${\lambda }_2$	MLE	1.0059	0.0044	0.0406	0.6175	1.3943	0.953	0.4946	1.5172	0.992
		LSE	0.9280	−0.0735	0.0124
	${\beta }_1$	MLE	0.7135	0.0135	0.0086	0.5323	0.8947	0.949	0.4750	0.9520	0.988
		LSE	0.7201	0.0201	0.0047
	${\beta }_2$	MLE	0.9285	0.0285	0.0225	0.6366	1.2204	0.948	0.5442	1.3128	0.992
		LSE	0.9569	0.0569	0.0041

,

Table 3. The results of MLEs, LSEs, interval estimates, and CPrs when ${\beta }_1=1,{\beta }_2=1.2$, ${\lambda }_1=0.7642$, ${\lambda }_2=1.1061,{\tau }_1=0.4,{\tau }_2=0.6,a=-0.7,b=4.6$.

n	$\mathbf{Par}$	$\mathbf{Method}$	$\mathbf{Mean}$	$\mathbf{Bias}$	$\mathbf{MSE}$	${\mathbf{LB}}_{95\%}$	${\mathbf{UB}}_{95\%}$	${\mathbf{CPr}}_{95\%}$	${\mathbf{LB}}_{99\%}$	${\mathbf{UB}}_{99\%}$	${\mathbf{CPr}}_{99\%}$
50	${\lambda }_1$	MLE	0.9009	0.1367	0.2990	0.2552	1.5466	0.940	0.0509	1.7509	0.987
		LSE	0.8342	0.0700	0.0793
	${\lambda }_2$	MLE	1.1245	0.0185	0.1404	0.4643	1.7848	0.938	0.2554	1.9937	0.978
		LSE	1.0138	−0.0923	0.0563
	${\beta }_1$	MLE	1.0870	0.0870	0.1224	0.5168	1.6572	0.944	0.3364	1.8376	0.982
		LSE	1.0527	0.0527	0.0511
	${\beta }_2$	MLE	1.3984	0.1984	0.2188	0.5822	2.2146	0.940	0.324	2.4727	0.989
		LSE	1.3959	0.1959	0.0055
100	${\lambda }_1$	MLE	0.8139	0.0497	0.0504	0.3951	1.2326	0.950	0.2627	1.3651	0.988
		LSE	0.8284	0.0643	0.0451
	${\lambda }_2$	MLE	1.1094	0.0033	0.0576	0.6578	1.5613	0.942	0.5152	1.7038	0.984
		LSE	1.0331	−0.0730	0.0343
	${\beta }_1$	MLE	1.0354	0.0354	0.0390	0.6529	1.4179	0.955	0.5319	1.5389	0.988
		LSE	1.0504	0.0504	0.0249
	${\beta }_2$	MLE	1.2926	0.0926	0.0779	0.7531	1.8323	0.951	0.5823	2.003	0.989
		LSE	1.3152	0.1152	0.0015
200	${\lambda }_1$	MLE	0.7859	0.0217	0.0208	0.4976	1.0741	0.957	0.4064	1.1653	0.994
		LSE	0.8156	0.0514	0.0219
	${\lambda }_2$	MLE	1.1142	0.0082	0.0290	0.7977	1.4307	0.947	0.6976	1.5309	0.987
		LSE	1.0520	−0.0541	0.0170
	${\beta }_1$	MLE	1.0168	0.0168	0.0184	0.7512	1.2823	0.963	0.6672	1.3664	0.993
		LSE	1.0339	0.0339	0.0122
	${\beta }_2$	MLE	1.2398	0.0398	0.0382	0.8716	1.6081	0.949	0.7551	1.7246	0.992
		LSE	1.2693	0.0693	0.0005

and

Table 4. The results of MLEs, LSEs, interval estimates, and CPrs when ${\beta }_1=0.8,{\beta }_2=1.2$, ${\lambda }_1=0.3679$, ${\lambda }_2=0.0802,{\tau }_1=1,{\tau }_2=2,a=0.5,b=0.3$.

n	$\mathbf{Par}$	$\mathbf{Method}$	$\mathbf{Mean}$	$\mathbf{Bias}$	$\mathbf{MSE}$	${\mathbf{LB}}_{95\%}$	${\mathbf{UB}}_{95\%}$	${\mathbf{CPr}}_{95\%}$	${\mathbf{LB}}_{99\%}$	${\mathbf{UB}}_{99\%}$	${\mathbf{CPr}}_{99\%}$
50	${\lambda }_1$	MLE	0.3705	0.0027	0.0048	0.2254	0.5157	0.977	0.1795	0.5616	0.996
		LSE	0.3761	0.0082	0.0061
	${\lambda }_2$	MLE	0.0295	−0.0508	0.5604	0	0.2688	0.981	0	0.3445	0.994
		LSE	0.0426	−0.0376	0.0001
	${\beta }_1$	MLE	0.8371	0.0371	0.0278	0.5147	1.1595	0.954	0.4127	1.2615	0.990
		LSE	0.8159	0.0159	0.0105
	${\beta }_2$	MLE	1.3156	0.1156	0.0826	0.7707	1.8606	0.871	0.5984	2.0329	0.936
		LSE	1.3586	0.1586	0.0133
100	${\lambda }_1$	MLE	0.3715	0.0036	0.0024	0.2717	0.4714	0.972	0.2401	0.5029	0.991
		LSE	0.3737	0.0058	0.0031
	${\lambda }_2$	MLE	0.0908	0.0106	0.0077	0	0.2173	0.969	0	0.2573	0.993
		LSE	0.0552	−0.0250	0.0001
	${\beta }_1$	MLE	0.8204	0.0204	0.0141	0.5984	1.0424	0.943	0.5282	1.1126	0.995
		LSE	0.8161	0.0161	0.0052
	${\beta }_2$	MLE	1.2495	0.0495	0.0363	0.8793	1.6197	0.923	0.7622	1.7368	0.969
		LSE	1.2895	0.0895	0.0035
200	${\lambda }_1$	MLE	0.3711	0.0032	0.0012	0.3016	0.4405	0.965	0.2796	0.4625	0.994
		LSE	0.3709	0.0030	0.0015
	${\lambda }_2$	MLE	0.0791	−0.0011	0.0019	0	0.1625	0.965	0	0.1889	0.997
		LSE	0.0635	−0.0167	0.0001
	${\beta }_1$	MLE	0.8058	0.0058	0.0060	0.6519	0.9596	0.956	0.6033	1.0083	0.993
		LSE	0.8101	0.0101	0.0026
	${\beta }_2$	MLE	1.2360	0.0360	0.0161	0.9859	1.4861	0.915	0.9068	1.5652	0.971
		LSE	1.2551	0.0551	0.0014

. The results include the mean, bias, mean square error (MSE), lower bounds ($L{B}_{95\%}$), upper bounds ($U{B}_{95\%}$), and coverage probabilities ($CP{r}_{95\%}$) of 95% confidence interval and lower bounds ($L{B}_{99\%}$), upper bounds ($U{B}_{99\%}$), and coverage probabilities ($CP{r}_{99\%}$) of 99% confidence intervals.

Table 1 shows the simulation results when the hazard functions under the two stresses are both bathtub-shaped with $n=50,n=100,n=200$.

Table 2 shows the simulation results when the distributions under the two stress levels are the same as those in Table 1, but the lag time (${\tau }_2-{\tau }_1$) is shortened.

Table 3 shows the simulation results when the hazard functions under the two different stresses both increase monotonically with $n=50,n=100,n=200$.

Table 4 shows the simulation results when the hazard function under the first stress level is bathtub-shaped and in the second stress level is monotonically increasing with $n=50,n=100,n=200$.

Based on Table 1, Table 2, Table 3 and Table 4, some conclusions are summarized as follows.

(1)

No matter which values the parameters take, the estimated values are close to the real values, and mostly the bias and mean square errors decrease with the increase in sample size, which shows that the two estimations are effective.

(2)

From the perspective of bias, the results of LSE are generally better than MLE when $n=50$; the results of MLE are generally better than LSE when $n=100$ and $n=200$. This means that LSE is preferred when the sample size is small, while MLE is preferred when the sample size is large.

(3)

Under different parameters, the mean square errors of LSEs are generally less than that of MLEs, and the advantage of LSE in the mean square errors is more obvious when the sample size n is small.

(4)

In terms of the asymptotic confidence intervals, generally, the coverage probabilities of the 95% are close to 95%, and the coverage probabilities of the 99% are close to 99%, which verifies the correctness of the methods. The coverage probabilities are closer to 1-$\alpha$ with the increase in the sample size, which means the asymptotic confidence intervals will be more precise when the sample size is larger.

(5)

In general, the estimations perform better when the hazard function under the first stress is bathtub-shaped and under the second stress is monotonically increasing. The coverage probabilities fit better when the risk function is monotonically increasing under both stress levels.

(6)

Comparing Table 1 and Table 2, it can be seen that when the lagged-effect time is shortened, the mean square errors of MLEs and LSEs both increase under the small sample size.

6. A Special Case

Since the parameter $\beta$ determines whether the shape of the hazard function is a bathtub shape or not and, in many cases, the stress does not change the shape of the hazard function, a special case will be discussed below.

When assuming that parameter ${\beta }_1$ is equal to parameter ${\beta }_2$ and denoting them as $\beta$, the model becomes the following form.

The hazard functions under the two stresses are:

$$h_1\left(t\right)=\begin{array}{c}\hfill {\lambda }_1\beta t^{\beta -1}{e}^{t^{\beta }},t>0,\end{array}$$

(35)

$$h_2\left(t\right)=\begin{array}{c}\hfill {\lambda }_2\beta t^{\beta -1}{e}^{t^{\beta }},t>0.\end{array}$$

(36)

Under the CRM, the hazard function is obtained by:

$$h\left(t\right)=\left\{\begin{array}{c}\hfill {\lambda }_1\beta t^{\beta -1}{e}^{t^{\beta }},0<t<{\tau }_1\\ \hfill a+bt,{\tau }_1\le t<{\tau }_2\\ \hfill {\lambda }_2\beta t^{\beta -1}{e}^{t^{\beta }},t\ge {\tau }_2\end{array}\right.$$

(37)

Replace parameters ${\lambda }_1$ and ${\lambda }_2$ with a and b, and the cumulative hazard function $H\left(t\right)$ under the model is:

$$H\left(t\right)=\left\{\begin{array}{c}\hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{t^{\beta }}-1),0<t<{\tau }_1\\ \\ \hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)+a(t-{\tau }_1)+\frac{b}{2}(t^2-{\tau }_1^2),{\tau }_1\le t<{\tau }_2\\ \\ \hfill \frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)+a({\tau }_2-\tau )+\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\\ \hfill +\frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{\beta {{\tau }_2}^{\beta -1}}(e^{t^{\beta }}-e^{{\tau }_2^{\beta }}),t\ge {\tau }_2\end{array}\right.$$

(38)

The survival function $S\left(t\right)$ under the model is:

$$S\left(t\right)=\left\{\begin{array}{c}\hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{t^{\beta }}-1)\right\},0<t<{\tau }_1\\ \\ \hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\\ \\ \hfill exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{{\beta }_1-1}}(e^{{\tau }_1^{\beta }}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\\ \hfill \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{\beta {{\tau }_2}^{\beta -1}}(e^{t^{\beta }}-e^{{\tau }_2^{\beta }})\right\},t\ge {\tau }_2\end{array}\right.$$

(39)

The probability density function $f\left(t\right)$ of the lifetime is as follows:

$$f\left(t\right)=\left\{\begin{array}{c}\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{{{\tau }_1}^{\beta -1}}{t}^{\beta -1}{e}^{t^{\beta }}exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{t^{\beta }}-1)\right\},0<t<{\tau }_1\hfill \\ \\ (a+bt)exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\hfill \\ \\ \frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{{{\tau }_2}^{\beta -1}}{t}^{\beta -1}{e}^{t^{\beta }}exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{{\beta }_1{{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\hfill \\ \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{\beta {{\tau }_2}^{\beta -1}}(e^{t^{\beta }}-e^{{\tau }_2^{\beta }})\right\},t\ge {\tau }_2\hfill \end{array}\right.$$

(40)

The corresponding cumulative distribution function $F\left(t\right)$ is given by:

$$F\left(t\right)=\left\{\begin{array}{c}\hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{t^{\beta }}-1)\right\},0<t<{\tau }_1\\ \\ \hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)-a(t-{\tau }_1)-\frac{b}{2}(t^2-{\tau }_1^2)\right\},{\tau }_1\le t<{\tau }_2\\ \\ \hfill 1-exp\left\{-\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}(e^{{\tau }_1^{\beta }}-1)-a({\tau }_2-{\tau }_1)-\frac{b}{2}({\tau }_2^2-{\tau }_1^2)\right\}\\ \hfill \times exp\left\{-\frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{\beta {{\tau }_2}^{\beta -1}}(e^{t^{\beta }}-e^{{\tau }_2^{\beta }})\right\},t\ge {\tau }_2\end{array}\right.$$

(41)

Accordingly, the log-likelihood function $l(\beta ,a,b)$ can be written as:

$$\begin{array}{cc}\hfill l(\beta ,a,b)=& n_1\left[ln(a+b{\tau }_1)-{{\tau }_1}^{\beta }\right]+\sum _{i=n_1+1}^{n_1+n_2}ln(a+b{t}_i)+n_3\left[ln(a+b{\tau }_2)-{{\tau }_2}^{\beta }\right]\hfill \\ \hfill & +(\beta -1)[\sum _{i=1}^{n_1}ln\left(\frac{t_i}{{\tau }_1}\right)+\sum _{i=n_1+n_2+1}^{n}ln\left(\frac{t_i}{{\tau }_1}\right)]+\sum _{i=1}^{n_1}{t_i}^{\beta }+\sum _{i=n_1+n_2+1}^n{t_i}^{\beta }\hfill \\ \hfill & -a[\sum _{i=n_1+1}^{n_1+n_2}(t_i-{\tau }_1)+({\tau }_2-{\tau }_1)n_3]-\frac{b}{2}[\sum _{i=n_1+1}^{n_1+n_2}(t_i^2-{\tau }_1^2)+({\tau }_2^2-{\tau }_1^2)n_3]\hfill \\ \hfill & -\frac{(a+b{\tau }_1)e^{-{\tau }_1^{\beta }}}{\beta {{\tau }_1}^{\beta -1}}[\sum _{i=1}^{n_1}(e^{t_i^{\beta }}-1)+(n_2+n_3)+(e^{{{\tau }_1}^{\beta }}-1)]\hfill \\ \hfill & -\frac{(a+b{\tau }_2)e^{-{\tau }_2^{\beta }}}{\beta {{\tau }_2}^{\beta -1}}\sum _{i=n_1+n_2+1}^n(e^{t_i^{\beta }}-e^{{\tau }_2^{\beta }})\hfill \end{array}$$

(42)

Other relevant parameter estimations can also be obtained. The corresponding methods are similar to those of previous sections and the specific steps are omitted.

7. Conclusions

In this paper, the parameter estimations and the statistical inference of the Chen distribution under the step-stress model with lagged effect are studied. Maximum likelihood estimation is used for point estimation, and the Newton–Raphson algorithm is used when solving the nonlinear equations. Based on the unique linear form of risk function under CRM, another point estimation is obtained based on the large sample theory and the least squares estimation method. Different from maximum likelihood estimation, it gives the specific expressions of a, b for the given ${\beta }_1,{\beta }_2$. Moreover, using the observed Fisher matrix and the asymptotic normality of the maximum likelihood estimators, a method to construct the asymptotic confidence interval and coverage probabilities is provided. The performance of those estimation methods is evaluated by Monte Carlo simulation. It can be seen from the simulation results that the accuracy of the two point estimations is different when parameters or sample sizes change, which may be due to distinct forms of the Chen distribution’s risk functions.

The bathtub-shaped hazard function of the Chen distribution is of great significance in real life. The step-stress model is practical in survival analysis and the lagged effect makes it more consistent with reality. This paper can also be further extended by considering competing risks or a censoring scheme.