1.
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 of stress, in order to obtain data about their failure times. We can call this kind of life test under tough conditions accelerated life testing (ALT), and the failure times of the test units from such ALT can be used to determine life characteristics under normal usage conditions.
There are many models to apply acceleration in the experiments such as Arrhenius model and inverse power model. Choosing the appropriate model of acceleration depends on the type of stress (voltage, pressure, or temperature) that the researcher wants to apply. Arrhenius model and inverse power model are usually used for the thermal and non-thermal stresses, respectively. The common kinds of ALTs are constant-stress ALT (CSALT), step-stress ALT (SSALT) and progressive-stress ALT (PSALT). In the CSALT the stress is constant during the whole time of the experiment, but the stress is increased gradually at specific times, in the SSALT. On the Contrast, the stress is increasing function of the time in the PSALT. ALTs can be divided based on the number of stress levels into two types: Simple and multiple level ALT. Simple ALT only has two levels of stress while, multiple ALT contains more than wo levels of stress.
Several authors have investigated the CSALT, Mohie El-Din et al. [32] considered the estimation problem of the CSALT for the extension of the exponential distribution under progressive type-II censoring. Mohie El-Din et al. [33] discussed the problem of obtaining the optimal CSALT designs for the Lindley distribution. The concern of estimation for CSALT for generalized half normal distribution under complete sampling was addressed by Wang and Shi [39]. Abd El-Raheem [6] developed CSALT's optimal designs for extension of the exponential distribution. The problem of the optimum configuration of CSALT to extension of the exponential distribution under censoring was considered by Abd El-Raheem [7]. For further references on the CSALT, see Abd El-Raheem [8,9,10].
The SSALTs are discussed by several authors see for example, Balakrishnan and Han [17], Ismail [24], Mohie El-Din et al. [30,31], Chandra and Khan [37] and Hafez et al. [22]. For more reading about PSALTs, we can refer to Abdel-Hamid and Al-Hussaini [3,2,4], AL-Hussaini et al. [12], Abdel-Hamid and Abushul [5], Mohie El-Din et al. [35,34] and Abd El-Raheem [11].
Censoring appears normally in reliability examinations. It happens when lifetimes are known only for a portion of the units under examination, the residual lifetimes being known only to exceed determined values. There two main types of censoring schemes which are type-I and type-II censoring schemes (CSs). The fundamental difference between type-I and type-II CSs is that the first depends on the time of ending the experiment and the second depends on the number of failures. The test units in these two CSs cannot be excluded from the test until the end the test. A progressive type-II (PT-II) censoring is an extension of type-II censoring. It encourages the experimenter to remove experimented units at different times throughout the test. For more reading about PT-II CS, we can refer to Balakrishnan and Aggarwala [16], Almetwally et al. [13], Alshenawy et al. [14] and El-Sherpieny et al. [21].
Due to the great importance of PT-II censoring in reliability experiments, many researchers addressed the issue of statistical inference of ALTs under PT-II censored data, we can refer to Abdel-Hamid [1], Abdel-Hamid and AL-Hussaini [2,4], Jaheen et al. [25] and Mohie El-Din et al. [29].
The motivation and contribution in this research is that we apply CSALT for the lifetime of units that follow the MKEx distribution under PT-II censored sample with binomial removal, also we estimate the parameters using Bayesian and classical methods. Moreover, we use a real life censored data of transforms insulation to illustrate the proposed methods.
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.
2.
Lifetime distribution and test assumptions
2.1. 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
The corresponding probability density function (PDF) of (2.1) is given by
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 nl units are exposed to an accelerated condition, so that the number of units under the lifetime experiment are ∑Ll=1nl=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 xl1:ml:nl, Rl1 of the nl−1 remaining units are randomly excluded from the test. In the same manner Rl2 of the surviving units, nl−Rl1−1, are randomly excluded from the test after the second failure xl2:ml:nl is detected. This mechanism continues until the failure of mthl occurs. The remaining surviving units Rlml=nl−ml−∑ml−1j=1Rlj are excluded from the test after the mthl 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 R1∼binomial(nl−ml,P), Rlj∼binomial(nl−ml−∑ml−1l=1Rlj,P),l=2,...,ml and Rlml=nl−ml−∑ml−1j=1Rlj. 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.
where ζ∈(−∞,∞) and β>0 are the unknown model parameters and ηl=η(ϕl) is an increasing function of ϕ.
(a) If η(ϕl)=log(ϕl), the model in (2.3) becomes the inverse power model.
(b) If η(ϕl)=1−ϕl, the model in (2.3) becomes Arrhenius model.
(c) If η(ϕl)=ϕl, the model in (2.3) becomes exponential model.
For more extensive reading about acceleration and its different models we can refer to the book of Nelson [36], specifically Chapter 2.
From Eq (2.3), we have
where θ=exp{β(η1−η0)}=σ1σ0>1, σ0>0 is the scale parameter of the MKEx distribution under usage conditions ϕ0, and
The transformation from the parameters (α,σl)=(α,ζ,β) to the new parameters (α,σ0,θ) is an one-to-one mapping. Since the Jacobian determinant from (α,ζ,β) to (α,σ0,θ) does not equal zero. Thus, the unknown parameters should be estimated are α, σ0 and θ.
3.
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 Rlj has a binomial distribution, then the PDF of Rlj, l=1,2,...,L is given as follows:
while, for j=2,3,...,ml−1:
where 0≤rlj≤nl−ml−∑l−1j=1rlj. Furthermore, we suppose that Rlj is independent of Xlj:ml:nl for all l, l=1,2,...,L. Therefore, the likelihood function α, σ0, θ under PT-II censoring with binomial removal is given by
since Xlj:ml:nl and Rlj for all l=1,2,...,L are independent, then the MLE of P can be derived by maximizing Pr(Rlj=rlj) directly. {Hence the MLE of P is given by
The likelihood function under PT-II censoring after estimating P, have the following form:
where τlj=xlj:ml:nl, Cl=nl(nl−1−Rl1)(nl−2−Rl1−Rl2)...(nl−ml+1−∑ml−1j=1Rlj).
The log-likelihood function after removing the normalizing constant, can be formed as in Eq (3.3)
By finding partial derivatives of ℓ with respect to the distribution parameters, we have
where ϵlj=σ0θzlτlj. 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.
4.
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
and
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. It's 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 MH algorithm produces a series of draws from MKEx distribution as follows:
1. initiate with α(0)=ˆα, σ(0)0=^σ0, θ(0)=ˆθ.
2. Set i=1.
3. Simulate α∗ from proposal distribution N(α(i−1),var(α(i−1))).
4. Evaluate the acceptance probability A(α(i−1)|α∗)=min[1,C∗(α∗|σ(i−1)0,θ(i−1))C∗(α(i−1)|σ(i−1)0,θ(i−1))].
5. Draw U∼U(0,1).
6. If U≤A(α(i−1)|α∗), put α(i)=α∗, else put α(i)=α(i−1).
7. Do the Steps 2-6 for σ0 and θ.
8. Put i=i+1.
9. Repeat Steps 3-8, N times to obtain (α(1),σ(1)0,θ(1)), ..., (α(N),σ(N)0,θ(N)).
Then, the BEs of u(α,σ0,θ) using MCMC under SE, and LINEX loss functions are respectively
where M is the burn-in period.
The conditional posterior distributions used in MH algorithm are given as follows:
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].
5.
Confidence intervals
In this section, the asymptotic, percentile Bootstrap and Bootstrap-t confidence intervals (CIs) for the unknown distribution parameters α,σ0,θ are obtained.
5.1. 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 Zq is the 100q−th percentile of a standard normal distribution.
5.2. 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.
1. Find the values of the MLE of MKEx distribution.
2. To find the estimates of the bootstrap, we must generate a bootstrap samples, (αb,σb0,θb), using MLEs of (α,σ0,θ).
3. Repeat step number (2) B times to have (αb(1),αb(2),…,αb(B)),(σb(1)0,σb(2)0,…σb(B)0) and (θb(1),θb(2),…,θb(B)).
4.Arrange (αb(1),αb(2),…,αb(B)),(σb(1)0,σb(2)0,…σb(B)0) and (θb(1),θb(2),…,θb(B)) from smallest to the biggest, (αb[1],αb[2],…,αb([B])),(σb[1]0,σb[2]0,…σb[B]0) and (θb[1],θb[2],…,θb[B]).
5. The two side 100(1−γ)% percentile bootstrap confidence interval for α,σ0 and θ are evaluated by [αb([Bγ2]),αb([B(1−γ2)])],[σb([Bγ2])0,σb([B(1−γ2)])0] and [θb([Bγ2]),θb([B(1−γ2)])].
The bootstrap-t CIs can be obtained in such a way.
1. Repeat the first two steps in the percentile bootstrap algorithm.
2. Evaluate T=ˆϑb−ˆϑ√V(ˆϑb), where ϑ is α,σ0 or θ, and V(ˆϑb) is asymptotic variances of ˆϑb.
3. Repeat the above two steps B times and rearrange (T(1),T(2),…,T(B)) in ascending order.
4. A two side 100(1−γ)% bootstrap-t CI for α,σ0 and θ are given by [α−Tb([Bγ2])α√V(αb),α+√V(αb)Tb([B(1−γ2)])α], [σ0−Tb([Bγ2])σ0√V(σb0),σ0+√V(σb0)Tb([B(1−γ2)])σ0] and [θ−Tb([Bγ2])θ√V(θb),θ+√V(θb)Tb([B(1−γ2)])θ]
6.
Simulation study
In this section, we conduct a Monte Carlo simulation to find the estimates of the distribution parameters with different sample sizes nl and different censoring schemes Rl 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, nl, different number of failures, ml, and different ratio of failures, rl, where rl=ml/nl. 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.
Lengths of asymptotic, percentile bootstrap, bootstrap-t, and credible CIs when L=2, α=2, σ0=1.5 and θ=2.
Lengths of asymptotic, percentile bootstrap, bootstrap-t, and credible CIs when L=4, α=0.8, σ0=1.5 and θ=2.
The following observations are conducted from the simulation study:
1. For fixed ml, P, ηl and true values of the parameters, the RABias, MSE and length of CIs decrease as nl 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.
6. The BEs under LINEX loss function with parameter c=0.5 have MSE less than corresponding estimates with c=−0.5, for all parameters.
7. Lengths of credible CIs are always shorter than the corresponding lengths of asymptotic CIs.
8. The percentile bootstrap CI has the shortest length among all considered CIs.
7.
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:
1. For ϕ1=35.4kv: n1=10, m1=8 and R1j=0, j=1,...,7, R18=2, with P=0.0277.
2.For ϕ2=42.4kv: n2=10, m2=9 and R2j=0, j=1,...,7, R28=1, R29=0, with P=0.0123.
3. For ϕ3=46.7kv: n3=10, m3=9 and R3j=0, j=1,...,8, R39=1, with P=0.0123.
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.
7.1. Modified KS algorithm for fitting progressive censored data
When the data is PT-II censored data, we have to use modified Kolmogorov-Smirnov (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 (R1,R2,...,Rm) from a distribution function F(t;θ), then the modified KS statistic for PT-II censored data is
where
and
where νi:m:n=E(Ui:m:n) is the expected value of the ith PT-II censored order statistic from the U(0,1) distribution, given by
and ui:m:n=F(ti:m:n;ˆθ) for i=1,2,...,m.
The following algorithm was proposed by Pakyari and Balakrishnan [38] to apply the KS test for PT-II censored data.
1.Find the MLE of the parameter θ, and calculate αi:m:n=F(ti:m:n;ˆθ) for i=1,2,...,m.
2.Evaluate yi:m:n=Φ−1(αi:m:n) for i=1,2,...,m.
3.Considering y1:m:n,y2:m:n,...,ym:m:n as a PT-II censored data from a normal distribution with mean μ and standard deviation σ, calculate the MLEs ˆμ and ˆσ.
4.Evaluate ui:m:n=Φ{(yi:m:n−ˆμ)/ˆσ} for i=1,2,...,m.
5.Evaluate Dm:n according to (7.1).
6.Evaluate the P-value of the test, and reject the null hypothesis H0 at significance level 0.05 if P-value less than 0.05.
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.
Figures 2-5 display the reliability function of transformers insulation under different stress levels.
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.
8.
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.
Acknowledgments
The authors are thankful for the Taif University researchers supporting project number TURSP-2020/160, Taif University, Taif, Saudi Arabia.
Conflict of interest
The authors declare no conflict of interest.