Robust inference for non-destructive one-shot device testing under step-stress model with exponential lifetimes

One-shot devices analysis involves an extreme case of interval censoring, wherein one can only know whether the failure time is either before or after the test time. Some kind of one-shot devices do not get destroyed when tested, and so can continue within the experiment, providing extra information for inference, if they did not fail before an inspection time. In addition, their reliability can be rapidly estimated via accelerated life tests (ALTs) by running the tests at varying and higher stress levels than working conditions. In particular, step-stress tests allow the experimenter to increase the stress levels at pre-fixed times gradually during the life-testing experiment. The cumulative exposure model is commonly assumed for step-stress models, relating the lifetime distribution of units at one stress level to the lifetime distributions at preceding stress levels. In this paper,vwe develop robust estimators and Z-type test statistics based on the density power divergence (DPD) for testing linear null hypothesis for non-destructive one-shot devices under the step-stress ALTs with exponential lifetime distribution. We study asymptotic and robustness properties of the estimators and test statistics, yielding point estimation and confidence intervals for different lifetime characteristic such as reliability, distribution quantiles and mean lifetime of the devices. A simulation study is carried out to assess the performance of the methods of inference developed here and some real-life data sets are analyzed finally for illustrative purpose.


Introduction
One-shot device testing is an increasingly important problem in the area of reliability.This involves an extreme case of interval censoring, wherein one only knows if the device works when it is tested.Most of the existing literature considers the case of "destructive" one-shot devices.This is the case when, once the device is used, it is either destroyed or must be rebuilt.Some typical examples are automobile air bags, fuel injectors, disposable napkins, missiles (Olwell and Sorell, 2001) and fire extinguishers and munition (Newby, 2008).Mainly motivated by the work of Fan et al. (2009), Balakrishnan and Ling (2012, 2013, 2014) developed efficient EM algorithms for the estimation of model parameters under the assumption of exponential, Weibull and gamma lifetime distributions, respectively.One may refer to the recent book by Balakrishnan et al. (2021) for a detailed review of all these works.Balakrishnan and Castilla (2021) recently developed results for the lognormal lifetime distribution.Some other related works about one-shot devices can be found in Mun et al. (2013), Sharma and Upadhyay (2018) and Zhu et al. (2021).However, the destructiveness assumption is not always necessary as in many experiments, the tested devices can be reused if has not failed during the test.We will refer to this type of devices as "non-destructive" one-shot devices.Their major advantage is that operating devices can continue in the experiment, providing extra information about their lifetime characteristics.Some typical examples are metal fatigue, thermal ageing of electrical insulation, spare wheels, safety valves, hot spare disks, electronics components, light bulbs, electric motors and stability of pharmaceuticals.
A common practice in reliability is to employ accelerated life tests (ALTs) to shorten the lifetime of a product by increasing some stress factors associated with it, such as temperature, pressure or humidity.This way, the experimental time and cost can be reduced.After suitable inference is developed, we can then extrapolate the results to normal operating conditions; see Meeter and Mecker (1994) and Meeker et al. (1998).There are different types of ALTs resulting in different statistical models.For example, constant-stress ALTs assume that each device is subject to only pre-specified stress levels, while step-stress ALTs apply stress to devices in such a way they will get changed at pre-specified times, and progressive-stress ALT continuously increases the stress level.The constantstress and step-stress ALTs have been widely studied for destructive one-shot devices; see Ling (2019), Lee and Bae (2020), Wu et al. (2020) and Ling and Hu (2020), among others.We focus here on the step-stress model for non-destructive one-shot devices.In particular, we adopt a parametric approach in which the lifetimes of the devices are assumed to follow exponential distribution.
While classical estimation methods are based on the maximum likelihood estimators (MLE), recent works have shown the advantage of using divergence-based methods in terms of robustness, with an unavoidable loss of efficiency in the case of uncontaminated data.Balakrishnan et al. (2019aBalakrishnan et al. ( , 2019bBalakrishnan et al. ( , 2020aBalakrishnan et al. ( , 2020bBalakrishnan et al. ( , 2021) ) developed robust estimation methods based on density power divergence (DPD) for the constant-stress model and then constructed robust test statistics for testing linear hypothesis.In this paper, we develop robust estimators and test statistics for non-destructive one-shot devices under the multiple step-stress model and exponential lifetimes, and illustrate their robustness features both theoretically and empirically.
The rest of the paper is organized as follows: Section 2 describes the multiple step-stress accelerated life test (SSALT) under exponential lifetimes and introduces the classical MLE for the SSALT model and Section 3 presents the minimum DPD and the minimum restricted DPD estimators under linear constraints, with the corresponding asymptotic results.In Section 4 the robustness of the proposed estimators is examined through their influence function analysis.Section 5 describes point estimation and confidence intervals for the reliability, distributional quantile and mean lifetime of the device based on minimum DPD estimators.In Section 6, robust test statistics are developed, including Z-type and Rao-type tests, and their asymptotic properties are examined, providing approximate power functions of the tests.Sections 7 and 8 empirically illustrate the performance of the proposed methods thorough an extensive simulation study and real data analysis respectively.Finally, Section 9 presents some concluding remarks.

Model formulation and the maximum likelihood estimator
Let us consider a SSALT with k ordered stress levels, x 1 < x 2 < • • • < x k , and N one-shot devices under test.At pre-fixed times τ i , called times of stress change, we increase the stress level from x i to x i+1 , for i = 1, ..., k − 1, and we denote τ k the time at which the experiment terminates.Let us also consider a sequence of length L of inspection times during the experiment, including the times of stress change τ i , i = 1, ..., k, where l i denotes the number of inspection times before the i-th stress change and L = l k .Under this set-up, the simple step-stress model corresponds to the case when k = 2, l 1 = 1 and l 2 = 2.This model has been widely studied in the literature; for example, Nelson (1980) discussed general cumulative exposure model, including the simple stress model, while Balakrishnan (2009) reviewed exact inferential procedures for exponential step-stress models.
We further assume that the lifetime, T, of a device follows an exponential distribution, under stress level x i , with failure rate λ i depending on the stress.The distribution of the lifetime of a device during the test is then formed under applying the cumulative exposure model, which relates the lifetime distribution of a device at one stress level to the distributions at preceding stress levels by assuming the residual life of that device depends only on the cumulative exposure it had experienced, with no memory of how this exposure was accumulated.Then, if G i (•) denotes the exponential lifetime distribution function at the i-th stress level, the distribution function of T, G T (t), is given by with for i = 1, ..., k − 1.For notational convenience, we set a −1 = τ −1 = 0.The corresponding density function of T is given by Although the distribution function is continuos in (0, ∞), the density function has k points of discontinuity at times of stress change.We further assume that at stress level x i , the rate parameter λ i of a device has a log-linear relationship with stress level given by where θ = (θ 0 , θ 1 ) ∈ R + × R = Θ is an unknown parameter vector of the model.Note that the mean lifetime of a device is inverse of the exponential parameter, and so it would decrease with an increase in the level.The log-linear relation in ( 4) is frequently assumed in accelerated life test models, as it can be shown to be equivalent to the well-known inverse power law model or the Arrhenius reaction rate model.Suppose n j failures of test devices are observed in the interval (t j−1 , t j ], j = 1, .., k, and for notational ease, let us denote n L+1 for the number of surviving devices at the end of the experiment.Then, the probability of failure of a device in the j-th interval is and the probability of survival at the end of the experiment is π L+1 (θ) = 1 − G T (t L ).Accordingly, a multinomial model with probability vector π(θ) = (π 1 (θ)), ..., π L+1 (θ)) T and N trials can be used to Remark 1 We could have alternatively derived the likelihood function of the model using binomial distribution for each interval, by using conditional probabilities of failure, given that the device did not fail in earlier time intervals.However, both approaches yield the same likelihood function.
Now, let p = (n 1 /N, ..., n L+1 /N ) be the empirical probability vector obtained from the observed data.Then, the Kullback-Leibler divergence between the empirical and theoretical probability vectors, p and π(θ), is given by .
It is straightforward to see that the Kullback-Leibler divergence is related to the log likelihood function in the form where the constant c = L+1 i=1 p i log( p i ) does not depend on θ.Hence, the MLE can equivalently be defined as θ From an asymptotic point of view, it is well-known that the MLE is a BAN (Best Asymptotically Normal) estimator, and it has therefore been widely used for the SSALT model.However, despite its high efficiency, the MLE lacks robustness as contaminated data could influence the parameter estimation considerably.In the next section, we present a robust family of estimators for the SSALT model based on the DPD.

Minimum density power divergence estimator
The density power divergence (DPD) family, introduced by Basu et al (1998), a rich class of densitybased divergences, produces robust estimators with relative small loss in efficiency.Given two density or mass functions, f θ and g, the DPD between them is defined as The parameter β, indexing the DPD divergence, controls the trade-off between efficiency and robustness.In fact, the DPD can be defined at β = 0 by taking continuous limits leading to the Kullback-Leibler divergence.Following the discussion in the last section, we consider the DPD between the empirical and theoretical probability vectors, p and π(θ), and correspondingly define the minimum DPD estimator (MDPPE) as Note that the value β = 0 corresponds to the MLE of θ.Hence, the proposed family could be considered as a generalization of the MLE with a tuning parameter β accounting for the compromise between efficiency and robustness.Moreover, the last term of each addend in (8) does not depend on the model parameter, and so it can be ignored in the minimization process.
The next result presents the estimating equations for the MDPDE.

Result 2
The estimating equations associated with the MDPDE for the SSALT model, under exponential lifetimes, satisfying the log-linear relation in (4), are given by where 0 2 is the 2-dimensional null vector, D π(θ) denotes a (L + 1) × (L + 1) diagonal matrix with diagonal entries π j (θ), j = 1, ..., L + 1, and W is a (L + 1) × 2 matrix with rows w j = z j − z j−1 , where z −1 = z L+1 = 0 and i is the stress level at which the units are tested after the j−th inspection time.
For the MLE, the estimating equations are obtained by deriving the Kullback-Leibler divergence given in (6), yielding Next we present the asymptotic distribution of the proposed estimator, for any positive value of β.
For β = 0, the Fisher information matrix associated with the SSALT model under exponential lifetimes coincides with the matrices J β (θ 0 ) and K β (θ 0 ), and so we obtain the asymptotic distribution of the MLE as a particular case, i.e., where Remark 4 As θ β is a consistent estimator of θ 0 , the asymptotic variances of θ β 0 and θ β 1 for β > 0, denoted by σ 2 (θ β 0 ) and σ 2 (θ β 1 ), respectively, can be estimated by the diagonal entries of ).Therefore, asymptotic confidence intervals for θ 0 and θ 1 , at confidence level (1 − α), are given by with z α/2 being the lower α/2-quantile of a standard normal distribution.Moreover, we have and consequently an associated ellipsoidal confidence region for θ = (θ 0 , θ 1 ) is given by Choosing c = χ 2 2,α , the 100(1 − α)-percentile of the chi-square distribution with 2 degrees of freedom, we have P θ (C α N,β ) tending to 1 − α as n → ∞.Hence, C α n,β represents an ellipsoidal confidence region for θ having limiting confidence coefficient 1 − α as n → ∞ (see Serfling (2009) for more details).The volume of the ellipsoidal region C α N,β , based on Cramer (1946), is given by Thus, a measure of the asymptotic relative efficiency of θ β , for β > 0, with respect to the MLE, θ 0 , is given by

Influence function of the MDPDE
The influence function (IF), first introduced by Hampel et al. ( 1986), plays a central role in the study of robustness properties of an estimator.Intuitively, it quantifies the impact of an infinitesimal perturbation in the true distribution underlying the data on the asymptotic value of the resulting parameter estimate.An estimator is said to be robust if its influence function is bounded.Mathematically, the IF of an estimator is computed in terms of its corresponding statistical functional.Let F θ and G be the assumed distribution of the model and the true density underlying the data, respectively.We use T (G) to denote the statistical functional associated with the estimator θ.Then, the IF of the estimator θ at a point t is computed as where G ε = (1 − ε)G + ε∆ t is the contaminated version of G, with ε being the contamination proportion, and ∆ t being the degenerate distribution at the contamination point t.An estimator is said to be robust if its IF at the model distribution F θ is bounded.For the SSALT model, we could consider only one cell contamination, and so the contamination point t would have all elements equal to zero except for only one component.Let us denote F θ for the assumed distribution of the multinomial model with mass function π(θ) given by the SSALT model with exponential lifetimes and G denote the true distribution underlying the data, with mass function g.We define the statistical functional T β (G) as the minimizer of the DPD between the two mass functions, π θ and g, given in (3).Then, an expression of the IF can be computed from (14) as stated in the following result.

Result 5
The IF of the MDPDE of the SSALT model, θ β , at a point contamination n and the assumed model distribution F θ 0 is given by Remark 6 The matrix J β (θ 0 ) is assumed to be bounded, and so the robustness of the estimators depends on the boundedness of the second factor of the IF, given by where z j is as defined in (10).All the terms in ( 16) are bounded for fixed stress levels and inspection times at any contamination point n, as the maximum value of its components is the number of trials N.Then, any of the proposed MDPDE for β ≥ 0 is robust against vertical outliers, including the MLE.Conversely, the IF boundedness is affected by leverage points, i.e. outlier inspection times or outlier stress levels.
Let us first consider the situation wherein an inspection time, t j , tends to infinity, for fixed j.We denote i the fixed stress level corresponding to the j−th inspection time.As the inspection times are ordered, there will be no more terms in the summation after the j-th term.Then, we can write All the terms depending on times before t j are bounded, and the values λ i (θ 0 ), a i−1 , a * i−1 and τ i−1 are positive constants.Then, taking limits as T j → ∞, we get lim Hence, the IF of the MDPDEs, for positive values of β, is bounded when any inspection time gets increased, whereas the IF of the MLE is unbounded for this class of leverage points.
Similarly, let us consider a stress level x i and let x i → ∞.We take t j such that t j = τ i , the time of stress change for the i-th stress level.Again, as the stress levels are ordered, we can consider that the devices at subsequent steps are subjected to the same stress x i .Then, we need to stablish the boundedness of all terms from j onwards.The lifetime rates are not constant since they depend on the stress level.Therefore, taking limits on (4), we have The limiting behaviour of the IF of θ 0 and θ 1 may be different, since the first only depends on the stress level at g T (T j ), whereas the IF of θ 1 includes a term in g T (T j )T j x i .Therefore, for θ 1 > 0, lim For θ 1 < 0, we must deal with the IF of θ 0 and the IF of θ 1 separately.For the IF of the first parameter θ 0 , we have whereas taking limits in the IF of the second parameter θ 1 , we obtain Thus, the IF of the proposed MDPDE is bounded for all β > 0 regardless of the sign of the true parameter value θ 1 .In contrast, the IF of the MLE of the parameter θ 0 is unbounded for positive true parameter values of θ 1 and so is the IF of the MLE of the parameter θ 1 .This means that the proposed estimators are also robust for all type of outliers, whereas the MLE lacks robustness against this "bad" leverage points.

Point estimation and confidence intervals of reliability and mean lifetime
One may be interested in studying the reliability of non-destructive one-shot devices or in estimating its expected lifetime.Technically, the reliability of a device is the probability that it will perform its intended function, under operating condition, for a specified period of time.Then, the reliability can be measured as the probability of survival until a pre-specified time under normal operating conditions.Following the notation in Section 2, the reliability or survival function of the lifetime T of a device is given by where a i−1 is as defined in (2) and R i (t), i = 1, ...k, is the reliability function at the i − th stress level, which depends in turn on the model parameters θ = (θ 0 , θ 1 ).Therefore, an estimated reliability at a certain time can be obtained from the above formula.For cumulative exposure model with exponential lifetime distributions, the reliability of the device at the inspection interval [τ i , τ i+1 ] (assuming that the stress level will not be increased) corresponds to the reliability function of a translated exponential distribution with parameter λ i , i = 1, ..., k.If the device is subjected to a constant stress level x i , then its reliability at time t can be computed as the reliability function R i (t).Let us denote x 0 for the stress level at normal operating conditions.Then, the reliability of the device is given by, for a fixed time t, R 0 (t) = exp (−λ 0 t) = exp (−θ 0 exp(θ Also, for planning purposes, one may need to estimate the time at which more than a certain percentage of devices are expected to fail under normal operating conditions.Mathematically, those times are the distribution quantiles, computed as the inverse distribution (or reliability) function, with 1 − α being the proportion of surviving units.Further, the mean lifetime of a device under an exponential lifetime distribution with parameter λ is E[T ] = 1/λ.From (4), the mean lifetime of the device depends on the stress level through a log-linear relationship as Hence, under normal operating conditions, the expected lifetime of the device is simply Given the MDPDEs of the model parameters, its is straightforward to obtain point estimate of the reliability at a mission time, estimate quantiles and mean lifetime of the devices under normal operating conditions by substituting the estimated parameters in ( 17)- (19), yielding the estimators , respectively.Additional interest may be on confidence intervals (CI) for such quantities.We first present the asymptotic distribution of the reliability, quantiles and mean lifetime estimators based on the MDPDEs, θ β , under normal operating conditions.These results can be obtained readily from the asymptotic distribution of the MDPDE by employing the Delta method.
Result 7 Let θ 0 be the true value of the parameter θ.Let θ β be the MDPDE, with tuning parameter β.Then, the asymptotic distribution of the estimated reliability at a mission time t, under normal operating conditions, based on the MDPDE θ where the matrices J β (θ 0 ) and K β (θ 0 ) are as defined in (12) and ∇h (θ Result 8 Under the same assumptions as in Result 7, the asymptotic distribution of the estimated (1 − α)-quantile, under normal operating conditions, based on the MDPDE θ where the matrices J β (θ 0 ) and K β (θ 0 ) are as defined in (12) and is the gradient of the function h 1 (θ) = 1 θ 0 exp(−θ 1 x 0 ).
Result 9 Under the same assumptions as in Result 7, the asymptotic distribution of the estimated mean lifetime, under normal operating conditions, based on the MDPDE θ , is given by where the matrices J β (θ 0 ) and K β (θ 0 ) are as defined in ( 12) and ∇h 2 (θ As θ β are consistent estimators, from the above results, we can easily obtain approximate twosided 100(1 − α)% CI for the reliability, (1 − α)-quantile and mean lifetime, under normal operating conditions, to be where σ( R β 0 (t)), σ( Q β 1−α ) and σ( E β T ) are as defined in Results 7-9, respectively.The above asymptotic confidence intervals are based on the asymptotic properties of the estimators and so they may be satisfactory only for large sample sizes.In small samples, we may have to truncate the confidence intervals as the mean lifetime and quantiles must be positive and the reliability should be between 0 and 1.In this regard, Viveros and Balakrishnan (1993) employed a logit transformation of the estimated reliability to obtain more accurate CIs based on the MLE.The transformed reliability is defined as where φ ∈ R. Thus, the range for this transformed reliability would not require truncation.The logit transformation is a natural choice when dealing with parameters that represent probabilities, since it results in R. Estimated values of the transformed reliabilities based on the MDPDEs, φ β , can be easily obtained by substituting the corresponding estimated reliabilities R β 0 (t) in (20), and their asymptotic distribution and CIs can be derived by using Delta method.Inverting such a CI, after some algebra, we obtain the asymptotic CI for the reliability as and σ(R 0 (t)) as defined in Result 7.
A similar idea can be applied for the quantiles and mean lifetimes in (18) and (19).As both quantities must be positive, the logarithm is a natural choice for transforming them to R. Transformed quantiles and mean lifetimes are then Again, using Delta method for deriving the asymptotic distributions, and then inverting the logarithmic transformations, we obtain CIs for Q 1−α and E T as , respectively, with σ(E T ) and σ(Q 1−α ) as defined in Results 8 and 9.

Robust tests of hypotheses
In this section, we consider linear hypothesis tests on the model parameter θ, of the form where m = (m 0 , m 1 ) T ∈ R 2 .In particular, the linear hypothesis with m T = (0, 1) and d = 0 would test if the stress level affects the lifetime of the one-shot devices or not.We present here the a testing procedure based on the MDPDE and then we study it robustness and asymptotic behaviour.Specifically, we define the Z-type statistics based on the MDPDE and then study theoretically its asymptotic distribution under the null and contiguous hypotheses and robustness properties Definition 10 The Z-type statistic based on the MDPDE θ β , for testing null hypothesis (22), is given by The asymptotic distribution of this statistic is given in the following result

Result 11
The asymptotic distribution of the Z-type statistic (23), under the null hypothesis (22), is a standard normal distribution.
Based on Result 11, for any β ≥ 0 and m ∈ R 2 , the critical region with significance level α for the hypothesis test with linear null hypothesis in (22), with level α is given by where z α/2 denotes the upper α/2-quantile of the standard normal distribution.

Remark 12
We can generalize the null hypothesis in (22) to with M being a r × 2 (r ≤ 2) matrix and d being a r−dimensional vector.Then, we can define the corresponding test statistic as Note that the matrix )M is symmetric, and so the statistic is well defined.It is not difficult to stablish that, under the generalized null hypothesis above, and therefore, Consequently, a critical region corresponding to the generalized null hypothesis is where χ 2 r,α denotes the upper α-quantile of a chi-square distribution with r degrees of freedom.
The robustness of the proposed Z-type test statistic can be established by its IF.The IF of a testing procedure at a contamination point n is defined as the Gateaux derivative of the functional, defining the test statistic at the contamination direction given by ∆ n .In the present context, the functional associated with the proposed Z-type test statistic, Z N ( θ β ), under the null hypothesis is given by Therefore, the IF of the proposed Z-type test statistic can be easily derived from the IF of the MDPDE, as The boundedness of the IF of the Z-type test statistic at a contamination point n and the true distribution F θ 0 can be discussed by the boundedness of the IF of the corresponding MDPDE, and thus, robust estimators results in robust test statistics.
On the other hand, we can obtain the asymptotic distribution of the Z-type test in (23) at a contiguous alternative hypothesis.Let θ L ∈ Θ \ Θ 0 be an alternative and take θ 0 as the closest element to the boundary of Θ 0 in the sense of Euclidean distance.We consider contiguous alternative hypothesis of the form with θ L = θ 0 + 1 √ N , for a fixed vector ∈ R 2 .Note that, defining * = m T , we have so that the contiguous hypothesis in ( 26) can be equivalently stated by the condition g(θ L ) = 1 √ N * .

Result 13
The asymptotic distribution of the Z-type statistic in (23), under the contiguous hypothesis (26), is a normal distribution, with mean (m From the above result, we can obtain an approximation for the power function of the test statistic in (22) at the contiguous hypothesis in (26), as It is clear that lim N →∞ β N (θ L ) = 1 and so the Z-type statistic is consistent in the sense of Fraser (1957).More generally, the following result provides an asymptotic approximation to the power function.
Result 14 Let θ * ∈ Θ be the true value of the parameter θ with m T θ * = d.Then, the approximate power function of the test statistic in ( 22) is given by where Φ(•) denotes the standard normal distribution function.

Simulation study
In this section, we examine the behaviour of the proposed robust MDPDEs, Z-type tests and Raotype tests.for the SSALT model with exponential lifetime distribution under different contamination scenarios.
For multinomial sampling, we must consider "outlying cells" instead of "outlying devices", see Balakrishnan et al. (2019a).Then, to introduce contamination in our context, we should increase (or decrease) the probability of failure in ( 5) for (at least) one interval (i.e., one cell).So, the probability of failure is switched in such contaminated cells as for some j = 2, ..., L, where θ = ( θ 0 , θ 1 ) is a contaminated parameter with θ0 ≤ θ 0 and θ1 ≤ θ 1 .
It is important to point out that, after the contamination of the probability of failure in a cell, the probability vector of the multinomial model must be normalized to add up to 1.

Minimum density power divergence estimators
Let us consider a 2-step stress ALT experiment with L = 11 inspection times and a total of N = 180 devices under test.At the beginning of the experiment, all the devices are subjected to a stress level x 1 = 35 until the first time of stress change τ 1 = 25.Then, the surviving units are subjected to an increased stress level, x 2 = 45, till the end of the experiment at τ 2 = 70.During the experiment, inspection is performed at a grid of inspection times containing the times of stress change, IT = (10,15,20,25,30,35,40,45,50,60,70).We set the true value of the true parameter θ 0 = (0.003, 0.03), and then generate data from the corresponding multinomial model described in Section 2, with exponential lifetimes.Moreover, we contaminate the data by increasing the probability of failure in the third interval as mentioned in (27).
In order to evaluate the performance of the proposed estimators, we calculate the root mean square error (RMSE) of the MDPDE for different values of β ∈ {0, 0.2, 0.4, 0.6, 0.8, 1}, including the MLE for β = 0. Further, to asses the efficiency loss of an estimator with respect to the MLE, we define a measure ρ( θ Then, ρ( θ β ) measures the efficiency loss of an estimator with respect to the MLE.Clearly, when ρ( θ β ) < 0, the MDPDE is more accurate than the MLE, and evidently ρ( In addition, we test different scenarios of contamination.In the first scenario, we generate an outlying cell in the third interval by decreasing the value of the first parameter, θ 0 , and in the second scenario we perform similarly, but decreasing the second parameter θ 1 .In both cases, the lifetime rate λ( θ), is decreased; the smaller is the contamination parameter, the greater is the contamination.
Figures 1 and 2 show the RMSE and the RMSE ratio, ρ, produced with different values of β and the two contamination scenarios determined from R = 1000 replications.In the left side plots, the contamination is introduced by decreasing θ 0 , yielding a contamination rate of ε = 1 − θ 0 θ 0 while on the right side plots they are computed by decreasing the second parameter θ 1 , and the corresponding contamination rate is then ε = 1 − θ 1 θ 1 .These results show the advantage of the proposed MDPDE in terms of robustness.The larger is the parameter β, the more robust is the corresponding estimator.As expected, in the absence of contamination, the MLE is the most efficient estimator, even though all proposed MDPDEs perform competitively in this uncontaminated scenario.On the other hand, when a "great" outlier cell is generated, the efficiency loss of MLE with respect to the proposed MDPDE is seen to be quite pronounced.More specifically, under the contamination rates greater than 20%, the MDPDEs outperform the MLE.

Z-type tests
We empirically examine the performance of the Z-type test statistic based on the MDPDE.We adopt again the 2-step stress ALT experiment with L = 11 inspection times and N = 180 devices described in Section 7.1, and we consider testing the hypothesis The true value of the parameter is set to be θ 0 = (0.003, 0.03) T so as to fit the null hypothesis.Then, the Z-type test statistic is defined using (23) with m = (0, 1) T and d = 0.03, and the critical region of the test is given by (24). Figure 3 shows the empirical level of the test against cell contamination, ε, with the two different contaminated scenarios considered in the last section, θ 0 -contaminated third cell (left) and θ 1 -contaminated third cell (right).The empirical level is computed as the proportion of rejected Z-type test statistic over R = 1000 replications of the model under the null hypothesis for a significance level of α = 0.05.The empirical level of the Z-type tests based on the MLE promptly grows when the contamination rate gets increased, whereas the empirical level of the Z-type test based on the MDPDE, with large values of β, remains low in heavily contaminated scenarios, highlighting its robustness property.Next, we examine the empirical power of the linear hypothesis test by considering a different true parameter value, θ = (0.003, 0.6) T , and the hypothesis test H 0 : θ 1 = 0.03 against the alternative H 1 : θ 1 = 0.03.Now, the true value of the parameter does not satisfy the null hypothesis, and the Z-type test statistic has its the empirical power under the two different contaminated scenarios as displayed in Figure 4, with R = 1000 replications.Again, the robustness of the Z-type test based on MDPDE becomes better when large values of β are used to construct the test statistics.

Choice of the tuning parameter
The tuning parameter β of the DPD loss function controls the trade-off between efficiency and robustness of the resulting MDPPE.Following the discussions in the preceding sections, larger values of β produce more robust but less efficient estimators.Therefore, the optimal value of β will be good to determine.From our empirical results, a moderately large value of β (over β = 0.4) is expected to provide robust estimators without a high loss of efficiency with respect to the MLE in the absence of contamination.Determining this optimal value for the best compromise is, therefore, of great practical interest.Optimal values of β will produce robust estimators without coming at the cost of a high efficiency loss.Then, a criterion measuring the efficiency loss in favour of robustness gain should be adopted.Warwick and Jones (2005) introduced an useful data-based procedure for the choice of the tuning parameter for the MDPDE.However, this method depends on the choice of a pilot estimator and Basak et al. (2021) improved the method by removing the dependency on an initial estimator.The approach of Warwick and Jones (2005) minimizes the asymptotic MSE of the MDPDE given by where θ P is a pilot estimator and Tr denote the trace of the matrix.Several proposals of this pilot estimator have been studied in the literature.However, the choice of the pilot significantly impact on the optimal tuning parameter, as it invariably draws the final estimator towards itself.To overcome this drawback, Basak et al. (2021) proposed an iterative algorithm that replaces in each step, the value of the pilot estimator by the estimator obtained with the optimal value of β until the optimal choice of the tuning parameter (or equivalently, the pilot estimator) gets stabilized.The process should be initialized with a suitable robust pilot estimator, but the final choice of β gets more pilot-independent.Basak et al. (2021) empirically showed that when the pilot estimators are within the MDPDE class, all robust pilots lead to the same iterated optimal choice, and moreover the performance of the algorithm improves even with pure data.These are all summarized in the following algorithm: Algorithm [Choice of the tuning parameter] 1. Fix the convergence rate ε and choose an initial pilot estimator θ P from the MDPDE family; 2. Update the optimal value of the tuning parameter, β * , using the minimum asymptotic MSE in (29); Let us consider the previous simulation set up with true parameter value θ = (0.003, 0.03) T , but now let us choose the optimal value of β according to the presented data-based procedure detailed in the above algorithm.We initialize the method with the MDPDE with different tuning parameters β p = 0, 0.5 and 1, yielding the pilot estimators θ 0 , θ 0.5 and θ

1
, respectively, and we fix the convergence rate to be ε = 0.001.The minimization of ( 29) is carried out over a grid search in [0, 1] of size 100.
Figure 6 shows the optimal values of the tuning parameter β against data contamination over R = 1000 repetitions.As expected, optimal values of β are greater with high contamination rates.Further, the choice of optimal β is almost entirely independent of the pilot estimator, and so the presented algorithm does not seem to get affected by this initial choice.Optimal values of the tuning parameter are generally larger when the contamination is introduced on the first parameter.Then, the model is more sensitive to contamination in such direction.Moreover, in Figure 7, the RMSE of the resulting (optimal) estimator is compared to the RMSE of the estimators with fixed values of β ∈ {0, 0.2, 0.4, 0.6, 0.8, 1} based on R = 1000 repetitions.As expected, the data-based method for choosing optimal β outperforms any of the methods with a pre-fixed value of β, since it adapts the tuning parameter value to the amount of contamination present in the data.However, it tends to be slightly conservative and selects insufficiently high values of β when a high contamination rate is introduced.In this section, we discuss two real-life applications of the MDPDE for the SSALT model, developed in the preceding sections.

Electronic components data
The first example was studied by Wang and Fei (2003) to get the reliability indices of a kind of electronic components at the normal temperature of x 0 = 25 • C. N = 100 items from a batch of products were randomly selected for a simple SSALT (k = 2), with two stress levels x 1 = 100 • C and x 2 = 150 • C. In the original experiment, the stress level rises when 30 products had failed and the test continues until 20 more products had failed, obtaining a total of 50 failures.Their failure times are as follows: • Failure times at the first stress level   For illustrating the performance of the MPDPE for the SSALT model with one-shot devices, we will assume that we only know how many devices had failed before certain pre-specified inspection times, t = 270, 430, 600, 910, 975, 1015, 1040, 1096.Additionally, the time of stress change, τ 1 , is pre-fixed at t = 910.
Table 1 shows the estimated model parameters with different values of the tuning parameter β, along with an approximate CI constructed from (13).The last row of the table contains the estimates with the optimal value of β obtained with the algorithm presented in Section 7.3.The algorithm is initialized with the pilot estimate θ 0.5 and the optimum value of β is reached at β = 0.027, implying moderately low contamination in the data.Further, to fairly analyze the role of each of the parameters in the model, we present the logarithm of the first parameter, θ 0 , so both model parameters are reported on the same scale.The mean lifetime of the electronic component at low temperatures is appreciably high, and increasing a degree on the temperature multiplies the lifetime of the device, according to all estimators, in approximately 0.97 times.Consequently, the increase of the temperature from x 1 = 100 • to x 2 = 150 • shortens the lifetime by more than 0.22 times.Robust methods tend to estimate with a higher value the first parameter, and consequently decrease the estimate of the second one.Conversely, the variance of the estimator increases with β, producing wider intervals.
Table 2 shows the mean lifetime (in hours) of the electronic components under different (constant) stress level, and their associated direct and transformed CI.Under normal operating conditions (x 0 = 25), the devices are expected to last for more than 6 hours, but their lifetime gets severely decreased when exposed to very high temperatures, which allows to infer about the reliability in a short period of time.It is interesting to note that direct CI of the mean lifetime under normal operating conditions gets truncated due to the positivity constraint.On the other hand, transformed CI is wider, and the right end point is quite far from than the the one obtained with direct CIs.This difference gets reduced at increased temperatures.
On the other hand, one may be interested in estimating the reliability of the devices when it is exposed to different constant temperatures.We fix a "mission time" at t = 600s and we report the estimated reliabilities and CIs under different stress levels in Table 3. Again, direct CIs had to be truncated under normal operating conditions so as to remain within the interval (0, 1).As expected, the reliability of the devices decreases when increasing the stress level, and here all estimates remain close for all values of the tuning parameter.Finally, one may be interested in determining the time at which 10% of the devices are expected to fail, under different (constant) temperatures.Table 4 presents the estimated 0.9−quantiles of the lifetime distribution (or equivalently the 0.1−quantiles of the reliability).Here, the direct CI of the estimated quantiles under normal operating conditions are again truncated, demonstrating again the drawback of the direct method, while transformed CI provides a good alternative for such intervals without the problem of constraints.

Light bulbs data
Zhu (2010) conducted an accelerated life testing experiment in the Quality and Reliability Engineering Laboratory of the Industrial and Systems Engineering Department of Rutgers University so as to examine the reliability of light bulbs.Two sets of 32 miniature light bulbs were placed in a temperature and humidity chamber where humidity was held constant, and the long term failure due to bulb filament fatigue was then studied.When the switch was turned on, full current suddenly flowed to the filament at the speed of light.This sudden massive vibration caused the filament to wildly bounce causing fatigue behaviour of the filament which resulted in breakage of the filament.Long Term Failure occurred when the filament eventually become so fatigued that its electrical resistance increased to the point that current would not flow.Each light bulb was connected with a resistor, across which Voltage was measured to monitor the status of the light bulbs.Normal operating condition of the light bulbs is 2V.To carry out the SSALT, they applied 2.25V for 96hr and then increased the voltage to 2.44V.The step-voltage test got stopped at 140hr.Failure times during the experiment are as follows: The remaining 11 light bulbs continued to provide light when the experiment was terminated.To illustrate the performance of the MDPDE of the SSALT model, we transformed the collected data into one-shot devices data with inspection times t = 25, 50, 96, 110, 120, 140.Table 5 shows the estimated values of the model parameters with different values of the tuning parameter β.Applying the databased choice of β described in Section 7.3, the optimum value is approximately β = 0.12, thus showing slightly higher contamination to be present in this data than in the last electronic component example.Results for this optimum value are presented in the last row of Table 5.As in the previous example, we report the values of log( θ 0 ) so both model parameters are in same scale.Unlike in the previous Appendix: Proofs of the main results

Proof of Result 2
Proof.The MDPDE is defined as the minimizer of the DPD-based loss in (8), and so must satisfy: Next, the derivative of the probability of success π j (θ) depends on the stress level at which the device is being tested.We denote x i for the stress level at which the units are tested after the τ i −th inspection time.Taking derivatives in (5), we get ∂π j (θ) ∂θ = ∂G T (t j ) ∂θ − G T (t j−1 ) ∂θ .

Figure 1 :
Figure 1: RMSE of different estimators against data contamination in R = 1000 replications

Figure 2 :
Figure 2: Efficiency loss with respect to the MLE against data contamination in R = 1000 replications θ 1 -contaminated cell

3 .
If the optimal estimate θ β * differs from the pilot estimator by less than the convergence rate: Stop; Else, replace the pilot estimator by the optimal θ β * and return to Step 2.

Table 6 :
Estimated mean lifetime (in minutes) and asymptotic confidence intervals for the light bulbs under different voltages.

Table 8 :
Estimated 90% quantile and asymptotic (direct and transformed) confidence intervals (in seconds) of the light bulbs under three constant stress levels.