Order restricted Bayesian inference for exponential simple step-stress model

ABSTRACT A step-stress model has received a considerable amount of attention in recent years. In the usual step-stress experiment, a stress level is allowed to increase at each step to get rapid failure of the experimental units. The expected lifetime of the experimental unit is shortened as the stress level increases. Although extensive amount of work has been done on step-stress models, not enough attention has been paid to analyze step-stress models incorporating this information. We consider a simple step-stress model and provide Bayesian inference of the unknown parameters under cumulative exposure model assumption. It is assumed that the lifetime of the experimental units are exponentially distributed with different scale parameters at different stress levels. It is further assumed that the stress level increases at each step, hence the expected lifetime decreases. We try to incorporate this restriction using the prior assumptions. It is observed that different censoring schemes can be incorporated very easily under a general setup. Monte Carlo simulations have been performed to see the effectiveness of the proposed method, and two datasets have been analyzed for illustrative purposes.


Introduction
In many reliability experiments, often an investigator has to wait a long period of time to observe failures. Accelerated Life Testing (ALT) experiment has been used quite often to observe early failures. In this setup, the experimental units are exposed to higher stress levels than usual to reduce the time to failure, hence to observe more failures within an affordable time. Data, collected in this method, needs to extrapolate to get back to the normal condition.
A particular case of the accelerated life-tests is step-stress life-test (SSLT), where the experimenter is allowed to change the stress levels during the life-testing experiment. In this case, a number of experimental units, say n, are placed on a test at an initial stress level s 1 and then the stress levels are changed to s 2 , s 3 , . . ., s m+1 at some prefixed times, say at τ 1 < τ 2 < · · · < τ m , respectively. A simple SSLT is a special case of SSLT, where only two stress levels are considered, and the stress level is changed from s 1 to s 2 at a prefixed time τ 1 . Moreover, to analyze such data we need a model that relates the distributions of lifetimes under different stress levels to that of lifetimes under the step-stress pattern. One such model is cumulative exposure model (CEM), first introduced by Seydyakin (1966), further studied by several authors, see, for example, Bagdonavicius (1978) and Nelson (1980). This model has been extensively discussed in the literature.
In this article, we consider a simple step-stress model, and it is assumed that the lifetime distributions at two different stress levels follow exponential distribution with mean lifetimes λ −1 1 and λ −1 2 , respectively, where λ 1 < λ 2 . Moreover, CEM is assumed. It may be mentioned that although, extensive amount of work has been done on step-stress models, not much attention has been paid to develop the order restricted inference. Balakrishnan et al. (2009a) considered the order restricted inference for exponential step-stress models when the data are Type-I or Type-II censored. They have mainly adopted the frequentist approach, and the maximum likelihood estimators (MLEs) of the unknown parameters are obtained using isotonic regression. It is observed that obtaining the exact joint distribution of the MLEs is not very easy, hence they derived the asymptotic distribution of the MLEs. Based on the asymptotic distribution, the asymptotic confidence intervals (CI) of the unknown parameters can be constructed. It is not immediate that how this method can be extended for other censoring schemes like hybrid and progressive censoring schemes.
In the life-testing experiment often the data are censored. The most popular censoring schemes are Type-I and Type-II censoring schemes. The hybrid censoring scheme (HCS) is a mixture of the Type-I and Type-II censoring schemes, and it was introduced by Epstein (1954). From now on, the HCS proposed by Epstein (1954) will be called as Type-I HCS. Similar to conventional Type-I censoring scheme, the main disadvantage of Type-I HCS is that almost all inferential results are based on the assumption that there are at least one failure. Moreover, there may be very few failures, hence the efficiency of the estimator might be very low. For these reasons, Childs et al. (2003) introduced Type-II HCS, which guarantees a minimum number of failures during the experiment. For an extensive survey of different hybrid censoring schemes, the readers are referred to Balakrishnan and Kundu (2013). Recently, progressive censoring scheme (PCS) has received significant attention in the statistical literature. The main advantage of progressive censoring schemes is that it is possible to remove experimental units during the experiment, even if they do not fail. For an exhaustive survey on PCSs, the readers are referred to the review article by Balakrishnan (2007). A brief review of the different censoring schemes is provided in Section 2.
Simple step-stress models under different censoring schemes are extensively studied based on the assumption that the lifetime of the experimental units follow exponential distributions with different scale parameters at different stress levels. From now on, we call them as exponential step-stress models. The simple exponential step-stress model under Type-I censoring is considered by Balakrishnan et al. (2009b). Balakrishnan et al. (2007) considered simple exponential step-stress model under the Type-II censoring scheme. Simple exponential stepstress models under HCS-I and HCS-II are considered by Balakrishnan and Xie (2007b) and Balakrishnan and Xie (2007a), respectively. In all these cases, the exact distributions of the unknown parameters are obtained, and they can be used to construct exact CIs. However, it is observed that the exact distribution and hence the construction of associated CI is quite complicated in all these cases. Moreover, all the inferential issues are obtained without the order restriction on the unknown parameters. It is clear that the ordered restricted inference will be quite complicated in the frequentist setup, see Balakrishnan et al. (2009a). It seems Bayesian analysis is a natural choice in these cases. Some work has been done on the Bayesian inference of the step-stress model, see, for example, Dorp et al. (1996), Lee and Pan (2011), Leu and Shen (2007), Fan et al. (2008), and Liu (2010). However, none of them dealt with the ordered restricted inference.
The main aim of this article is to consider the order restricted Bayesian inference of the unknown parameters of a simple exponential step-stress model under different censoring schemes. We have assumed fairly flexible priors on the unknown parameters. It is observed that in all the cases the Bayes estimates of the unknown parameters cannot be obtained in an explicit form. We propose to use importance sampling technique to compute Bayes estimate (BE) and also to construct associated credible interval (CRI). We also discuss construction of credible set for model parameters. Extensive Monte Carlo simulations are performed to see the effectiveness of the proposed method, and the performances are quite satisfactory. The analyses of two datasets have been performed for illustrative purposes.
The rest of the article is organized as follows. In Section 2, we briefly discuss different censoring schemes and available data. Model assumptions and the prior information of the unknown parameters are considered in Section 3. In Section 4, the maximum likelihood estimation of model parameters is briefly discussed under the order restriction for Type-I censored data. In Section 5, we provide the posterior analysis for different censoring schemes under the order restriction. Monte Carlo simulation results are presented in Section 6.1 Data analyses have been provided in Section 6.2 Finally, we conclude the article in Section 7.

Different censoring schemes and available data
A total of n units is placed on a simple SSLT experiment. The stress level is changed from s 1 to s 2 at a prefixed time τ 1 , and τ 2 > τ 1 is another prefixed time. The positive integer r ≤ n is also prefixed. The role of r and τ 2 will be clear later. Let the ordered lifetimes of the items be denoted by t 1:n < · · · < t n:n . Now we briefly describe different censoring schemes, and available data in each case.
Type-II progressive censoring scheme. In this case, it is assumed that n experimental units are put in a life test. R 1 , . . . , R m are m prefixed non-negative integers such that At the time of the first failure, say t 1:n , R 1 units are chosen at random from the remaining (n − 1) units and they are removed from the experiment. Similarly, at the time of the second failure, say t 2:n , R 2 units are chosen at random from the remaining (n − R 1 − 2) surviving units and they are removed from the test, and so on. Finally at the time of the mth failure, say t m:n , the rest of the n − m − m−1 j=1 R j = R m units are removed and the experiment is stopped. In this case, the available data are of the form (a) {τ 1 < t 1:n < · · · < t m:n }, (b) {t 1:n < · · · < t n 1 :n < τ 1 < t n 1 +1:n < · · · < t m:n }, (c) {t 1:n < · · · < t m:n < τ 1 }.
In Case (b), n 1 is number of failures at the stress level s 1 .

Model assumption and prior information
We consider a simple SSLT, where n identical units are placed on a life testing experiment at the initial stress level s 1 . The stress level is increased to a higher level s 2 at a prefixed time τ 1 . It is assumed that the lifetimes of the experimental units are independently and exponentially distributed random variables with different scale parameters at different stress levels. Probability density function (PDF) and the cumulative distribution function (CDF) of the lifetime under stress level s i for i = 1, 2 are given by and respectively. Let us assume that the stress level is changed from s 1 to s 2 at the time point τ 1 . It is further assumed that the failure time data come from a CEM, hence, it has the following CDF; The corresponding PDF is given by For developing the Bayesian inference, we need to assume some priors on the unknown parameters. We want the prior assumption on λ 1 and λ 2 , so that it maintains the order restriction, namely, λ 1 < λ 2 . We take the following priors on λ 1 and λ 2 . It is assumed that λ 2 has a Gamma(a, b) distribution with a > 0 and b > 0, i.e., it has the following PDF: Moreover, λ 1 = α λ 2 and α has a beta distribution with parameters c > 0 and d > 0, i.e., the PDF of α is given by and the distribution of α is independent of λ 2 . Therefore, the joint prior of (λ 1 , λ 2 ) can be written as As the joint prior on (λ 1 , λ 2 ) is little complicated, a gray-scale plot is provided in Fig. 1 for different values of hyper-parameters. In the plot, black color represents the maximum value of the density function, whereas white color represents the minimum value (which is zero) of the density function. We have taken b = 1.0 only, as different values of b only effects the spread of the density function keeping the shape fixed.

Maximum likelihood estimator under Type-I censoring scheme
In this section, we present maximum likelihood estimation of the scale parameters under the restriction λ 1 ≤ λ 2 , when data are Type-I censored. Let n * 1 and n * 2 denote the number of failures before the time τ 1 and between τ 1 and τ 2 , respectively. They can be zero also. Let τ * and n * denote the termination time of the experiment and total number of failures observed before τ * , respectively. Note that τ * and n * depend on the censoring scheme. In case of Type-I censoring scheme, τ * = τ 2 . For Case (a): n * 1 = 0, n * 2 = n 2 > 0, Case (b): n * 1 = n 1 > 0, n * 2 = n 2 > 0, Case (c): n * 1 = n 1 > 0, n * 2 = 0. In all the cases, n * = n * 1 + n * 2 . Based on the observations from a simple SSLT under Type-I censoring scheme, the likelihood can be written as where Note that d 1 and d 2 are the total time elapsed by all the units at stress level s 1 and s 2 , respectively. The unrestricted MLEs of λ 1 and λ 2 are given by Clearly, if λ * 1 ≤ λ * 2 , MLEs of the scale parameters under the restriction λ 1 ≤ λ 2 are given by As l 1 (λ 1 , λ 2 | Data) is unimodal function, if λ * 1 > λ * 2 , maximization of l 1 (λ 1 , λ 2 | Data) under the order restriction λ 1 ≤ λ 2 is equivalent to maximization of l 1 (λ 1 , λ 2 | Data) under λ 1 = λ 2 , and hence, in this case the MLEs of the scale parameters under the restriction λ 1 ≤ λ 2 are given by

Type-I censoring scheme
Based on the likelihood function in (8), priors π 1 (·) and π 2 (·) mentioned in Section 3, posterior density function of (α, λ 2 ) becomes The right-hand side of (9) is integrable if n * 1 + c > 0 and n * + a > 0. Bayes estimate of some function of α and λ 2 , say g(α, λ 2 ), with respect to the squared error loss function, is posterior expectation of g(α, λ 2 ), i.e., Unfortunately, the close form of (10) cannot be obtained in most of the cases. One may use numerical techniques to compute (10). Alternatively, other approximations, like Lindey's approximation, can be used to compute (10). However, CRI for a parametric function cannot be constructed by these numerical methods. Hence, we propose to use importance sampling to compute BE as well as to construct CRI of a parametric function in this article. Note that for 0 < α < 1 and λ 2 > 0, l 2 (α, λ 2 | Data) can be expressed as where The proportionality constant, c 1 , for the posterior distribution of α given in (12) can be found using numerical techniques. However, generation from (12) is not a trivial issue. Hence, we propose to use the importance sampling (see Algorithm 5.1) to compute the BE and as well as to construct CRI of g(α, λ 2 ) noting the following representation of l 2 (α, λ 2 | Data). For 0 < α < 1 and λ 2 > 0 where Algorithm 5.1.
Step 3. Find the integer M γ such that Step 4. Construct the HPD credible set for (α, λ 2 ) as Similar methodology can be applied for other censoring schemes, and we briefly mention all the cases in the subsequent subsections for completeness purposes.

Type-II censoring scheme
Based on the observed sample, the likelihood function is given in (8), where τ * = t r:n , in Case (a), n * 1 = 0, n * 2 = r, in Case (b), n * 1 = n 1 , n * 2 = r − n 1 and in Case (c), n * 1 = r, n * 2 = 0. d 1 and d 2 have the same expression as given in case of Type-I censoring.

Progressive Type-II censoring scheme
With the observed Progressive Type-II censoring data, the likelihood function is given by (8), where for Case (a), n * 1 = 0, n * 2 = m, for Case (b), n * 1 = n 1 , n * 2 = m − n 1 and for Case (c) n * 1 = m, n * 2 = 0. For all the cases τ * = t m:n , d 1 = N * 1 k=1 (R k + 1)t k:n + (n − N * 1 − N * k=1 R k )τ 1 and d 2 = m k=N * 1 +1 (R k + 1)(t k:n − τ 1 ). In all the above cases, the likelihood function are in the same form as Type-I censoring scheme and hence, the posterior density is also in the same form as given in (9). In all these cases, computation of the BE and construction of the associated CRI for some function of α and λ 2 can be done exactly along the same line. One can also construct credible set for (λ 1 , λ 2 ) following the same methodology.

Simulation results
In this section, we present some simulation results to see how the BE works for different sample sizes and for different values of τ 1 and τ 2 . Along with the coverage percentage (CP) and average length (AL) of symmetric CRI, HPD CRI, same of bootstrap CI is also presented for a comparison purpose. Here, we choose λ 1 = 1/12 0.083 and λ 2 = 1/4.5 0.222. We also choose a = 0.001, b = 0.001, c = 1 and d = 1, i.e., the noninformative prior and hence, the comparison with MLE is meaningful. All the results are based on 5000 simulations and symmetric CRI, HPD CRI, and bootstrap CI for same n, τ 1 , τ 2 , r are reported in Tables 2, 5 , 7, 9, and 11. Interested readers are referred to supplementary file for more simulation results. We also report the CP of the HPD credible set for (α, λ 2 ) in Table 3 using Algorithm 5.2 for the same parametric values when the data are Type-I censored.
The following points are quite clear from the simulation results. MSE of estimator of λ 1 decreases as τ 1 increases when other parameters are held constant. MSE of estimator of λ 2 decreases as τ 2 increases, whereas it increases as τ 1 increases. Also MSEs of all the estimators decrease as n increases. MSE of MLE is close to that of BE for λ 1 , but MSE of MLE is higher than MSE of BE for λ 2 . This difference decreases as expected number of failures at second stress level increases.
The performance of CI and CRIs of λ 1 are not satisfactory for all small sample sizes, but they are quite satisfactory for moderate and large sample sizes. We note that average length of HPD CRI of λ 1 decreases as τ 1 or n increases, keeping the other fixed. HPD CRI of λ 1 performs well compared to the symmetric CRI and bootstrap CI with respect to CP of the respective intervals, though AL of HPD CRI is larger than that of the bootstrap CI for Type-I, Type-II, and hybrid Type-I censoring schemes, but smaller than that of symmetric CRI. For moderate and large sample sizes, AL of symmetric CRI, HPD CRI, and bootstrap CI are very close to each other for all the censoring schemes.
The performance of the HPD CRI is not so satisfactory for λ 2 with respect to CP. However, CP of symmetric CRI and bootstrap CI is close to nominal level when expected number of failures at the first stress level is large. For small sample size performance of CRIs as well as bootstrap CI for λ 2 not at all satisfactory, especially when expected number of failures at first stress level is small. ALs of CRIs and bootstrap CI of λ 2 decrease as expected number of failures at the second stress level or n increases, keeping the other parameter constant. Also note that ALs of bootstrap CI and HPD CRI of the same parameter increase as τ 1 increases.

Data analysis
Example 1. Here, we consider the data (see Table 12) presented by Xiong (1998) to illustrate the methods of estimation discussed previously. This is actually a Type-II censored data from a simple step stress life experiment, where n = 20 units are placed on the test, the data are right censored at 16th failure time and stress changing time is τ 1 = 5. The average lifetimes are 3.52 and 7.70 at first and second stress levels, whereas standard deviations are 1.05 and 1.72. Balakrishnan et al. (2009b) used these data for illustrative example for the step-stress model under Type-I censoring scheme choosing τ 2 = 7, 8, 9, and 12. They reported the MLE of 1/λ 1 and 1/λ 2 and associated CIs for above mentioned τ 1 and τ 2 's under the assumption that the data are coming from exponential CEM. Choosing a = 0.001, b = 0.001, c = 1, d = 1 and M = 8000, the BE and Bayesian CRIs for 1/λ 1 and 1/λ 2 are reported in Tables 13 and 14,  respectively. We also obtain HPD credible set for (λ 1 , λ 2 ) and is given by where c 0 , n 2 , and D 2 depend on τ 2 and are presented in Table 15. Figure 2 shows the plot of the HPD credible set of (λ 1 , λ 2 ) for different values of τ 2 .
Example 2. Next, we consider the data (see Table 16) presented by Balakrishnan et al. (2009b). Choice made by Balakrishnan et al. (2009b) were n = 35, λ 1 = e −3.5 , λ 2 = e −2.0 , and τ 1 = 8. This is a complete dataset. Average lifetimes are 4.45 and 15.06 at the first and second stress levels, respectively. Standard deviations are 1.85 and 5.41, respectively. To make it a Type-I censored data one may take any τ 2 greater than 8. Balakrishnan et al. (2009b) took different choices for τ 2 , viz., 16, 20, and 24. Along these choices of τ 2 , we take τ 2 = 12 also. Assuming that the data are coming from the exponential CEM under Type-I censoring, they presented  Figure . Credible set of (λ 1 , λ 2 ) for the data in Table . MLE and associated CIs of 1/λ 1 and 1/λ 2 . Here, we present the BE of 1/λ 1 and 1/λ 2 and associated CRIs for the above-mentioned values of τ 2 under the same assumption. The results are presented in Tables 17 and 18 with a = 0.001, b = 0.001, c = 1, d = 1, and M = 8000. Like the previous example, we also obtain HPD credible set for (λ 1 , λ 2 ) and is given by C γ = (αλ 2 , λ 2 ) ∈ R 2 : c 0 (n 2 +8) α 8 λ n 2 +7 2 e 251.60α+D 2 ≥ c γ ,   Figure . Credible set of (λ 1 , λ 2 ) for the data in Table . where c 0 , n 2 , and D 2 depend on τ 2 and are presented in Table 19. Figure 3 shows the plot of the HPD credible set of (λ 1 , λ 2 ) for different values of τ 2 .

Conclusion
We have considered the Bayesian estimation of the unknown parameters in a simple SSLT under the restriction λ 1 < λ 2 and under different censoring schemes. The analysis is performed under exponentially distributed lifetimes and under CEM assumption. We have taken mainly the squared error loss function, though other loss functions can also be handled in a very similar way. We have seen that the BE of some parametric function under the square error loss function does not exist in a close form in most of the cases. Algorithms based on importance sampling are proposed to compute BE and to construct CRIs and credible set.
We have done a simulation study to judge the performance of the procedures described. We also considered two datasets to illustrate the estimation procedures. We have noticed that the performance of BE and CRIs for λ 1 is quite satisfactory at least for moderate and large sample sizes. It is also noticed that the performance of BE and CRI for λ 1 is better than that of MLE and bootstrap CI for small sample size and for the small values of τ 1 . For moderate and large sample sizes the performance of BE and CRI is quite close to that of MLE and bootstrap CI. We have also noticed that the performance of BE, CRI, MLE, and bootstrap CI of λ 2 are not satisfactory for small values of n and small values of expected number of failures at the second stress level. However, BE and CRI work quite well for moderate or large sample sizes and for large values of expected number of failures at the second stress level and the performance of BE and CRI is close to that of MLE and bootstrap CI of λ 2 in these cases. It is also noticed that HPD CRI works well for λ 1 , where symmetric CRI works well for λ 2 . Therefore, we recommend to use HPD CRI for λ 1 and symmetric CRI for λ 2 . Note that one may generate α from other distributions having support on (0,1) instead of uniform distribution as mentioned in Algorithm 5.1. A right-truncated gamma distribution and the distribution obtained by spline fitting to the posterior distribution of α have been tried. However, no significant improvement has been noticed. However, the results are quite prior dependent. The choice of priors is an important issue, which has not been pursued here and more research is needed in this direction.