Statistical Inference and Optimal Design of Accelerated Life Testing for the Chen Distribution under Progressive Type-II Censoring

: This paper discusses statistical inference and optimal design of constant-stress accelerated life testing for the Chen distribution under progressive Type-II censoring. The scale parameter of the life distribution is assumed to be a logarithmic linear function of the stress level. The maximum likelihood estimates of the parameters are obtained. Then, the observed Fisher information matrix is derived and utilized to construct asymptotic conﬁdence intervals. Meanwhile, the parametric bootstrap methods are provided for the interval estimation. In addition, the Bayes estimates under the squared error loss function are obtained by applying the Tierney and Kadane technique and Lindley’s approximation. As for the optimal design, D- and A-optimality criteria are considered to determine the optimal transformed stress level. Finally, the simulation is carried out to demonstrate the proposed estimation techniques and the optimal criteria, and a real data set is discussed.


Introduction
In factory production, life testing is of great significance to evaluate product quality. Conventional life testing is to mainly observe the failure time information of all test units under normal usage in order to reveal the life characteristics. However, products are produced more and more reliably currently, so it is difficult or even impossible for them to fail under normal usage. To reduce the cost and time, it is necessary to make the test units fail more quickly. A common way is to perform life testing at higher stress levels than normal, which is called accelerated life testing (ALT). Then, the failure time information obtained is analyzed to extrapolate the life characteristics under the normal stress level. One type of classifying ALT is the way stress changes with time. ALTs are mainly classified into constant-stress ALT (CSALT), step-stress ALT (SSALT) and progressive-stress ALT (PSALT). The different focus of the three types is the relationship between stress loading and testing time. For CSALT, the stress is fixed with time; for SSALT, the stress is increased step by step at specific moments; for PSALT, stress is increased continuously with time. Another type of classifying ALTs is by the number of stress levels. They are mainly divided into simple ALT and multiple ALT. The simple ALT has only two stress levels while the multiple ALT has more than two stress levels.
Some of the earlier works on ALT are mainly based on McCool [1], Nelson [2] and Miller and Nelson [3]. Gradually, the study of ALT has attracted the interest of more and more researchers. For the studies of CSALT, Ismail [4] discussed the estimation of parameters of Weibull distribution in the constant-stress partial ALT under hybrid censoring. Statistical inferences on the acceleration factor and the distribution parameters were made utilizing maximum likelihood and percentile bootstrap. Mohie EI-Din et al. [5] derived the optimal plans for k-level CSALT for Lindley failure data and compared the Considering that ALT is relatively rarely used for the Chen distribution which has a bathtub-shaped hazard function and that optimal design studies are less concerned with solving for optimal transformed stress levels, we apply the CSALT model to test units whose lifetimes follow the Chen distribution in the presence of progressive Type-II censoring in this paper. To better estimate the lifetime characteristics of products, we discuss a variety of parameter and interval estimation methods. The maximum likelihood estimates of the parameters for the Chen distribution are obtained under progressive censoring. The Bayes estimates of the parameters are derived by both the Tierney and Kadane technique and Lindley's approximation. We also discuss three types of confidence intervals which are the asymptotic confidence interval, the bootstrap percentile confidence interval, and the bootstrap-t confidence interval. In addition, to better estimate the life parameters of products under the normal stress level, two optimal criteria are utilized to determine and specify the optimal transformed stress level of CSALT. To measure the performance of the estimates, we conduct some simulations.
The rest of this paper is organized as follows. In Section 2, we discussed some properties of the Chen distribution. In Section 3, we propose the progressive censoring scheme and necessary model assumptions. Then, the maximum likelihood estimates are obtained in Section 4. In Section 5, we derive the Fisher information matrix and the variance-covariance matrix. The asymptotic confidence intervals are also constructed. In Section 6, we present the parametric bootstrap confidence intervals. In Section 7, we discuss the Bayesian estimation and utilize the Tierney and Kadane technique and Lindley's approximation to obtain the approximate estimates. In Section 8, we apply two optimal criteria to determine the optimal transformed stress levels. In Section 9, a numerical study is conducted to illustrate the methods proposed and examine their performance. In addition, a real data set is discussed in Section 10. Finally, we make conclusions.

The Chen Distribution
The probability density function (pdf), the cumulative distribution function (cdf), and the hazard function (hf) for the Chen distribution are expressed respectively as follows: where α is the scale parameter and β is the shape parameter. When β is less than one, the hazard function has a bathtub shape, while the shape is increasing when β ≥ 1. Different plots of pdf and hf are shown in Figures 1-6. The shape of pdf or hf varies for different parameter values. Assume that β is fixed. If 0 < β < 1: when 0 < α < 1, pdf is monotonically decreasing or decreasing, then increasing and finally decreasing, and the tail is thinner with a larger value of α; when α ≥ 1, pdf is monotonically decreasing, and the decrease is sharper with a larger value of α. If β = 1: when 0 < α < 1, the pdf is unimodal; when α ≥ 1, the pdf is monotonically decreasing and the decrease is sharper with a larger value of α. If β > 1: the pdf is unimodal no matter what α value is taken, and the peak value is larger with a larger value of α. From Figures 5 and 6, it can be seen that the hazard rate is larger with a larger value of α if 0 < β < 1; the rising rate of hf is higher with a larger value of α if β ≥ 1.  When β is fixed at some value less than 1, the relationship between the mean of the Chen distribution and α is shown in Figure 7. The mean decreases when the value of α increases. Additionally, letting the cdf of the Chen distribution equal 1 2 , we can solve for the median x m = ln 1 + 1 α ln 2 1 β . Similarly, when β is fixed and α increases, the median decreases. In addition, the Chen distribution has a great deal of flexibility. Let X be a nonnegative random variable and follow the Chen distribution. Following Ahmed et al. [21], Y = e X β − 1 follows the exponential distribution and Y = e X β − 1 1 θ follows the Weibull distribution.
Furthermore, the Chen distribution reduces to Gompertz distribution when β = 1.

Model Assumptions
For CSALT, assume a total of n units under the test and k ≥ 2 stress levels, s 1 < · · · < s k . At each stress level s i , i = 1, 2, . . . , k, n i identical units are assigned for the test, such that ∑ k i=1 n i = n. The lifetime of each unit X ij , j = 1, 2, . . . , n i , follows the Chen distribution. s 0 is the normal stress level. The progressive Type-II censoring scheme in CSALT is briefly described as follows: For i = 1, 2, . . . , k, when the first failure occurs at time t i1:m i :n i , R i1 units are randomly withdrawn from the remaining n i − 1 units whose lifetime X ij > t i1:m i :n i . Then, when the second failure occurs at time t i2:m i :n i , R i2 units are randomly withdrawn from the remaining n i − 2 − R i1 units whose lifetime X ij > t i2:m i :n i . Similarly, when the m ith (1 ≤ m i ≤ n i ) failure occurs at a random time of t im i :m i :n i , all remaining R im i = n i − m i − ∑ m i −1 j=1 R ij units are withdrawn and the test is terminated. The censoring scheme is shown in Figure 8. Under stress level s i , i = 1, 2, . . . k, R i1 , R i2 , · · · , R im i is the progressive censoring scheme and t ij:m i :n i is the observed progressive censored data that t i1:m i :n i < · · · < t im i :m i :n i . Obviously, the sample is complete when R i1 , R i2 , · · · , R im i = (0, 0, · · · , 0) and the sample is Type-II censored when R i1 , R i2 , · · · , R im i = (0, · · · , 0, n i − m i ). t km k :m k :n k R km k t Figure 8. For n i assigned units, progressive Type-II censoring scheme R i1 , R i2 , · · · , R im i at stress level s i , i = 1, 2, . . . k, in constant-stress accelerated life testing.
The following model assumptions are used throughout the whole paper: (1) Under stress level s i , i = 1, 2, . . . , k, the lifetimes of the units follow (2) The relationship between the stress level s i and the lifetime distribution parameter α i is where ϕ i = ϕ(s i ) is an incremental function of s i . b 0 and b > 0 are unknown parameters. To simplify the expression, according to the model (5), the scale parameter α i can be presented as where α 0 is the scale parameter for the Chen distribution under normal stress level s 0 ,

Maximum Likelihood Estimation
Based on the above assumptions and notations, when the stress level is s i , the likelihood function of parameters α 0 , λ and β is given by where j=1 R ij and t ij = t ij:m i :n i . Then, the likelihood function is given by Correspondingly, the log-likelihood function ignoring the constants is given by The partial derivatives of l(α 0 , λ, β) with respect to α 0 , λ and β are given by ∂l ∂λ ∂l ∂β We can solve the maximum likelihood estimates (MLE) α 0 ,λ,β for α 0 , λ and β by making the functions (10)-(12) equal to 0 simultaneously. Unfortunately, explicit solutions cannot be obtained and numerical methods such as the Newton-Raphson method can be utilized to calculate the estimates. According to MacDonald [22], for a constrained optimization problem, one can enforce the constraint by means of a transformation or use of a constrained optimizer. Considering the three parameters are constrained (α > 0, λ > 1, β > 0), we use the optim function in R software with method L-BFGS-B and default settings. L-BFGS-B is an application of the quasi-Newton method.

Fisher Information Matrix
Let Θ = (α 0 , λ, β) and T be the observed failure data. For a large-size sample, the inverse of the Fisher information matrix (FIM) can be regarded as the approximation of the variance-covariance matrix. First, FIM is given as follows by evaluating the negative of the expected values of the second-order derivatives and mixed partial derivatives of log-likelihood function in (9) with respect to α 0 , λ and β, Due to the complexity of the pdf and the integral formula, in practice, the observed FIM which utilizes the MLEs of parametersΘ = α 0 ,λ,β is used to approximate the exact values, where According to the observed FIM in (14), the asymptotic variance-covariance matrix of Θ can be constructed as O −1 T (Θ). Then, on the basis of the respective asymptotic variances and the asymptotic normality ofΘ, the 100(1 − τ)% asymptotic confidence intervals can be given by where θ is α 0 , λ or β, and Z r is the 100r-th percentile of the standard normal distribution.

Parametric Bootstrap Intervals
The asymptotic confidence interval based on the normal approximation which we mentioned in Section 5 has a better performance for the large sample size. In this section, two parametric bootstrap methods are discussed to find more accurate approximate confidence intervals for unknown parameters α 0 , λ and β when the sample size is not large. The algorithms are as follows.
By following the above steps, the 100(1 − τ)% bootstrap percentile confidence interval of θ is given by

The Bootstrap-t Confidence Interval
-Steps 1 and 2 are the same as those in the algorithm of the bootstrap percentile confidence interval. - Step 3: Compute the MLEθ * using the progressive censored samples from Step 2 where θ is α 0 , λ or β . Then, obtain the statistic Step 4: Repeat Step 2 and 3, B times.

Bayesian Estimation
The aim of this section is to obtain the Bayes estimates (BE) of unknown parameters α 0 , λ and β, utilizing progressive censored samples of CSALT. First, considering α 0 > 0, λ > 1 and β > 0, the prior distributions of α 0 and β are assumed to be the gamma distributions and the prior distribution of λ is assumed to be the left truncated gamma distribution as where all the hyperparameters a i and b i , (i = 1, 2, 3) are considered to be known and non-negative. When the hyperparameters all tend to zero, the BEs are obtained in the case of non-informative prior. Because of the independence of α 0 , λ and β, their joint prior distribution is given by Denote the observed failure data as T. Then using both the likelihood function (8) and the joint prior distribution (25), the posterior distribution is given by Let g(α 0 , λ, β) be an arbitrary function of α 0 , λ and β. The BE of g(α 0 , λ, β) under squared error loss function is given bŷ For the ease of calculation, the BE of g(α 0 , λ, β) is rewritten aŝ where l(α 0 , λ, β) = ln(L(α 0 , λ, β)) and ρ(α 0 , λ, β) = ln(π(α 0 , λ, β)).
Considering the complexity of the likelihood function (8) and the ratio of two multiple integrals,ĝ BE (α 0 , λ, β) cannot be given by an analytical form. Therefore, it is necessary to employ some approximation methods to obtain the approximate posterior expectations. Among others, the methods suggested by Tierney and Kadane [23] and Lindley [24] are commonly used and perform well. In the next subsections, the Tierney and Kadane technique and Lindley's approximation are applied.

Optimal Design
In addition to the study of parameter estimation, this paper also includes the problem of selecting the optimal transformed stress level h * i (i = 1, . . . , k) of CSALT for the Chen distribution under progressive Type-II censoring. In CSALT, test units are assigned to run under different stress levels to better estimate the life characteristics of the units under the normal stress level. For different combinations of stress levels, the performance of parameter estimation will be different. Therefore, to improve the accuracy of product life characteristic estimation, the selection of stress levels is very important.
According to Nelson [2], the best test plan has two stress levels to produce the most precise estimation of life distribution parameters (minimum variance) under design stress levels. Therefore, we consider and discuss the case k = 2 when solving the optimal transformed stress level. Since the minimum transformed stress level is constant at h 1 = 1, we only need to solve the optimal stress level of h 2 . Here, we investigate this problem according to the following two optimal criteria.

D-Optimality
D-optimality is very common in designing ALT, which maximizes the determinant of FIM. The determinant of FIM is known to be proportional to the determinant of the inverse of the asymptotic variance-covariance matrix, which is proportional to the total volume of the joint confidence interval of the parameter (α 0 , λ, β). To estimate the life parameters of the unit under the normal stress level with greater accuracy, when the progressive censoring scheme is determined, the smaller the total volume of the joint confidence region of the parameter (α 0 , λ, β), the better the estimation. Therefore, maximizing the determinant of the FIM within a suitable range can achieve this aim. Here, we continue to use the observed FIM O T (Θ) derived from the preceding section. The criterion can be expressed as follows max |O T (α 0 , λ, β)|. (52)

A-Optimality
Similarly, according to the relationship between the asymptotic variance-covariance matrix and the joint confidence region of the parameter (α 0 , λ, β), we can also minimize the trace of the inverse of FIM to obtain the optimal solution. It provides a global indicator to measure the overall change of MLE. We still use observed FIM O T (Θ) to obtain the optimal transformed stress to estimate the life parameters of the unit under the normal stress level more accurately when the progressive censoring scheme is determined. The criterion can be expressed as follows

Simulation Study
In this section, a numerical study is carried out with the Monte Carlo simulation using R software to assess the performance of the methods for point and interval estimation discussed in the preceding sections. Initially, we set the true values of the parameters to α 0 = 0.45, λ = 2.0 and β = 0.7. We also set the CSALT with two stress levels that h = (1, 8.0) and k = 2. Under both stress levels s 1 and s 2 , suppose that the number of assigned test units n i , (i = 1, 2), the number of failed test units m i , and the censoring schemes R i do not depend on i. We set three different sample sizes (30, 15), (45, 25) and (60, 35), and consider three different censoring schemes in each sample size. Among the censoring schemes, (0 * 7, 15, 0 * 7) denotes that no unit is withdrawn for the first 7 times, then 15 units are withdrawn at one time, and no unit is withdrawn for the last 7 times. It is similar to other notations. By utilizing the algorithm from Balakrishnan and Sandhu [26], the progressive Type-II censored samples are generated.
For all point and interval estimation methods, the simulation is repeated 5000 times. Then, we take the expected values (EV) as the final point estimates and take the mean values of lower bounds (LB) and upper bounds (UB) as the final interval estimates. Additionally, we calculate the mean square error (MSE) for point estimation and the average interval length (AIL) and the coverage probability (COVP) for interval estimation to compare the performance of the estimation from various methods.
According to the log-likelihood function (9), we use the progressively censored samples to obtain the MLEs of the parameters in Table 1. In computing BEs, we set (a 1 , a 2 , a 3 ) = (0.1, 0.1, 0.1) and (b 1 , b 2 Tables 2 and 3, respectively. Further, Table 4 shows the corresponding 95% asymptotic confidence intervals based on the asymptotic property of MLEs and the approximate variance of the parameters in the matrix O −1 T (Θ). For the interval estimation, we also apply the parametric bootstrap method to find the 95% bootstrap percentile confidence interval and the 95% bootstrap-t confidence interval that are shown in Tables 5 and 6, where we set B = 1000. Finally, utilizing D-and A-optimal cri-teria, we obtain the average optimal transformed stress level for each progressive censoring scheme that are shown in Tables 7 and 8.  From Tables 1-8, the following conclusions can be drawn. (1) According to Tables 1-6, the parameter estimation under different progressive censoring schemes all roughly follows the pattern of better performance with increasing sample size. Different censoring schemes have different performance in different estimation methods. As the censoring scheme of withdrawing one surviving unit for each unit observed to fail from the beginning of the test (for example, (1 * 15)), the performance in point estimation is poor, especially when the sample size is relatively small and the estimation bias is relatively large. However, this censoring scheme has no obvious difference from other schemes in interval estimation.
(2) According to the results of the point estimates in Tables 1-3, the estimated values are all close to the true values and, roughly, the MSE decreases as the sample size increases, so all the point estimation methods mentioned are valid. Based on the MSEs of the parameter estimators, the three methods perform similarly when the progressive censoring scheme is determined. For the estimation of α 0 , the BE obtained from the Lindley's approximation performs better than BE obtained from the Tierney and Kadane technique. Additionally, BE obtained from the Tierney and Kadane technique performs slightly better than MLE. For λ estimation, BE obtained from the Lindley's approximation performs best, and the BE obtained from the Tierney and Kadane technique and MLE perform similarly. For β estimation, BE obtained from the Lindley's approximation performs best, followed by BE obtained from Lindley's approximation and MLE. Overall, BE obtained from the Lindley's approximation performs well for all three parameters.
(3) According to the results of the interval estimation in Tables 4-6, generally speaking, the asymptotic confidence intervals perform best. The bootstrap percentile confidence intervals and the bootstrap-t confidence intervals perform similarly when the progressive censoring scheme is determined. For all three interval estimation methods, the AIL of λ is longer than those of α 0 and β. The COVP of the intervals of α 0 outperforms those of both λ and β in parametric bootstrap intervals.
(4) According to Tables 7 and 8 for the optimal transformed stress level, the values based on A-optimality are smaller than those based on D-optimality. For a similar censoring scheme, the results based on D-optimality become larger as the sample size increases. However, the results based on A-optimality are the opposite. In general, the optimal transformed stress levels do not differ significantly in each of the optimal criteria.      Table 7. The average optimal transformed stress level h * k using D-optimality for different sample sizes (n i , m i ) and different censoring schemes R i .

Real Data Analysis
A real data set on the breakdown times (in minutes) of a type of insulating fluid from Chapter Three of Nelson [2] is considered to illustrate the methods. The normal stress level is s 0 = 20 kilovolt (kV). In CSALT, the complete failure data under both stress levels are presented in Table 9. To ensure the fit of the data to the Chen distribution, we calculate the Kolmogorov-Smirnov statistics and p-values, also shown in Table 9. Furthermore, the plots of the empirical cumulative distribution function and the fitted cumulative distribution function are shown in Figure 9. It is clear that the data set fits the model well.
Based on the chemical reaction mechanism leading to the failure, the inverse power model is a suitable acceleration model. Therefore, we let φ(s i ) = ln(s i ), (i = 1, 2). The value of h is then calculated as h = (1, 1.44966). Under both stress levels s 1 = 30 kV and s 2 = 36 kV, the number of assigned test units n i , i = 1, 2, the number of failed test units m i , the progressive Type-II censoring schemes R i and the censored data are shown in Table 10. The point estimates are shown in Table 11 and the intervals estimates are shown in Table 12. The number below each interval in Table 12 is the length of that interval. For Bayesian estimation, the results are calculated with non-informative priors in which the hyperparameters all tend to zero.
From Tables 11 and 12, the following conclusions can be drawn: (1) The point estimates of α 0 and β obtained from the three different point estimation methods are similar, but the point estimates of λ are slightly different. For λ estimation, the MLE is maximal, followed by the BE obtained from the Lindley's approximation and the BE obtained from the Tierney and Kadane technique; (2) In terms of AIL, the results obtained by the three interval estimation methods do not differ significantly. The bootstrap-t confidence intervals outperform slightly the asymptotic confidence intervals and the bootstrap percentile confidence intervals. This is consistent with the pattern that parametric bootstrap intervals outperform asymptotic intervals when the sample size is small.
(3) For the bounds of the interval estimates of α 0 , the bootstrap percentile confidence intervals and the bootstrap-t confidence intervals are similar. For the interval estimates of λ, the upper and lower bounds of the asymptotic confidence intervals are significantly larger than those of the bootstrap percentile confidence intervals and the bootstrap-t confidence intervals. The bounds of the bootstrap-t confidence intervals are slightly larger than those of the bootstrap percentile confidence intervals. For the bounds on the interval estimates of β, the results obtained by the three methods are similar.

Conclusions
This paper considers statistical inference and optimal design of CSALT for the Chen distribution under progressive Type-II censoring. Specifically, maximum likelihood estimates are obtained utilizing the Newton-Raphson algorithm. Bayesian estimation under the squared error loss function is also considered and estimated by the Tierney and Kadane technique and Lindley's approximation. At the same time, we establish confidence intervals for the lifetime distribution parameters and compare them. We construct the asymptotic intervals based on the observed Fisher information matrix and calculate the bootstrap percentile confidence intervals and the bootstrap-t confidence intervals to address the problem of small sample sizes. In addition, two optimal criteria are applied to determine the optimal transformed stress level that would provide more accurate estimates of the life parameters under the normal stress level. The performance of these estimation methods and the optimal criteria are compared using Monte Carlo simulation methods.
The Chen distribution possessing a bathtub-shaped hazard function is of great significance and practical importance in life experiments. Accelerated life tests and progressive censoring also play important roles. Further research into these models has great potential for estimating lifetime characteristics. This paper can also be extended by considering other distributions or other censoring schemes.