# A Robust Bayesian approach for reliability prognosis of one-shot devices
under cumulative risk model

Shanya Baghel  Shuvashree Mondal11footnotemark: 1<EMAIL_ADDRESS>
###### Abstract

Reliability prognosis of one-shot devices is drawing increasing attention
because of its wide applicability. The present study aims to determine the
lifetime prognosis of highly durable one-shot device units under step-stress
accelerated life testing (SSALT) experiment applying a cumulative risk model
(CRM). In an SSALT experiment, CRM retains the continuity of hazard function
by allowing the lag period before effects of stress change emerges. In an
analysis of such lifetime data, plentiful datasets might have outliers where
conventional methods like maximum likelihood estimation or likelihood-based
Bayesian estimation frequently fail. This work develops a robust estimation
method based on density power divergence in classical and Bayesian frameworks.
The hypothesis is tested by implementing the Bayes factor based on a
robustified posterior. In Bayesian estimation, we exploit Hamiltonian Monte
Carlo, which has certain advantages over the conventional Metropolis-Hastings
algorithms. Further, the influence functions are examined to evaluate the
robust behaviour of the estimators and the Bayes factor. Finally, the
analytical development is validated through a simulation study and a real data
analysis.

###### keywords:

Bayes factor , cumulative risk model , density power divergence , influence
function , one-shot device , robust Bayes estimation

[inst1]organization=Department of Mathematics and Computing, Indian Institute
of Technology (ISM) Dhanbad,postcode=826004, state=Jharkhand, country=India

label1label1footnotetext: Tel number: (0326)2235458

Step-stress cumulative risk model

Robust estimation in classical and Bayesian frameworks

Testing of hypothesis based on robust Bayes factor

Influence function analysis

##  1 Introduction

The one-shot device is a frequently employed term to describe a system or
component that perishes upon accomplishing its intended task. Classic examples
that satisfy assumptions of one-shot devices are automobile airbags, fuses,
electro explosives, pharmaceutical stability, safety valves, light bulbs and
many such examples. For testing of such devices, the observation is mostly
restricted to recording if device failure occurs before or after a specified
inspection time, which leads to the study of dichotomous data only. For highly
durable products, typically accelerated life and degradation testing (ALT and
ADT) are executed for reliability analysis within limited resources. Such
testing employs higher-than-use-condition stress for the analysis, and results
are then extrapolated to the normal-use conditions. There are plenty of
studies devoted to ALT and ADT [1, 2, 3, 4, 5, 6]. The present study
incorporates the reliability prognosis under ALT and it is found in the
literature that several authors indulged in the analysis of one-shot device
testing under constant stress ALT [7, 8, 9, 10, 11]. However, step-stress ALT
(SSALT) for one-shot device testing has garnered the attention of researchers
[12, 13, 14, 15] recently as it can produce a large number of failures in a
short span of time. In SSALT, stress increases step-wise over prefixed time
points.

For SSALT, a connection model is needed to relate lifetime distributions at
different levels. The cumulative exposure model (CEM)[12, 13, 14, 15, 16] is
widely used for SSALT. However, a change in stress level in this model is
instantaneous, leading to discontinuity in hazard function at the stress
change point. To address this shortcoming, Van Dorp and Mazzuchi [17] proposed
a model based on hazard rate function, which was subsequently referred to as
cumulative risk model (CRM) by Kannan et al.[18]. This model removes
discontinuity by allowing the lag period before effects of stress change
emerges. Although various authors have studied CRM in the past [19, 20, 21,
22], the SSALT experiment for one-shot device data under CRM is yet to be
explored. The present study focuses towards this direction where failure of a
one-shot device is subject to SSALT experiment under CRM with interval
monitoring over intermediate inspection time points. The lifetime of one-shot
device is assumed to follow the well-known standard family of Lehman
distributions.

When data comes with potential outliers, conventional maximum likelihood
estimation (MLE) may fail to provide robust results as it lacks stability.
Basu et al. [23] proposed a robust estimation method based on density power
divergence (DPD) measure to tackle this. It was further employed in robust
inferential study of one-shot device testing data by several authors [24, 25,
26, 27, 28]. The present study incorporates a DPD-based estimation procedure
for inference. Moreover, the tuning parameter in DPD measure plays a pivotal
role as it balances robustness and efficiency. Castilla and Chocano [29]
focused on the choice of optimal tuning parameter specifically for one-shot
device testing. In this study, we have suggested an approach that does not
require any iterative computational procedure.

The Bayesian approach comes into picture with the availability of prior
knowledge about data. The conventional Bayes estimation based on likelihood-
based posterior is still quite popular [30, 31, 32, 33] but lacks robustness.
Consequently, several authors divulge into the development of robust Bayesian
inference [34, 35, 36]. The present work incorporates robust Bayes estimation
developed by Ghosh and Basu [37] where DPD measure has substituted likelihood
in the posterior density function. Another intriguing feature of the Bayesian
framework is that testing of hypothesis can be processed by using the Bayes
factor introduced by Jeffreys [38] and extensively exploited by many authors
later on.

This work presents a robust estimation both in classical and Bayesian
frameworks where the estimation procedure relies on density power divergence
(DPD) measure between empirical and assumed lifetime distribution. The
asymptotic distribution of DPD estimates is derived and optimization of
optimal tuning parameter is designed. Another essential contribution of this
work is to develop the robust testing of the hypothesis by applying a Bayes
factor derived from the robustified posterior. Further, the influence
functions are derived and examined thoroughly to assess robust behaviour of
the point estimators in both classical and Bayesian frameworks.

In prior selection of the Bayesian framework, Normal and Dirichlet prior are
considered. Under both prior assumptions for the considered model, a closed
form of posterior cannot be obtained. For estimation purposes, one can rely on
frequently used Markov Chain Monte Carlo methods like Gibbs sampler and
Metropolis-Hastings algorithms [39, 40]. However, they behave inefficiently in
exploring the target distribution in the presence of high dimensional or
highly correlated variables [41]. To deal with it, the Hamiltonian Monte Carlo
(HMC) approach can be implemented to generate posterior samples. The HMC
algorithm is helpful for complex models as it delivers more accurate findings
and flexible settings [42, 43]. HMC was first brought by Neal [44] to the
applications of statistics. For an in-depth explanation of HMC, one can refer
to Neal et al. [45], Monnahan et al. [46] and the references therein.

The rest of the article proceeds as follows. Section 2 focuses on building a
cumulative risk model with computation of maximum likelihood estimate. Robust
estimation methods are discussed in Sections 3 and 4. In Section 5 testing of
hypothesis based on robust Bayes factor is developed. Section 6 studies the
property of robustness through influence function. Next, Section 7 discusses
the model with two special cases of the Lehman family of distributions.
Finally, Sections 8 and 9 contain the simulation study and data analysis
respectively. Concluding remarks are given in the Section 10.

##  2 Cumulative risk step-stress model

This section discusses experimental setup for analysing one-shot device
testing data under the cumulative risk step-stress model.

###  2.1 Model set-up

Consider a sample of n one-shot devices put to step-stress accelerated life
testing (SSALT) experiment with $k$ stress levels denoted by
$x_{i}\,;\,i=1,2,\dots,k.$ Let us denote $\tau_{i}$ as pre-specified time
point where the stress level changes from $x_{i}$ to $x_{i+1}$ and
$\tau_{0}=0\,;\,\tau_{k}=\infty.$ At each stress level, lifetime distribution
of one-shot device is assumed to follow the Lehman family of distributions
with different shape and scale parameters whose cumulative distribution
function (cdf) and probability density function (pdf) are defined as

$\displaystyle F_{i}(t)$
$\displaystyle=1-\exp\big{\\{}-\lambda_{i}Q(t;\gamma_{i})\big{\\}}.$ (1)
$\displaystyle f_{i}(t)$
$\displaystyle=\lambda_{i}Q^{{}^{\prime}}(t;\gamma_{i})\exp\big{\\{}-\lambda_{i}Q(t;\gamma_{i})\big{\\}}.$

where $\lambda_{i}>0$ and $\gamma_{i}>0$ are shape or scale parameters and
$Q(t;\gamma_{i})$ is strictly increasing function of $t$. Various lifetime
distributions such as exponential, Weibull and Gompertz are members of this
family for $Q(t;\gamma_{i})=t,t^{\gamma_{i}},(e^{\gamma_{i}t}-1)$
respectively. Here, in equation (1), parameter $\lambda_{i}$ is related to
stress factor in log-linear form as

$\lambda_{i}=\exp(c_{0}+c_{1}x_{i})\;;\;i=1,2,\dots,k.$ (2)

Thus the set of model parameters to be estimated is denoted by

$\bm{\theta}=\\{(c_{j},\gamma_{i})^{T}\,;\,j=0,1\,;\,i=1,2,\dots,k\\}.$

###  2.2 Cumulative risk model

Introduced by Van Dorp and Mazzuchi [17] and Kannan et al.[18], the cumulative
risk model (CRM) allows for an induction of lag period $\delta$ before the
stress change effects emerge. In CRM, the hazard function is linearly modelled
in the interval $(\tau_{i},\tau_{i}+\delta)$. Thus, piece-wise hazard rate
function for Lehman family of distribution takes the following form

$h(t)=\begin{cases}\lambda_{1}Q^{{}^{\prime}}(t;\gamma_{1})&;\;0<t\leq\tau_{1}.\\\
a_{i-1}+b_{i-1}t&;\;\tau_{i-1}<t\leq\tau_{i-1}+\delta\;,\;i=2,3,\dots,k-1.\\\
\lambda_{i}Q^{{}^{\prime}}(t;\gamma_{i})&;\;\tau_{i-1}+\delta<t\leq\tau_{i}\;;\;i=2,3,\dots,k-1.\\\
\lambda_{k}Q^{{}^{\prime}}(t;\gamma_{k})&;\;\tau_{k-1}<t<\infty.\end{cases}$
(3)

To ensure continuity of $h(t)$ in equation (3), $a_{i-1}$ and $b_{i-1}$ must
satisfy

$\displaystyle a_{i-1}+b_{i-1}\tau_{i-1}$
$\displaystyle=\lambda_{i-1}Q^{{}^{\prime}}(\tau_{i-1};\gamma_{i-1}).$ (4)
$\displaystyle a_{i-1}+b_{i-1}(\tau_{i-1}+\delta)$
$\displaystyle=\lambda_{i}Q^{{}^{\prime}}(\tau_{i-1}+\delta;\gamma_{i-1}).$

By solving equation (4), we obtain

$\displaystyle a_{i-1}$
$\displaystyle=\frac{1}{\delta}\bigg{\\{}(\delta+\tau_{i-1})\lambda_{i-1}Q^{{}^{\prime}}(\tau_{i-1};\gamma_{i-1})-\tau_{i-1}\lambda_{i}Q^{{}^{\prime}}(\tau_{i-1}+\delta;\gamma_{i})\bigg{\\}}.$
$\displaystyle b_{i-1}$
$\displaystyle=\frac{1}{\delta}\bigg{\\{}\lambda_{i}Q^{{}^{\prime}}(\tau_{i-1}+\delta;\gamma_{i})-\lambda_{i-1}Q^{{}^{\prime}}(\tau_{i-1};\gamma_{i-1})\bigg{\\}}.$

Therefore, Survival function $S(t)=e^{-\int_{0}^{t}h(x)dx}$ is obtained as

$S(t)=\begin{cases}\exp\Big{\\{}-\lambda_{1}Q(t;\gamma_{1})\Big{\\}}\quad;\;0<t\leq\tau_{1}.\\\
\exp\Big{\\{}-D^{(\delta)}(t;\gamma_{i-1,i})\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{i-1}Q(\tau_{i-1};\gamma_{i-1})+\sum\limits_{l=1}^{i-2}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{i-1}<t\leq\tau_{i-1}+\delta\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{i}Q(t;\gamma_{i})\Big{\\}}\exp\bigg{\\{}-\sum\limits_{l=1}^{i-1}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\bigg{\\}}\\\
\hskip 110.96556pt;\;\tau_{i-1}+\delta<t\leq\tau_{i}\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{k}Q(t;\gamma_{k})\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{k-1}Q(\tau_{k-1};\gamma_{k-1})+\sum\limits_{i=1}^{k-2}E^{(\delta)}(\tau_{i};\gamma_{i+1,i})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{k-1}<t<\infty.\end{cases},$ (5)

where,

$\displaystyle D^{(\delta)}(t;\gamma_{i-1,i})$
$\displaystyle=\frac{(t-\tau_{i-1})^{2}}{2\delta}\Bigg{[}\Big{\\{}2\delta(t-\tau_{i-1})^{-1}-1\Big{\\}}\lambda_{i-1}Q^{{}^{\prime}}(\tau_{i-1};\gamma_{i-1})+\lambda_{i}Q^{{}^{\prime}}(t;\gamma_{i})\Bigg{]}.$
$\displaystyle E^{(\delta)}(\tau_{l};\gamma_{l+1,l})$
$\displaystyle=\lambda_{l}Q(\tau_{l};\gamma_{l})-\lambda_{l+1}Q(\tau_{l}+\delta;\gamma_{l+1})$
$\displaystyle\qquad+\frac{\delta}{2}\Big{\\{}\lambda_{l}Q^{{}^{\prime}}(\tau_{l};\gamma_{l})+\lambda_{l+1}Q^{{}^{\prime}}(\tau_{l}+\delta;\gamma_{l+1})\Big{\\}}.$

###  2.3 CRM under SSALT with interval monitoring

As considered earlier, n one-shot devices are exposed to CRM SSALT experiment
inspected at pre-fixed time points with termination of the experiment at
$\tau_{k}$. Let $q_{i}$ be the number of inspection time points at stress
level $x_{i}$ and $\tau_{i,m}$ be $mth$ inspection time point at $ith$ stress
level with
$\tau_{i,q_{i}}=\tau_{i}\,;\,i=1,2,\dots,k\,;\,m=1,2,\dots,q_{i}\,;\,\tau_{0}=0$.
Figure (1) depicts the layout of CRM SSALT experiment with interval monitoring
and intermediate inspection time points (IMIIP).

Figure 1: SSALT under CRM with IMIIP

Let us denote $n_{im}\,;\,i=1,2,\dots,k\,;\,m=1,2,\dots,q_{i}$ as number of
observed failures in the interval $(\tau_{i(m-1)},\tau_{im}]$. Then,
$n_{i}=\sum_{i=1}^{q_{i}}n_{im}$ is total number of failures at $ith$ stress
level and the total number of observed failures is thus given by
$n_{f}=\sum_{i=1}^{k}n_{i}$. Hence, $n_{s}=n-n_{f}$ is the number of survived
units after time point $\tau_{k}.$ If $T$ is lifetime of a one-shot device,
then failure and survival probabilities using equation (5) are given as

$\displaystyle p_{i1}$
$\displaystyle=P(\tau_{i-1}<T\leq\tau_{i,1})=\int_{\tau_{i-1}}^{\tau_{i-1}+\delta}f(x)\,dx+\int_{\tau_{i-1}+\delta}^{\tau_{i,1}}f(x)\,dx.$
$\displaystyle=G^{(\delta)}(\tau_{i-1}+\delta;\gamma_{i-1,i})+\exp\left\\{-\sum_{l=1}^{i-1}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\right\\}G^{(1)}(\tau_{i-1,i};\gamma_{i}).$
(6) $\displaystyle p_{im}$
$\displaystyle=P(\tau_{i,m-1}<T\leq\tau_{i,m})\;;\;m=2,3,\dots,q_{i}.$
$\displaystyle=\exp\left\\{-\sum_{l=1}^{i-1}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\right\\}G^{(m)}(\tau_{i};\gamma_{i}).$
(7) $\displaystyle p_{s}$
$\displaystyle=P(T>\tau_{k})=\exp\left\\{-\sum_{i=1}^{k-1}E^{(\delta)}(\tau_{i};\gamma_{i+1,i})\right\\}\exp\Big{\\{}-\lambda_{k}Q(\tau_{k};\gamma_{k})\Big{\\}},$
(8)

where,

$\displaystyle G^{(\delta)}(\tau_{i-1}+\delta;\gamma_{i-1,i})$
$\displaystyle=\exp\bigg{[}-\Big{\\{}\lambda_{i-1}Q(\tau_{i-1};\gamma_{i-1})+\sum\limits_{l=1}^{i-2}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\Big{\\}}\bigg{]}$
$\displaystyle\quad\qquad\bigg{[}1-\exp\Big{\\{}-D^{(\delta)}(\tau_{i-1}+\delta;\gamma_{i-1,i})\Big{\\}}\bigg{]}.$
$\displaystyle G^{(m)}(\tau_{i};\gamma_{i})$
$\displaystyle=\exp\Big{\\{}-\lambda_{i}Q(\tau_{i,m-1};\gamma_{i})\Big{\\}}-\exp\Big{\\{}-\lambda_{i}Q(\tau_{i,m};\gamma_{i})\Big{\\}}.$
$\displaystyle G^{(1)}(\tau_{i-1,i};\gamma_{i})$
$\displaystyle=\exp\Big{\\{}-\lambda_{i}Q(\tau_{i-1}+\delta;\gamma_{i})\Big{\\}}-\exp\Big{\\{}-\lambda_{i}Q(\tau_{i,1};\gamma_{i})\Big{\\}}.$

The likelihood function based on observed failure count data is given by

$L(\bm{\theta})\propto\left\\{\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}(p_{im})^{n_{im}}\right\\}\bigg{\\{}p_{s}^{n_{s}}\bigg{\\}}.$

The log-likelihood function is obtained as

$ln\,L(\bm{\theta})\propto\left(\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}n_{im}\,ln\,p_{im}\right)+\bigg{(}n_{s}\,ln\,p_{s}\bigg{)}.$

Therefore, maximum likelihood estimate (MLE) can be obtained as

$\hat{\bm{\theta}}=arg\mathop{max}_{\bm{\theta}}\,ln\,L(\bm{\theta})\;;\;\sum\limits_{i=1}^{k}\sum\limits_{m=1}^{q_{i}}n_{im}>0.$

######  Result 1.

The system of likelihood equations is derived as

$\sum_{i=1}^{k}\left[\frac{n_{i1}p^{{}^{\prime}(\bm{\theta})}_{i1}}{p_{i1}}+\sum_{m=2}^{q_{i}}\frac{n_{im}A^{(i)}_{1}(\bm{\theta})}{G^{(m)}(\tau_{i};\gamma_{i})}\right]=n_{s}A^{(k)}_{2}(\bm{\theta}),$

where description of expressions is given in the appendix A.1.

The presence of outliers in data makes it challenging for MLE to produce a
valid estimate. A robust estimation method helps improve the estimation
procedure in such situations. The following section confides in studying power
divergence-based estimation method.

##  3 Density power divergence method of estimation

Developed by Basu et al.[23] with the assumptions stated in their study,
density power divergence measure (DPD) between any two probability densities,
say, $f_{1}$ and $f_{2}$ are defined as

$D_{\alpha}(f_{2},f_{1})=\int\left\\{f^{1+\alpha}_{1}(t)-\left(1+\frac{1}{\alpha}\right)f_{2}(t)f^{\alpha}_{1}(t)+\left(\frac{1}{\alpha}\right)g^{1+\alpha}(t)\right\\}\,dt\;;\;\alpha>0,$

where $\alpha$ is termed as the tuning parameter. If $\alpha\to 0$, DPD
measure tends to Kullback-Leibler divergence measure.

Here, the DPD measure is computed using empirical and theoretical failure and
survival probabilities for a one-shot device testing unit. For CRM under SSALT
with IMIIP, the empirical failure and survival probabilities are given as
$\left(\frac{n_{im}}{n},\frac{n_{s}}{n}\right)$ where
$\,;\,i=1,2,\dots,k\,;\,m=1,2,\dots,q_{i}.$ Then, DPD measure is described as

$\displaystyle
D_{\alpha}(\bm{\theta})=\left\\{p^{\alpha+1}_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}p^{\alpha+1}_{im}\right\\}$
$\displaystyle-\left(1+\frac{1}{\alpha}\right)\left\\{\frac{n_{s}}{n}p^{\alpha}_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\frac{n_{im}}{n}p_{im}^{\alpha}\right\\}$
$\displaystyle\frac{1}{\alpha}\left\\{\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left(\frac{n_{s}}{n}\right)^{\alpha+1}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left(\frac{n_{im}}{n}\right)^{\alpha+1}\right\\}.$
(9)

The minimum DPD estimator (MDPDE) can be obtained as

$\hat{\bm{\theta}}_{\alpha}=arg\mathop{min}_{\bm{\theta}}D_{\alpha}(\bm{\theta}).$
(10)

The tuning parameter $\alpha$ plays a crucial role in estimation as it
balances efficiency and robustness.

######  Result 2.

The system of equations for estimating the parameter under study
$\bm{\theta}{\,=\,}\\{(c_{j},\gamma_{i})^{T}\,;\,j{\,=\,}0,1\,;\,i{\,=\,}1,2,\dots,k\\}$
is obtained as follows.

$\sum_{i=1}^{k}\Bigg{[}\left(p_{i1}-\frac{n_{i1}}{n}\right)p^{\alpha-1}_{i1}p^{{}^{\prime}(\bm{\theta)}}_{i1}+\sum_{m=2}^{q_{i}}\left(p_{im}-\frac{n_{im}}{n}\right)\frac{p^{\alpha}_{im}A^{(i)}_{1}(\bm{\theta})}{G^{(m)}(\tau_{i};\gamma_{i})}\Bigg{]}=\left(p_{s}-\frac{n_{s}}{n}\right)p^{\alpha}_{s}A^{(k)}_{2}(\bm{\theta}).$

It is evident from Results (1) and (2) that explicit forms of MLE and MDPDE
cannot be obtained from the given set-up. One can rely on the iterative Co-
ordinate Descent method to obtain estimates. The steps of implementing this
method are presented in Algorithm (1).

Algorithm 1 Co-ordinate Descent

  * •

Initialize the starting value
$\bm{\theta}^{(0)}=(c^{(0)}_{j},\gamma^{(0)}_{i});j=0,1;i=1,2,\dots,k.$

  * •

For $t=0,1,\dots;$ the parameter estimates at $t+1^{th}$ iteration are

$\displaystyle c^{(t+1)}_{0}$ $\displaystyle=c^{(t)}_{0}-h\frac{\partial
D(c^{(t)}_{j},\gamma^{(t)}_{i})}{\partial c_{0}}.$ $\displaystyle
c^{(t+1)}_{1}$ $\displaystyle=c^{(t)}_{1}-h\frac{\partial
D(c^{(t+1)}_{0},c^{(t)}_{1},\gamma^{(t)}_{i})}{\partial c_{1}}.$
$\displaystyle\gamma^{(t+1)}_{1}$
$\displaystyle=\gamma^{(t)}_{1}-h\frac{\partial
D(c^{(t+1)}_{j},\gamma^{(t)}_{i})}{\partial\gamma_{1}}\;;\;j=0,1.$
$\displaystyle\dots$ $\displaystyle\qquad\dots\qquad\dots$
$\displaystyle\gamma^{(t+1)}_{k}$
$\displaystyle=\gamma^{(t)}_{k}-h\frac{\partial
D(c^{(t+1)}_{j},\gamma^{(t+1)}_{i},\gamma^{(t)}_{k})}{\partial\gamma_{k}}\;;\;i=1,2,\dots,k-1,$

where, $D(.)=-ln\,L(\bm{\theta})$ for MLE and $D(.)=D_{\alpha}(\bm{\theta})$
for DPDE and $h$ is the learning rate.

  * •

Continue the process until {$max|\bm{\theta}^{(t+1)}-\bm{\theta}^{(t)}|<c$}
where $c$ is a predefined threshold value.

###  3.1 Asymptotic Property

Due to non-linearity, the MDPD estimators (MDPDE) do not follow tractable
finite sample property. Thus, a general approch is to study the asymptotic
properties. In this section, asymptotic distribution of the MDPDE
${\hat{\bm{\theta}}}_{\alpha}$ is derived with the help of the procedure
followed by Calvino et al. [47].

######  Result 3.

Let $\bm{\theta}_{0}$ be true value of parameter $\bm{\theta}.$ The asymptotic
distribution of MDPD estimator ${\hat{\bm{\theta}}}_{\alpha}$ is given by

$\sqrt{n}({\hat{\bm{\theta}}}_{\alpha}-{\bm{\theta}_{0}})\xrightarrow[n\to\infty]{\mathscr{L}}N\Big{(}\bm{0}_{k+2},J^{-1}_{\alpha}({\bm{\theta}_{0}})K_{\alpha}({\bm{\theta}_{0}})J^{-1}_{\alpha}({\bm{\theta}_{0}})\Big{)},$

where, $J_{\alpha}({\bm{\theta}_{0}})$ and $K_{\alpha}({\bm{\theta}_{0}})$ are
defined in Appendix.

###### Proof.

Given in the appendix A.2. ∎

###  3.2 Optimal choice of tuning parameter

The optimization of tuning parameter $\alpha$ in DPD estimation is crucial as
it strikes a balance between robustness and efficiency. An approach is
suggested here to find optimal tuning parameter, which does not require any
iterative computational procedure. The aim is to find optimum $\alpha$, which
increases robustness and precision of the estimation. Thus, the objective
function is defined as

$\Phi_{\alpha}(\hat{\bm{\theta}})=C_{1}\,D_{\alpha}({\hat{\bm{\theta}}})+C_{2}\,tr\Big{(}J^{-1}_{\alpha}({\bm{\theta}_{0}})K_{\alpha}({\bm{\theta}_{0}})J^{-1}_{\alpha}({\bm{\theta}_{0}})\Big{)},$
(11)

where, $D_{\alpha}(\hat{\bm{\theta}})$ is DPD measure of the estimates and
$C_{1},C_{2}$ are predefined positive weight values with $C_{1}+C_{2}=1$. The
optimal tuning parameter will minimize $\Phi_{\alpha}(\hat{\bm{\theta}})$.

##  4 Robust Bayes method of estimation

Conventional Bayesian inference based on the likelihood based posterior may
not provide valid estimated values in presence of outliers in the data set. In
this section, a robust Bayesian inference is studied based on density power
divergence measure. As suggested by Ghosh and Basu [37], a maximizer equation
based on DPD measure for one-shot device testing unit under SSALT CRM with
IMIIP is presented as

$B_{\alpha}(\bm{\theta})=\frac{1}{\alpha}\left\\{\frac{n_{s}}{n}p^{\alpha}_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\frac{n_{im}}{n}p^{\alpha}_{im}\right\\}-\frac{1}{\alpha+1}\left\\{p_{s}^{\alpha+1}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}p^{\alpha+1}_{im}\right\\},$
(12)

where, MDPDE with $\alpha>0$ is the maximizer of $B_{\alpha}(\bm{\theta})$.
Therefore, robust posterior density, a pseudo posterior, can be defined as

$\pi_{\alpha}(\bm{\theta}|data)=\frac{\exp(B_{\alpha}(\bm{\theta}))\pi(\bm{\theta})}{\int\exp(B_{\alpha}(\bm{\theta}))\pi(\bm{\theta})d\bm{\theta}}.$
(13)

Here, $\pi_{\alpha}(\bm{\theta}|data)$ is the proper density for $\alpha>0.$
For $\alpha\to 0$, robust pseudo posterior will converge to conventional
likelihood-based posterior density. For any loss function $Loss(.,.)$, robust
Bayes estimator (RBE) can be obtained as

$arg\min_{t}\int
Loss(\bm{\theta},t)\pi_{\alpha}(\bm{\theta}|data)d\bm{\theta}.$

Particularly, for squared error loss function, robust Bayes estimator can be
derived as

$\hat{\bm{\theta}}^{(b)}_{\alpha}=\int\bm{\theta}\pi_{\alpha}(\bm{\theta}|data)d\bm{\theta}.$
(14)

###  4.1 Choice of priors

In Bayesian inference, choice of prior determines the estimation process. As
considered by Fan et al. [48], we have taken prior information on $p_{im}$
instead of model parameters $\bm{\theta}.$ To avoid zero-frequency situation,
we follow the idea of Lee and Morris [49] and modify empirical probabilities
as

$(\tilde{p}_{s},\tilde{p}_{im})=\left(\frac{n_{s}+1}{n+k\sum_{i=1}^{k}q_{i}+1},\frac{n_{im}+1}{n+k\sum_{i=1}^{k}q_{i}+1}\right),$
(15)

where, $i=1,2,\dots,k\;;\;m=1,2,\dots,q_{i}.$

####  4.1.1 Normal prior based on data

Assume $e_{im}$ is the error representing difference between empirical and
true failure probabilities. Therefore, it can be expressed that

$\tilde{p}_{im}=p_{im}+e_{im}\;;\;i=1,2,\dots,k\;;\;m=1,2,\dots,q_{i},$ (16)

where, the error $e_{im}$ are assumed to be independent $N(0,\sigma^{2})$
variables. The conditional likelihood function as prior distribution of
$\bm{\theta}$ given $\sigma^{2}$ can be obtained by

$L(\bm{\theta}|\sigma^{2})\propto\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}\frac{1}{\sigma\sqrt{2\pi}}\exp\left\\{\frac{1}{2\sigma^{2}}(p_{im}-\tilde{p}_{im})^{2}\right\\},$

and $\pi(\sigma^{2})\propto\frac{1}{\sigma^{2}}$ is the non-informative prior
of $\sigma^{2}.$ The joint prior density of $\bm{\theta}$ can be obtained as

$\pi^{(Nor)}(\bm{\theta})\propto\int_{0}^{\infty}L(\bm{\theta}|\sigma^{2})\pi(\sigma^{2})\,d\bm{\theta}\propto\left\\{\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}(p_{im}-\tilde{p}_{im})^{2}\right\\}^{-\sum_{i=1}^{k}q_{i}/2}.$
(17)

Thus, posterior density would be given as

$\pi_{\alpha}^{(Nor)}(\bm{\theta}|data)\propto\exp(B_{\alpha}(\bm{\theta}))\left\\{\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}(p_{im}-\tilde{p}_{im})^{2}\right\\}^{-\sum_{i=1}^{k}q_{i}/2}.$
(18)

####  4.1.2 Dirlichet prior based on data

Beta prior is a natural choice if a parameter can be interpreted as a
probability. Extending this idea, a Dirichlet prior is considered for the
failure and survival probabilities as

$\pi^{(Dir)}(\bm{\theta})=\frac{p^{\beta_{s}-1}_{s}\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}p^{\beta_{im}-1}_{im}}{Beta(\bm{\beta})},$
(19)

where, $\beta_{s},\beta_{im}>0$ for $i=1,2,\dots,k\;;\;m=1,2,\dots,q_{i}$ and

$Beta(\bm{\beta})=\frac{\Gamma\beta_{s}\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}\Gamma\beta_{im}}{\Gamma(\beta_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\beta_{im})}.$

The hyper-parameters $\bm{\beta}$ are chosen such that

$\displaystyle E(p_{im})$
$\displaystyle=\frac{\beta_{im}}{\beta_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}}=\tilde{p}_{im}\;,\;E(p_{s})=\frac{\beta_{s}}{\beta_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}}=\tilde{p}_{s}$
(20) $\displaystyle Var(p_{s})$
$\displaystyle=\frac{\beta_{s}\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}\beta_{im}}{\left(\beta_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\beta_{im}\right)^{2}\left(\beta_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\beta_{im}+1\right)}=\sigma_{(p)}^{2}.$
(21)

where, $\sigma_{(p)}^{2}$ is assumed to be a prefix quantity. The estimates of
hyper-parameters can be obtained by equations (20) and (21) as

$\displaystyle\hat{\beta}_{im}$
$\displaystyle=\tilde{p}_{im}\left\\{\frac{\tilde{p}_{s}(1-\tilde{p})}{\sigma_{(p)}^{2}}-1\right\\}\;;\;\hat{\beta}_{s}=\left\\{\frac{\tilde{p}_{s}(1-\tilde{p})}{\sigma_{(p)}^{2}}-1\right\\}-\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\beta_{im}.$
(22)

Therefore, the joint posterior density can be obtained as

$\pi_{\alpha}^{(Dir)}(\bm{\theta}|data)\propto\exp(B_{\beta}(\bm{\theta}))\left\\{p^{\hat{\beta}_{s}-1}_{s}\prod_{i=1}^{k}\prod_{m=1}^{q_{i}}p^{\hat{\beta}_{im}-1}_{im}\right\\}.$
(23)

Under both of the prior assumptions, Bayes estimate cannot be obtained in
closed form. In such situations, Monte Carlo Markov Chain simulation methods
can be used to approximate the Bayes estimates.

###  4.2 Hamiltonian Monte Carlo

As mentioned in the introduction, widely used methods like Gibbs sampler and
Metropolis-Hastings (MH) algorithm struggle with high dimensional or highly
correlated variables. To avoid this, Hamiltonian Monte Carlo (HMC) approach is
implemented to generate the posterior samples as it delivers more accurate
findings and flexible settings for complex models. As well explained by Thach
and Bris [50], Thomas and Tu [51], brief description of HMC proceeds as
follows. HMC is based on Hamiltonian dynamics, which describes motion of an
object in terms of its location $\bm{\theta}$ and momentum $\bm{\phi}$ at any
specific time $t$. For each location and momentum taken by any object, there
is an associated potential energy $U(\bm{\theta})$ and kinetic energy
$K(\bm{\phi})$ respectively. Thus, total energy of the system can be described
as the Hamiltonian function,
$H(\bm{\theta},\bm{\phi})=U(\bm{\theta})+K(\bm{\phi}),$ where,
$H(\bm{\theta},\bm{\phi})$ is a constant.

The objective is to generate $\bm{\theta}$ from a posterior density
$\pi_{\alpha}(\bm{\theta})$. In HMC, $U(\bm{\theta})$ is taken as
$U(\bm{\theta})=-\log\,\pi_{\alpha}(\bm{\theta})$ and assume $\bm{\phi}\sim
N_{k+2}(\bm{0},\bm{M}),$ where, covariance matrix $\bm{M}$ is user defined.
Then, we get

$H(\bm{\theta},\bm{\phi})=-\log\,\pi_{\alpha}(\bm{\theta})+\frac{1}{2}\bm{\phi}^{T}\bm{M}^{-1}\bm{\phi}.$
(24)

Over time, $t$, HMC is governed by Hamiltonian equations

$\displaystyle\frac{d\bm{\phi}}{dt}$ $\displaystyle=-\frac{\partial
H(\bm{\theta},\bm{\phi})}{\partial\bm{\theta}}=-\frac{\partial
U(\bm{\theta})}{\partial\bm{\theta}}=\nabla_{\theta}\log\,\pi_{\alpha}(\bm{\theta}),$
(25) $\displaystyle\frac{d\bm{\theta}}{dt}$ $\displaystyle=\frac{\partial
H(\bm{\theta},\bm{\phi})}{\partial\bm{\phi}}=\frac{\partial
K(\bm{\phi})}{\partial\bm{\phi}}=\bm{M}^{-1}\bm{\phi},$

where, $\nabla_{\theta}\log\,\pi_{\alpha}(\bm{\theta})$ is gradient of log
posterior density. A solution to the equations (25) determines trajectory
through which $\bm{\theta}$ values can be sampled. The leapfrog method is a
good approximation to the solutions of Hamiltonian equations. A leapfrog step
updates $\bm{\theta}$ and $\bm{\phi}$ with the step size $\epsilon$ over time
$t$ as

$\displaystyle\bm{\phi}\left(t+\frac{\epsilon}{2}\right)$
$\displaystyle=\bm{\phi}(t)+\frac{\epsilon}{2}\nabla_{\theta}\log\,\pi_{\alpha}\big{(}\bm{\theta}(t)\big{)},$
(26) $\displaystyle\bm{\theta}(t+\epsilon)$
$\displaystyle=\bm{\theta}(t)+\epsilon\,\bm{M}^{-1}\bm{\phi}\left(t+\frac{\epsilon}{2}\right),$
$\displaystyle\bm{\phi}(t+\epsilon)$
$\displaystyle=\bm{\phi}\left(t+\frac{\epsilon}{2}\right)+\frac{\epsilon}{2}\nabla_{\theta}\log\,\pi_{\alpha}\big{(}\bm{\theta}(t+\epsilon)\big{)}.$

For HMC, $L$ denotes the number of leapfrog steps required. The steps are
given in the Algorithm (2). The negation of momentum $\bm{\phi}$ is taken
after final leapfrog step to ensure the reversibility of HMC transition.
Since, the leapfrog solution is an approximation, the acceptance rate of HMC
proposals is less than 100% but higher than MH algorithm. After the accept
(reject) step, $\phi^{(t)}$ instantly gets updated at the beginning of
subsequent iteration.

Algorithm 2 Hamiltonian Monte Carlo

  * •

Initialize the position state $\bm{\theta}^{(0)}.$

  * For $t=1,2,\dots,N$

  * •

Sample $\bm{\phi}^{(t)}\sim N(\bm{0},\bm{M}).$

  * •

Run leapfrog as given in equation (26) starting at
$(\bm{\theta}^{(t)},\bm{\phi}^{(t)})$ for $L$ step with step size $\epsilon$
to produce proposed state $(\bm{\theta}^{*},\bm{\phi}^{*})$.

  * •

Compute acceptance probability
$acc=min\Big{\\{}1,\exp\big{(}U(\bm{\theta}^{(t)})-U(\bm{\theta}^{*})+K(\bm{\phi}^{(t)})\big{)}-K(\bm{\phi}^{*})\big{)}\Big{\\}}.$

  * •

Generate a random number $u\sim U(0,1)$ and set
$\bm{\theta}^{(t+1)}=\begin{cases}\bm{\theta}^{*}&;\;u\leq acc.\\\
\bm{\theta}^{(t)}&;\;\text{otherwise}.\end{cases}$

For estimation purpose, $m^{{}^{\prime}}$ chains of $N$ values each are
generated through HMC, and the first $N^{{}^{\prime}}$ values are discarded as
burn-in period in each chain. The Bayes estimator based on squared error loss
can be approximated as

$\bm{\hat{\theta}}=\frac{1}{m^{{}^{\prime}}(N-N^{{}^{\prime}})}\sum_{l=1}^{m^{{}^{\prime}}}\sum_{t=N^{{}^{\prime}}+1}^{N}\bm{\theta}^{(t)}_{l}.$

where $\bm{\theta}^{(t)}_{l}$ is the value of $t^{th}$ iteration for $l^{th}$
chain.

##  5 Testing of hypothesis based on robust Bayes factor

Validating whether available data supports hypothesis of interest is essential
for inferential study. For datasets with outliers, robust testing of
hypothesis is pertinent. This section develops robust testing of hypothesis
based on Bayes factor inspired by the procedure followed by Ghosh et al. [52].
For parameter $\bm{\theta}=(c_{0},c_{1},\gamma_{i}\;;\;i=1,2,\dots,k),$
consider the vector-valued function
$fn:\mathbb{R}^{k+2}\xrightarrow[]{}\mathbb{R}^{w}.$ The null and alternative
hypotheses are given as

$\bm{H}_{0}:\bm{\theta}\in\bm{\Theta}_{0}\quad\text{against}\quad\bm{H}_{1}:\bm{\theta}\in\bm{\Theta}_{1},$

where,
$\bm{\Theta}_{0}=\\{\bm{\theta}\in\bm{\Theta}_{0}:fn(\bm{\theta})=\bm{0}_{w}\\}$
and $\bm{\Theta}_{1}=\\{\bm{\theta}\notin\bm{\Theta}_{0}\\}$. Let $\rho_{0}$
and $1-\rho_{0}$ be prior probabilities under $\bm{\Theta}_{0}$ and
$\bm{\Theta}_{1}$ respectively. Let $\pi_{j}(\bm{\theta})$ be prior density of
$\bm{\theta}$ under $\bm{\Theta}_{j}$ such that,
$\int_{\bm{\Theta}_{j}}\pi_{j}(\bm{\theta})d\bm{\theta}=1\,;\,j=0,1$. Then,
the prior can be expressed as

$\pi(\bm{\theta})=\rho_{0}\pi_{0}(\bm{\theta})I\\{\bm{\theta}\in\bm{\Theta}_{0}\\}+(1-\rho_{0})\pi_{1}(\bm{\theta})I\\{\bm{\theta}\in\bm{\Theta}_{1}\\}.$

Therefore, posterior probabilities under $\bm{\Theta}_{0}$ and
$\bm{\Theta}_{1}$ are

$\displaystyle P_{\pi_{\alpha}}(\bm{\theta}\in\bm{\Theta}_{0}|data)$
$\displaystyle=\frac{\rho_{0}}{M_{\alpha}(\pi)}\int_{\bm{\Theta}_{0}}\exp(B_{\alpha}(\bm{\theta}))\pi_{0}(\bm{\theta})\,d\bm{\theta}.$
$\displaystyle P_{\pi_{\alpha}}(\bm{\theta}\in\bm{\Theta}_{1}|data)$
$\displaystyle=\frac{(1-\rho_{0})}{M_{\alpha}(\pi)}\int_{\bm{\Theta}_{1}}\exp(B_{\alpha}(\bm{\theta}))\pi_{1}(\bm{\theta})\,d\bm{\theta},$

where, $M_{\alpha}(\pi)$ is the marginal density expressed as

$M_{\alpha}(\pi)=\rho_{0}\int_{\bm{\Theta}_{0}}\exp(B_{\alpha}(\bm{\theta}))\pi_{0}(\bm{\theta})d\bm{\theta}+(1-\rho_{0})\int_{\bm{\Theta}_{1}}\exp(B_{\alpha}(\bm{\theta}))\pi_{1}(\bm{\theta})\,d\bm{\theta}.$

The posterior odds ratio of $H_{0}$ relative to $H_{1}$ is given as

$\frac{P_{\pi_{\alpha}}(\bm{\theta}\in\bm{\Theta}_{0}|data)}{P_{\pi_{\alpha}}(\bm{\theta}\in\bm{\Theta}_{1}|data)}=\left(\frac{\rho_{0}}{1-\rho_{0}}\right)BF_{01},$
(27)

where, $BF_{01}$ is the Bayes factor given as

$BF_{01}=\frac{\int_{\bm{\Theta}_{0}}\exp(B_{\alpha}(\bm{\theta}))\,d\bm{\theta}}{\int_{\bm{\Theta}_{1}}\exp(B_{\alpha}(\bm{\theta}))\,d\bm{\theta}}.$
(28)

The Bayes factor measures strength of evidence the data offers supporting one
hypothesis over another; Jeffreys [38] suggested a scale to interpret the
Bayes factor and Kass and Raftery [53] simplified it further, which is given
in Table 1.

Table 1: Interpretation of Bayes factor [53]. BF01 | Support for $\bm{H}_{0}$
---|---
$<1$ | Negative
1 to 3 | Not worth more than a bare mention
3 to 20 | Positive
20 to 150 | Strong
$>150$ | Very Strong

##  6 Property of robustenss

This section includes robustness analysis through influence function (IF).
Suppose, for a true distribution $F_{\bm{\theta}}$, functional of any
estimator is denoted by $T_{\alpha}(F_{\bm{\theta}}).$ Then, influence
function is defined as

$IF(t;T_{\alpha},F_{\bm{\theta}})=\lim_{\epsilon\to
0}\frac{T_{\alpha}(U_{\varepsilon})-T_{\alpha}(F_{\bm{\theta}})}{\epsilon}=\left.\frac{\partial(T_{\alpha}(U_{\epsilon}))}{\partial\epsilon}\right|_{\epsilon\to
0^{+}}.$ (29)

Here, $U_{\epsilon}=(1-\epsilon)F_{\bm{\theta}}+\epsilon\Delta_{t}$ is the
contaminated model where $\epsilon,$ $(0<\epsilon<1)$ is the proportion of
contamination and $\Delta_{t}$ denotes the degenerate distribution at point
$t$.

###  6.1 Influence function of MDPDE

Let $F_{\bm{\theta}}$ be the true distribution from where data is generated.
If $T_{\alpha}(F_{\bm{\theta}})$ is statistical functional of MDPDE
$\hat{\bm{\theta}}_{\alpha}$, $T_{\alpha}(F_{\bm{\theta}})$ will be the value
of $\bm{\theta}$ which will minimize,

$p^{\alpha+1}_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}p^{\alpha+1}_{im}-\left(1+\frac{1}{\alpha}\right)\left\\{\left(\int_{I_{s}}dF_{\bm{\theta}}\right)p^{\alpha}_{s}+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left(\int_{I_{im}}dF_{\bm{\theta}}\right)p^{\alpha}_{im}\right\\},$
(30)

where $I_{s}=\\{t:t>\tau_{k}\\}$ and
$I_{im}=\\{t:\tau_{i,m-1}<t\leq\tau_{i,m}\\}$.

######  Result 4.

The influence function of $\hat{\bm{\theta}}_{\alpha}$ under SSALT CRM with
IMIIP for one-shot device testing units is given as

$\displaystyle IF(t;T_{\alpha},F_{\bm{\theta}})=J^{-1}_{\alpha}(\bm{\theta})$
$\displaystyle\left[\left\\{\Delta^{(I_{s})}_{t}-p_{s}\right\\}p^{\alpha-1}_{s}\frac{\partial(p_{s})}{\partial\bm{\theta}}\right.$
$\displaystyle\qquad\left.+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left\\{\Delta^{(I_{im})}_{t}-p_{im}\right\\}p^{\alpha-1}_{im}\frac{\partial(p_{im})}{\partial\bm{\theta}}\right].$
(31)

Here, $\Delta^{(I)}_{t}=\begin{cases}1\quad\text{if}\ t\in I\\\
0\quad\text{otherwise}\\\ \end{cases}.$

###### Proof.

Replacing $F_{\bm{\theta}}$ by contaminated model
$U_{\epsilon}=(1-\epsilon)F_{\bm{\theta}}+\epsilon\Delta_{t}$ in equation (30)
and differentiating it with respect to $\epsilon$ and putting $\epsilon\to
0^{+}$ the desired result can be obtained. ∎

###  6.2 Influence function of RBE

To study robustness through IF [37], Bayes functional of
$\hat{\bm{\theta}}^{(b)}_{\alpha}$ under squared error loss function is given
as

$T^{(b)}_{\alpha}(F_{\bm{\theta}})=\frac{\int\bm{\theta}\exp\big{\\{}B_{\alpha}(\bm{\theta};F_{\bm{\theta}})\big{\\}}\pi(\bm{\theta})\,d\bm{\theta}}{\int\exp\big{\\{}B_{\alpha}(\bm{\theta};F_{\bm{\theta}})\big{\\}}\pi(\bm{\theta})\,d\bm{\theta}},$
(32)

where,

$\displaystyle B_{\alpha}(\bm{\theta};F_{\bm{\theta}})$
$\displaystyle=\frac{1}{\alpha}\left\\{\left(\int_{I_{s}}dF_{\bm{\theta}}\right)p^{\alpha}_{s}(\bm{\theta})+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left(\int_{I_{im}}dF_{\bm{\theta}}\right)p^{\alpha}_{im}(\bm{\theta})\right\\}$
$\displaystyle\qquad\qquad-\frac{1}{\alpha+1}\left\\{p^{\alpha+1}_{s}(\bm{\theta})+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}p^{\alpha+1}_{im}(\bm{\theta})\right\\}$

######  Result 5.

The influence function of Bayes estimator $\hat{\bm{\theta}}^{(b)}_{\alpha}$
under SSALT CRM with IMIIP for one-shot device testing units is given by

$IF(t;T^{(b)}_{\alpha},F_{\bm{\theta}})=Cov_{(p)}(\bm{\theta},X_{\alpha}(\bm{\theta};t,f_{\bm{\theta}})).$

###### Proof.

Given in the appendix A.3. ∎

###  6.3 Influence function of Bayes factor

Here, robustness property of the Bayes factor is examined by deriving its IF
when null hypothesis is true. Let $F_{\theta_{0}}$ be the true distribution
under the null hypothesis $H_{0}:\bm{\theta}\in\bm{\Theta}_{0}$ and therefore
functional related to the Bayes factor can be defined as

$T^{(\alpha)}_{\bm{\Theta}}(F_{\bm{\theta}_{0}})=\frac{\int_{\bm{\Theta}_{0}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0};F_{\bm{\theta}_{0}})\big{\\}}\pi_{0}(\bm{\theta})\,d\bm{\theta}}{\int_{\bm{\Theta}_{1}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1};F_{\bm{\theta}_{1}})\big{\\}}\pi_{1}(\bm{\theta})\,d\bm{\theta}}.$
(33)

Here,
$B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{j};F_{\bm{\theta}_{0}})\,;\,j=0,1$, is
expressed as

$\displaystyle B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{j};F_{\bm{\theta}_{0}})$
$\displaystyle=\frac{1}{\alpha}\left\\{\left(\int_{I_{s}}dF_{\bm{\theta}_{0}}\right)p^{\alpha}_{s}(\bm{\theta}\in\bm{\Theta}_{j})+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left(\int_{I_{im}}dF_{\bm{\theta}_{0}}\right)p^{\alpha}_{im}(\bm{\theta}\in\bm{\Theta}_{j})\right\\}$
$\displaystyle\qquad\qquad-\frac{1}{\alpha+1}\left\\{p^{\alpha+1}_{s}(\bm{\theta}\in\bm{\Theta}_{j})+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}p^{\alpha+1}_{im}(\bm{\theta}\in\bm{\Theta}_{j})\right\\}.$

Let contamination in the true distribution $F_{\theta_{0}}$ under
$H_{0}:\bm{\theta}\in\bm{\Theta}_{0}$ be given as
$U_{\epsilon}=(1-\epsilon)F_{\bm{\theta}_{0}}+\epsilon\Delta_{t}$. The IF with
respect to contamination $U_{\epsilon}$ is obtained as

$IF(t;T^{(\alpha)}_{\bm{\Theta}},F_{\bm{\theta}_{0}})=\left.\frac{\partial(T^{(\alpha)}_{\bm{\Theta}}(U_{\epsilon}))}{\partial\epsilon}\right|_{\epsilon\to
0^{+}}.$

The following results provide explicit expression of IF under the given setup.

######  Result 6.

The influence function of Bayes factor BF01 under SSALT CRM with IMIIP for
one-shot device testing units is obtained as

$IF(t;T^{(\alpha)}_{\bm{\Theta}},F_{\bm{\theta}_{0}})=Y_{\alpha}(\bm{\Theta})\bigg{\\{}E\Big{[}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\Big{]}-E\Big{[}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\Big{]}\bigg{\\}}.$

###### Proof.

Given in the appendix A.4. ∎

The maximum value of IF shows the degree of bias resulting from contamination.
Therefore, the smaller the value of IF, the more robust the estimator or Bayes
factor. Also, for all the estimators and Bayes factor, IF is a bounded
function of $t.$

##  7 Standard Lehman Family of Distributions: Special cases

As discussed in section 2, standard Lehman family of distributions is a
generalized distribution family. The members of this family can emerge from
the different choices of $Q(t;\gamma_{i})$. Two special cases of the family
are discussed here.

###  7.1 Special case 1: Weibull lifetime distribution

Weibull distribution is one of the most widely used lifetime distributions for
inferential analysis. It is frequently used for one-shot device testing [54,
55]. It emerges from Lehman family of distributions if we put

$Q(t;\gamma_{i})=t^{\gamma_{i}}\;;\;\gamma_{i}>0.$

The cdf and pdf from equation (1) takes the form

$\displaystyle F^{(1)}_{i}(t)$
$\displaystyle=1-\exp\big{\\{}-\lambda_{i}t^{\gamma_{i}}\big{\\}}\;;\;\gamma_{i}\,,\lambda_{i}>0.$
$\displaystyle f^{(1)}_{i}(t)$
$\displaystyle=\lambda_{i}\gamma_{i}t^{\gamma_{i}-1}\exp\big{\\{}-\lambda_{i}t^{\gamma_{i}}\big{\\}}.$

Therefore, survival function under cumulative risk step-stress model is given
as

$S_{1}(t)=\begin{cases}\exp\Big{\\{}-\lambda_{1}t^{\gamma_{1}}\Big{\\}}\quad;\;0<t\leq\tau_{1}.\\\
\exp\Big{\\{}-D_{1}^{(\delta)}(t;\gamma_{i-1,i})\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{i-1}\tau^{\gamma_{i-1}}_{i-1}+\sum\limits_{l=1}^{i-2}E_{1}^{(\delta)}(\tau_{l};\gamma_{l+1,l})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{i-1}<t\leq\tau_{i-1}+\delta\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{i}t^{\gamma_{i}}\Big{\\}}\exp\bigg{\\{}-\sum\limits_{l=1}^{i-1}E_{1}^{(\delta)}(\tau_{l};\gamma_{l+1,l})\bigg{\\}}\\\
\hskip 110.96556pt;\;\tau_{i-1}+\delta<t\leq\tau_{i}\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{k}t^{\gamma_{k}}\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{k-1}\tau_{k-1}^{\gamma_{k}-1}+\sum\limits_{i=1}^{k-2}E_{1}^{(\delta)}(\tau_{i};\gamma_{i+1,i})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{k-1}<t<\infty,\end{cases}$

where,

$\displaystyle D_{1}^{(\delta)}(t;\gamma_{i-1,i})$
$\displaystyle=\frac{(t-\tau_{i-1})^{2}}{2\delta}\Bigg{[}\Big{\\{}2\delta(t-\tau_{i-1})^{-1}-1\Big{\\}}\lambda_{i-1}\gamma_{i-1}\tau^{\gamma_{i-1}-1}_{i-1}+\lambda_{i}\gamma_{i}t^{\gamma_{i}-1}\Bigg{]}.$
$\displaystyle E_{1}^{(\delta)}(\tau_{l};\gamma_{l+1,l})$
$\displaystyle=\lambda_{l}\tau^{\gamma_{l}}_{l}-\lambda_{l+1}(\tau_{l}+\delta)^{\gamma_{l}+1}+\frac{\delta}{2}\Big{\\{}\lambda_{l}\gamma_{l}\tau^{\gamma_{l}-1}_{l}+\lambda_{l+1}\gamma_{l+1}(\tau_{l}+\delta)^{\gamma_{l+1}-1}\Big{\\}}.$

Rest of the analysis is similar to the developments of previous sections.

###  7.2 Special case 2: Gompertz lifetime distribution

Though Gompertz lifetime distribution has yet not been seen for one-shot
device testing data analysis, there is enough literature to show that it is
used as a lifetime distribution to fulfil different objectives [56, 57]. It
emerges from Lehman family of distributions if we put

$Q(t;\gamma_{i})=(e^{\gamma_{i}t}-1)\;;\;\gamma_{i}>0.$

The cdf and pdf from equation (1) take the form

$\displaystyle F^{(2)}_{i}(t)$
$\displaystyle=1-\exp\big{\\{}-\lambda_{i}(e^{\gamma_{i}t}-1)\big{\\}}\;;\;\gamma_{i}\,,\lambda_{i}>0.$
$\displaystyle f^{(2)}_{i}(t)$
$\displaystyle=\lambda_{i}\gamma_{i}\exp\big{\\{}\gamma_{i}t-\lambda_{i}(e^{\gamma_{i}t}-1)\big{\\}}.$

Therefore, survival function under the cumulative risk step-stress model is

$S_{2}(t)=\begin{cases}\exp\Big{\\{}-\lambda_{1}(e^{\gamma_{1}t}-1)\Big{\\}}\quad;\;0<t\leq\tau_{1}.\\\
\exp\Big{\\{}-D_{2}^{(\delta)}(t;\gamma_{i-1,i})\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{i-1}(e^{{\gamma_{i-1}}\tau_{i-1}}-1)+\sum\limits_{l=1}^{i-2}E_{2}^{(\delta)}(\tau_{l};\gamma_{l+1,l})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{i-1}<t\leq\tau_{i-1}+\delta\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{i}(e^{{\gamma_{i}}t}-1)\Big{\\}}\exp\bigg{\\{}-\sum\limits_{l=1}^{i-1}E_{2}^{(\delta)}(\tau_{l};\gamma_{l+1,l})\bigg{\\}}\\\
\hskip 110.96556pt;\;\tau_{i-1}+\delta<t\leq\tau_{i}\;;\;i=2,3,\dots,k-1.\\\
\exp\Big{\\{}-\lambda_{k}(e^{{\gamma_{k}}t}-1)\Big{\\}}\exp\bigg{[}-\Big{\\{}\lambda_{k-1}(e^{{\gamma_{k-1}}\tau_{k-1}}-1)+\sum\limits_{i=1}^{k-2}E_{2}^{(\delta)}(\tau_{i};\gamma_{i+1,i})\Big{\\}}\bigg{]}\\\
\hskip 110.96556pt;\;\tau_{k-1}<t<\infty,\end{cases}$

where,

$\displaystyle D_{2}^{(\delta)}(t;\gamma_{i-1,i})$
$\displaystyle=\frac{(t-\tau_{i-1})^{2}}{2\delta}\Bigg{[}\Big{\\{}2\delta(t-\tau_{i-1})^{-1}-1\Big{\\}}\lambda_{i-1}\gamma_{i-1}e^{\gamma_{i-1}\tau_{i-1}}+\lambda_{i}\gamma_{i}e^{\gamma_{i}t}\Bigg{]}.$
$\displaystyle E_{2}^{(\delta)}(\tau_{l};\gamma_{l+1,l})$
$\displaystyle=\lambda_{l}(e^{\gamma_{l}\tau_{l}}-1)-\lambda_{l+1}(e^{\gamma_{l+1}(\tau_{l}+\delta)}-1)+\frac{\delta}{2}\Big{\\{}\lambda_{l}\gamma_{l}e^{\gamma_{l}\tau_{l}}+\lambda_{l+1}\gamma_{l+1}e^{\gamma_{l+1}(\tau_{l}+\delta)}\Big{\\}}.$

In next sections, performance of the theoretical results developed in previous
sections has been assessed numerically through simulation experiments and data
analysis under Weibull and Gompertz lifetime distributions.

##  8 Simulation study

For simulation analysis, $100$ one-shot device testing units are put into
three-step-stress ALT with $m=9$ inspection time points under cumulative risk
model having lag period $\delta=0.001.$ Stress levels and stress change time
points are taken as
$(x_{1}{=}1.5,\tau_{1}{=}3;\,x_{2}{=}2.5,\tau_{2}{=}6;\,x_{3}{=}3.5,\tau_{3}{=}9).$
The experiment is terminated at $\tau_{3}=9.$ The intermediate inspection
times, including stress change time points, are $(1,2,3,4,5,6,7,8,9).$

###  8.1 Special case 1: Weibull lifetime distribution

To generate data from Weibull lifetime distribution under the given set-up,
true model parameters are set as $\bm{\theta}=$
$(c_{0}=0.03,\,c_{1}=-0.08,\,g_{1}=0.1,\,g_{2}=0.2,\,g_{3}=0.3)^{{}^{\prime}}.$
To study robustness, a contaminated version of the true distribution is taken
by setting
$(c_{0}=0.03+0.003,\,c_{1}=-0.08+0.003,\,g_{1}=0.1-0.001,\,g_{2}=0.2-0.02,\,g_{3}=0.3-0.001)$
and then data are generated from it. Robustness can be studied through bias of
the estimates. Hence, mean absolute bias (MAB) and mean square error (MSE) are
obtained through Monte Carlo simulation based on 1000 generations. The
asymptotic behaviour of the estimates is also studied by calculating $95\%$
coverage probability (CP) and average width (AW). For classical inference,
maximum likelihood estimate (MLE) and minimum density power divergence
estimate (MDPDE) are obtained through the use of a Coordinate-descent method
given in the algorithm 1. The outcomes are reported in the table 2.

Table 2: Outcomes of MLE and MDPDE with Weibull lifetime distribution.

| Pure data | Contamination
---|---|---
| MAB | MSE | CP | AL | MAB | MSE | CP | AL
$\bm{\hat{c}}_{0}$ |  |  |  |  |  |  |  |
MLE | 0.006210 | 0.000352 | 94.7 | 0.093012 | 0.060324 | 0.008665 | 78.7 | 0.443638
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.035460 | 0.004293 | 90.8 | 0.289933 | 0.038359 | 0.004233 | 89.2 | 0.329038
$\bm{0.4}$ | 0.015613 | 0.001126 | 97.8 | 0.266549 | 0.017686 | 0.001219 | 96.9 | 0.235131
$\bm{0.6}$ | 0.005463 | 0.000076 | 100.0 | 0.143675 | 0.005963 | 0.000098 | 100.0 | 0.168868
$\bm{0.8}$ | 0.004224 | 0.000017 | 100.0 | 0.162807 | 0.004230 | 0.000017 | 100.0 | 0.162813
$\bm{1.0}$ | 0.004152 | 0.000017 | 100.0 | 0.223977 | 0.004155 | 0.000017 | 100.0 | 0.223983
$\bm{\hat{c}}_{1}$ |  |  |  |  |  |  |  |
MLE | 0.001421 | 0.000848 | 91.8 | 0.096852 | 0.020130 | 0.007350 | 63.9 | 0.122242
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.016174 | 0.003425 | 80.1 | 0.099060 | 0.016878 | 0.003245 | 78.6 | 0.102229
$\bm{0.4}$ | 0.010992 | 0.000814 | 94.4 | 0.089141 | 0.011625 | 0.000896 | 93.4 | 0.088933
$\bm{0.6}$ | 0.006206 | 0.000079 | 99.7 | 0.080844 | 0.006496 | 0.000094 | 99.4 | 0.081741
$\bm{0.8}$ | 0.005130 | 0.000026 | 100.0 | 0.080211 | 0.005133 | 0.000026 | 100.0 | 0.080215
$\bm{1.0}$ | 0.005083 | 0.000025 | 100.0 | 0.080777 | 0.005086 | 0.000025 | 100.0 | 0.080780
$\bm{\hat{g}}_{1}$ |  |  |  |  |  |  |  |
MLE | 0.003564 | 0.000545 | 95.8 | 0.184265 | 0.027673 | 0.001669 | 92.1 | 0.162307
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.035498 | 0.001682 | 93.2 | 0.157878 | 0.035690 | 0.001756 | 92.8 | 0.157747
$\bm{0.4}$ | 0.040026 | 0.001716 | 96.9 | 0.158761 | 0.040559 | 0.001782 | 96.0 | 0.157029
$\bm{0.6}$ | 0.040292 | 0.001631 | 99.7 | 0.163775 | 0.040520 | 0.001652 | 99.8 | 0.164146
$\bm{0.8}$ | 0.040006 | 0.001600 | 100.0 | 0.172950 | 0.040013 | 0.001601 | 100.0 | 0.172942
$\bm{1.0}$ | 0.040006 | 0.001600 | 100.0 | 0.185714 | 0.040010 | 0.001600 | 100.0 | 0.185710
$\bm{\hat{g}}_{2}$ |  |  |  |  |  |  |  |
MLE | 0.010504 | 0.003447 | 96.3 | 0.324029 | 0.072695 | 0.006452 | 99.9 | 0.406256
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.087867 | 0.008145 | 100.0 | 0.400015 | 0.085221 | 0.007781 | 99.9 | 0.395964
$\bm{0.4}$ | 0.095773 | 0.009283 | 100.0 | 0.414747 | 0.094511 | 0.009080 | 100.0 | 0.407944
$\bm{0.6}$ | 0.099468 | 0.009903 | 100.0 | 0.424099 | 0.099218 | 0.009856 | 100.0 | 0.424193
$\bm{0.8}$ | 0.099958 | 0.009991 | 100.0 | 0.430833 | 0.099951 | 0.009990 | 100.0 | 0.430805
$\bm{1.0}$ | 0.099977 | 0.009995 | 100.0 | 0.437649 | 0.099973 | 0.009994 | 100.0 | 0.437633
$\bm{\hat{g}}_{3}$ |  |  |  |  |  |  |  |
MLE | 0.031766 | 0.023953 | 95.5 | 0.132358 | 0.069585 | 0.006410 | 88.0 | 0.595815
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.059301 | 0.004117 | 94.3 | 0.326769 | 0.061222 | 0.004384 | 93.9 | 0.362652
$\bm{0.4}$ | 0.053443 | 0.002990 | 98.6 | 0.254988 | 0.054399 | 0.003119 | 98.7 | 0.241889
$\bm{0.6}$ | 0.050519 | 0.002563 | 100.0 | 0.164869 | 0.050755 | 0.0025914 | 100.0 | 0.179242
$\bm{0.8}$ | 0.050043 | 0.002504 | 100.0 | 0.155494 | 0.050049 | 0.002504 | 100.0 | 0.155503
$\bm{1.0}$ | 0.050025 | 0.002502 | 100.0 | 0.162836 | 0.050029 | 0.002502 | 100.0 | 0.162842

Table 3: Outcomes of BE and RBE with Weibull lifetime distribution.

| Normal Prior | Dirichlet Prior
---|---|---
| Pure data | Contamination | Pure | Contamination
| MAB | MSE | MAB | MSE | MAB | MSE | MAB | MSE
$\bm{\hat{c}}_{0}$ |  |  |  |  |  |  |  |
BE | 0.000996 | 0.000001 | 0.003993 | 0.000006 | 0.001028 | 0.000002 | 0.002731 | 0.000006
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.002006 | 0.000004 | 0.002080 | 0.000004 | 0.002025 | 0.000004 | 0.001994 | 0.000004
$\bm{0.4}$ | 0.001999 | 0.000001 | 0.002013 | 0.000003 | 0.001992 | 0.000008 | 0.002027 | 0.000004
$\bm{0.6}$ | 0.001992 | 0.000003 | 0.002693 | 0.000005 | 0.001295 | 0.000009 | 0.002007 | 0.000004
$\bm{0.8}$ | 0.001990 | 0.000008 | 0.002124 | 0.000009 | 0.001882 | 0.000004 | 0.001997 | 0.000004
$\bm{1.0}$ | 0.003983 | 0.000016 | 0.002334 | 0.000012 | 0.001681 | 0.000001 | 0.002024 | 0.000004
$\bm{\hat{c}}_{1}$ |  |  |  |  |  |  |  |
BE | 0.000998 | 0.000001 | 0.002924 | 0.000003 | 0.000989 | 0.000001 | 0.002715 | 0.000004
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.001993 | 0.000004 | 0.002324 | 0.000004 | 0.002835 | 0.000003 | 0.001997 | 0.000004
$\bm{0.4}$ | 0.002021 | 0.000002 | 0.002147 | 0.000004 | 0.002003 | 0.000002 | 0.001970 | 0.000004
$\bm{0.6}$ | 0.001945 | 0.000009 | 0.002254 | 0.000012 | 0.001991 | 0.000004 | 0.002009 | 0.000004
$\bm{0.8}$ | 0.001998 | 0.000002 | 0.002165 | 0.000008 | 0.002032 | 0.000003 | 0.001988 | 0.000004
$\bm{1.0}$ | 0.005010 | 0.000025 | 0.002270 | 0.000028 | 0.002280 | 0.000005 | 0.001965 | 0.000004
$\bm{\hat{g}}_{1}$ |  |  |  |  |  |  |  |
BE | 0.010002 | 0.000100 | 0.049913 | 0.000468 | 0.009198 | 0.000168 | 0.009983 | 0.00019
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.019998 | 0.000422 | 0.019989 | 0.000399 | 0.009796 | 0.000099 | 0.010001 | 0.000100
$\bm{0.4}$ | 0.020004 | 0.000445 | 0.020011 | 0.000476 | 0.009768 | 0.000098 | 0.010015 | 0.000100
$\bm{0.6}$ | 0.019966 | 0.000398 | 0.022219 | 0.000411 | 0.009844 | 0.000096 | 0.010019 | 0.000100
$\bm{0.8}$ | 0.019983 | 0.000399 | 0.021889 | 0.000414 | 0.009933 | 0.000114 | 0.009268 | 0.000093
$\bm{1.0}$ | 0.040001 | 0.001600 | 0.019994 | 0.000452 | 0.010013 | 0.000154 | 0.009488 | 0.000090
$\bm{\hat{g}}_{2}$ |  |  |  |  |  |  |  |
BE | 0.009964 | 0.000099 | 0.029805 | 0.000610 | 0.010009 | 0.000100 | 0.019718 | 0.000131
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.020003 | 0.000442 | 0.019982 | 0.000399 | 0.009966 | 0.000099 | 0.009886 | 0.000100
$\bm{0.4}$ | 0.019991 | 0.000399 | 0.019976 | 0.000435 | 0.009988 | 0.000148 | 0.010006 | 0.000100
$\bm{0.6}$ | 0.020021 | 0.000401 | 0.020239 | 0.000421 | 0.009981 | 0.000099 | 0.009010 | 0.000099
$\bm{0.8}$ | 0.019991 | 0.000399 | 0.020248 | 0.000432 | 0.010027 | 0.000181 | 0.010023 | 0.000100
$\bm{1.0}$ | 0.049974 | 0.002497 | 0.020535 | 0.000487 | 0.010004 | 0.000137 | 0.010008 | 0.000100
$\bm{\hat{g}}_{3}$ |  |  |  |  |  |  |  |
BE | 0.009968 | 0.000099 | 0.04927 | 0.000622 | 0.010008 | 0.000100 | 0.013388 | 0.000100
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.019986 | 0.000399 | 0.026911 | 0.000453 | 0.009992 | 0.000142 | 0.010015 | 0.000100
$\bm{0.4}$ | 0.019996 | 0.000411 | 0.020265 | 0.000426 | 0.009982 | 0.000099 | 0.010008 | 0.000100
$\bm{0.6}$ | 0.019993 | 0.000435 | 0.023388 | 0.000468 | 0.010013 | 0.000155 | 0.009937 | 0.000099
$\bm{0.8}$ | 0.020010 | 0.000470 | 0.020140 | 0.000832 | 0.009977 | 0.000099 | 0.009975 | 0.000097
$\bm{1.0}$ | 0.049987 | 0.002498 | 0.020555 | 0.000426 | 0.010032 | 0.000192 | 0.009996 | 0.000100

(a) MLE and DPDE

(b) MLE and DPDE

(c) BE and RBE (Normal prior)

(d) BE and RBE (Normal prior)

(e) BE and RBE (Dirichlet prior)

(f) BE and RBE (Dirichlet prior)

Figure 2: MAB (left) and MSE (right) for reliability estimates with Weibull
lifetime distribution under contamination.

For Bayesian inference, Bayes estimate (BE) and robust Bayes estimate (RBE)
are obtained by using Hamiltonian Monte Carlo (HMC) given in the algorithm 2.
For smooth running of the HMC algorithm, we consider step size $\epsilon=0.01$
and number of steps $L=10.$ Define $v=(0.01,0.02,0.02,0.02,0.02)$, andn
$\bm{M}$ is taken as diagonal matrix whose diagonal elements are $1/v$. Two
chains of $N=1200$ values are generated through HMC, and first
$N^{{}^{\prime}}=200$ values are discarded as burn-in period. The MAB and MSE
for BE and RBE under normal and Dirichlet prior are reported in table 3. For
Dirichlet prior $\sigma^{2}_{(p)}=0.03$ is taken.

After observing table 2, it is evident that MAB and MSE of MLE are lesser than
that of MDPDE in pure data. The CP is also closer to the nominal level for MLE
with smaller AL than that of MDPDE when data is pure. However, when
contamination is incorporated, MDPDE performs better than MLE in terms of MAB
and MSE. The CP is now relatively closer to the nominal level for MDPDE,
resulting in a smaller AL than that for MLE. Further observation shows that
from pure to contaminated data, MAB and MSE for MLE are increased, but there
is not much change in the values of MAB and MSE for MDPDE. This behaviour of
MDPDE reflects its robustness over MLE. Similarly, in table 3, close
observation indicates that MAB of RBE is lesser than that of BE under normal
and Dirichlet prior for contamination scheme and vice-versa for pure data.
Under this setup, both the priors are equally suitable as there is no
significant difference in the estimated values under both the prior
assumptions. Also, change of MAB is not very significant for RBE from pure to
contaminated scheme as compared to BE. Thus, it can be interpreted that RBE is
robust. If we compare tables 2 and 3, we observe that BE and RBE have the
least MAB among the four methods of estimation in pure and contaminated data,
respectively. Further, to illustrate robustness graphically, MAB and MSE of
reliability estimates for parameters in the contaminated settings are plotted
in Figure 2. The superiority of robust estimates under contamination is
clearly visible from these figures.

###  8.2 Special case 2: Gompertz lifetime distribution

Data from Gompertz lifetime distribution is generated by taking the true model
parameters as $\bm{\theta}=$
$(c_{0}=0.8,\,c_{1}=-0.1,\,g_{1}=0.1,\,g_{2}=0.2,\,g_{3}=0.2)^{{}^{\prime}}.$
The data is contaminated by using the contaminated version of the true
distribution by taking
$(c_{0}=0.8-0.1,\,c_{1}=-0.1+0.002,\,g_{1}=0.1+0.002,\,g_{2}=0.2-0.03,\,g_{3}=0.2-0.01).$
MAB, MSE, 95% CP and AL for MLE and MDPDE are reported in table 4. In Bayesian
inference, for HMC implementation, we have set $\epsilon=0.001$, $L=5$ and
$v=(0.004,0.006,0.01,0.01,0.01)$ and $\sigma^{2}_{(p)}=0.03.$ The MAB and MSE
for BE and RBE under normal prior and Dirichlet prior are presented in table
3.

Table 4: Outcomes of MLE and MDPDE with Gompertz lifetime distribution.

| Pure data | Contamination
---|---|---
| MAB | MSE | CP | AL | MAB | MSE | CP | AL
$\bm{\hat{c}}_{0}$ |  |  |  |  |  |  |  |
MLE | 0.008307 | 0.000120 | 97.1 | 0.346066 | 0.008621 | 0.000105 | 89.1 | 2.167666
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.009956 | 0.000013 | 100.0 | 2.545396 | 0.009960 | 0.000022 | 100.0 | 3.116749
$\bm{0.4}$ | 0.009995 | 0.000090 | 100.0 | 0.319118 | 0.009981 | 0.000062 | 100.0 | 1.038814
$\bm{0.6}$ | 0.009995 | 0.000073 | 100.0 | 0.262667 | 0.009990 | 0.000090 | 100.0 | 0.661661
$\bm{0.8}$ | 0.009996 | 0.000024 | 100.0 | 0.264314 | 0.009994 | 0.000094 | 100.0 | 0.556532
$\bm{1.0}$ | 0.009998 | 0.000036 | 100.0 | 0.314978 | 0.009996 | 0.000096 | 100.0 | 0.557034
$\bm{\hat{c}}_{1}$ |  |  |  |  |  |  |  |
MLE | 0.013544 | 0.001083 | 95.1 | 0.043646 | 0.016713 | 0.001241 | 93.9 | 0.044877
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.011467 | 0.000131 | 100.0 | 0.064463 | 0.011987 | 0.000144 | 100.0 | 0.066855
$\bm{0.4}$ | 0.010368 | 0.000107 | 100.0 | 0.042604 | 0.010972 | 0.000120 | 100.0 | 0.055235
$\bm{0.6}$ | 0.010312 | 0.000106 | 100.0 | 0.044559 | 0.010512 | 0.000110 | 100.0 | 0.058068
$\bm{0.8}$ | 0.010238 | 0.000104 | 100.0 | 0.052028 | 0.010285 | 0.000105 | 100.0 | 0.067381
$\bm{1.0}$ | 0.010186 | 0.000103 | 100.0 | 0.081008 | 0.010175 | 0.000103 | 100.0 | 0.088434
$\bm{\hat{g}}_{1}$ |  |  |  |  |  |  |  |
MLE | 0.010096 | 0.000215 | 96.3 | 0.041026 | 0.020060 | 0.000407 | 96.2 | 0.142888
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.017631 | 0.000311 | 100.0 | 0.165632 | 0.017539 | 0.000308 | 100.0 | 0.194334
$\bm{0.4}$ | 0.017225 | 0.000296 | 100.0 | 0.045792 | 0.017260 | 0.000298 | 100.0 | 0.076492
$\bm{0.6}$ | 0.017423 | 0.000303 | 100.0 | 0.053297 | 0.017441 | 0.000304 | 100.0 | 0.067158
$\bm{0.8}$ | 0.017757 | 0.000315 | 100.0 | 0.067310 | 0.017748 | 0.000315 | 100.0 | 0.075996
$\bm{1.0}$ | 0.018074 | 0.000326 | 100.0 | 0.088888 | 0.018061 | 0.000326 | 100.0 | 0.094482
$\bm{\hat{g}}_{2}$ |  |  |  |  |  |  |  |
MLE | 0.010081 | 0.000213 | 93.2 | 0.037677 | 0.042079 | 0.001969 | 88.3 | 0.138608
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.053141 | 0.002825 | 100.0 | 0.180930 | 0.051715 | 0.002677 | 100.0 | 0.214539
$\bm{0.4}$ | 0.020974 | 0.000440 | 100.0 | 0.045709 | 0.046895 | 0.002199 | 100.0 | 0.090937
$\bm{0.6}$ | 0.015966 | 0.000254 | 100.0 | 0.051827 | 0.042957 | 0.001845 | 100.0 | 0.074955
$\bm{0.8}$ | 0.011947 | 0.000142 | 100.0 | 0.065408 | 0.039628 | 0.001570 | 100.0 | 0.077528
$\bm{1.0}$ | 0.018591 | 0.000730 | 100.0 | 0.089864 | 0.036800 | 0.001354 | 100.0 | 0.089440
$\bm{\hat{g}}_{3}$ |  |  |  |  |  |  |  |
MLE | 0.001093 | 0.000273 | 94.4 | 0.040282 | 0.016640 | 0.001500 | 88.5 | 0.119585
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.001758 | 0.000003 | 100.0 | 0.143980 | 0.002555 | 0.000006 | 100.0 | 0.173854
$\bm{0.4}$ | 0.001239 | 0.000001 | 100.0 | 0.050538 | 0.000973 | 0.000001 | 100.0 | 0.078295
$\bm{0.6}$ | 0.000327 | 0.000002 | 100.0 | 0.061725 | 0.000503 | 0.000003 | 100.0 | 0.069509
$\bm{0.8}$ | 0.001740 | 0.000003 | 100.0 | 0.082397 | 0.001815 | 0.000003 | 100.0 | 0.077250
$\bm{1.0}$ | 0.002961 | 0.000008 | 100.0 | 0.117836 | 0.002966 | 0.000008 | 100.0 | 0.095132

Table 5: Outcomes of BE and RBE with Gompertz lifetime distribution.

| Normal Prior | Dirichlet Prior
---|---|---
| Pure data | Contamination | Pure data | Contamination
| MAB | MSE | MAB | MSE | MAB | MSE | MAB | MSE
$\bm{\hat{c}}_{0}$ |  |  |  |  |  |  |  |
BE | 0.004850 | 0.000024 | 0.005855 | 0.011024 | 0.004580 | 0.000024 | 0.104974 | 0.011024
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.005006 | 0.000028 | 0.004854 | 0.000026 | 0.005027 | 0.000025 | 0.049801 | 0.002495
$\bm{0.4}$ | 0.004751 | 0.000029 | 0.004802 | 0.000024 | 0.005165 | 0.000027 | 0.050071 | 0.002500
$\bm{0.6}$ | 0.005110 | 0.000025 | 0.005277 | 0.000029 | 0.004966 | 0.000028 | 0.050034 | 0.002572
$\bm{0.8}$ | 0.005153 | 0.000026 | 0.004941 | 0.000023 | 0.005211 | 0.000025 | 0.049457 | 0.002459
$\bm{1.0}$ | 0.004755 | 0.000027 | 0.004877 | 0.000024 | 0.004883 | 0.000026 | 0.050016 | 0.002532
$\bm{\hat{c}}_{1}$ |  |  |  |  |  |  |  |
BE | 0.004004 | 0.000025 | 0.005991 | 0.000048 | 0.004000 | 0.000025 | 0.008762 | 0.000048
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.004981 | 0.000025 | 0.004668 | 0.000022 | 0.004960 | 0.000023 | 0.005439 | 0.000034
$\bm{0.4}$ | 0.004781 | 0.000028 | 0.005578 | 0.000021 | 0.004978 | 0.000029 | 0.008086 | 0.000033
$\bm{0.6}$ | 0.004895 | 0.000029 | 0.004976 | 0.000020 | 0.005123 | 0.000027 | 0.007000 | 0.000013
$\bm{0.8}$ | 0.004802 | 0.000026 | 0.004477 | 0.000024 | 0.005377 | 0.000028 | 0.007000 | 0.000040
$\bm{1.0}$ | 0.005051 | 0.000026 | 0.004944 | 0.000027 | 0.004989 | 0.000024 | 0.007768 | 0.000018
$\bm{\hat{g}}_{1}$ |  |  |  |  |  |  |  |
BE | 0.004262 | 0.000024 | 0.006163 | 0.000090 | 0.004285 | 0.000024 | 0.002999 | 0.000008
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.004832 | 0.000028 | 0.004585 | 0.000021 | 0.004965 | 0.000027 | 0.003430 | 0.000009
$\bm{0.4}$ | 0.005041 | 0.000025 | 0.004782 | 0.000028 | 0.004959 | 0.000026 | 0.002874 | 0.000003
$\bm{0.6}$ | 0.005179 | 0.000028 | 0.005342 | 0.000025 | 0.005207 | 0.000025 | 0.003688 | 0.000009
$\bm{0.8}$ | 0.004996 | 0.000025 | 0.005066 | 0.000027 | 0.004999 | 0.000026 | 0.003026 | 0.000002
$\bm{1.0}$ | 0.004455 | 0.000026 | 0.004714 | 0.000024 | 0.005260 | 0.000026 | 0.002735 | 0.000004
$\bm{\hat{g}}_{2}$ |  |  |  |  |  |  |  |
BE | 0.004323 | 0.000026 | 0.005168 | 0.000624 | 0.004228 | 0.000022 | 0.064324 | 0.000982
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.005277 | 0.000021 | 0.004565 | 0.000024 | 0.005080 | 0.000028 | 0.025000 | 0.000625
$\bm{0.4}$ | 0.005071 | 0.000025 | 0.005184 | 0.000025 | 0.004986 | 0.000029 | 0.024302 | 0.000624
$\bm{0.6}$ | 0.005257 | 0.000026 | 0.005125 | 0.000022 | 0.005016 | 0.000025 | 0.024686 | 0.000637
$\bm{0.8}$ | 0.004956 | 0.000024 | 0.004647 | 0.000021 | 0.005164 | 0.00008 | 0.026864 | 0.000635
$\bm{1.0}$ | 0.004890 | 0.000024 | 0.004962 | 0.000024 | 0.004999 | 0.000027 | 0.025000 | 0.000625
$\bm{\hat{g}}_{3}$ |  |  |  |  |  |  |  |
BE | 0.004387 | 0.000020 | 0.014719 | 0.000224 | 0.004447 | 0.000024 | 0.019004 | 0.000501
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.005253 | 0.000025 | 0.004759 | 0.000025 | 0.005405 | 0.000025 | 0.014799 | 0.000240
$\bm{0.4}$ | 0.004955 | 0.000024 | 0.004707 | 0.000022 | 0.004756 | 0.000027 | 0.015160 | 0.000225
$\bm{0.6}$ | 0.004894 | 0.000024 | 0.004414 | 0.000024 | 0.005371 | 0.000027 | 0.015168 | 0.000282
$\bm{0.8}$ | 0.004950 | 0.000024 | 0.004873 | 0.000028 | 0.005413 | 0.000028 | 0.015406 | 0.000222
$\bm{1.0}$ | 0.005205 | 0.000025 | 0.005174 | 0.000025 | 0.004939 | 0.000029 | 0.015126 | 0.000284

(a) MLE and DPDE

(b) MLE and DPDE

(c) BE and RBE (Normal prior)

(d) BE and RBE (Normal prior)

(e) BE and RBE (Dirichlet prior)

(f) BE and RBE (Dirichlet prior)

Figure 3: MAB (left) and MSE (right) for reliability estimates with Gompertz
lifetime distribution under contamination.

The close observation demonstrates that MDPDE has smaller MAB and MSE than MLE
under contamination in table 4. By careful examination of table 5, it is
observed that RBE has less MAB than BE under both the priors when outliers are
present in the data. It is also pointed out that normal prior is the better
choice than Dirichlet prior as smaller MAB is observed under normal prior. The
comparison of tables 4 and 5 shows that RBE provides the least MAB among all
four methods of estimation under contamination. The better performance of
MDPDE and RBE than MLE and BE under contamination reflects their robustness,
which can also be visible from figure 3.

##  9 Data analysis

For practical implementation of the results previously discussed, a dataset
examining reliability of light bulbs is taken here from the experimental study
conducted by Zhu [58]. Balakrishnan et al. [14] used this data for one-shot
device testing under step-stress ALT. In light bulb experiment $n=64$
miniature light bulbs are put into a two-step SSALT experiment where voltage
$x_{1}=2.25V$ is applied for up to $\tau_{1}=96$ hours and then voltage is
increased to $x_{2}=2.44V$. The stopping time of experiment is $\tau_{2}=140$
hours. The failure times of light bulbs within 140 hours are recorded as
follows.

12.07, 19.5, 22.1, 23.11, 24, 25.1, 26.9, 36.64, 44.1, 46.3, 54, 58.09, 64.17,
72.25, 86.9, 90.09, 91.22, 102.1, 105.1, 109.2, 114.4, 117.9, 121.9, 122.5,
123.6, 126.5, 130.1, 14, 17.95, 24, 26.46, 26.58, 28.06, 34, 36.13, 40.85,
41.11, 42.63, 52.51, 62.68, 73.13, 83.63, 91.56, 94.38, 97.71, 101.53, 105.11,
112.11, 119.58, 120.2, 126.95, 129.25, 136.31.

$n_{s}=11$ lightbulbs survived after termination of the experiment.
Intermediate inspection times are taken as $\tau=(32,64,96,111,126,140)$ and
lag period $\delta=0.001$ is considered. To test if standard Lehman family of
distribution is fitted to data for the given model, a bootstrap-based goodness
of fit test is performed, and an approximated p-value is obtained. The
distance-based test statistic to conduct testing is defined as

$TS=\left|\frac{n_{s}-\hat{n}_{s}}{\hat{n}_{s}}\right|+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left|\frac{n_{im}-\hat{n}_{im}}{\hat{n}_{im}}\right|.$
(34)

Here, $\hat{n}_{im}$ and $\hat{n}_{s}$ are expected number of failures and
survivals obtained through MLE. For deriving MLE and MDPDE, coordinate descent
method is used. For deriving BE and RBE, HMC is used where we consider
$\epsilon=0.001$, $L=10$, $v=(0.01,0.01,0.01,0.01)$ and $M=1/v.$ For Dirichlet
prior $\sigma^{2}_{(p)}=0.05$ is taken.

Table 6: Information regarding the analysis under special cases from Lehman
family of distributions for the data.

Information | Weibull | Gompertz
---|---|---
Initial values | $(-0.9,-0.06,0.2,0.7)$ | $(0.03,0.22,0.08,0.07)$
Test statistic | $TS_{1}=7.702629$ | $TS_{2}=11.06017$
p-value | $0.687$ | $0.252$

For computational convenience, failure times and inspection time points are
multiplied by 0.2 and 0.1 for Weibull and Gompertz lifetime distributions,
respectively. Extensive grid search provides the initial values
$(c_{0},c_{1},g_{1},g_{2})$ for both lifetime distributions which are given in
table 6. A bootstrap-based goodness of fit test is conducted with the test
statistic given in the equation (34). The value of test statistics and
corresponding p-value for both lifetime distributions are reported in table 6.
The significant p-values indicate the suitability of both lifetime
distributions to the data.

Table 7: Optimal tuning parameter with Weibull
$\big{(}\Phi^{(1)}_{\alpha}(\hat{\bm{\theta}})\big{)}$ and Gompertz
$\big{(}\Phi^{(2)}_{\alpha}(\hat{\bm{\theta}})\big{)}$ lifetime distribution
for the data.

$\bm{\alpha}$ | $\Phi^{(1)}_{\alpha}(\hat{\bm{\theta}})$ | $\Phi^{(2)}_{\alpha}(\hat{\bm{\theta}})$ | $\bm{\alpha}$ | $\Phi^{(1)}_{\alpha}(\hat{\bm{\theta}})$ | $\Phi^{(2)}_{\alpha}(\hat{\bm{\theta}})$
---|---|---|---|---|---
0.10 | 0.471993 | 0.400612 | 0.15 | 0.488222 | 0.377846
0.20 | 0.507091 | 0.354708 | $\bm{0.25}$ | 0.528757 | $\bm{0.332835}$
0.30 | 0.544979 | 0.385612 | 0.35 | 0.563394 | 0.385350
0.40 | 0.575179 | 0.387732 | 0.45 | 0.571144 | 0.393123
0.50 | 0.544086 | 0.402022 | 0.55 | 0.506264 | 0.415094
0.60 | 0.468243 | 0.433223 | $\bm{0.65}$ | $\bm{0.443078}$ | 0.457578
0.70 | 0.462979 | 0.489712 | 0.75 | 0.485322 | 0.531703
0.80 | 0.510201 | 0.586361 | 0.85 | 0.537709 | 0.657534
0.90 | 0.567938 | 0.750579 | 0.95 | 0.600977 | 0.873099
1.00 | 0.636911 | 1.036144 |  |  |

The optimal choice of tuning parameter is obtained using method described in
subsection (3.2), and value of objective function
$\big{(}\Phi^{(1)}_{\alpha}(\hat{\bm{\theta}})\big{)}$ by equation (11) for a
grid of tuning parameters are reported in table 7. It is observed that
$\alpha_{opt}=0.65$ is the optimal value of tuning parameter. Similarly, for
Gompertz lifetime distribution optimal choice of tuning parameter observed
from table 7 is $\alpha_{opt}=0.25$. The estimates derived from MLE and MDPDE
with $95\%$ asymptotic confidence intervals (CI) are reported in tables 8. The
estimates derived from BE and RBE with $95\%$ highest posterior density
credible interval ((HPD CRI) are reported in tables 9 and 10. The bootstrap
bias (BT Bias) and root mean square of error (RMSE) of the four estimation
methods are given in tables 11 and 12. From these tables, we observe the
smaller magnitude of BT bias and RMSE for MDPDE and RBE in classical and
Bayesian frameworks respectively.

Table 8: Parameter estimates (95% Asymp. CI) with Weibull and Gompertz
lifetime distribution for the data.

| $\bm{\hat{c}}_{0}$ | $\bm{\hat{c}}_{1}$ | $\bm{\hat{g}}_{1}$ | $\bm{\hat{g}}_{2}$
---|---|---|---|---
| Est.(CI) | Est.(CI) | Est.(CI) | Est.(CI)
| Weibull lifetime distribution
MLE | -0.896599 | -0.019708 | 0.227660 | 0.672829
| (-0.9935,-0.7996) | (-0.1594,0.1200) | (0.1305,0.3247) | (0.6213,0.7243)
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | -0.895879 | -0.012961 | 0.218379 | 0.660822
| (-1.0031,-0.7885) | (-0.1836,0.1576) | (0.1082,0.3285) | (0.6002,0.7214)
$\bm{0.4}$ | -0.896026 | -0.016272 | 0.208129 | 0.659650
| (-1.0060,-0.7859) | (-0.2094,0.1768) | (0.0882,0.3280) | (0.5924,0.7268)
$\bm{0.6}$ | -0.899175 | -0.051281 | 0.200906 | 0.692237
| (-0.9916,-0.8067) | (-0.2128,0.1103) | (0.0974,0.3043) | (0.6345,0.7499)
$\bm{0.8}$ | -0.899919 | -0.059150 | 0.200011 | 0.699200
| (-0.9929,-0.8069) | (-0.2336,0.1153) | (0.0871,0.3128) | (0.6372,0.7611)
$\bm{1.0}$ | -0.899937 | -0.059344 | 0.199956 | 0.699341
| (-0.9974,-0.8023) | (-0.2611,0.1424) | (0.0691,0.3308) | (0.6288,0.7698)
$\bm{0.65}$ | -0.899901 | -0.058965 | 0.200091 | 0.699082
| (-0.9895,-0.8102) | (-0.2161,0.0982) | (0.0988,0.3013) | (0.6426,0.7555)
| Gompertz lifetime distribution
MLE | 0.029733 | 0.219952 | 0.080961 | 0.071048
| (-0.0217,0.0812) | ( 0.0579,0.3819) | (0.0436,0.1182) | (0.0389,0.1031)
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.028974 | 0.219194 | 0.081802 | 0.073907
| (-0.0324,0.0903) | (0.0745,0.3638) | (0.0480,0.1155) | (0.0453,0.1024)
$\bm{0.4}$ | 0.029931 | 0.219863 | 0.079909 | 0.070284
| (-0.0453,0.1052) | (0.0528,0.3868) | (0.0449,0.1148) | (0.0393,0.1012)
$\bm{0.6}$ | 0.029969 | 0.219882 | 0.079790 | 0.070148
| ( -0.0715,0.1314) | (0.0423,0.3973) | (0.0448,0.1147) | (0.0388,0.1014)
$\bm{0.8}$ | 0.029988 | 0.219904 | 0.079761 | 0.070077
| (-0.1198,0.1798) | (0.0208,0.4189) | (0.0441,0.1154) | (0.0377,0.1024)
$\bm{1.0}$ | 0.029998 | 0.219924 | 0.079771 | 0.070040
| (-0.2116,0.2716) | (-0.0225,0.4623) | (0.0421,0.1174) | (0.0353,0.1047)
$\bm{0.25}$ | 0.028474 | 0.218576 | 0.081885 | 0.075741
| (-0.0367,0.0937) | (0.0823,0.3548) | (0.0495,0.1142) | (0.0490,0.1024)

Table 9: Bayes estimates (95% HPD CRI) of parameters with Weibull lifetime
distribution for the data.

| $\bm{\hat{c}}_{0}$ | $\bm{\hat{c}}_{1}$ | $\bm{\hat{g}}_{1}$ | $\bm{\hat{g}}_{2}$
---|---|---|---|---
| Est.(CRI) | Est.(CRI) | Est.(CRI) | Est.(CRI)
| Normal Prior
BE | -0.899879 | -0.059033 | 0.200273 | 0.699513
| (-0.9016,-0.8976) | (-0.0614,-0.0572) | (0.1984,0.2025) | (0.6979,0.7010)
RBE with $\bm{\alpha}$
$\bm{0.2}$ | -0.899399 | -0.061377 | 0.198872 | 0.701296
| (-0.9011,-0.8973) | (-0.0630,-0.0596) | (0.1969,0.2008) | (0.6998,0.7027)
$\bm{0.4}$ | -0.899050 | -0.060713 | 0.199531 | 0.700412
| (-0.9011,-0.8974) | (-0.0623,-0.0588) | (0.1974,0.2012) | (0.6988,0.7017)
$\bm{0.6}$ | -0.899727 | -0.059471 | 0.200222 | 0.699272
| (-0.9014,-0.8979) | (-0.0616,-0.0577) | (0.1984,0.2020) | (0.6979, 0.7009)
$\bm{0.8}$ | -0.900680 | -0.058280 | 0.199315 | 0.699322
| (-0.9026,-0.8988) | (-0.0596,-0.0559) | (0.1973,0.2015) | (0.6978, 0.7008)
$\bm{1.0}$ | -0.898687 | -0.060482 | 0.201138 | 0.700334
| (-0.9004,-0.8965) | (-0.0623,-0.0582) | (0.1991,0.2027) | (0.6987,0.7016)
$\bm{0.65}$ | -0.901501 | -0.060721 | 0.200756 | 0.699890
| (-0.9031,-0.8996) | (-0.0624,-0.0587) | (0.1989,0.2027) | (0.6985,0.7018)
| Dirichlet Prior
BE | -0.901220 | -0.060709 | 0.198343 | 0.700696
| (-0.9032,-0.8992) | (-0.0627,-0.0587) | (0.1964,0.2004) | (0.6993,0.7025)
RBE with $\bm{\alpha}$
$\bm{0.2}$ | -0.900204 | -0.059956 | 0.198501 | 0.698962
| (-0.9019,-0.8979) | (-0.0619,-0.0580) | (0.1964,0.2005) | (0.6974,0.7005)
$\bm{0.4}$ | -0.899018 | -0.060378 | 0.201762 | 0.701691
| (-0.9011,-0.8972) | ( -0.0625,-0.0586) | (0.1996,0.2036) | (0.6998,0.7031)
$\bm{0.6}$ | -0.898608 | -0.058784 | 0.199480 | 0.699089
| (-0.9004,-0.8963) | (-0.0607,-0.0570) | (0.1975,0.2013) | (0.6971,0.7005)
$\bm{0.8}$ | -0.899173 | -0.061371 | 0.199810 | 0.699252
| (-0.9011,-0.8968) | (-0.0635,-0.0591) | (0.1977,0.2015) | (0.6976,0.7009)
$\bm{1.0}$ | -0.9004 | -0.0588 | 0.1984 | 0.6991
| (-0.9025,-0.8983) | (-0.0603,-0.0563) | (0.1967,0.2006) | (0.6976,0.7008)
$\bm{0.65}$ | -0.901899 | -0.058795 | 0.198142 | 0.698615
| (-0.9040,-0.9002) | (-0.0607,-0.0570) | (0.1964,0.1998) | (0.6968,0.6998)

Table 10: Bayes estimates (95% HPD CRI) of parameters with Gompertz lifetime
distribution for the data.

| $\bm{\hat{c}}_{0}$ | $\bm{\hat{c}}_{1}$ | $\bm{\hat{g}}_{1}$ | $\bm{\hat{g}}_{2}$
---|---|---|---|---
| Est.(CRI) | Est.(CRI) | Est.(CRI) | Est.(CRI)
| Normal Prior
BE | 0.028192 | 0.218093 | 0.079698 | 0.070048
| (0.0264,0.0303) | (0.2157,0.2196) | (0.0784,0.0815) | (0.0680,0.0720)
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.030546 | 0.219769 | 0.078377 | 0.071961
| (0.0290,0.0327) | (0.2179,0.2217) | (0.0767,0.0798) | (0.0702,0.0739)
$\bm{0.4}$ | 0.030583 | 0.220174 | 0.080330 | 0.068060
| (0.0282,0.0322) | (0.2183,0.2223) | (0.0789,0.0823) | (0.0660,0.0701)
$\bm{0.6}$ | 0.029263 | 0.218703 | 0.079530 | 0.068788
| (0.0271,0.0309) | (0.2169,0.2204) | (0.0777,0.0810) | (0.0667,0.0705)
$\bm{0.8}$ | 0.030077 | 0.218702 | 0.078272 | 0.070884
| (0.0284,0.0318) | (0.2168,0.2206) | (0.0769,0.0798) | (0.0693,0.0730)
$\bm{1.0}$ | 0.031034 | 0.218872 | 0.081396 | 0.069204
| (0.0291,0.0331) | (0.2166,0.2203) | (0.0796,0.0827) | (0.0673,0.0709)
$\bm{0.25}$ | 0.031885 | 0.221929 | 0.079725 | 0.071999
| (0.0299,0.0336) | (0.2202,0.2240) | (0.0783,0.0811) | (0.0702,0.0740)
| Dirichlet Prior
BE | 0.028054 | 0.218051 | 0.078103 | 0.068342
| (0.0262,0.0301) | (0.2163,0.2195) | (0.0767,0.0796) | (0.0666,0.0700)
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.030799 | 0.220753 | 0.079553 | 0.071053
| (0.0291,0.0326) | (0.2186,0.2223) | (0.0779,0.0809) | (0.0694,0.0730)
$\bm{0.4}$ | 0.030076 | 0.218702 | 0.078272 | 0.070884
| (0.0284,0.0318) | (0.2168,0.2206) | (0.0769,0.0798) | (0.0693,0.0730)
$\bm{0.6}$ | 0.031885 | 0.221929 | 0.079725 | 0.071998
| (0.0299,0.0336) | (0.2202,0.2240) | (0.0783,0.0811) | (0.0702,0.0740)
$\bm{0.8}$ | 0.029923 | 0.219409 | 0.081313 | 0.068159
| (0.0280,0.0317) | (0.2173,0.2208) | (0.0800,0.0828) | (0.0659,0.0697)
$\bm{1.0}$ | 0.030372 | 0.220228 | 0.078431 | 0.071512
| (0.0283,0.0320) | (0.2187,0.2222) | (0.0770,0.0804) | (0.0698,0.0733)
$\bm{0.25}$ | 0.030065 | 0.219617 | 0.081656 | 0.071046
| (0.0282,0.0319) | (0.2177,0.2218) | (0.0798,0.0834) | (0.0693,0.0731)

Table 11: BT bias and RMSE of the estimates with Weibull lifetime distribution
for the data.

| $\bm{\hat{c}}_{0}$ | $\bm{\hat{c}}_{1}$ | $\bm{\hat{g}}_{1}$ | $\bm{\hat{g}}_{2}$
---|---|---|---|---
| BT Bias | RMSE | BT Bias | RMSE | BT Bias | RMSE | BT Bias | RMSE
MLE | 0.000362 | 0.009260 | 0.024155 | 0.089807 | -0.059080 | 0.094638 | -0.118604 | 0.220863
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | -0.000365 | 0.004297 | 0.009252 | 0.044963 | -0.027467 | 0.043435 | -0.056180 | 0.103130
$\bm{0.4}$ | -0.000251 | 0.000438 | 0.003091 | 0.004230 | -0.009852 | 0.014875 | -0.021145 | 0.030138
$\bm{0.6}$ | -0.000087 | 0.000216 | 0.001635 | 0.002156 | -0.003843 | 0.006427 | -0.008969 | 0.013393
$\bm{0.8}$ | -0.000022 | 0.000046 | 0.001126 | 0.001155 | -0.001713 | 0.001836 | -0.004558 | 0.004697
$\bm{1.0}$ | -0.000005 | 0.000021 | 0.000963 | 0.000980 | -0.001129 | 0.001152 | -0.003269 | 0.003276
$\bm{0.65}$ | -0.000067 | 0.000163 | 0.001434 | 0.001564 | -0.003115 | 0.004905 | -0.007440 | 0.010374
| Normal Prior
BE | -0.000018 | 0.001148 | -0.000100 | 0.001170 | 0.000023 | 0.001179 | 0.000047 | 0.001138
RBE with $\bm{\alpha}$
$\bm{0.2}$ | -0.000032 | 0.001184 | -0.000023 | 0.001153 | 0.000069 | 0.001159 | -0.000032 | 0.001158
$\bm{0.4}$ | -0.000007 | 0.001140 | -0.000036 | 0.001169 | -0.000017 | 0.001160 | -0.000092 | 0.001170
$\bm{0.6}$ | 0.000023 | 0.001129 | 0.000081 | 0.001166 | 0.000009 | 0.001144 | -0.000003 | 0.001146
$\bm{0.8}$ | 0.000009 | 0.001148 | -0.000042 | 0.001147 | -0.000034 | 0.001163 | -0.000030 | 0.001145
$\bm{1.0}$ | 0.000100 | 0.001142 | -0.000016 | 0.001152 | -0.000029 | 0.001169 | 0.000077 | 0.001167
$\bm{0.65}$ | 0.000014 | 0.001176 | 0.000021 | 0.001173 | 0.000044 | 0.001127 | 0.000021 | 0.001158
| Dirichlet Prior
BE | -0.000025 | 0.001184 | 0.000095 | 0.001132 | 0.000029 | 0.001119 | -0.000051 | 0.0011402
RBE with $\bm{\alpha}$
$\bm{0.2}$ | -0.000043 | 0.001178 | 0.000019 | 0.001149 | -0.000028 | 0.001152 | 0.000001 | 0.001169
$\bm{0.4}$ | -0.000128 | 0.001177 | 0.000040 | 0.001149 | 0.000141 | 0.001141 | 0.000126 | 0.001134
$\bm{0.6}$ | -0.000092 | 0.001133 | 0.000039 | 0.001130 | -0.000030 | 0.001131 | 0.000015 | 0.001165
$\bm{0.8}$ | 0.000098 | 0.001137 | -0.000042 | 0.001141 | 0.000046 | 0.001156 | -0.000006 | 0.001161
$\bm{1.0}$ | 0.000055 | 0.001151 | 0.000024 | 0.001136 | -0.000007 | 0.001166 | -0.000094 | 0.001148
$\bm{0.65}$ | 0.000058 | 0.001131 | -0.000065 | 0.001161 | -0.000091 | 0.001185 | 0.000007 | 0.001146

Table 12: BT bias and RMSE of the estimates with Gompertz lifetime
distribution for the data.

| $\bm{\hat{c}}_{0}$ | $\bm{\hat{c}}_{1}$ | $\bm{\hat{g}}_{1}$ | $\bm{\hat{g}}_{2}$
---|---|---|---|---
| BT Bias | RMSE | BT Bias | RMSE | BT Bias | RMSE | BT Bias | RMSE
MLE | 0.001979 | 0.004073 | 0.008310 | 0.016801 | 0.010338 | 0.014504 | 0.000887 | 0.004937
MDPDE with $\bm{\alpha}$
$\bm{0.2}$ | 0.000280 | 0.001093 | 0.001181 | 0.004475 | 0.001951 | 0.004056 | 0.000353 | 0.001364
$\bm{0.4}$ | 0.000120 | 0.000628 | 0.000461 | 0.002398 | 0.000630 | 0.001311 | 0.000032 | 0.000798
$\bm{0.6}$ | 0.000017 | 0.000032 | 0.000062 | 0.000119 | 0.000236 | 0.000261 | 0.000065 | 0.000095
$\bm{0.8}$ | 0.000009 | 0.000018 | 0.000031 | 0.000063 | 0.000123 | 0.000151 | 0.000035 | 0.000051
$\bm{1.0}$ | 0.000004 | 0.000011 | 0.000015 | 0.000035 | 0.000064 | 0.000096 | 0.000018 | 0.000027
$\bm{0.25}$ | 0.000060 | 0.000395 | 0.000217 | 0.001448 | 0.001007 | 0.002140 | 0.000354 | 0.000740
| Normal Prior
BE | -0.000074 | 0.001164 | -0.000026 | 0.001182 | -0.000030 | 0.001180 | -0.000080 | 0.001187
RBE with $\bm{\alpha}$
$\bm{0.2}$ | -0.000006 | 0.001130 | -0.000006 | 0.001152 | -0.000046 | 0.001135 | 0.000012 | 0.001163
$\bm{0.4}$ | -0.000025 | 0.001188 | -0.000081 | 0.001150 | -0.000026 | 0.001111 | 0.000047 | 0.001156
$\bm{0.6}$ | 0.000094 | 0.001118 | -0.000070 | 0.001146 | -0.000029 | 0.001144 | 0.000039 | 0.001163
$\bm{0.8}$ | 0.000017 | 0.001164 | 0.000054 | 0.001136 | -0.000078 | 0.001172 | 0.000015 | 0.001137
$\bm{1.0}$ | -0.000080 | 0.001154 | 0.000068 | 0.001158 | 0.000045 | 0.001119 | 0.000024 | 0.001170
$\bm{0.25}$ | 0.000053 | 0.001125 | 0.000051 | 0.001130 | -0.000028 | 0.001132 | 0.0000205 | 0.001184
| Dirichlet Prior
BE | 0.000049 | 0.001144 | 0.000058 | 0.001160 | -0.000023 | 0.001162 | -0.000096 | 0.001196
RBE with $\bm{\alpha}$
$\bm{0.2}$ | 0.000047 | 0.001195 | -0.000036 | 0.001152 | 0.000013 | 0.001184 | -0.000049 | 0.001164
$\bm{0.4}$ | 0.000017 | 0.001164 | 0.000054 | 0.001136 | -0.000078 | 0.001172 | 0.000015 | 0.001137
$\bm{0.6}$ | 0.000053 | 0.001125 | 0.000051 | 0.001130 | -0.000028 | 0.001132 | 0.000020 | 0.001184
$\bm{0.8}$ | 0.000092 | 0.001137 | -0.000046 | 0.001163 | -0.0000074 | 0.001183 | 0.000061 | 0.001174
$\bm{1.0}$ | -0.000015 | 0.001136 | 0.000026 | 0.001116 | -0.000015 | 0.001100 | 0.000030 | 0.001137
$\bm{0.25}$ | -0.000031 | 0.001165 | 0.000113 | 0.001133 | 0.000069 | 0.001154 | 0.000003 | 0.001192

For testing of hypothesis, robust Bayes factor is used as test statistic for
any particular hypothesis for the given data. Let us define a simple null
hypothesis against an alternative hypothesis as

$\bm{H}_{0}:\bm{\theta}=\bm{\theta}_{0}\quad\text{against}\quad\bm{H}_{1}:\bm{\theta}\neq\bm{\theta}_{0}.$

A continuous prior density would lead to zero prior probability to test
$\bm{H}_{0}$. Therefore, it is suggestive to take an
$\varepsilon$-neighbourhood (spherical) around $\bm{\theta}_{0}$ and assign
prior probability $\rho_{0}$ under $\bm{H}_{0}$. The empirical prior and
posterior probabilities are calculated to obtain empirical Bayes factor. From
equation (27), the Bayes factor can be calculated using relation

$\text{Posterior odds}=\text{Prior odds}\times\text{Bayes factor}.$

Here, the simple null hypothesis under Weibull lifetime distribution is taken
as $\bm{\theta}^{(1)}_{0}=(-0.09,-0.06,0.2,0.7)^{{}^{\prime}}$ and
$\varepsilon=0.007.$ The values of empirical Bayes factor (BF01) for different
tuning parameters under Normal prior and Dirichlet prior are reported in table
13. The interpretation of Bayes factor values (BF01) can be made based on the
scale given in table 1. Since BF01 value lies in $20$ to $150$ under both
Normal and Dirichlet prior, support for $\bm{H}_{0}$ is strong. For Gompertz
lifetime distribution
$\bm{\theta}^{(2)}_{0}=(0.03,0.22,0.08,0.07)^{{}^{\prime}}$ and
$\varepsilon=0.0028$ is taken. Bayes factor values are provided in table 14.
From the interpretation given in table 1, support for $\bm{H}_{0}$ under
Normal prior is positive, whereas support for $\bm{H}_{0}$ under Dirichlet
prior is strong.

Table 13: Empirical value of Bayes factor with Weibull lifetime distribution. Tuning | Prior | Posterior | Bayes Factor
---|---|---|---
Parameter | odds | odds | BF01
Normal Prior
$\bm{0.2}$ | 0.231481 | 12.299990 | 53.135980
$\bm{0.4}$ | 09.783787 | 42.265970
$\bm{0.6}$ | 09.783787 | 42.265970
$\bm{0.8}$ | 08.279072 | 35.765600
$\bm{1.0}$ | 10.083330 | 43.560080
$\bm{0.65}$ | 15.624990 | 67.499980
Dirichlet Prior
$\bm{0.2}$ | 0.147126 | 09.499999 | 64.570320
$\bm{0.4}$ | 12.758620 | 86.718770
$\bm{0.6}$ | 08.975003 | 61.001990
$\bm{0.8}$ | 08.975003 | 61.001990
$\bm{1.0}$ | 12.758620 | 86.718770
$\bm{0.65}$ | 08.499996 | 57.773420
Table 14: Empirical value of Bayes factor with Gompertz lifetime distribution. Tuning | Prior | Posterior | Bayes Factor
---|---|---|---
Parameter | odds | odds | BF01
Normal Prior
$\bm{0.2}$ | 0.417613 | 5.835616 | 13.97374
$\bm{0.4}$ | 5.397437 | 12.92449
$\bm{0.6}$ | 4.870588 | 11.66292
$\bm{0.8}$ | 6.447763 | 15.43956
$\bm{1.0}$ | 6.920635 | 16.57189
$\bm{0.25}$ | 5.835616 | 13.97374
Dirichlet Prior
$\bm{0.2}$ | 0.128959 | 4.670455 | 36.21659
$\bm{0.4}$ | 5.480520 | 42.49816
$\bm{0.6}$ | 6.231882 | 48.32452
$\bm{0.8}$ | 5.835616 | 45.25172
$\bm{1.0}$ | 4.606742 | 35.72253
$\bm{0.25}$ | 6.128569 | 47.52339

##  10 Conclusion

This work has proposed a robust estimation in classical and Bayesian
frameworks under step-stress accelerated life testing experiment with
cumulative risk model where lifetime of one-shot device units comes from a
standard Lehman family of distributions. An intensive simulation study has
demonstrated robustness of minimum density power divergence estimation (MDPDE)
and robust Bayes estimation (RBE) over conventional maximum likelihood
estimation (MLE) and Bayes estimation (BE). The asymptotic property of MDPDE
is discussed, and an iteration-free method to find optimal tuning parameter is
suggested. Robust testing of hypotheses exploiting the Bayes factor has also
been conducted. Further, influence function to measure robustness analytically
has been derived. Finally, a data analysis has been conducted to establish
utility of the theoretical results developed in this work.

This work can be extended to the non-parametric approach for inferential
analysis. The step-stress model can be reanalysed under a competing risk set-
up. The missing cause of failure analysis can also be conducted. Efforts in
this direction are in the pipeline and we will report these findings soon.

## CRediT authorship contribution statement

Shanya Baghel: Conceptualization, Formal analysis, Methodology, Software,
Validation, Visualization, Writing - original draft.
Shuvashree Mondal: Conceptualization, Methodology, Supervision, Validation,
Visualization, Writing - review $\&$ editing.

## Declaration of competing interest

The authors declare no competing interests.

## Funding

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

##  Appendix A Proofs and description

###  A.1 Description of expressions in Result 1

$\displaystyle
A^{(i)}_{1}(\bm{\theta})=G^{{}^{\prime}(m)}_{\bm{\theta}}(\tau_{i};\gamma_{i})-G^{(m)}(\tau_{i};\gamma_{i})\sum_{l=1}^{i-1}E^{{}^{\prime}(\delta)}_{\bm{\theta}}(\tau_{l};\gamma_{l+1,l}).$
$\displaystyle
A^{(k)}_{2}(c_{j})=x^{(c_{j})}_{k}\lambda_{k}Q(\tau_{k};\gamma_{k})+\sum_{i=1}^{k-1}E^{{}^{\prime}(\delta)}_{c_{j}}(\tau_{i};\gamma_{i+1,i})\;;\;x^{(c_{0})}_{k}=1,x^{(c_{1})}_{k}=x_{k}.$
$\displaystyle
A^{(k)}_{2}(\gamma_{i})=\lambda_{i}\bigg{[}\Big{\\{}Q^{{}^{\prime}}_{\gamma_{i}}(\tau_{i};\gamma_{i})-Q^{{}^{\prime}}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i})\Big{\\}}+\frac{\delta}{2}\Big{\\{}Q^{{}^{\prime\prime}}_{\gamma_{i}}(\tau_{i};\gamma_{i})+Q^{{}^{\prime\prime}}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i})\Big{\\}}\bigg{]}.$
$\displaystyle
p^{{}^{\prime}(\bm{\theta})}_{i1}=\exp\bigg{[}-\Big{\\{}\lambda_{i-1}Q(\tau_{i-1};\gamma_{i-1})+\sum_{l=1}^{i-2}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})+D^{(\delta)}(\tau_{i-1}+\delta;\gamma_{i-1,i})\Big{\\}}\bigg{]}$
$\displaystyle\qquad\qquad
D^{{}^{\prime}(\delta)}_{\bm{\theta}}(\tau_{i-1}+\delta;\gamma_{i-1,i})+A^{(1)}_{3}(\bm{\theta})\exp\Big{\\{}-\sum_{l=1}^{i-1}E^{(\delta)}(\tau_{l};\gamma_{l+1,l})\Big{\\}}.$
$\displaystyle
E^{{}^{\prime}(\delta)}_{c_{j}}(\tau_{l};\gamma_{l+1,l})=x^{(c_{j})}_{l}\lambda_{l}\left\\{Q(\tau_{l};\gamma_{l})+\frac{\delta}{2}Q^{{}^{\prime}}(\tau_{l};\gamma_{l})\right\\}$
$\displaystyle\hskip
142.26378pt+x^{(c_{j})}_{l+1}\lambda_{l}\left\\{\frac{\delta}{2}Q^{{}^{\prime}}(\tau_{l}+\delta;\gamma_{l+1})-Q(\tau_{l}+\delta;\gamma_{l+1})+\right\\}.$
$\displaystyle
E^{{}^{\prime}(\delta)}_{\gamma_{i}}(\tau_{i-1};\gamma_{i,i-1})=\lambda_{i}\left\\{\frac{\delta}{2}Q^{{}^{\prime\prime}}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i})-Q^{{}^{\prime}}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i})\right\\}.$
$\displaystyle
A^{(1)}_{3}(\bm{\theta})=G^{{}^{\prime}(1)}_{\bm{\theta}}(\tau_{i-1,i};\gamma_{i})-G^{(1)}(\tau_{i-1,i};\gamma_{i})\sum_{l=1}^{i-1}E^{{}^{\prime}(\delta)}_{\bm{\theta}}(\tau_{l};\gamma_{l+1,l}).$
$\displaystyle
D^{{}^{\prime}(\delta)}_{c_{j}}(\tau_{i-1}+\delta;\gamma_{i-1,i})=\frac{\delta}{2}\Big{\\{}x^{(c_{j})}_{i-1}\lambda_{i-1}Q^{{}^{\prime}}(\tau_{i-1};\gamma_{i-1})+x^{(c_{j})}_{i}\lambda_{i}Q^{{}^{\prime}}(\tau_{i-1}+\delta;\gamma_{i})\Big{\\}}.$
$\displaystyle
D^{{}^{\prime}(\delta)}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i-1,i})=\frac{\delta}{2}\lambda_{i}Q^{{}^{\prime\prime}}_{\gamma_{i}}(\tau_{i-1}+\delta;\gamma_{i}).$

###  A.2 Proof of Result 3

For $M=k\sum_{i=1}^{k}q_{i}+1$, define $p_{im}=P_{r}$, $n_{im}=N_{r}$,
$p_{s}=P_{M}$ and $n_{s}=N_{M}$ where $r=q_{i}(i-1)+m$ for
$i=1,2,\dots,k\;;\;m=1,2,\dots,q_{i}$. Then, DPD measure ignoring terms
independent of parameters is

$D^{(M)}_{\alpha}(\bm{\theta})=\sum_{r=1}^{M}P^{\alpha+1}_{r}-\left(\frac{\alpha+1}{\alpha}\right)\sum_{r=1}^{M}\frac{N_{r}}{n}P^{\alpha}_{r}.$

Define, $X_{u}=(X_{u1},X_{u2},\dots,X_{uM})\sim
MN(1,P_{1},P_{2},\dots,P_{M})$. Therefore, $N_{r}=\sum_{u=1}^{n}X_{ur}$ and
$D^{(M)}_{\alpha}(\bm{\theta})$ is rewritten as

$D^{(M)}_{\alpha}(\bm{\theta})=\frac{1}{n}\sum_{u=1}^{n}V_{\alpha}(X_{u},\bm{\theta}),$

where,
$V_{\alpha}(X_{u},\bm{\theta})=\sum\limits_{r=1}^{M}\left\\{P^{\alpha+1}_{r}-\left(\frac{\alpha+1}{\alpha}\right)X_{ur}P^{\alpha}_{r}\right\\}$.
Denoting $\frac{\partial(D_{\alpha}^{(M)}(\bm{\theta}))}{\partial\bm{\theta}}$
$=D_{n}^{(\alpha)}(\bm{\theta})$, we get
$E\left[\frac{\partial}{\partial\bm{\theta}}V_{\alpha}(X_{u},\bm{\theta})\right]=0$
and

$\displaystyle
Var\left[\frac{\partial}{\partial\bm{\theta}}V_{\alpha}(X_{u},\bm{\theta})\right]$
$\displaystyle=(\alpha+1)^{2}\left\\{\sum_{r=1}^{M}P^{2\alpha-1}_{r}(1-P_{r})\left(\frac{\partial(P_{r})}{\partial\bm{\theta}}\right)^{2}\right\\}$
$\displaystyle\qquad-2\sum_{1\leq(h_{1},h_{2})\leq
M}P^{\alpha}_{r_{1}}P^{\alpha}_{r_{2}}\frac{\partial(P_{r_{1}})}{\partial\bm{\theta}}\frac{\partial(P_{r_{2}})}{\partial\bm{\theta}},$

and

$\displaystyle
Cov\left[\frac{\partial\big{\\{}V_{\alpha}(X_{u};\bm{\theta})\big{\\}}}{\partial\bm{\theta}_{h_{1}}},\frac{\partial\big{\\{}V_{\alpha}(X_{u};\bm{\theta})\big{\\}}}{\partial\bm{\theta}_{h_{1}}}\right]$
$\displaystyle{=}(\alpha+1)^{2}\left\\{\sum_{r=1}^{M}P^{2\alpha-1}_{r}(1-P_{r})\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}\right.$
$\displaystyle\Bigg{.}{-}2\sum_{1\leq(h_{1},h_{2})\leq
M}P^{\alpha}_{r_{1}}P^{\alpha}_{r_{2}}\frac{\partial(P_{r_{1}})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r_{2}})}{\partial\bm{\theta}_{h_{2}}}\Bigg{\\}};h_{1}{\neq}h_{2}.$

Define matrix $K_{\alpha}(\bm{\theta})$ such that

$\displaystyle K_{\alpha}(\bm{\theta})_{h,h}$
$\displaystyle=\frac{1}{(\alpha+1)^{2}}Var\left[\frac{\partial}{\partial\bm{\theta}}V_{\alpha}(X_{u},\bm{\theta})\right].$
$\displaystyle K_{\alpha}(\bm{\theta})_{h_{1},h_{2}}$
$\displaystyle=\frac{1}{(\alpha+1)^{2}}Cov\left[\frac{\partial\big{\\{}V_{\alpha}(X_{u};\bm{\theta})\big{\\}}}{\partial\bm{\theta}_{h_{1}}},\frac{\partial\big{\\{}V_{\alpha}(X_{u};\bm{\theta})\big{\\}}}{\partial\bm{\theta}_{h_{1}}}\right].$

Denote
$T_{\alpha}=-\sqrt{n}D_{n}^{(\alpha)}(\bm{\theta})=-\sqrt{n}\left.\frac{\partial(D_{\alpha}^{(M)})(\bm{\theta})}{\partial\bm{\theta}}\right|_{\bm{\theta}=\bm{\theta}_{0}}.$

Then, $T_{\alpha}\sim
N\bigg{(}\bm{0}_{k+2},(\alpha+1)^{2}K_{\alpha}(\bm{\theta_{0}})\bigg{)}$ (by
CLT). For

$\displaystyle\frac{\partial(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{2}}}$
$\displaystyle=\frac{\partial^{2}V_{\alpha}(X_{u},\bm{\theta})}{\partial
h_{1}\partial h_{2}}$ $\displaystyle=(\alpha+1)\sum_{r=1}^{M}\left[\alpha
P^{\alpha-1}_{r}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}+P^{\alpha}_{r}\frac{\partial^{2}(P_{r})}{\partial\bm{\theta}_{h1}\partial\bm{\theta}_{h2}}-(\alpha-1)X_{u_{r}}P^{\alpha-2}_{r}\right.$
$\displaystyle\hskip
71.13188pt\left.\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}-X_{ur}P^{(\alpha-1)}_{r}\frac{\partial^{2}(P_{r})}{\partial\bm{\theta}_{h1}\partial\bm{\theta}_{h2}}\right].$

Since $\frac{1}{n}\sum_{u=1}^{n}X_{ur}\xrightarrow[\text{}]{\text{ p}}P_{r}.$
Therefore,
$\frac{\partial(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{2}}}\xrightarrow[\text{}]{\text{
p}}(\alpha+1)\sum_{r=1}^{M}\left(P^{\alpha-1}_{r}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}\right).$
Consider $\bm{\theta}_{0}$ to be the true value of parameters, then by Taylor
series expansion ignoring higher-order terms

$\displaystyle D^{(\alpha)}_{n}(\bm{\theta})$
$\displaystyle=D^{(\alpha)}_{n}(\bm{\theta}_{0})+\left.\sum_{h_{1}=1}^{k+2}\frac{\partial(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}}\right|_{\bm{\theta=\theta_{0}}}(\hat{\bm{\theta}}_{h_{1}}-\bm{\theta}_{h_{1}0})$
$\displaystyle\hskip
71.13188pt+\left.\frac{1}{2}\sum_{h_{1}=1}^{k+2}\sum_{h_{2}=1}^{k+2}\frac{\partial^{2}(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}\partial\bm{\theta}_{h_{2}}}\right|_{\bm{\theta=\theta_{0}}}(\hat{\bm{\theta}}_{h_{1}}-\bm{\theta}_{h_{1}0})(\hat{\bm{\theta}}_{h_{2}}-\bm{\theta}_{h_{2}0}).$

As $D^{(\alpha)}_{n}(\hat{\bm{\theta}}_{\alpha})=0,$ therefore,

$\displaystyle-\sqrt{n}D^{(\alpha)}_{n}(\bm{\theta}_{0})$
$\displaystyle=\sqrt{n}\sum_{h_{1}=1}^{k+2}\left[\left.\frac{\partial(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}}\right|_{\bm{\theta=\theta_{0}}}\right.$
$\displaystyle\quad\Bigg{.}+\frac{1}{2}\left.\sum_{h_{2}=1}^{k+2}\frac{\partial^{2}(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}\partial\bm{\theta}_{h_{2}}}\right|_{\bm{\theta=\theta_{0}}}(\hat{\bm{\theta}}_{\alpha
h_{2}}-\bm{\theta}_{h_{2}0})\Bigg{]}(\hat{\bm{\theta}}_{\alpha
h_{1}}-\bm{\theta}_{h_{1}0}).$

Define,
$H_{(h_{1},h_{2})}=\Bigg{[}\left.\frac{\partial(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}}\right|_{\bm{\theta=\theta_{0}}}+\frac{1}{2}\left.\sum_{h_{2}=1}^{k+2}\frac{\partial^{2}(D^{(\alpha)}_{n}(\bm{\theta}))}{\partial\bm{\theta}_{h_{1}}\partial\bm{\theta}_{h_{2}}}\right|_{\bm{\theta=\theta_{0}}}(\hat{\bm{\theta}}_{\alpha
h_{2}}-\bm{\theta}_{h_{2}0})\Bigg{]},$ and it can be easily seen that,
$H_{(h_{1},h_{2})}\xrightarrow[\text{}]{\text{
p}}(\alpha+1)\sum_{r=1}^{M}\left(P^{\alpha-1}_{r}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}\right).$
Where we define $H_{\alpha}$ to be the $(k+2)\times(k+2)$ matrix with
$(h_{1},h_{2})th$ element $H_{(h_{1},h_{2})}.$ Thus

$J_{\alpha}(\bm{\theta}_{0})=\Bigg{[}\left(\sum_{r=1}^{M}P_{r}^{\alpha-1}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{1}}}\frac{\partial(P_{r})}{\partial\bm{\theta}_{h_{2}}}\right)_{h_{1},h_{2}}\Bigg{]}.$

Define, $Z_{h_{1}}=\sqrt{n}(\hat{\bm{\theta}}_{\alpha
h_{1}}-\bm{\theta}_{h_{1}0})$. Then
$-\sqrt{n}D^{(\alpha)}_{n}(\bm{\theta}_{0})=Z_{h_{1}}H_{(h_{1},h_{2})}$ which
implies $T_{\alpha}=Z_{\alpha}H_{\alpha}$. Hence
$Z_{\alpha}=H_{\alpha}^{-1}T_{\alpha}$, which results in

$\sqrt{n}(\hat{\bm{\theta}}_{\alpha}-\bm{\theta}_{0})=Z_{\alpha}\sim
N(\bm{0}_{k+2},J^{-1}_{\alpha}({\bm{\theta}_{0}})K_{\alpha}({\bm{\theta}_{0}})J^{-1}_{\alpha}({\bm{\theta}_{0}})).$

###  A.3 Proof of Result 5

$\text{Denote,}\quad
T^{(b)}_{\alpha}(F_{\bm{\theta}})=\frac{\int\bm{\theta}\exp\big{\\{}B_{\alpha}(\bm{\theta};F_{\bm{\theta}})\big{\\}}\pi(\bm{\theta})\,d\bm{\theta}}{\int\exp\big{\\{}B_{\alpha}(\bm{\theta};F_{\bm{\theta}})\big{\\}}\pi(\bm{\theta})\,d\bm{\theta}}=\frac{T_{1}(F_{\bm{\theta}})}{T_{2}(F_{\bm{\theta}})}.$
$\displaystyle\text{Then,}\quad IF(t;T_{\alpha}^{(b)},F_{\bm{\theta}})$
$\displaystyle=\left.\frac{\partial}{\partial\epsilon}T_{\alpha}^{(b)}(U_{\epsilon})\right|_{\epsilon\to
0^{+}}$
$\displaystyle=\left.\frac{T_{2}(U_{\epsilon})\frac{\partial}{\partial\epsilon}T_{1}(U_{\epsilon})-T_{1}(U_{\epsilon})\frac{\partial}{\partial\epsilon}T_{2}(U_{\epsilon})}{\left\\{T_{2}(U_{\epsilon})\right\\}^{2}}\right|_{\epsilon\to
0^{+}}$
$\displaystyle=\frac{\int\bm{\theta}X_{\alpha}(\bm{\theta};t,f_{\bm{\theta}})\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}{\int\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}$
$\displaystyle\qquad-\left[\frac{\int\bm{\theta}\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}{\int\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}\right.$
$\displaystyle\qquad\qquad\left.\times\frac{\int
X_{\alpha}(\bm{\theta};t,f_{\bm{\theta}})\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}{\int\exp\left\\{B_{\alpha}(\bm{\theta})\right\\}\pi(\bm{\theta})d\bm{\theta}}\right]$
$\displaystyle=Cov_{(p)}\left(\bm{\theta},X_{\alpha}(\bm{\theta};t,f_{\bm{\theta}})\right),$

where, $Cov_{(p)}()$ is the covariance for posterior distribution and

$\displaystyle X_{\alpha}(\bm{\theta};t,f_{\bm{\theta}})$
$\displaystyle=\left.\frac{\partial(B_{\alpha}(\bm{\theta};F_{\bm{\theta}}))}{\partial\epsilon}\right|_{\epsilon\to
0^{+}}$
$\displaystyle=\frac{1}{\alpha}\left[\left\\{\Delta^{(I_{s})}_{t}-p_{s}(\bm{\theta})\right\\}p^{\alpha}_{s}(\bm{\theta})+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left\\{\Delta^{(I_{im})}_{t}-p_{im}(\bm{\theta})\right\\}p^{\alpha}_{im}(\bm{\theta})\right].$
$\displaystyle\Delta^{(I)}_{t}$ $\displaystyle=\begin{cases}1\quad\text{if}\
t\in I\\\ 0\quad\text{otherwise}\\\ \end{cases}.$

###  A.4 Proof of Result 6

$\displaystyle\text{Denote,}\quad
T^{(\alpha)}_{\bm{\Theta}}(F_{\bm{\theta}_{0}})=\frac{\int_{\bm{\Theta}_{0}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0};F_{\bm{\theta}_{0}})\big{\\}}\pi_{0}(\bm{\theta})\,d\bm{\theta}}{\int_{\bm{\Theta}_{1}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1};F_{\bm{\theta}_{1}})\big{\\}}\pi_{1}(\bm{\theta})\,d\bm{\theta}}=\frac{T_{0}(\bm{\theta}\in\bm{\Theta}_{0})}{T_{1}(\bm{\theta}\in\bm{\Theta}_{1})}.$
$\displaystyle\text{Then,}\;IF(t;T^{(\alpha)}_{\bm{\Theta}},F_{\bm{\theta}_{0}})$
$\displaystyle=\left.\frac{\partial(T^{(\alpha)}_{\bm{\Theta}}(U_{\epsilon}))}{\partial\epsilon}\right|_{\epsilon\to
0^{+}}.$
$\displaystyle=\left[\frac{\int_{\bm{\Theta}_{0}}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\big{\\}}\pi_{0}(\bm{\theta})\,d\bm{\theta}}{\int_{\bm{\Theta}_{0}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\big{\\}}\pi_{0}(\bm{\theta})\,d\bm{\theta}}\times
Y_{\alpha}(\bm{\Theta})\right]$
$\displaystyle\qquad-\left[Y_{\alpha}(\bm{\Theta})\times\frac{\int_{\bm{\Theta}_{1}}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\big{\\}}\pi_{1}(\bm{\theta})\,d\bm{\theta}}{\int_{\bm{\Theta}_{1}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\big{\\}}\pi_{1}(\bm{\theta})\,d\bm{\theta}}\right]$
$\displaystyle=Y_{\alpha}(\bm{\Theta})\bigg{\\{}E\Big{[}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\Big{]}-E\Big{[}X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\Big{]}\bigg{\\}},$

where,

$\displaystyle Y_{\alpha}(\bm{\Theta})$
$\displaystyle=\frac{\int_{\bm{\Theta}_{0}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{0})\big{\\}}\pi_{0}(\bm{\theta})\,d\bm{\theta}}{\int_{\bm{\Theta}_{1}}\exp\big{\\{}B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{1})\big{\\}}\pi_{1}(\bm{\theta})\,d\bm{\theta}}.$
$\displaystyle X_{\alpha}(\bm{\theta}\in\bm{\Theta}_{j})$
$\displaystyle=\left.\frac{\partial(B_{\alpha}(\bm{\theta}\in\bm{\Theta}_{j};F_{\bm{\theta}_{0}}))}{\partial\epsilon}\right|_{\epsilon\to
0^{+}}\;;\;j=0,1.$
$\displaystyle=\frac{1}{\alpha}\left[\left\\{\Delta^{(I_{s})}_{t}-p_{s}(\bm{\theta}_{0})\right\\}p^{\alpha}_{s}(\bm{\theta}\in\bm{\Theta}_{j})\right.$
$\displaystyle\qquad\qquad\left.+\sum_{i=1}^{k}\sum_{m=1}^{q_{i}}\left\\{\Delta^{(I_{im})}_{t}-p_{im}(\bm{\theta}_{0})\right\\}p^{\alpha}_{im}(\bm{\theta}\in\bm{\Theta}_{j})\right].$

## References

  * [1] K. Song, L. Cui, Fiducial inference-based failure mechanism consistency analysis for accelerated life and degradation tests, Applied Mathematical Modelling 105 (2022) 340–354.  doi:https://doi.org/10.1016/j.apm.2021.12.048.
  * [2] M. Khakifirooz, M. Fathi, I.-C. Lee, S.-T. Tseng, Neural ordinary differential equation for sequential optimal design of fatigue test under accelerated life test analysis, Reliability Engineering & System Safety 235 (2023) 109242.
  * [3] M. Kateri, N. I. Nikolov, Product of spacings estimation in step-stress accelerated life testing: An alternative to maximum likelihood, IEEE Transactions on Reliability (2024).  doi:10.1109/TR.2024.3369977.
  * [4] C. Chen, B. Li, J. Guo, Z. Liu, B. Qi, C. Hua, Bearing life prediction method based on the improved fides reliability model, Reliability Engineering & System Safety 227 (2022) 108746.  doi:https://doi.org/10.1016/j.ress.2022.108746.
  * [5] W. Yan, X. Xu, D. Bigaud, W. Cao, Optimal design of step-stress accelerated degradation tests based on the tweedie exponential dispersion process, Reliability Engineering & System Safety 230 (2023) 108917.  doi:https://doi.org/10.1016/j.ress.2022.108917.
  * [6] Y. Li, H. Gao, H. Chen, C. Liu, Z. Yang, E. Zio, Accelerated degradation testing for lifetime analysis considering random effects and the influence of stress and measurement errors, Reliability Engineering & System Safety (2024) 110101doi:https://doi.org/10.1016/j.ress.2024.110101.
  * [7] N. Balakrishnan, M. H. Ling, Gamma lifetimes and one-shot device testing analysis, Reliability Engineering & System Safety 126 (2014) 54–64.  doi:https://doi.org/10.1016/j.ress.2014.01.009.
  * [8] N. Balakrishnan, H. Y. So, M. H. Ling, Em algorithm for one-shot device testing with competing risks under exponential distribution, Reliability Engineering & System Safety 137 (2015) 129–140.  doi:https://doi.org/10.1016/j.ress.2014.12.014.
  * [9] S.-J. Wu, C.-C. Hsu, S.-R. Huang, Optimal designs and reliability sampling plans for one-shot devices with cost considerations, Reliability Engineering & System Safety 197 (2020) 106795.  doi:https://doi.org/10.1016/j.ress.2020.106795.
  * [10] X. Zhu, K. Liu, M. He, N. Balakrishnan, Reliability estimation for one-shot devices under cyclic accelerated life-testing, Reliability Engineering & System Safety 212 (2021) 107595.  doi:https://doi.org/10.1016/j.ress.2021.107595.
  * [11] X. Zhu, N. Balakrishnan, One-shot device test data analysis using non-parametric and semi-parametric inferential methods and applications, Reliability Engineering & System Safety 221 (2022) 108319.  doi:https://doi.org/10.1016/j.ress.2022.108319.
  * [12] M. H. Ling, Optimal design of simple step-stress accelerated life tests for one-shot devices under exponential distributions, Probability in the Engineering and Informational Sciences 33 (1) (2019) 121–135.  doi:https://doi.org/10.1017/S0269964818000049.
  * [13] M. H. Ling, X. Hu, Optimal design of simple step-stress accelerated life tests for one-shot devices under weibull distributions, Reliability Engineering & System Safety 193 (2020) 106630.  doi:https://doi.org/10.1016/j.ress.2019.106630.
  * [14] N. Balakrishnan, E. Castilla, M. Jaenada, L. Pardo, Robust inference for nondestructive one-shot device testing under step-stress model with exponential lifetimes, Quality and Reliability Engineering International 39 (4) (2023) 1192–1222.  doi:https://doi.org/10.1002/qre.3287.
  * [15] N. Balakrishnan, M. Jaenada, L. Pardo, Non-destructive one-shot device test under step-stress experiment with lognormal lifetime distribution, Journal of Computational and Applied Mathematics 437 (2024) 115483.  doi:https://doi.org/10.1016/j.cam.2023.115483.
  * [16] N. Balakrishnan, M. Jaenada, L. Pardo, Step-stress tests for interval-censored data under gamma lifetime distribution, Quality Engineering 36 (1) (2024) 3–20.  doi:https://doi.org/10.1080/08982112.2023.2199826.
  * [17] J. René Van Dorp, T. A. Mazzuchi, A general bayes exponential inference model for accelerated life testing, Journal of Statistical Planning and Inference 119 (1) (2004) 55–74.  doi:https://doi.org/10.1016/S0378-3758(02)00411-1.
  * [18] N. Kannan, D. Kundu, N. Balakrishnan, Survival models for step-stress experiments with lagged effects, Advances in Degradation Modeling: Applications to Reliability, Survival Analysis, and Finance (2010) 355–369.
  * [19] J.-y. Yao, R.-m. Luo, Step-stress accelerated degradation test model of storage life based on lagged effect for electronic products, in: The 19th International Conference on Industrial Engineering and Engineering Management: Engineering Management, Springer, 2013, pp. 541–550.  doi:10.1007/978-3-642-38433-2-59.
  * [20] J. Beltrami, Weibull lagged effect step-stress model with competing risks, Communications in Statistics-Theory and Methods 46 (11) (2017) 5419–5442.  doi:https://doi.org/10.1080/03610926.2015.1102283.
  * [21] N. Kannan, D. Kundu, Simple step-stress models with a cure fraction, Brazilian Journal of Probability and Statistics 34 (1) (2020) 2–17.  doi:https://www.jstor.org/stable/26924216.
  * [22] Y. Qiao, W. Gui, Inference for cumulative risk model under step-stress experiments and its application in nanocrystalline data, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability 237 (1) (2023) 195–209.  doi:https://doi.org/10.1177/1748006X211072643.
  * [23] A. Basu, I. R. Harris, N. L. Hjort, M. Jones, Robust and efficient estimation by minimising a density power divergence, Biometrika 85 (3) (1998) 549–559.  doi:https://doi.org/10.1093/biomet/85.3.549.
  * [24] N. Balakrishnan, E. Castilla, N. Martín, L. Pardo, Robust estimators for one-shot device testing data under gamma lifetime model with an application to a tumor toxicological data, Metrika 82 (8) (2019) 991–1019.  doi:https://doi.org/10.1007/s00184-019-00718-5.
  * [25] N. Balakrishnan, E. Castilla, N. Martín, L. Pardo, Robust inference for one-shot device testing data under exponential lifetime model with multiple stresses, Quality and Reliability Engineering International 36 (6) (2020) 1916–1930.  doi:https://doi.org/10.1002/qre.2665.
  * [26] E. M. Castilla Gonzalez, Robust statistical inference for one-shot devices based on divergences, Universidad Complutense de Madrid, New York, 2021.
  * [27] N. Balakrishnan, M. H. Ling, H. Y. So, Accelerated life testing of one-shot devices: Data collection and analysis, John Wiley & Sons, 2021.
  * [28] S. Baghel, S. Mondal, Analysis of one-shot device testing data under logistic-exponential lifetime distribution with an application to seer gallbladder cancer data, Applied Mathematical Modelling 126 (2024) 159–184.  doi:https://doi.org/10.1016/j.apm.2023.10.037.
  * [29] E. Castilla, P. J. Chocano, On the choice of the optimal tuning parameter in robust one-shot device testing analysis, Trends in Mathematical, Information and Data Sciences: A Tribute to Leandro Pardo (2022) 169–180doi:10.1007/978-3-031-04137-2-16.
  * [30] S. Shuto, T. Amemiya, Sequential bayesian inference for weibull distribution parameters with initial hyperparameter optimization for system reliability estimation, Reliability Engineering & System Safety 224 (2022) 108516.  doi:https://doi.org/10.1016/j.ress.2022.108516.
  * [31] C. Zhang, L. Wang, X. Bai, J. Huang, Bayesian reliability analysis for copula based step-stress partially accelerated dependent competing risks model, Reliability Engineering & System Safety 227 (2022) 108718.  doi:https://doi.org/10.1016/j.ress.2022.108718.
  * [32] Z. Saberzadeh, M. Razmkhah, M. Amini, Bayesian reliability analysis of complex k-out-of-n:$\ell$ systems under degradation performance, Reliability Engineering & System Safety 231 (2023) 109020.  doi:https://doi.org/10.1016/j.ress.2022.109020.
  * [33] F. Hong, P. Wei, J. Fu, M. Beer, A sequential sampling-based bayesian numerical method for reliability-based design optimization, Reliability Engineering & System Safety 244 (2024) 109939.  doi:https://doi.org/10.1016/j.ress.2024.109939.
  * [34] R. Giacomini, T. Kitagawa, M. Read, Robust bayesian inference in proxy svars, Journal of Econometrics 228 (1) (2022) 107–126.  doi:https://doi.org/10.1016/j.jeconom.2021.02.003.
  * [35] L. Ouyang, S. Zhu, K. Ye, C. Park, M. Wang, Robust bayesian hierarchical modeling and inference using scale mixtures of normal distributions, IISE Transactions 54 (7) (2022) 659–671.  doi:https://doi.org/10.1080/24725854.2021.1912440.
  * [36] H. J. Park, N. H. Kim, J.-H. Choi, A robust health prediction using bayesian approach guided by physical constraints, Reliability Engineering & System Safety 244 (2024) 109954.  doi:https://doi.org/10.1016/j.ress.2024.109954.
  * [37] A. Ghosh, A. Basu, Robust bayes estimation using the density power divergence, Annals of the Institute of Statistical Mathematics 68 (2016) 413–437.  doi:https://doi.org/10.1007/s10463-014-0499-0.
  * [38] H. Jeffreys, The theory of probability, OuP Oxford, 1998.
  * [39] J. Chan, I. Papaioannou, D. Straub, An adaptive subset simulation algorithm for system reliability analysis with discontinuous limit states, Reliability Engineering & System Safety 225 (2022) 108607.  doi:https://doi.org/10.1016/j.ress.2022.108607.
  * [40] B. Abba, H. Wang, H. S. Bakouch, A reliability and survival model for one and two failure modes system with applications to complete and censored datasets, Reliability Engineering & System Safety 223 (2022) 108460.
  * [41] T. T. Thach, R. Bris, Reparameterized weibull distribution: a bayes study using hamiltonian monte carlo, in: Proceedings of the 29th European safety and reliability conference. Singapore: Research Publishing, 2019, pp. 997–1004.  doi:10.3850/978-981-11-2724-3-0494-cd.
  * [42] T. Thanh Thach, R. Briš, An additive chen-weibull distribution and its applications in reliability modeling, Quality and Reliability Engineering International 37 (1) (2021) 352–373.  doi:https://doi.org/10.1002/qre.2740.
  * [43] B. Abba, H. Wang, A new failure times model for one and two failure modes system: A bayesian study with hamiltonian monte carlo simulation, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability (2023) 1748006X221146367doi:https://doi.org/10.1177/1748006X221146367.
  * [44] R. M. Neal, Bayesian learning for neural networks, Springer New York, NY, 1996.  doi:https://doi.org/10.1007/978-1-4612-0745-0.
  * [45] R. M. Neal, et al., Mcmc using hamiltonian dynamics, Handbook of markov chain monte carlo 2 (11) (2011) 2.
  * [46] C. C. Monnahan, J. T. Thorson, T. A. Branch, Faster estimation of bayesian models in ecology using hamiltonian monte carlo, Methods in Ecology and Evolution 8 (3) (2017) 339–348.  doi:https://doi.org/10.1111/2041-210X.12681.
  * [47] A. Calvino, N. Martin, L. Pardo, Robustness of minimum density power divergence estimators and wald-type test statistics in loglinear models with multinomial sampling, Journal of Computational and Applied Mathematics 386 (2021) 113214.  doi:https://doi.org/10.1016/j.cam.2020.113214.
  * [48] T.-H. Fan, N. Balakrishnan, C.-C. Chang, The bayesian approach for highly reliable electro-explosive devices using one-shot device testing, Journal of Statistical Computation and Simulation 79 (9) (2009) 1143–1154.  doi:https://doi.org/10.1080/00949650802142592.
  * [49] H. L. Lee, C. Morris A., A multinomial logit model for the spatial distribution of hospital utilization, Journal of Business & Economic Statistics 3 (2) (1985) 159–168.  doi:https://doi.org/10.1080/07350015.1985.10509445.
  * [50] T. T. Thach, R. Bris, Improved new modified weibull distribution: A bayes study using hamiltonian monte carlo simulation, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability 234 (3) (2020) 496–511.  doi:https://doi.org/10.1177/1748006X19896740.
  * [51] S. Thomas, W. Tu, Learning hamiltonian monte carlo in r, The American Statistician 75 (4) (2021) 403–413.  doi:https://doi.org/10.1080/00031305.2020.1865198.
  * [52] J. K. Ghosh, M. Delampady, T. Samanta, An introduction to Bayesian analysis: theory and methods, Vol. 725, Springer, 2006.
  * [53] R. E. Kass, A. E. Raftery, Bayes factors, Journal of the american statistical association 90 (430) (1995) 773–795.  doi:https://doi.org/10.1080/01621459.1995.10476572.
  * [54] N. Balakrishnan, E. Castilla, N. Martín, L. Pardo, Robust inference for one-shot device testing data under weibull lifetime model, IEEE transactions on Reliability 69 (3) (2019) 937–953.  doi:10.1109/TR.2019.2954385.
  * [55] N. Balakrishnan, M. Jaenada, L. Pardo, Non-destructive one-shot device testing under step-stress model with weibull lifetime distributions, arXiv preprint arXiv:2208.02674 (2022).  doi:https://doi.org/10.48550/arXiv.2208.02674.
  * [56] A. M. Almarashi, G. A. Abd-Elmougod, Accelerated competing risks model from gompertz lifetime distributions with type-ii censoring scheme, Thermal Science 24 (Suppl. 1) (2020) 165–175.  doi:https://doi.org/10.2298/TSCI20S1165A.
  * [57] M. A. Milevsky, M. A. Milevsky, Life and death in continuous time: Gompertz 101, Retirement Income Recipes in R: From Ruin Probabilities to Intelligent Drawdowns (2020) 131–157doi:https://doi.org/10.1007/978-3-030-51434-1-8.
  * [58] Y. Zhu, Optimal design and equivalency of accelerated life testing plans, Rutgers The State University of New Jersey, School of Graduate Studies, 2010.  doi:https://www.proquest.com/dissertations-theses/optimal-design-equivalency-accelerated-life/docview/847547908/se-2?accountid=49706.