Accelerated life tests for modified Kies exponential lifetime distribution: binomial removal, transformers turn insulation application and numerical results

Abstract: This paper is concerned with statistical inference of multiple constant-stress testing for progressive type-II censored data with binomial removal. The failure times of the test units are assumed to be independent and follow the modified Kies exponential (MKEx) distribution. The maximum likelihood method as well as Bayes method are used to derive both point and interval estimates of the parameters. Furthermore, a real data application for transformers turn insulation is used to illustrate the proposed methods. Moreover, this real data set is used to show that MKEx distribution can be a possible alternative model to the exponential, generalized exponential and Weibull distributions. Finally, simulation studies are carried out to investigate the accuracy of the different estimation methods.


Introduction
The purpose of a life testing is to analyze failure times of test units that obtained under normal operating conditions. As it is well known, the modern products are designed to last longer, so for these products, collecting failure data under ordinary circumstances is entirely difficult or even impractical. In such situations, items should be exposed to stress higher than the manufacture levels The paper is organized as follows: The MKEx distribution and test assumptions for CSALT are discussed in Section 2. The maximum likelihood estimation (MLE) for the model parameters are provided in Section 3. Bayes estimates (BEs) under different loss functions are discussed in Section 4. In Section 5, asymptotic and bootstrap confidence intervals (CIs) are presented. Section 6 contains the simulation numerical results. Section 7 contains the transformers turn insulation application. Section 8 contains the conclusion of the paper and the major findings in this research.

MKEx distribution
In 2013, Kumar and Dharmaja [26] introduced and studied the reduced Kies distribution. In many references, the reduced Kies distribution is known as modified Kies (MK) distribution. Kumar and Dharmaja [26] find out that the MK distribution can perform better than common lifetime distributions such as Weibull distribution and some of its extensions in modelling lifetime data. Kumar and Dharmaja [27] presented the exponentiated MK distribution and introduced its statistical properties. Dey et al. [20] estimated the distribution parameters of MK distribution under progressive type-II censoring and introduced the recurrence relations for the moments of the MK distribution. Based on MK distribution and the T − X family, Al-Babtain et al. [15] introduced a new lifetime distribution, they called it modified Kies exponential (MKEx) distribution. Aljohani et al. [40] discussed the parameter estimation of the MKEx distribution using the MLE method, based on ranked set sampling. It has bathtub shape, increasing and decreasing failure rate. Furthermore, It has the ability to model negatively and positively skewed data. Moreover, it has a closed form cumulative distribution function (CDF) and very easy to handle which make the distribution is candidate to use in different fields such as life testing, reliability, biomedical studies and survival analysis. In this context, Al-Babtain et al. [15] used two different types of real data to show that this distribution may be a good alternative to many popular distributions such as Weibull, Marshall-Olkin exponential, Kumaraswamy exponential, beta exponential, gamma, and exponentiated exponential distribution. The CDF of the MKEx distribution is F(x; α, σ) = 1 − e −(e σx −1) α , x > 0, α, σ > 0. (2.1) The corresponding probability density function (PDF) of (2.1) is given by f (x; α, σ) = ασ e ασx 1 − e −σx α−1 e −(e σx −1) α , x > 0, α, σ > 0. (2.2)

Multiple CSALT's assumptions
In this subsection, we introduce the assumption of CSALT under PT-II CS with binomial removal. Suppose an ALT contains number of stress levels L ≥ 2 such that the stress is arranged ascendingly where φ 1 < φ 2 < ... < φ L , within the level l, l = 1, 2, ..., L, identical n l units are exposed to an accelerated condition, so that the number of units under the lifetime experiment are L l=1 n l = n, where n is the whole number of tested items in the test. In each stress level, φ l , l = 1, 2, ..., L, at the time of the first failure x l1:m l :n l , R l1 of the n l − 1 remaining units are randomly excluded from the test. In the same manner R l2 of the surviving units, n l − R l1 − 1, are randomly excluded from the test after the second failure x l2:m l :n l is detected. This mechanism continues until the failure of m th l occurs. The remaining surviving units R lm l = n l − m l − m l −1 j=1 R l j are excluded from the test after the m th l occurs, and the test is terminated. Suppose that the elimination of an individual unit from the test is independent of the others but with the same probability of removal P. Then, the number of units withdrawn at each failure time has a binomial distribution. That is R 1 ∼ binomial(n l − m l , P), R l j ∼ binomial(n l − m l − m l −1 l=1 R l j , P), l = 2, ..., m l and R lm l = n l − m l − m l −1 j=1 R l j . In this context, the assumptions of multiple CSALT are as follows: 1. In each stress level φ l , the lifetime of the experimental units follow MKEx(α, σ l ) distribution. 2. The scale parameter in each stress level σ l and the stress level φ l is linked by the following relation. log(σ l ) = ζ + βη l , l = 0, 1, ..., L, where ζ ∈ (−∞, ∞) and β > 0 are the unknown model parameters and η l = η(φ l ) is an increasing function of φ. For more extensive reading about acceleration and its different models we can refer to the book of Nelson [36], specifically Chapter 2.

Estimation via maximum likelihood method
In this section, the classical estimates of the parameters of MKEx distribution under PT-II CS with binomial removal are obtained. As we mentioned later that the R l j has a binomial distribution, then the PDF of R l j , l = 1, 2, ..., L is given as follows: while, for j = 2, 3, ..., m l − 1: Pr(R l j = r l j |R (l−1) j ) = n l − m l − l−1 j=1 r l j r l j P r l j (1 − P) (n l −m l − l j=1 r l j ) , where 0 ≤ r l j ≤ n l − m l − l−1 j=1 r l j . Furthermore, we suppose that R l j is independent of X l j:m l :n l for all l, l = 1, 2, ..., L. Therefore, the likelihood function α, σ 0 , θ under PT-II censoring with binomial removal is given by since X l j:m l :n l and R l j for all l = 1, 2, ..., L are independent, then the MLE of P can be derived by maximizing Pr(R l j = r l j ) directly. Hence the MLE of P is given bŷ The likelihood function under PT-II censoring after estimating P, have the following form: where τ l j = x l j:m l :n l , The log-likelihood function after removing the normalizing constant, can be formed as in Eq (3.3) 1 + R l j e σ 0 θ z l τ l j − 1 α .

(3.3)
By finding partial derivatives of with respect to the distribution parameters, we have where l j = σ 0 θ z l τ l j . The MLE of (α, σ 0 , θ) is (α,σ 0 ,θ), which may be derived simultaneously by solving Eqs (3.4)-(3.6). Regrettably, solving these equations will be very hard, so we have to use numerical techniques like the Newton-Raphson method.

Bayesian estimation
This section includes the BEs of α, σ 0 and θ. We assume that α, σ 0 and θ are independent and have gamma priors. Gamma prior for the acceleration factor θ > 1 was first considered by DeGroot and Goel [19]. They stated that in most problems of accelerated life testing the accelerating factor θ will be greater than 1. However, in order to not restrict the applicability of the acceleration model we shall consider prior distributions for θ that assign positive density to all positive values of θ. If the experimenter is almost certain that θ > 1, then he can choose a gamma prior distribution that assigns a suitably small probability to the interval 0 < θ < 1. For more details about this point, see Section 3 of DeGroot and Goel [19]. The gamma priors for distribution parameters are as follows: The joint prior of α, σ 0 and θ is obtained as To determine suitable and superior values of the hyper-parameters of the independent joint prior, we can use estimate and variance-covariance matrix of MLE method. By equating mean and variance of gamma priors as the following equations, the estimated hyper-parameters can be computed as where, B is the number of iteration andΩ 1 =α,Ω 2 =σ 0 , andΩ 3 =θ . By multiplying Eq (3.2) by (4.4), and making some simplifications the posterior distribution C * (α, σ 0 , θ) is formed as follows The BEs of u(Θ) = u(α, σ 0 , θ) using squared error (SE) and LINEX loss functions are as follows It is obvious that both BEs of u(α, σ 0 , θ) in (4.8) and (4.9) are considered as the division of more than one integration over each other. As we know multiple integrals is very tough to be solved analytically or even mathematically by hand. Therefore, we have to use the Markov Chain Monte Carlo (MCMC) technique to find an approximate value of integrals. An important methods of the MCMC technique, is the Metropolis-Hastings (MH) algorithm, some times they call it the random walk algorithm. Its similar to acceptance-rejection sampling, the MH algorithm consider for each iteration of the algorithm, a candidate value can be generated from normal proposal distribution.
The conditional posterior distributions used in MH algorithm are given as follows: e ασ 0 θ z l τ l j 1 − e −σ 0 θ z l τ l j α−1 e −(1+R l j) e σ 0 θ z l τ l j −1 α , (4.10) and The Bayesian estimates have CIs which are called the credible intervals or some times we call it the highest posterior density (HPD) intervals, for more information see, Chen and Shao [18]. They performed a technique that was used extensively to generate the HPD intervals of unknown parameters of the distribution. In this method, samples are drawn with the proposed MH algorithm that are used to generate estimates, for the HPD algorithm see Chen and Shao [18].

Confidence intervals
In this section, the asymptotic, percentile Bootstrap and Bootstrap-t confidence intervals (CIs) for the unknown distribution parameters α, σ 0 , θ are obtained.

Asymptotic confidence intervals
Asymptotic CI is the most popular approach to establish the approximate confidence limits for parameters, we use the MLE to obtain the observed Fisher information matrix I(Ω), which consists of of the negative second derivative of the natural logarithm of the likelihood function evaluated at Ω = (α,σ 0 ,θ), where and by inverseing this matrix we can find the asymptotic variance-covariance matrix. Now we can find the asymptotic variance-covariance matrix of the parameter vector Ω is V(Ω) = I −1 (Ω). So the 100(1 − γ)% asymptotic confidence intervals for parameters α, σ 0 and θ can be established as follows: where ϑ is α, σ 0 or θ, and Z q is the 100q − th percentile of a standard normal distribution.

Bootstrap confidence interval
In this subsection, two parametric bootstrap methods: Percentile bootstrap (B-P) and the bootstrap-t (B-T) [23] are considered to obtain CIs for α, σ 0 , and θ. The percentile bootstrap CIs can be obtained in such a way.

Repeat step number
). 5. The two side 100(1 − γ)% percentile bootstrap confidence interval for α, σ 0 and θ are evaluated The bootstrap-t CIs can be obtained in such a way.
1. Repeat the first two steps in the percentile bootstrap algorithm.

4.
A two side 100(1 − γ)% bootstrap-t CI for α, σ 0 and θ are given by In this section, we conduct a Monte Carlo simulation to find the estimates of the distribution parameters with different sample sizes n l and different censoring schemes R l with different probability of binomial removal P. We compare the performance of the MLEs and the BEs under different loss functions in terms of their relative absolute bias (RABias), and the mean square error (MSE). Furthermore, we estimate the length of asymptotic CI (L.CI), length of B-P CI (L.BP), length of B-T CI (L.BT) and for Bayesian estimation method we estimate length of credible interval (L.Cr). In the simulation, we consider different sample sizes, n l , different number of failures, m l , and different ratio of failures, r l , where r l = m l /n l . In addition, probability of binomial removal P is considered to be 0.35, and 0.8 for each stress level l, l = 1, 2, ..., L. The true values of the parameters used in the simulation study are (α = 2, σ 0 = 1.5, θ = 2) and (α = 0.8, σ 0 = 1.5, θ = 2). The simulation study is done using 10, 000 iterations, and the average of the results of these iterations are tabulated in Tables  1-9.
The MCMC iterations and the kernel histograms of the posterior samples of the parameters α, σ 0 , and θ are plotted in Figure 1. Figure 1 shows that the MCMC samples are well mixed and stationary achieved. Also, it indicates that posterior distributions of the three parameters are symmetric. Table 1. RABias and MSE of the MLE for binomial removal parameter P based on CSALT. The following observations are conducted from the simulation study: 1. For fixed m l , P, η l and true values of the parameters, the RABias, MSE and length of CIs decrease as n l increases. 2. The best method of estimation is the Bayesian estimation according to the the values of the MSE. 3. The BEs under SE loss function are better than the corresponding estimates under LINEX loss function with c = −0.5 for all parameters. 4. BEs under LINEX loss function for α and σ 0 with c = 0.5 are better than the corresponding estimates under SE loss function for all cases. 5. BEs under LINEX loss function for θ, with c = 0.5 is better than SE loss function when α is less than 1, while, when α is greater than 1, BEs under SE loss function are better.     Table 6.  Table 7.  Table 8.  Table 9.

Application of transformer insulation
In this section, a real data set is analyzed to illustrate the proposed methods in the previous sections. Furthermore, this data set is used to show that the MKEx distribution can be a possible alternative to widely known distributions such as exponential distribution, generalized exponential distribution and Weibull distribution.
Nelson [36] presented in Chapter three of his book the results of a constant-stress accelerated life test of a transformer insulation. The test consisted of three levels of constant voltage, which are respectively 35.4kv, 42.4kv and 46.7kv with normal voltage is 14.4kv. The results of such test are presented in Table  10. In this table, the sign "+" refers to the censored data. For each level of this test, we suggested using the following progressives CSs as follows: In the following subsection, we explain how to perform a goodness of fit test for the data in Table 10 and the proposed MKEx distribution.

Modified KS algorithm for fitting progressive censored data
When the data is PT-II censored data, we have to use modified KolmogorovSmirnov (KS) goodness of fit test. The modified Kolmogorov-Smirnov statistic for PT-II censored data was originally introduced by Pakyari and Balakrishnan [38]. This algorithm is based on several steps, first, find the estimates of the parameters for the proposed distribution and next transforming the data to normality, then testing the goodness of fit of the transformed data to normality. Let τ 1:m:n < τ 2:m:n < ... < τ m:m:n be a PT-II censored sample with CS (R 1 , R 2 , ..., R m ) from a distribution function F(t; θ), then the modified KS statistic for PT-II censored data is The following algorithm was proposed by Pakyari and Balakrishnan [38] to apply the KS test for PT-II censored data. Table 11 contains the MLEs and BEs using SE and LINEX loss functions of α, σ 0 and θ for the real data set. Based on the results of MLEs, the summary of results of the modified KS test for MKEx distribution are displayed in Table 12.  Table 13 contains the value of the life characteristics σ l , l = 0, 1, 2, 3 for each stress level, also it contains the mean time to failure (MTTF) for each stress level.     From Figures 2-5, we can note that with the increase in the stress value, the reliability function tends to zero faster.
The remainder of this section deals with the comparison of the proposed model, MKEx, generalized exponential (GE), Weibull, and exponential distributions. The comparison is performed using the real data shown in Table 10 and the PT-II CSs previously presented in this section. To clarify which of these distributions is more suitable for the data in Table 10, the parameters for the four distributions are estimated and the Akaike information criteria (AIC) is calculated for each distribution. The results of these calculations are summarized in Table 14. Furthermore, the results of modified KS test for the GE, Weibull, and exponential distributions are presented in Table 15. From Tables 14 and 15, we can see that the MKEx distribution provides a better fit to the given data compared to exponential, GE, and Weibull regarding AIC.

Conclusions
This paper discussed the statistical inference of CSALT under PT-II censoring with binomial removal when the lifetimes of test units follow the MKEx distribution. In this context, we obtained the point and interval estimates for the unknown parameters using both classical and Bayesian methods. We concluded that the Bayesian method was better than the classical method according to MSE and relative absolute bias of the estimates. Regarding the interval estimates, we noted that the percentile bootstrap interval was the best one according to the shortness of the interval length.
Furthermore, An application about the insulation of transformers was discussed and used to illustrate the theoretical results. Moreover, the data of insulation of transformers was used to show that the suggested model, MKEx, can be a possible alternative to some well known distributions.