# One-Shot Device Testing Data Analysis under Logistic-Exponential Lifetimes
with an Application to SEER Gallbladder Cancer Data

Shanya Baghela Shuvashree Mondalb<EMAIL_ADDRESS>
###### Abstract

In the literature, the reliability analysis of one-shot devices is found under
accelerated life testing in the presence of various stress factors. The
application of one-shot devices can be extended to the bio-medical field,
where we often evidence that inflicted with a certain disease, survival time
would be under different stress factors like environmental stress, co-
morbidity, the severity of disease etc. This work is concerned with a one-shot
device data analysis and applies it to SEER Gallbladder cancer data. The two-
parameter logistic exponential distribution is applied as a lifetime
distribution. For robust parameter estimation, weighted minimum density power
divergence estimators (WMDPDE) is obtained along with the conventional maximum
likelihood estimators (MLE). The asymptotic behaviour of the WMDPDE and the
robust test statistic based on the density power divergence measure are also
studied. The performances of estimators are evaluated through extensive
simulation experiments. Later those developments are applied to SEER
Gallbladder cancer data. Citing the importance of knowing exactly when to
inspect the one-shot devices put to the test, a search for optimum inspection
times is performed. This optimization is designed to minimize a defined cost
function which strikes a trade-off between the precision of the estimation and
experimental cost. The search is accomplished through the population-based
heuristic optimization method Genetic Algorithm.

###### keywords:

Density Power Divergence Estimator, Genetic Algorithm , Kullback-Leibler
divergence, Logistic-Exponential distribution, One-Shot Devices.

###### MSC:

62F10, 62F12, 62NO2.

11footnotemark: 1

organization=Department of Mathematics and Computing, Indian Institute of
Technology (Indian School of Mines) Dhanbad, postcode=826004, state=Jharkhand,
country=India

organization=Department of Mathematics and Computing, Indian Institute of
Technology (Indian School of Mines) Dhanbad, postcode=826004, state=Jharkhand,
country=India

A Robust inference based on weighted density power divergence

Elaborated study of the asymptotic property

Testing of hypothesis based on a robust statistic

Search for optimum inspection times using Genetic algorithm

Application to SEER Gallbladder Cancer Data

##  1 Introduction

Applications of one-shot devices prevail in a broad spectrum of life. One-shot
devices stay in a torpid state until activated and are immediately destroyed
after the operation. Amidst the COVID-19 pandemic, extensive use of PPE kits,
RT-PCR and Rapid Antigen testing kits etc., is evidenced, which are only for
one-time use. Apart from these, the sugar level testing strips, pregnancy
testing kits, fire extinguishers, fuses, fuel injectors, missiles, nuclear
weapons, and space probes all come into the category of one-shot devices. The
reliability analysis of such devices becomes quite challenging as the
defectiveness of the devices can be discovered only after being tested, and
the exact failure times of such devices cannot be recorded accurately in most
instances. Hence it can only be observed whether the device’s failure occurs
before or after an inspection time, which yields dichotomous data. With highly
reliable devices, typically accelerated life tests (ALTs) are executed to
conduct the reliability testing of these devices under a limited time and
experimental budget. With high-stress factors under ALT, reliability is
estimated, which is extrapolated to real-life operating situations. Readers
may refer to Wang [1], Xu et al. [2], Roy [3], Srivastava and Jain [4] for the
widespread use of ALT in reliability analysis. The application of one-shot
devices can be extended in the bio-medical field, where diseases like cancer
and pulmonary infection are often exhibited with different stress factors like
environmental stress, lifestyle stress, co-morbidity and severity of disease
etc.

Reviewing the literature for one-shot devices shows that Exponential, Gamma,
and Weibull distributions are commonly used as lifetime distributions.
Balakrishnan et al. [5] developed an EM algorithm for point estimation under
ALT for an exponentially distributed lifetime of one-shot devices. They also
compared it with the Bayesian approach using normal prior developed by Fan et
al. [6]. Balakrishnan et al. [7] studied the reliability of one-shot devices
in the presence of competing risk factors. Balakrishnan and Ling [8] applied
gamma distributions as a lifetime of one-shot devices and provided an
inference study. Balakrishnan and Castilla [9] studied the reliability of one-
shot devices under log-normal distribution and provided an EM algorithm for
estimation purposes under constant stress accelerated life tests.

Introduced by Lam and Leemis [10] as a survival distribution, a two-parameter
Logistic-Exponential (LE) distribution is chosen in the present work for
depicting the lifetime of the one-shot unit model because of its high
flexibility regarding the shape of the hazard function. It can be observed
that it is the only two parameter distribution which demonstrates five
different hazard rate shapes namely, constant, increasing, decreasing, bath-
tub and upside-down-bath-tub [11]. Figure 1 illustrates the various hazard
rate shapes of LE distribution for different values of shape and scale
parameters. As a result, LE distribution can be used as an alternative to some
well-known two-parameter models like gamma, log-normal, Weibull, exponentiated
exponential, inverse Gaussian, Birnbaum–Saunders distributions [12] and in
terms of model fitting, it might work nicely than the stated distributions in
some situations. In the present work, to incorporate the stress factors, the
shape and scale parameters of the two-parameter Logistic-Exponential
distribution are linked with the stress factors through a log-linear function.

(a) $\alpha=1,\beta=0.5$

(b) $\alpha=0.8,\beta=0.7$

(c) $\alpha=4,\beta=0.2$

(d) $\alpha=0.6,\beta=5$

(e) $\alpha=1.5,\beta=4$

(f) $\alpha=2,\beta=0.2$

Figure 1: Hazard function of LE distribution for different choices of shape
and scale parameters

For estimation of the lifetime distribution, though maximum likelihood
estimation (MLE) is a prevalent and common method because of its well-known
properties such as asymptotic efficiency, consistency, sufficiency, and
invariance, it may not possess the property of robustness. Hence, in this work
a robust estimator based on density power divergence measure is used for
estimating the lifetime distribution as an alternative. Minimum density power
divergence (MDPD) estimator was first proposed by Basu et.al [13]. In this
work, certain weights are imposed on the density power divergence measure and
the corresponding estimator is termed as weighted minimum density power
divergence estimator (WMDPDE). The literature shows that Balakrishnan et
al.[14, 15, 16] observed the robustness of WMDPD estimators for analysing one-
shot device testing data under ALT based on gamma, exponential and Weibull
lifetime distributions.

In this work, the asymptotic distribution of the WMDPD estimator is derived.
Inspired by the idea of Basu et al. [17], a testing of hypothesis based on
WDPD is also developed in which robust WDPD-based test statistic is presented
and an approximation to its power function is determined. The WMDPDEs and the
test statistic based on the WDPD measure satisfy robustness without
compromising efficiency.

Another essential contribution of this work is in the design aspect of the
life testing experiment, which consists of finding the optimal inspection
times. Optimal inspection times are the set of inspection times that will
optimize certain cost functions, which may be based on maximizing the
precision of the estimation, minimizing the experimental cost, or both. In the
literature, it is found that Wu et al. [18] utilized the D-optimality
criterion with cost constraints to determine the optimal inspection times.
Ling and Hu [19] minimized the asymptotic variance of MLE under normal
operating conditions for Weibull distributions concerning sample allocations
and inspection times for determining the optimal designs of simple ALT for
one-shot devices. Balakrishnan and Castilla [9] studied optimal design under
constant stress accelerated life tests for one-shot devices with budget
constraints.

Further, a cost function is proposed based on minimizing the determinants of
the asymptotic covariance matrix of the WMDPD estimator and the number of
expected failures. Population-based heuristic optimization method Genetic
Algorithm (GA) in continuous design space is used in finding optimal
inspection times. The immense popularity of GA lies on the merits of its easy
understanding of the concept, ability to converge global or near global optima
reasonably well, to avoid the trap of local optima, to work without
computation of complex mathematical derivatives, to deal with various
constraints, to use probabilistic transition rules, to show efficiency even
when many parameters are involved.

The rest of the work proceeds as follows: Section (2) comprises model
description and likelihood function computation. Density power divergence
(DPD) measure with the weighted estimating equations, its asymptotic
properties and hypothesis testing, and the power function are discussed in
Section (3). Section (4) is comprised of an extensive simulation study to
assess the performances of the derived estimators. The cost function is
defined in section (5), and the search for the optimum inspection times
applying GA is provided. In section (6), outcomes of previous sections are
applied to SEER Gallbladder Cancer Data [20]. The study is concluded with a
summary discussion in section (7) along with some future insights.

##  2 Model Description

In the context of life testing experiment of the one-shot devices, the devices
are set down to $I$ observation groups where the number of devices in $i$th
group is $k_{i}.$ These devices  are subject to $J$ type of stress factors,
quantified by $x_{ij};j=0,1,2,...,J,$ where $x_{i0}=1.$ In the $i$th
observation group, the devices are inspected at a fixed time point say
$\tau_{i}$, and the number of devices fail before $\tau_{i}$, denoted by
$n_{i},$ are recorded for $i=1,\ldots,I.$ Therefore, the observed data can be
represented in tabular form as given in Table (1).

Table 1: Model layout and the data representation |  |  |  | Covariates
---|---|---|---|---
Groups | Inspection | Mice/ | Deaths/ | Stress 1 | … | Stress J
| Times | Devices | Failures |  |  |
1 | $\tau_{1}$ | $k_{1}$ | $n_{1}$ | $x_{11}$ | … | $x_{1J}$
2 | $\tau_{2}$ | $k_{2}$ | $n_{2}$ | $x_{21}$ | … | $x_{2J}$
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮
I | $\tau_{I}$ | $k_{I}$ | $n_{I}$ | $x_{I1}$ | … | $x_{IJ}$

In this work, the lifetime $T$ of a device is assumed to follow two parameters
Logistic-Exponential distribution with shape parameter $\alpha$ and scale
parameter $\lambda$. The cumulative distribution function and probability
density function of $T$ are given as follows:

$\displaystyle F(t)=$ $\displaystyle 1-(1+(e^{\lambda
t}-1)^{\alpha})^{-1}\;;t>0,\ \alpha,\lambda>0,$ $\displaystyle f(t)=$
$\displaystyle\frac{\alpha\lambda e^{\lambda t}(e^{\lambda
t}-1)^{\alpha-1}}{(1+(e^{\lambda t}-1)^{\alpha})^{2}};\;t>0,\
\alpha,\lambda>0.$

Within each observation group, shape and scale parameters are assumed to have
log-linear link with the stress factors as follows:

$\displaystyle\alpha_{i}$
$\displaystyle=\exp\left\\{\sum_{j=0}^{J}a_{j}x_{ij}\right\\}$
$\displaystyle\lambda_{i}$
$\displaystyle=\exp\left\\{\sum_{j=0}^{J}b_{j}x_{ij}\right\\}$

where $x_{i0}=1$ for $i=1,\ldots,I.$ Denote,
$\bm{\theta}=\\{a_{j},b_{j};j=0,...,J\\}$ for the model parameters to be
estimated.

The likelihood function on the basis of given observed data can be obtained as

$\displaystyle L(\bm{\theta})\propto$
$\displaystyle\prod_{i=1}^{I}F(\tau_{i};x_{ij},\bm{\theta})^{n_{i}}(1-F(\tau_{i};x_{ij},\bm{\theta}))^{k_{i}-n_{i}}$
$\displaystyle=$
$\displaystyle\prod_{i=1}^{I}\left[1-\left\\{1+\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}\right\\}^{-1}\right]^{n_{i}}\left[\left\\{1+\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}\right\\}^{-1}\right]^{k_{i}-n_{i}}.$
(1)

Therefore, the log-likelihood function without the normalizing constant can be
obtained as,

$\displaystyle ln\;L(\bm{\theta})\propto$
$\displaystyle\sum_{i=1}^{I}n_{i}\;ln\left[\frac{\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}}{1+\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}}\right]-(k_{i}-n_{i})\;ln\left[1+\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}\right]$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{I}\left[n_{i}\alpha_{i}\;ln\left(e^{\lambda_{i}\tau_{i}}-1\right)-k_{i}\;ln\left\\{1+\left(e^{\lambda_{i}\tau_{i}}-1\right)^{\alpha_{i}}\right\\}\right].$
(2)

Hence the MLE of $\bm{\theta},$ say
$\;\hat{\bm{\theta}}=\\{\hat{a}_{j},\hat{b}_{j};j=0,...,J\\}$ would be derived
as

$\hat{\bm{\theta}}=arg\mathop{max}_{\bm{\theta}}lnL(\bm{\theta})$ (3)

provided $\sum_{i=1}^{I}n_{i}>0.$

######  Theorem 1.

To obtain the MLE, the set of estimating equations is given as follows,
for the parameter $a_{j}$,

$\sum_{i=1}^{I}\alpha_{i}x_{ij}\;ln(e^{\lambda_{i}\tau_{i}}-1)\left[n_{i}-\frac{k_{i}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}}{\left\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\right\\}}\right]=0$
(4)

and for the parameter $b_{j}$,

$\sum_{i=1}^{I}\frac{\hat{\alpha}_{i}\lambda_{i}\tau_{i}x_{ij}e^{\lambda_{i}\tau_{i}}}{(e^{\lambda_{i}\tau_{i}}-1)}\left[n_{i}-\frac{k_{i}(e^{\lambda_{i}\tau_{i}}-1)^{\hat{\alpha}_{i}}}{\left\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\hat{\alpha}_{i}}\right\\}}\right]=0.$
(5)

Since MLE is not robust, we discuss a robust method of estimation in the
following section.

##  3 Density Power Divergence (DPD) Measure

The minimum divergence estimation method, namely density power divergence for
robust parameter estimation, was developed by Basu et al. [13]. Density power
divergence measure between any two probability distributions say, F and G on
the same variable with density respectively, f and g can be obtained as:

$\displaystyle
D_{\beta}(g,f)=\int\left\\{f^{\beta+1}(x)-\frac{\beta+1}{\beta}g(x)f^{\beta}(x)+\frac{1}{\beta}g^{\beta+1}(x)\right\\}dx$
(6)

where $\beta$ $(0\leq\beta\leq 1)$ is called the tuning parameter.

It is simple to show that, $\quad\mathop{lim}_{\beta\to
0^{+}}D_{\beta}\left(g,f\right)=D_{KL}\left(g|f\right).$
Here $D_{KL}(g|f)$ is the Kullback-Leibler (KL) divergence measure which is
the measure of the information lost when a probability distribution F is used
to approximate another probability distribution G on the same variable.

Based on the dichotomous data for one-shot devices, the DPD measure between
empirical probability distribution,
$\pi_{i}=\Big{(}\frac{n_{i}}{k_{i}},1-\frac{n_{i}}{k_{i}}\Big{)}$ and
theoretical probability distribution,
$p_{i}=\Big{(}1-\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{-1}\Big{.}$,
$\Big{.}\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{-1}\Big{)}$ for the
$ith$ group (for $i=1,\ldots,I$) can be obtained as

$\displaystyle D_{\beta}(p_{i},\pi_{i})=$
$\displaystyle\left[\frac{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}(\beta+1)}}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{\beta+1}}\right]-\frac{\beta+1}{\beta}\left[\left(\frac{n_{i}}{k_{i}}\right)\right.\left(\frac{(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}\beta}}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{\beta}}\right)$
$\displaystyle\left.+\left(1-\frac{n_{i}}{k_{i}}\right)\left(\frac{1}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{\beta}}\right)\right]$
$\displaystyle+\frac{1}{\beta}\left[\left(\frac{n_{i}}{k_{i}}\right)^{\beta+1}+\left(1-\frac{n_{i}}{k_{i}}\right)^{\beta+1}\right].$
(7)

Considering all I groups and imposing weights proportional to the group size,
the weighted DPD measure with the weights $w_{i}=\frac{k_{i}}{K}$ for
$i=1,\ldots,I$ where $K=k_{1}+k_{2}+...+k_{I}$ is given as:

$\displaystyle D^{w}_{\beta}(\bm{\theta})=$
$\displaystyle\sum_{i=1}^{I}\frac{k_{i}}{K}D_{\beta}(p_{i},\pi_{i})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{I}\frac{k_{i}}{K}\left\\{\frac{1}{1+\left(e^{\lambda(\tau_{i})}-1\right)^{\alpha_{i}}}\right\\}^{\beta}\left[\frac{\left\\{1+\left(e^{\lambda_{i}(\tau_{i})}-1\right)^{\alpha_{i}(\beta+1)}\right\\}}{\left\\{1+\left(e^{\lambda_{i}(\tau_{i})}-1\right)^{\alpha_{i}}\right\\}}\right.$
$\displaystyle\qquad\qquad\left.-\frac{\beta+1}{\beta}\left\\{\frac{n_{i}}{k_{i}}\left(e^{\lambda_{i}(\tau_{i})}-1\right)^{\alpha_{i}\beta}+\left(1-\frac{n_{i}}{k_{i}}\right)\right\\}\right]$
$\displaystyle\qquad\qquad+\frac{1}{\beta}\left\\{\left(\frac{n_{i}}{k_{i}}\right)^{\beta+1}+\left(1-\frac{n_{i}}{k_{i}}\right)^{\beta+1}\right\\}$
(8)

The relationship between the likelihood function and the weighted DPD measure
can be established as follows.

$\displaystyle\mathop{lim}_{\beta\to 0^{+}}D^{w}_{\beta}(\bm{\theta})=$
$\displaystyle C-\frac{1}{K}\;lnL(\bm{\theta})$

where $C$ being a constant which is $\bm{\theta}$ independent.
Weighted Minimum density power divergence estimators (WMDPDE) for
$\bm{\theta}$ can be defined as follows.

$\displaystyle\hat{\bm{\theta}}_{\beta}=$ $\displaystyle
arg\mathop{min}_{\bm{\theta}}D^{w}_{\beta}(\bm{\theta})$ (9)

In parameter estimation, the tuning parameter $\beta$ plays an important role.
It brings a balance between the efficiency and robustness of the WMDPD
estimator.

######  Theorem 2.

The set of estimating equations to minimize WMDPD measure is given as follows.

$\displaystyle\sum_{i=1}^{I}k_{i}$
$\displaystyle\left[\left\\{\left(1-(1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}})^{-1}\right)^{\beta-1}+(1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}})^{-(\beta-1)}\right\\}\right.$
$\displaystyle\left.\left\\{\left(1-(1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}})^{-1}\right)-\frac{n_{i}}{k_{i}}\right\\}\right]\left[\frac{\partial
F(\tau_{i};x_{ij},\bm{\theta})}{\partial\bm{\theta}}\right]=0_{2J}.$ (10)

with
$\bm{\theta}=(a_{1},..,a_{j},..,a_{J},b_{1},..,b_{j},..,b_{J})^{{}^{\prime}}$
and

$\displaystyle\frac{\partial F(\tau_{i};x_{ij},\bm{\theta})}{\partial a_{j}}$
$\displaystyle=\frac{\alpha_{i}x_{ij}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\,ln(e^{\lambda_{i}\tau_{i}}-1)}{\left\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\right\\}^{2}}$
$\displaystyle\frac{\partial F(\tau_{i};x_{ij},\bm{\theta})}{\partial b_{j}}$
$\displaystyle=\frac{\alpha_{i}\lambda_{i}\tau_{i}\,x_{ij}\,e^{\lambda_{i}\tau_{i}}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}-1}}{\left\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\right\\}^{2}}$

For $\beta$ tends to 0, WMDPDE converges to MLE.

###  3.1 Asymptotic Property

Here, we present the asymptotic distribution of the weighted minimum density
power divergence estimator under the logistic exponential distribution based
on the failure count data.

######  Theorem 3.

Suppose $\bm{\theta}_{0}$ is the true value of the parameter $\bm{\theta},$
when $k_{i}\rightarrow\infty$ and $\lim_{k_{i}\to\infty}\frac{k_{i}}{K}$
finite for all $i=1,\ldots,I$

$\displaystyle\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})\xrightarrow[K\to\infty]{\mathscr{L}}N\left(0_{2J},J_{\beta}^{-1}(\bm{\theta}_{0})K_{\beta}(\bm{\theta}_{0})J_{\beta}^{-1}(\bm{\theta}_{0})\right)$
(11)

######  Proof 1.

Given in Appendix B.

###  3.2 Testing of hypothesis based on weighted Density Power Divergence
(WDPD)

Testing of statistical hypothesis is of fundamental importance in inferential
analysis. Due to the lack of robustness of MLE, tests based on MLE may not be
able to perform well in the presence of outliers in the data set. Here we
present a Robust WDPD-based test statistic for testing a simple null
hypothesis based on the idea of Basu et al. [17]. The null and alternative
hypotheses are given as follows.

$H_{0}:\bm{\theta}=\bm{\theta}_{0}\qquad\text{against}\qquad
H_{1}:\bm{\theta}\neq\bm{\theta}_{0}$

where $\bm{\theta}=\\{a_{j},b_{j};j=1,2,\dots,J\\}$.
The robust WDPD-based test statistic is defined as,

$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})=2KD^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$
(12)

where,

$\displaystyle D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left[\\{F^{\beta+1}_{i}(\bm{\theta}_{0})+\bar{F}^{\beta+1}_{i}(\bm{\theta}_{0})\\}-\frac{\beta+1}{\beta}\left\\{F_{i}(\hat{\bm{\theta}}_{\beta})F^{\beta}_{i}(\bm{\theta}_{0})\right.\right.$
$\displaystyle\qquad\left.\left.+\bar{F}_{i}(\hat{\bm{\theta}}_{\beta})\bar{F}^{\beta}(\bm{\theta}_{0})\right\\}+\frac{1}{\beta}\left\\{F_{i}^{\beta+1}(\hat{\bm{\theta}}_{\beta})+\bar{F}_{i}^{\beta+1}(\hat{\bm{\theta}}_{\beta})\right\\}\right].$

Applying Taylor series expansion of second order around
$\bm{\theta}=\bm{\theta}_{0}$ at $\bm{\theta}=\hat{\bm{\theta}}_{\beta},$
$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ can be expressed
as,

$\displaystyle\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$
$\displaystyle=2K\left[D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})+(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})^{T}\nabla
D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})\right.$
$\displaystyle\qquad\qquad+\left.\frac{1}{2}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})^{T}\nabla^{2}D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})\right].$
When $H_{0}$ is true,
$\displaystyle\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$
$\displaystyle=\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})^{T}\nabla^{2}D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})$
(13)
$\displaystyle\qquad(\text{since}\;\;D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})=0\,,\;\nabla
D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})=0)$

For further development, the following result is presented.

######  Result 1.

Let $\bm{X}\sim N_{p}(0,\Sigma)$ and $\bm{A}$ be a p$\times$p real symmetric
matrix then $\bm{X^{{}^{\prime}}AX}$ can be expressed as

$\displaystyle\bm{X^{{}^{\prime}}AX}$
$\displaystyle=\bm{X^{{}^{\prime}}\Sigma^{-1/2}\Sigma^{1/2}A\Sigma^{1/2}\Sigma^{-1/2}X}$
$\displaystyle=\bm{X^{{}^{\prime}}\Sigma^{-1/2}P\Lambda
P^{{}^{\prime}}\Sigma^{-1/2}X}$ $\displaystyle=\bm{W^{{}^{\prime}}\Lambda W}.$

Here, $\bm{P}$ is an orthogonal matrix containing the eigen vectors of
$\bm{\Sigma^{1/2}A\Sigma^{1/2}}.$ $\bm{\Lambda}$ is a diagonal matrix with
diagonal elements $\lambda_{l}$ which are eigen values of
$\bm{\Sigma^{1/2}A\Sigma^{1/2}}$ and
$\bm{W}=\bm{P^{{}^{\prime}}\Sigma^{-1/2}X}$.

As $\bm{W\sim N(0,I)},$ therefore, the distribution of
$\bm{X^{{}^{\prime}}AX}$ is same as the distribution of
$\sum_{l=1}^{r}\lambda_{l}W^{2}_{l}$, where $W_{1},\ldots,W_{r}$ are the
independent standard normal variables. Here,
$r=rank\bm{(\Sigma^{1/2}A\Sigma^{1/2})}=rank\bm{(A\Sigma)}$ and $\lambda_{i}s$
are the non-zero eigenvalues of $\bm{A\Sigma}$.

Using the above result, the asymptotic distribution of test statistic
$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ can be described
by $\sum_{l=1}^{r}\lambda^{(\beta)}_{l}W^{2}_{l}$, where
$\lambda^{(\beta)}_{l}s$ are the non-zero eigen values of
$\bm{(\nabla^{2}D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})\Sigma_{\beta}(\bm{\theta}_{0})})$
and $\bm{\Sigma_{\beta}(\bm{\theta}_{0})}$=
$\bm{J^{-1}_{\beta}(\bm{\theta}_{0})K_{\beta}(\bm{\theta}_{0})}$
$\bm{J^{-1}_{\beta}(\bm{\theta}_{0})}$ with
$r=rank\bm{(\nabla^{2}D^{w}_{\beta}}$
$\bm{(\bm{\theta}_{0},\bm{\theta}_{0})\Sigma_{\beta}(\bm{\theta}_{0}))}$.
Further, let us define,

$\Lambda^{*}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})=\frac{\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})}{\lambda^{(\beta)}_{max}}\leq\sum_{l=1}^{r}W^{2}_{l}$

where $\lambda^{(\beta)}_{max}=max(\lambda^{(\beta)}_{l};l=1,2,\dots,r)$. As,
$\sum_{l=1}^{r}W^{2}_{l}\sim\chi^{2}_{(r)}$, $H_{0}$ is rejected when
$\Lambda^{*}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})\geq\chi^{2}_{(r,1-\alpha)}$
where $\chi^{2}_{(r,1-\alpha)}$ is the upper $(1-\alpha)$ quantile point of
$\chi^{2}_{(r)}$.

###  3.3 Power Function of WDPD Based Test

In this section, we study the power function of WDPD based test. The
alternative hypothesis is set as
$H_{1,K}:\bm{\theta}_{K}=\bm{\theta}_{0}+K^{-1/2}\mathbf{d}$, where
$\mathbf{\bm{\theta}}_{K}\in\mathbf{\bm{\theta}}\subset\bm{\mathbb{R}^{2J}}$
and $\mathbf{d}$ is a fixed vector in $\bm{\mathbb{R}^{2J}}$. Further
$\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}_{0})$ can be written as,

$\displaystyle\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}_{0})$
$\displaystyle=\sqrt{K}(\bm{\hat{\theta}_{\beta}}-\bm{\theta}_{K})+\bm{d}$
$\displaystyle\qquad\text{Under }H_{1,K},$
$\displaystyle\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}_{K})$
$\displaystyle\xrightarrow[K\to\infty]{\mathscr{L}}N\left(0_{\bm{2J}},\bm{\Sigma}_{\beta}(\bm{\theta}_{K})\right)\quad\text{and
therefore}$ $\displaystyle\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}_{0})$
$\displaystyle\xrightarrow[K\to\infty]{\mathscr{L}}N\left(\bm{d},\bm{\Sigma}_{\beta}(\bm{\theta}_{K})\right)$

The establishment of the asymptotic distribution of
$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ under $H_{1,K}$
is based on the Corollary 2.2 of Dik and De Gunst [21] which is given as
follows.

######  Result 2.

Let $\bm{X\sim N_{p}(\mu,\Sigma)}$, $\bm{A}$ be a real-symmetric non-negative
definite matrix of order p, $r=rank(\bm{\Sigma A\Sigma}),r\geq 1$ and
$\lambda_{1},\dots,\lambda_{r}$ be the positive eigenvalues of $\bm{A\Sigma}$.
Then the distribution of $\bm{X^{{}^{\prime}}AX}$ is equivalent to the
distribution of $\sum_{l=1}^{r}\lambda_{l}(W_{l}+v_{l})^{2}+\bm{\Psi}$, where
$W_{1},\dots,W_{r}$ are independent standard normal variables and
$\bm{v}=\bm{\Lambda^{-1}P^{{}^{\prime}}S^{{}^{\prime}}A\mu}$,
$\bm{\Psi=\mu^{{}^{\prime}}A\mu-v^{{}^{\prime}}\Lambda v}$. Here, $\bm{S}$ be
any square root of $\bm{\Sigma}$ and $\bm{S^{{}^{\prime}}AS=P\Lambda
P^{{}^{\prime}}}$ where $\bm{\Lambda}$ is a diagonal matrix with diagonal
elements $\lambda_{1},\dots,\lambda_{r}$ which are positive eigenvalues of
$\bm{S^{{}^{\prime}}AS}$ and $\bm{P}$ is an orthogonal matrix with column
vectors being the eigenvectors of the corresponding eigenvalues.

Using the above result, the asymptotic distribution of
$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ under $H_{1,K}$
is equivalent with the distribution of
$\sum_{l=1}^{r}\lambda^{(\beta)}_{l}(\bm{\theta}_{0})(W_{l}+v_{l})^{2}+\bm{\Psi}$
where,
$\lambda^{(\beta)}_{1}(\bm{\theta_{0}}),\dots,\lambda^{(\beta)}_{r}(\bm{\theta_{0}})$
are the positive eigenvalues of
$\bm{({\nabla^{2}D^{w}_{\beta}(\theta_{0},\theta_{0})\Sigma_{\beta}(\theta_{0})})}$,

$\displaystyle W_{l}\mathop{\sim}N(0,1)\ \text{independently for}\
l=1,\ldots,r,$
$\displaystyle\bm{v}=\bm{\Lambda^{-1}P^{{}^{\prime}}S^{{}^{\prime}}\nabla^{2}D^{w}_{\beta}(\theta_{0},\theta_{0})d}\
\text{where}\ \bm{v}=(v_{1},\dots,v_{r})^{{}^{\prime}},$
$\displaystyle\bm{\Psi}=\bm{d^{{}^{\prime}}\nabla^{2}D^{w}_{\beta}(\theta_{0},\theta_{0})d-v^{{}^{\prime}}\Lambda
v},$ $\displaystyle\bm{S}\;\text{is square root
of}\;\bm{\Sigma_{\beta}(\theta_{0})},$
$\displaystyle\bm{\Lambda}=diag(\lambda^{(\beta)}_{1}(\bm{\theta}_{0}),\dots,\lambda^{(\beta)}_{r}(\bm{\theta}_{0})),$
$\displaystyle\lambda^{(\beta)}_{l}(\bm{\theta}_{0})\ \text{for}\
l=1,\ldots,r,\ \text{are the positive eigenvalues of}\
\bm{(S^{{}^{\prime}}\nabla^{2}D^{w}_{\beta}(\theta_{0},\theta_{0})S)}$
$\displaystyle\text{and}\ \bm{P}\ \text{is an orthogonal matrix with column
vectors being the eigenvectors of}$ the corresponding eigenvalues.

Though the asymptotic distribution of
$\Lambda_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ under $H_{1,K}$
using the above method is very informative, yet it does not help determine the
power function due to its complex nature. The alternative approach for
determining the power function approximation is given as follows.
The first order Taylor expansion of
$D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})$ under
$\bm{\theta^{*}}$, $\bm{\theta^{*}}\neq\bm{\theta}_{0}$ is given as,

$D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})=D^{w}_{\beta}(\bm{\theta^{*}},\bm{\theta}_{0})+\bm{A^{{}^{\prime}}_{\beta}}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}^{*}),\,\text{where
}\,\bm{A_{\beta}}=\nabla D^{w}_{\beta}(\bm{\theta}^{*},\bm{\theta}_{0})$

As

$\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}^{*})\xrightarrow[K\to\infty]{\mathscr{L}}N\left(0_{\bm{2J}},\bm{\Sigma}_{\beta}(\bm{\theta}^{*})\right)$

Then,$\;\sqrt{K}(D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})-D^{w}_{\beta}(\bm{\theta}^{*},\bm{\theta}_{0}))$
and
$\bm{A^{{}^{\prime}}_{\beta}}\sqrt{K}(\bm{\hat{\theta}}_{\beta}-\bm{\theta}^{*})$
have the same asymptotic distribution and

$\sqrt{K}(D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})-D^{w}_{\beta}({\bm{\theta}}^{*},\bm{\theta}_{0}))\xrightarrow[K\to\infty]{\mathscr{L}}N\left(0_{\bm{2J}},\bm{\Sigma}^{*}_{\beta}(\bm{\theta}^{*})\right)$

where,
$\bm{\Sigma}^{*}_{\beta}(\bm{\theta}^{*})=\bm{A}^{{}^{\prime}}_{\beta}\bm{\Sigma}_{\beta}(\bm{\theta}^{*})\bm{A}_{\beta}$.

Therefore the power function can be obtained as

$\displaystyle\bm{\pi}^{(\beta)}_{K,\alpha}(\bm{\theta}^{*})$
$\displaystyle=Pr[2KD^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})>c^{(\beta)}_{\alpha}]$
$\displaystyle=Pr\left[\frac{\sqrt{K}(D^{w}_{\beta}(\hat{\bm{\theta}}_{\beta},\bm{\theta}_{0})-D^{w}_{\beta}(\bm{\theta}^{*},\bm{\theta}_{0}))}{\bm{\Sigma}^{*}_{\beta}(\bm{\theta}^{*})}\right.$
$\displaystyle\qquad\qquad\left.>\frac{\sqrt{K}}{\bm{\Sigma}^{*}_{\beta}(\bm{\theta}^{*})}\left(\frac{c^{(\beta)}_{\alpha}}{2K}-D^{w}_{\beta}(\bm{\theta}^{*},\bm{\theta}_{0})\right)\right]$
$\displaystyle\bm{\pi}^{(\beta)}_{K,\alpha}(\theta^{*})$
$\displaystyle=1-\Phi\left[\frac{\sqrt{K}}{\bm{\Sigma}^{*}_{\beta}(\bm{\theta}^{*})}\left(\frac{c^{(\beta)}_{\alpha}}{2K}-D^{w}_{\beta}(\bm{\theta}^{*},\bm{\theta}_{0})\right)\right]$
(14)

where, $\Phi(x)$ is the standard normal distribution function and
$c^{(\beta)}_{\alpha}$ is $(1-\alpha)$ percentile of the distribution of
$\Lambda_{\beta}(\hat{\bm{\theta}},\bm{\theta}_{0})$ under $H_{0}$.

##  4 Simulation Experiment

In this section, a simulated environment is created using the Markov Chain
Monte Carlo simulations based on 1000 generations and performances of MLE and
WMDPDEs are observed.

The lifetimes of the one-shot devices are considered to follow the Logistic-
Exponential distribution. An accelerated life test is conducted under three
testing groups, with three inspection times. A different number of devices are
put to the test in each group which are subjected to two stress factors. The
layout summaries are given in Table (2).

Table 2: Layout of one-shot device ALT design for simulation |  |  |  | Covariates
---|---|---|---|---
Groups | Inspection Times | Devices | Failures | Stress 1 | Stress 2
1 | $1.00$ | $\color[rgb]{0,0,0}{15}$ | $n_{1}$ | $0.2$ | $0.4$
2 | $1.00$ | $\color[rgb]{0,0,0}{20}$ | $n_{2}$ | $0.3$ | $0.6$
3 | $1.00$ | $\color[rgb]{0,0,0}{25}$ | $n_{3}$ | $0.4$ | $0.8$

The three sets of model parameters are taken so that the results can be
analysed with at least three different viewpoints and the performances of the
estimators can be assessed. The generated data has been contaminated in three
ways, and the results are analysed. Different sets of model parameters for
pure ALT data and contamination scheme are given in Table (3).

Table 3: Model parameters to generate data S.No. | $\;\mathbf{a_{1}}\;$ | $\;\mathbf{a_{2}}\;$ | $\;\mathbf{b_{1}}\;$ | $\;\mathbf{b_{2}}\;$
---|---|---|---|---
Pure Data |  |  |  |
$\bm{\theta}_{1}$ | $0.2$ | $-0.6$ | $-0.2$ | $0.4$
$\bm{\theta}_{2}$ | $0.4$ | $0.3$ | $-0.1$ | $-0.2$
$\bm{\theta}_{3}$ | $-0.06$ | $-0.06$ | $0.4$ | $-0.1$
Contamination |  |  |  |
$\bm{\theta}_{1}$ | $a_{1}+0.02$ | $a_{2}-0.09$ | $b_{1}-0.07$ | $b_{2}-0.01$
$\bm{\theta}_{2}$ | $a_{1}-0.07$ | $a_{2}+0.06$ | $b_{1}-0.05$ | $b_{2}+0.03$
$\bm{\theta}_{3}$ | $a_{1}+0.00$ | $a_{2}-0.07$ | $b_{1}+0.07$ | $b_{2}+0.05$

For the given data, the Coordinate-Descent method is used to obtain MLEs and
WMDPDEs. The steps of the Coordinate-Descent algorithm are given in the
Algorithm (1).

Algorithm 1 Coordinate-Descent Algorithm

  * •

Choose initial value of $\bm{\theta}=(a_{1},a_{2},b_{1},b_{2})$ say
$\bm{\theta}_{0}=(a^{(0)}_{1},a^{(0)}_{2},b^{(0)}_{1},b^{(0)}_{2})$.

  * •

Let at $t$th iteration for $t=0,1,2,\ldots,$ the estimate of $\bm{\theta}$ be
$\bm{\theta}^{t}=(a^{(t)}_{1},a^{(t)}_{2},b^{(t)}_{1},b^{(t)}_{2})$. Then at
$t+1$th iteration, the estimate can be derived as,
$a^{(t+1)}_{1}=a^{(t)}_{1}-h\frac{\partial
H(a^{(t)}_{1},a^{(t)}_{2},b^{(t)}_{1},b^{(t)}_{2})}{\partial a_{1}}\\\
a^{(t+1)}_{2}=a^{(t)}_{2}-h\frac{\partial
H(a^{(t+1)}_{1},a^{(t)}_{2},b^{(t)}_{1},b^{(t)}_{2})}{\partial a_{2}}\\\
b^{(t+1)}_{1}=b^{(t)}_{1}-h\frac{\partial
H(a^{(t+1)}_{1},a^{(t+1)}_{2},b^{(t)}_{1},b^{(t)}_{2})}{\partial b_{1}}\\\
b^{(t+1)}_{2}=b^{(t)}_{2}-h\frac{\partial
H(a^{(t+1)}_{1},a^{(t+1)}_{2},b^{(t+1)}_{1},b^{(t)}_{2})}{\partial b_{2}}\\\ $
where $H=-\ln L(\bm{\theta})$ in case of MLE and
$H=D^{w}_{\beta}(\bm{\theta})$ in case of WMDPDE where $h$ is the learning
rate. The value of the learning rate is chosen as $h=0.01.$

  * •

The process continues until
$max(|a^{(t+1)}_{1}-a^{(t)}_{1}|,|a^{(t+1)}_{2}-a^{(t)}_{2}|,|b^{(t+1)}_{1}-b^{(t)}_{1}|,|b^{(t+1)}_{2}-b^{(t)}_{2}|)$
is less than some pre-specified threshold value and if it is satisfied the
final estimate is obtained as
$\bm{\theta}^{t+1}=(a^{(t+1)}_{1},a^{(t+1)}_{2},b^{(t+1)}_{1},b^{(t+1)}_{2}).$

In Tables (4)-(7), bias and mean-square-error (MSE) of MLE and WMDPDE, as well
as estimated reliability at time point $t=2.5$ for pure data and contamination
settings, are provided. It can be seen that the bias and MSE of MLE are lower
than those of WMDPDE in the context of pure data. Thus, MLE beats WMDPDE in
the pure data scheme. But in case of contamination in data, bias and MSE of
WMDPDE are lower than those of MLE. The increase of bias of WMDPDE from pure
data to contaminated data is lesser than that of MLE. This indicates that
WMDPDE is a robust estimator without compromising efficiency much. The
graphical representation of bias of reliability estimates in pure data and
contaminated data settings at different time points is given in Figure 2. It
can be visually observed that in the pure data scheme, MLE performs better
than WMDPDE and in the contaminated data setting WMDPDE outperforms MLE in
terms of bias. Due to the superiority of WMDPDE numerically proven here, these
estimators are the preferable choice over MLEs when data is contaminated.

Table 4: Bias of Estimates of parameters in pure and contaminated data setting

| Pure Data | Contaminated Data
---|---|---
| $\mathbf{a_{1}}$ | $\mathbf{a_{2}}$ | $\mathbf{b_{1}}$ | $\mathbf{b_{2}}$ | $\mathbf{a_{1}}$ | $\mathbf{a_{2}}$ | $\mathbf{b_{1}}$ | $\mathbf{b_{2}}$
$\bm{\theta=\theta}_{1}$ |  |  |  |  |  |  |  |
MLE | -0.0073 | -0.0231 | -0.0245 | -0.0010 | -0.0390 | -0.0878 | -0.0720 | -0.0639
$\bm{\beta{=}0.2}$ | -0.0100 | -0.0303 | -0.0362 | -0.0020 | -0.0230 | -0.0561 | -0.0621 | -0.0441
$\bm{\beta{=}0.4}$ | -0.0137 | -0.0379 | -0.0402 | -0.0055 | -0.0249 | -0.0598 | -0.0646 | -0.0492
$\bm{\beta{=}0.6}$ | -0.0162 | -0.0425 | -0.0432 | -0.0069 | -0.0235 | -0.0570 | -0.0639 | -0.0479
$\bm{\beta{=}0.8}$ | -0.0139 | -0.0378 | -0.0414 | -0.0039 | -0.0211 | -0.0523 | -0.0631 | -0.0462
$\bm{\beta{=}1.0}$ | -0.0140 | -0.0381 | -0.0448 | -0.0099 | -0.0194 | -0.0489 | -0.0627 | -0.0454
$\bm{\theta=\theta}_{2}$ |  |  |  |  |  |  |  |
MLE | -0.0061 | 0.0029 | 0.0081 | -0.0419 | -0.0221 | 0.0398 | 0.0278 | -0.0582
$\bm{\beta{=}0.2}$ | -0.0147 | 0.0215 | 0.0152 | -0.0462 | -0.0104 | 0.0251 | 0.0159 | -0.0472
$\bm{\beta{=}0.4}$ | -0.0170 | 0.0258 | 0.0145 | -0.0495 | -0.0177 | 0.0245 | 0.0184 | -0.0407
$\bm{\beta{=}0.6}$ | -0.0209 | 0.0181 | 0.0213 | -0.0354 | -0.0266 | 0.0201 | 0.0235 | -0.0311
$\bm{\beta{=}0.8}$ | -0.0229 | 0.0140 | 0.0248 | -0.0416 | -0.0241 | 0.0196 | 0.0262 | -0.0262
$\bm{\beta{=}1.0}$ | -0.0221 | 0.0156 | 0.0227 | -0.0242 | -0.0229 | 0.0141 | 0.0279 | -0.0231
$\bm{\theta=\theta}_{3}$ |  |  |  |  |  |  |  |
MLE | -0.0010 | -0.0140 | 0.0301 | -0.0130 | -0.0171 | -0.0634 | 0.0717 | 0.0622
$\bm{\beta{=}0.2}$ | 0.0047 | -0.0204 | 0.0367 | -0.0060 | 0.0047 | -0.0203 | 0.0490 | 0.0179
$\bm{\beta{=}0.4}$ | 0.0019 | -0.0260 | 0.0307 | -0.0177 | 0.0023 | -0.0124 | 0.0499 | 0.0199
$\bm{\beta{=}0.6}$ | 0.0031 | -0.0236 | 0.0330 | -0.0133 | 0.0077 | -0.0145 | 0.0482 | 0.0163
$\bm{\beta{=}0.8}$ | 0.0027 | -0.0184 | 0.0371 | -0.0110 | 0.0034 | -0.0101 | 0.0452 | 0.0103
$\bm{\beta{=}1.0}$ | 0.0066 | -0.0167 | 0.0327 | -0.0140 | 0.0073 | -0.0113 | 0.0447 | 0.0113

Table 5: MSE of Estimates of parameters in pure and contaminated data setting

| Pure Data | Contaminated Data
---|---|---
| $\mathbf{a_{1}}$ | $\mathbf{a_{2}}$ | $\mathbf{b_{1}}$ | $\mathbf{b_{2}}$ | $\mathbf{a_{1}}$ | $\mathbf{a_{2}}$ | $\mathbf{b_{1}}$ | $\mathbf{b_{2}}$
$\bm{\theta=\theta}_{1}$ |  |  |  |  |  |  |  |
MLE | 0.0080 | 0.0345 | 0.0164 | 0.0721 | 0.0345 | 0.1383 | 0.0456 | 0.1626
$\bm{\beta{=}0.2}$ | 0.0087 | 0.0351 | 0.0249 | 0.0881 | 0.0122 | 0.0498 | 0.0222 | 0.0734
$\bm{\beta{=}0.4}$ | 0.0165 | 0.0714 | 0.0174 | 0.0930 | 0.0110 | 0.0451 | 0.0296 | 0.1024
$\bm{\beta{=}0.6}$ | 0.0174 | 0.0704 | 0.0265 | 0.0972 | 0.0106 | 0.0433 | 0.0259 | 0.0882
$\bm{\beta{=}0.8}$ | 0.0146 | 0.0592 | 0.0210 | 0.0639 | 0.0094 | 0.0385 | 0.0204 | 0.0667
$\bm{\beta{=}1.0}$ | 0.0094 | 0.0382 | 0.0181 | 0.0635 | 0.0033 | 0.0141 | 0.0221 | 0.0738
$\bm{\theta=\theta}_{2}$ |  |  |  |  |  |  |  |
MLE | 0.0600 | 0.3946 | 0.0287 | 0.1177 | 0.0717 | 0.3467 | 0.0229 | 0.1815
$\bm{\beta{=}0.2}$ | 0.0476 | 0.1829 | 0.0251 | 0.0975 | 0.0525 | 0.1919 | 0.0321 | 0.1236
$\bm{\beta{=}0.4}$ | 0.0372 | 0.1482 | 0.0306 | 0.1097 | 0.0212 | 0.0842 | 0.0217 | 0.0841
$\bm{\beta{=}0.6}$ | 0.0231 | 0.0908 | 0.0213 | 0.0818 | 0.0406 | 0.1672 | 0.0170 | 0.0641
$\bm{\beta{=}0.8}$ | 0.0321 | 0.1264 | 0.0167 | 0.0693 | 0.0211 | 0.0798 | 0.0187 | 0.0708
$\bm{\beta{=}1.0}$ | 0.0136 | 0.0529 | 0.0152 | 0.0523 | 0.0121 | 0.0467 | 0.0134 | 0.0494
$\bm{\theta=\theta}_{3}$ |  |  |  |  |  |  |  |
MLE | 0.0170 | 0.0472 | 0.0291 | 0.0828 | 0.0392 | 0.1560 | 0.0496 | 0.1742
$\bm{\beta{=}0.2}$ | 0.0050 | 0.0201 | 0.0140 | 0.0790 | 0.0106 | 0.0423 | 0.0193 | 0.0649
$\bm{\beta{=}0.4}$ | 0.0088 | 0.0360 | 0.0173 | 0.0637 | 0.0175 | 0.0498 | 0.0199 | 0.0681
$\bm{\beta{=}0.6}$ | 0.0099 | 0.0400 | 0.0177 | 0.0650 | 0.0084 | 0.0334 | 0.0167 | 0.0564
$\bm{\beta{=}0.8}$ | 0.0079 | 0.0208 | 0.0121 | 0.0485 | 0.0146 | 0.0321 | 0.0115 | 0.0368
$\bm{\beta{=}1.0}$ | 0.0034 | 0.0138 | 0.0098 | 0.0342 | 0.0105 | 0.0269 | 0.0136 | 0.0493

Table 6: Bias of Reliability estimates of parameters in pure and contaminated
data setting

| Pure Data | Contaminated Data
---|---|---
| Group 1 | Group 2 | Group 3 | Group 1 | Group 2 | Group 3
$\bm{\theta=\theta}_{1}$ |  |  |  |  |  |
MLE | 0.0057 | 0.0099 | 0.0107 | 0.0219 | 0.0338 | 0.0423
$\bm{\beta{=}0.2}$ | 0.0087 | 0.0141 | 0.0181 | 0.0143 | 0.0221 | 0.0296
$\bm{\beta{=}0.4}$ | 0.0079 | 0.0125 | 0.0165 | 0.0135 | 0.0215 | 0.0287
$\bm{\beta{=}0.6}$ | 0.0089 | 0.0139 | 0.0183 | 0.0126 | 0.0203 | 0.0278
$\bm{\beta{=}0.8}$ | 0.0079 | 0.0133 | 0.0177 | 0.0128 | 0.0199 | 0.0264
$\bm{\beta{=}1.0}$ | 0.0086 | 0.0131 | 0.0174 | 0.0138 | 0.0221 | 0.0303
$\bm{\theta=\theta}_{2}$ |  |  |  |  |  |
MLE | 0.0017 | 0.0038 | 0.0057 | 0.0067 | 0.0095 | 0.0114
$\bm{\beta{=}0.2}$ | 0.0054 | 0.0075 | 0.0092 | 0.0047 | 0.0061 | 0.0075
$\bm{\beta{=}0.4}$ | 0.0049 | 0.0062 | 0.0075 | 0.0045 | 0.0060 | 0.0075
$\bm{\beta{=}0.6}$ | 0.0059 | 0.0076 | 0.0093 | 0.0050 | 0.0062 | 0.0076
$\bm{\beta{=}0.8}$ | 0.0054 | 0.0069 | 0.0083 | 0.0058 | 0.0071 | 0.0084
$\bm{\beta{=}1.0}$ | 0.0045 | 0.0059 | 0.0073 | 0.0040 | 0.0049 | 0.0060
$\bm{\theta=\theta}_{3}$ |  |  |  |  |  |
MLE | 0.0049 | 0.0092 | 0.0135 | 0.0076 | 0.0188 | 0.0217
$\bm{\beta{=}0.2}$ | 0.0027 | 0.0060 | 0.0093 | 0.0020 | 0.0063 | 0.0087
$\bm{\beta{=}0.4}$ | 0.0059 | 0.0102 | 0.0133 | 0.0021 | 0.0048 | 0.0064
$\bm{\beta{=}0.6}$ | 0.0072 | 0.0136 | 0.0174 | 0.0018 | 0.0039 | 0.0051
$\bm{\beta{=}0.8}$ | 0.0041 | 0.0081 | 0.0113 | -0.0010 | 0.0057 | 0.0054
$\bm{\beta{=}1.0}$ | 0.0025 | 0.0049 | 0.0075 | -0.0022 | -0.0020 | -0.0014

Table 7: MSE of Reliability estimates of parameters in pure and contaminated
data setting

| Pure Data | Contaminated Data
---|---|---
| Group 1 | Group 2 | Group 3 | Group 1 | Group 2 | Group 3
$\bm{\theta=\theta}_{1}$ |  |  |  |  |  |
MLE | 0.0015 | 0.0039 | 0.0064 | 0.0053 | 0.0114 | 0.0147
$\bm{\beta{=}0.2}$ | 0.0024 | 0.0057 | 0.0082 | 0.0023 | 0.0044 | 0.0074
$\bm{\beta{=}0.4}$ | 0.0021 | 0.0045 | 0.0071 | 0.0023 | 0.0053 | 0.0083
$\bm{\beta{=}0.6}$ | 0.0022 | 0.0045 | 0.0070 | 0.0017 | 0.0041 | 0.0069
$\bm{\beta{=}0.8}$ | 0.0018 | 0.0046 | 0.0072 | 0.0023 | 0.0048 | 0.0074
$\bm{\beta{=}1.0}$ | 0.0022 | 0.0040 | 0.0066 | 0.0015 | 0.0036 | 0.0064
$\bm{\theta=\theta}_{2}$ |  |  |  |  |  |
MLE | 0.0028 | 0.0043 | 0.0048 | 0.0048 | 0.0070 | 0.0078
$\bm{\beta{=}0.2}$ | 0.0033 | 0.0048 | 0.0053 | 0.0035 | 0.0043 | 0.0047
$\bm{\beta{=}0.4}$ | 0.0033 | 0.0041 | 0.0045 | 0.0032 | 0.0043 | 0.0049
$\bm{\beta{=}0.6}$ | 0.0033 | 0.0044 | 0.0050 | 0.0028 | 0.0035 | 0.0039
$\bm{\beta{=}0.8}$ | 0.0026 | 0.0033 | 0.0036 | 0.0031 | 0.0038 | 0.0041
$\bm{\beta{=}1.0}$ | 0.0023 | 0.0031 | 0.0036 | 0.0024 | 0.0032 | 0.0036
$\bm{\theta=\theta}_{3}$ |  |  |  |  |  |
MLE | 0.0012 | 0.0032 | 0.0052 | 0.0042 | 0.0171 | 0.0205
$\bm{\beta{=}0.2}$ | 0.0009 | 0.0028 | 0.0048 | 0.0018 | 0.0064 | 0.0087
$\bm{\beta{=}0.4}$ | 0.0021 | 0.0053 | 0.0067 | 0.0019 | 0.0051 | 0.0062
$\bm{\beta{=}0.6}$ | 0.0022 | 0.0076 | 0.0092 | 0.0019 | 0.0046 | 0.0056
$\bm{\beta{=}0.8}$ | 0.0013 | 0.0039 | 0.0054 | 0.0010 | 0.0034 | 0.0058
$\bm{\beta{=}1.0}$ | 0.0010 | 0.0023 | 0.0035 | 0.0007 | 0.0019 | 0.0029

(a) Group 1

(b) Group 1

(c) Group 2

(d) Group 2

(e) Group 3

(f) Group 3

Figure 2: Bias of reliability estimates in pure (left) and contaminated
(right) data set.

##  5 Optimum Inspection Time

One-shot devices stay in a torpid state until used. Hence the behaviour of
such devices is revealed only when these are put to the test. The aim is to
find a set of inspection times which will increase the precision of the
estimator and reduce the experimental cost due to the destruction of
experimental units. Therefore, the cost function can be defined as follows.

$Cost=C_{1}\times|V|+C_{2}\times E(\sum_{i=1}^{I}n_{i})$ (15)

Here $|V|$ is the determinant value of the covariance matrix based on the
asymptotic distribution of WMDPDE, $E(\sum_{i=1}^{I}n_{i})$ is the expected
total number of failures from all I different groups. Cost $C_{1}$ is imposed
on $|V|$, which reflects the precision cost of the estimator and $C_{2}$ is
the cost per failure of experimental units. Given the values of group sizes,
the optimal inspection times will minimize the cost function defined in (15).
This cost function brings a trade-off between the precision of the estimation
and the experimental cost.

###  5.1 Search Algorithm

Here, the objective is to find the set of inspection times which will minimize
the above-defined cost. We apply population-based heuristic Natural
optimization method, namely Genetic Algorithm (GA), which was developed by
Holland [22] and popularized by Goldberg [23]. GA is an optimization method
that produces generations of the population through selection, crossover and
mutation. It applies to discrete as well as continuous variables and can deal
with a large number of variables. It can perform optimization with extremely
complex cost surfaces and is capable of discovering global optimums avoiding
the trap of local optima. GA is inherently parallel, modifiable, easily
distributed and adaptable to different problems [24]. GA works well with
numerically generated data, experimental data, or analytical functions [25].
Some of the applications of GA can be seen in Babaei and Sheidaii [26], Liou
et al. [27], Soolaki et al. [28], Lopez et al. [29] and references therein.

To search the optimal inspection times through GA, the strategy of keeping the
best or the elite observations of the current population in the new population
is followed; readers may see Chakraborty and Chaudhuri [30] for detailed
discussion. Thus traits of the observations with the best cost, i.e. minimum
cost, are retained from generation to generation. The steps of GA using
elitism is as follows.

  1. i)

Set the initial population size $N_{pop}$ where each element in the population
is a set of inspection times with dimension $N_{var}=I$.

  2. ii)

Obtain the cost for each set of inspection times and select the parents for
pairing.

  3. iii)

The parents are selected using a specified procedure based on cost. Among the
various selection methods, the random rank selection procedure and tournament
selection procedure are used, which are explained below.

  4. iv)

Each set of parents would produce two offspring. Odd numbered element is
chosen as the mother, and an even-numbered element is chosen as the father.
The offsprings are produced as:

$\text{Offspring 1}=b\times\text{mother}+(1-b)\times\text{father}$
$\text{Offspring 2}=(1-b)\times\text{mother}+b\times\text{father}$

where $b$ is a random uniform number. Thus a new population is formed
containing the parents and offspring.

  5. v)

Now, the population is mutated with the mutation rate $mr$. So the number of
mutations is $mn=mr\times N_{pop}\times N_{var}$. The values in randomly
selected $mn$ positions out of $N_{pop}\times N_{var}$ positions are replaced
with the randomly generated new values. Hence a mutated population is
obtained, which would function as the initial population for the next
generation.

  6. vi)

The process is repeated up to, say, B generations. For each generation, a
local minimum is obtained, i.e. the set of inspection times with the least
cost. Among these B local minima, the global minimum is obtained: the set of
inspection times with minimum cost based on all the generated populations.
This global minimum provides the set of optimum inspection times.

####  5.1.1 Random Rank Selection Procedure

     * •

Rank the observations concerning the cost in increasing order.

     * •

Define $N_{keep}=0.5\times N_{pop}$.

     * •

Keep first $N_{keep}$ observations and discard the rest.

     * •

Use a permutation of ranks and rearrange the $N_{keep}$ size data concerning
this permutation.

     * •

The resultant observations are chosen as the parents who are to be paired.

####  5.1.2 Tournament Selection Procedure

     * •

Randomly pick a pair of observations such that this pair is not selected
again.

     * •

The observation having lesser cost in the pair is kept as a parent, and the
other is discarded. Though the same pair can not be chosen again, the
discarded observation or the chosen observation in the pair can be chosen
again for pairing.

     * •

Repeat this process up to $N_{keep}$ times.

     * •

The resultant observations are chosen as the parents who are to be paired.

Table 8: Model Parameters S.No. | $\;\mathbf{a_{1}}\;$ | $\;\mathbf{a_{2}}\;$ | $\;\mathbf{b_{1}}\;$ | $\;\mathbf{b_{2}}\;$
---|---|---|---|---
$\bm{\theta}_{1}$ | $0.2$ | $-0.6$ | $-0.2$ | $0.4$
$\bm{\theta}_{2}$ | $-0.3$ | $-0.1$ | $-0.2$ | $0.1$

Here, a numerical experiment is conducted to derive optimal inspection times
implementing GA under the following setups. The chosen model parameters are
given in Table (8) with tuning parameter $\beta=0.1$, and the layout is given
in Table (9). Costs are set as $C_{1}=C_{2}=0.5$ where $N_{pop}=12,$
$N_{var}=5.$ GA continues till $B=250$ generations. In the crossover, the
random number $b$ follows $U(0,1),$ where in the mutation stage, the mutation
rate is set as $mr=0.20,$ and the value in the selected position is replaced
by a new random number from $U(0,10)$.

Table 9: Layout for the implementation of GA |  |  | Covariates
---|---|---|---
Groups | Devices | Failures | Stress 1 | Stress 2
1 | $8$ | $n_{1}$ | $20$ | $40$
2 | $12$ | $n_{2}$ | $30$ | $60$
3 | $16$ | $n_{3}$ | $40$ | $80$
4 | $20$ | $n_{4}$ | $30$ | $40$
5 | $24$ | $n_{5}$ | $20$ | $50$

The Optimum inspection times $\tau_{1},\tau_{2},\tau_{3},\tau_{4},\tau_{5}$ by
the random rank selection and tournament selection procedure are given in the
Table (10). Though the two selection procedures yield very different sets of
inspection times yet the resultant costs for both of them are very similar.
Thus these two selection procedures can be used alternatively for Genetic
Algorithm.

Table 10: Optimum Inspection Times

Selection Method | $\mathbf{\tau_{1}}$ | $\mathbf{\tau_{2}}$ | $\mathbf{\tau_{3}}$ | $\mathbf{\tau_{4}}$ | $\mathbf{\tau_{5}}$ | Cost
---|---|---|---|---|---|---
$\bm{\theta}=\bm{\theta}_{1}$ |  |  |  |  |  |
Random rank selection | 6.13382 | 0.98725 | 2.69014 | 8.90894 | 4.05132 | 2.46078
Tournament selection | 8.34024 | 1.56001 | 7.02301 | 0.98418 | 0.21274 | 2.32247
$\bm{\theta}=\bm{\theta}_{2}$ |  |  |  |  |  |
Random rank selection | 9.18581 | 2.63731 | 9.62269 | 0.05032 | 4.53493 | 1.80780
Tournament selection | 4.29600 | 9.97924 | 0.67811 | 8.40763 | 0.02557 | 1.95393

##  6 Application to Real Dataset

Actual data is adopted from the database named “Incidence$\cdot$SEER Research
Data, 8 Registries, Nov 2021 Sub (1975-2019)”[20] which is recorded by
National Cancer Institute (Surveillance, Epidemiology, and End Results Program
(SEER)). We have extracted data from patients between the periods 2016-2019
who are diagnosed with gallbladder cancer. The patients of three age groups
(40-49), (50-59), and (60-69) are considered for the study. Their median age
at diagnosis ($x_{i1}$) and the size of the tumour ($x_{i2}$) are taken as two
stress factors/covariates. The median age is divided by 10 for the sake of
computational ease. The size of the tumour is indicated below,

$x_{i2}=\begin{cases}1,&\text{if tumor size}\leq 50\,mm\\\
2,&\text{otherwise}\end{cases}$

Those patients are observed over a number of months, and the number of
failures was reported within follow-up times where death due to gallbladder
cancer indicates failure. It is also observed from the data that, there are
some patients who died within a month of diagnosis of disease and could not
contribute to the study. Those patients can be considered outliers. The
description of data is given in Table (11).

Table 11: Layout of Dataset extracted from [20] |  |  |  | Covariates
---|---|---|---|---
Groups | Observation | Patients | Death | Median | Tumor
| months | count | count | age | size
1 | 12 | 17 | 12 | 46.40 | 1
2 | 12 | 12 | 10 | 46.40 | 2
3 | 12 | 41 | 32 | 56.02 | 1
4 | 12 | 48 | 42 | 56.02 | 2
5 | 12 | 95 | 61 | 64.96 | 1
6 | 12 | 74 | 64 | 64.96 | 2

To check if a lifetime distribution can be fitted to the data, bootstrap
testing is conducted based on the distance-based test statistic proposed by
Balakrishnan and Ling [31]. The test statistic is given as,

$D=max_{i}|n_{i}-e_{i}|$

where $n_{i}$ is the observed failure, and $e_{i}$ is the expected failure for
$i=1,\ldots,6$. The MLE is used to estimate the expected failures,
$e_{i}=k_{i}\times F(\tau_{i};x_{ij},\bm{\theta})$. The bootstrap (BT) method
is utilized to approximate the p-value of the test based on 1000 generated
samples. The test statistic and corresponding p-value for Weibull, Gamma and
Logistic-Exponential distribution are given in Table (12). The approximate
p-value from this table suggests the suitability of the Logistic-Exponential
model to the given data. Therefore, we assume that the lifetime of gallbladder
cancer patients follows the LE distribution.

MLE and WMDPDE of the model parameters are obtained, and bootstrap estimates
of bias and root-mean-square-error generated from 1000 bootstrap samples are
computed. The coordinate descent algorithm’s initial parameter value is set as
(-4.46, 0.08, -0.21, 0.34), which is chosen through an extensive grid search
process. The outcomes of the MLE and WMDPDE are reported in Tables (13) and
(14). The graphical representation of the bias of reliability estimates at
different time points is given in Figure 3. It is observed that WMDPDE have
lower bias and RMSE than those of MLE. Since the data is contaminated, the
superiority of WMDPDE in this situation validates the objective of our study.

Table 12: Test statistic and corresponding p-value calculation Lifetime Distribution | Test Statistic | p-value
---|---|---
Weibull | 12.6873 | 0.004
Gamma | 33.5532 | 0.005
Logistic Exponential | 06.0709 | 0.192
Table 13: Estimates of parameters for real dataset

| MLE | $\bm{\beta=0.2}$ | $\bm{\beta=0.4}$ | $\bm{\beta=0.6}$ | $\bm{\beta=0.8}$ | $\bm{\beta=1.0}$
---|---|---|---|---|---|---
$\hat{\bm{a}}_{1}$ | -4.463718 | -4.462089 | -4.461218 | -4.460754 | -4.460504 | -4.460369
Bias | 0.002955 | 0.002513 | -0.000923 | -0.000666 | 0.000471 | -0.000070
RMSE | 0.059291 | 0.053992 | 0.044467 | 0.029429 | 0.022231 | 0.012952
$\hat{\bm{a}}_{2}$ | 0.080093 | 0.080022 | 0.080026 | 0.080026 | 0.080023 | 0.080020
Bias | 0.000089 | 0.000058 | -0.000023 | -0.000014 | 0.000007 | -0.000001
RMSE | 0.001739 | 0.001613 | 0.001240 | 0.000816 | 0.000608 | 0.000347
$\hat{\bm{b}}_{1}$ | -0.210451 | -0.210285 | -0.210189 | -0.210132 | -0.210098 | -0.210077
Bias | 0.000333 | 0.000176 | -0.000014 | -0.000048 | 0.000032 | -0.000005
RMSE | 0.004724 | 0.004219 | 0.003337 | 0.002143 | 0.001630 | 0.000902
$\hat{\bm{b}}_{2}$ | 0.339126 | 0.339590 | 0.339825 | 0.339939 | 0.339991 | 0.340012
Bias | 0.000906 | 0.000457 | -0.000069 | -0.000102 | 0.000063 | -0.000015
RMSE | 0.014378 | 0.012907 | 0.009650 | 0.006169 | 0.004589 | 0.002562

Table 14: Reliability estimates of parameters for real dataset

| MLE | $\bm{\beta=0.2}$ | $\bm{\beta=0.4}$ | $\bm{\beta=0.6}$ | $\bm{\beta=0.8}$ | $\bm{\beta=1.0}$
---|---|---|---|---|---|---
Group 1 |  |  |  |  |  |
Bias | 0.000818 | 0.000049 | 0.000341 | 0.000212 | 0.000008 | 0.000051
RMSE | 0.009974 | 0.008298 | 0.006667 | 0.004400 | 0.003277 | 0.001902
Group 2 |  |  |  |  |  |
Bias | 0.000991 | 0.000299 | 0.000410 | 0.000221 | 0.000037 | 0.000042
RMSE | 0.008364 | 0.006788 | 0.005318 | 0.003529 | 0.002575 | 0.001501
Group 3 |  |  |  |  |  |
Bias | 0.000763 | -0.000118 | 0.000309 | 0.000216 | -0.000014 | 0.000058
RMSE | 0.012092 | 0.010133 | 0.008180 | 0.005371 | 0.004052 | 0.002333
Group 4 |  |  |  |  |  |
Bias | 0.001228 | 0.0001184 | 0.000508 | 0.000305 | 0.000008 | 0.000064
RMSE | 0.015134 | 0.012468 | 0.009902 | 0.006535 | 0.004852 | 0.002809
Group 5 |  |  |  |  |  |
Bias | 0.000569 | -0.000213 | 0.000222 | 0.000174 | -0.000028 | 0.000050
RMSE | 0.011115 | 0.009351 | 0.007570 | 0.004949 | 0.003769 | 0.002156
Group 6 |  |  |  |  |  |
Bias | 0.000940 | -0.000202 | 0.000379 | 0.000267 | -0.000038 | 0.000062
RMSE | 0.017221 | 0.014293 | 0.011427 | 0.007481 | 0.005637 | 0.003233

(a) Group 1

(b) Group 2

(c) Group 3

(d) Group 4

(e) Group 5

(f) Group 6

Figure 3: Bias of reliability estimates for the real dataset.

For robust WDPD-based testing, the null and alternative hypothesis is taken as
follows.

$H_{0}:\bm{\theta}=\bm{\theta}_{0}\qquad\text{against}\qquad
H_{1}:\bm{\theta}\neq\bm{\theta}_{0}$

where, $\bm{\theta}_{0}=(-4.46,0.08,-0.21,0.34)^{{}^{\prime}}$. The values of
WDPD-based test statistic are given in Table (15). The
$rank(\nabla^{2}D^{w}_{\beta}(\bm{\theta}_{0},\bm{\theta}_{0})\Sigma_{\beta}(\bm{\theta}_{0}))=4$
and $\chi^{2}_{(4,0.95)}=0.711$. From the values of the test statistic in
Table (15), it is seen that in all cases $\chi^{2}_{(4,0.95)}$ value is
greater than the value of the observed test statistic. Hence the null
hypothesis fails to be rejected.

Table 15: WDPD-based test statistic Tuning parameter | Test Statistic
---|---
$\beta=0.2$ | 0.000998
$\beta=0.4$ | 0.000887
$\beta=0.6$ | 0.000791
$\beta=0.8$ | 0.000712
$\beta=1.0$ | 0.000648

##  7 Discussion and Conclusion

This work proposed a robust method of estimation, the robust weighted minimum
density power divergence estimators (WMDPDE), for point estimation of a
Logistic-Exponential lifetime distribution based on one-shot device failure
data. The derived estimators have been shown numerically to be more robust
than the conventional Maximum Likelihood Estimation (MLE) when dealing with
contaminated data. The asymptotic properties of the WMDPDEs and the density
power divergence-based test statistic were provided. Furthermore, a search for
optimal inspection times was conducted using a Genetic Algorithm based on a
cost function that balances the precision of the estimates and experimental
cost. Finally, the proposed method was applied to a real data set where the
Logistic-Exponential distribution was found to be a good fit.

In future we can incorporate several inspection times per group where failures
are subject to multiple competing risks. To conduct the experiment in a
limited time and budget, we may apply progressive censoring in life-testing
experiments of one-shot devices. Works in those directions are in the
pipeline, and we are optimistic about reporting those findings soon.

## Acknowledgements

The authors would like to thank two unknown reviewers and also the associate
editor for many constructive suggestions which have helped to improve the
paper significantly.

## Declarations of 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.

## Data availability

Surveillance, Epidemiology, and End Results (SEER) Program
(www.seer.cancer.gov) SEER*Stat Database: Incidence - SEER Research Data, 8
Registries, Nov 2021 Sub (1975-2019) - Linked To County Attributes - Time
Dependent (1990-2019) Income/Rurality, 1969-2020 Counties, National Cancer
Institute, DCCPS, Surveillance Research Program, released April 2022, based on
the November 2021 submission.

##  Appendix A Elements of Hessian Matrix for MLE

$\displaystyle\frac{\partial^{2}\big{(}ln\,L(\bm{\theta})\big{)}}{\partial
a_{j_{1}}\partial a_{j_{2}}}$
$\displaystyle=\sum_{i=1}^{I}\frac{\alpha_{i}x_{ij_{1}}x_{ij_{2}}\,ln\,(e^{\lambda_{i}\tau_{i}}-1)}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{2}}\Big{[}k_{i}\alpha_{i}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\,\,ln\,(e^{\lambda_{i}\tau_{i}}-1)\Big{.}$
$\displaystyle\qquad\quad\Big{.}+\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}\\{n_{i}+(n_{i}-k_{i})(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}\Big{]}$
(16)
$\displaystyle\frac{\partial^{2}\big{(}ln\,L(\bm{\theta})\big{)}}{\partial
b_{j_{1}}\partial b_{j_{2}}}$
$\displaystyle=\sum_{i=1}^{I}\frac{\alpha_{i}\lambda_{i}\tau_{i}x_{ij_{1}}x_{ij_{2}}e^{\lambda_{i}\tau_{i}}(e^{\lambda_{i}\tau_{i}}-1)^{-2}}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{2}}\Big{[}k_{i}\alpha_{i}\lambda_{i}\tau_{i}e^{\lambda_{i}\tau_{i}}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha}\Big{.}$
$\displaystyle\quad\qquad+\\{e^{\lambda_{i}\tau_{i}}-\tau_{i}\lambda_{i}-1\\}\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}$
$\displaystyle\quad\qquad\qquad\Big{.}\\{n_{i}+(n_{i}-k_{i})(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}\Big{]}$
(17)
$\displaystyle\frac{\partial^{2}\big{(}ln\,L(\bm{\theta})\big{)}}{\partial
a_{j_{1}}\partial b_{j_{2}}}$
$\displaystyle=\sum_{i=1}^{I}\frac{\alpha_{i}\lambda_{i}\tau_{i}x_{ij_{1}}x_{ij_{2}}e^{\lambda_{i}\tau_{i}}(e^{\lambda_{i}\tau_{i}}-1)^{-1}}{\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}^{2}}\Big{[}k_{i}\alpha_{i}(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\Big{.}$
$\displaystyle\quad\qquad
ln\,(e^{\lambda_{i}\tau_{i}}-1)+\\{1+(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}$
$\displaystyle\quad\qquad\qquad\Big{.}\\{n_{i}+(n_{i}-k_{i})(e^{\lambda_{i}\tau_{i}}-1)^{\alpha_{i}}\\}\Big{]}$
(18)

where, $j_{1},j_{2}=1,\dots,J$.

##  Appendix B Proof of Theorem (3).

Based on [32], the proof of the theorem carries on as follows.

WDPD measure keeping only the terms involving parameters can be written as,

$\displaystyle D_{K}$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left[P^{\beta+1}_{i1}+\bar{P}^{\beta+1}_{i1}-\frac{\beta+1}{\beta}\left\\{\frac{n_{i}}{k_{i}}P^{\beta}_{i1}+\left(\frac{k_{i}-n_{i}}{k_{i}}\right)\bar{P}^{\beta}_{i1}\right\\}\right]$
$\displaystyle\text{where,}\;P_{i1}$
$\displaystyle=F(\tau_{i};x_{ij},\bm{\theta})\;,\;\bar{P}_{i1}=1-P_{i1}=\bar{F}(\tau_{i};x_{ij},\bm{\theta})$

Let us define, $X_{ui}\sim Bin(1,P_{i1})$, then,
$\sum_{ui=1}^{K_{i}}X_{ui}=n_{i}.$ Therefore,

$\displaystyle D_{K}$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left[\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}(P^{\beta+1}_{i1}+\bar{P}^{\beta+1}_{i1})\right.$
$\displaystyle\qquad\qquad\qquad\left.-\frac{\beta+1}{\beta}\frac{1}{k_{i}}\left\\{\sum_{u_{i}=1}^{k_{i}}X_{u_{i}}P^{\beta}_{i1}+\sum_{u_{i}=1}^{k_{i}}(1-X_{u_{i}})\bar{P}^{\beta}_{i1}\right\\}\right]$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left[\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}\left\\{(P^{\beta+1}_{i1}+\bar{P}^{\beta+1}_{i1})\right.\right.$
$\displaystyle\qquad\qquad\qquad\left.\left.-\frac{\beta+1}{\beta}\left(X_{u_{i}}P^{\beta}_{i1}+(1-X_{u_{i}})\bar{P}^{\beta}_{i1}\right)\right\\}\right]$
$\displaystyle\implies D_{K}$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left[\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}H_{i\beta}(\bm{\theta})\right]$

where,

$\displaystyle H_{i\beta}(\bm{\theta})$
$\displaystyle=\left\\{P^{\beta+1}_{i1}+\bar{P}^{\beta+1}_{i1}-\frac{\beta+1}{\beta}\left(X_{u_{i}}(P^{\beta}_{i1}-\bar{P}^{\beta}_{i1})+\bar{P}^{\beta}_{i1}\right)\right\\}$

For the $i^{th}$ group,

$\displaystyle
Y_{i}=\frac{\partial}{\partial\bm{\theta}}H_{i\beta}(\bm{\theta})$
$\displaystyle=(\beta+1)P^{\beta}_{i1}\frac{\partial(P_{i1})}{\partial\bm{\theta}}-(\beta+1)\bar{P}^{\beta}_{i1}\frac{\partial(P_{i1})}{\partial\bm{\theta}}-\frac{\beta+1}{\beta}\bigg{[}X_{ui}\bigg{.}$
$\displaystyle\qquad\qquad\left.\left.\left(\beta
P^{\beta-1}_{i1}\frac{\partial(P_{i1})}{\partial\bm{\theta}}\right.+\beta\bar{P}^{\beta-1}_{i1}\frac{\partial(P_{i1})}{\partial\bm{\theta}}\right)-\beta\bar{P}^{\beta-1}_{i1}\frac{\partial(P_{i1})}{\partial\bm{\theta}}\right]$
$\displaystyle=(\beta+1)\left[P^{\beta}_{i1}-\bar{P}^{\beta}_{i1}-X_{ui}\left(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1}\right)+\bar{P}^{\beta-1}_{i1}\right]\frac{\partial(P_{i1})}{\partial\bm{\theta}}$
$\displaystyle=(\beta+1)\left[(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})(P_{i1}-X_{ui})\right]\frac{\partial}{\partial\bm{\theta}}P_{i1}.$
$\displaystyle\text{Therefore,}\;\;Y_{is}$
$\displaystyle=(\beta+1)\left[(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})(P_{i1}-X_{ui})\right]\frac{\partial}{\partial\bm{\theta}_{s}}P_{i1}$
$\displaystyle\noindent\text{where}\qquad Y_{i}$
$\displaystyle=(Y_{i1},Y_{i2},\dots,Y_{is},\dots,Y_{2J})^{{}^{\prime}}.$
$\displaystyle\text{Here, }\quad E(Y_{i})=0,$ $\displaystyle Var(Y_{is})$
$\displaystyle=(\beta+1)^{2}\left[(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s}}\right]^{2}Var(X_{u_{i}})$
$\displaystyle=(\beta+1)^{2}\left[(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s}}\right]^{2}P_{i1}\bar{P}_{i1}$

and,

$\displaystyle Cov(Y_{is_{1}},Y_{is_{2}})$
$\displaystyle=\left[(\beta+1)(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})\right]^{2}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{1}}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{2}}}Var(X_{u_{i}})$
$\displaystyle=\left[(\beta+1)(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})\right]^{2}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{1}}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{2}}}P_{i1}\bar{P}_{i1}$

Therefore,

$\displaystyle
Cov\left(\frac{\partial}{\partial\bm{\theta}}H_{i\beta}(\bm{\theta})\right)$
$\displaystyle=(\beta+1)^{2}K_{i\beta}(\bm{\theta})$ $\displaystyle
K_{i\beta}(\bm{\theta})_{ss}$
$\displaystyle=\frac{Var(Y_{is})}{(\beta+1)^{2}}\quad\text{(variance term)}$
$\displaystyle K_{i\beta}(\bm{\theta})_{s_{1}s_{2}}$
$\displaystyle=\frac{Cov(Y_{is_{1}},Y_{is_{2}})}{(\beta+1)^{2}}\quad\text{(Covariance
term)}$

Applying Central Limit Theorem (CLT), when $k_{i}\to\infty$ for $i=1,\dots,I,$

$\displaystyle\sqrt{k_{i}}\left(\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}\frac{\partial}{\partial\bm{\theta}}H_{i\beta}(\bm{\theta})\right)\sim
N\left(0,(\beta+1)^{2}K_{i\beta}(\bm{\theta})\right)$

We denote,

$\displaystyle\nabla D_{K}(\bm{\theta})$
$\displaystyle=(D_{K,1}(\bm{\theta}),\dots,D_{K,s}(\bm{\theta}),\dots,D_{K,2J}(\bm{\theta}))^{{}^{\prime}}\;\text{where,}\;\;D_{K,s}(\bm{\theta})=\frac{\partial}{\partial\bm{\theta}_{s}}D_{K}$
$\displaystyle\text{define,}\;\;T_{\beta}$ $\displaystyle=-\sqrt{K}\nabla
D_{K}(\bm{\theta})=-\sqrt{K}\sum_{i=1}^{I}\frac{k_{i}}{K}\left(\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}\frac{\partial}{\partial\bm{\theta}}H_{i\beta}(\bm{\theta})\right)$
$\displaystyle\implies T_{\beta}$ $\displaystyle\sim
N\left(0,(\beta+1)^{2}\sum_{i=1}^{I}\frac{k_{i}}{K}K_{i\beta}(\bm{\theta})\right)\,,\,$
(19) $\displaystyle\text{where}\;\frac{k_{i}}{K}\;\text{is finite
when}\;k_{i}\xrightarrow[]{}\infty$

Now,

$\displaystyle\frac{\partial(D_{K,s_{1}}(\bm{\theta}))}{\partial\bm{\theta}_{s_{2}}}$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}\left(\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}\frac{\partial^{2}H_{i\beta(\bm{\theta})}}{\partial\bm{\theta}_{s_{2}}\partial\bm{\theta}_{s_{1}}}\right)$

where,

$\displaystyle\frac{\partial^{2}H_{i\beta}(\bm{\theta})}{\partial\bm{\theta}_{s_{2}}\partial\bm{\theta}_{s_{1}}}$
$\displaystyle=(\beta+1)\left[\left\\{(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})(P_{i1}-X_{ui})\frac{\partial^{2}(P_{i1})}{\partial\bm{\theta}_{s_{2}}\bm{\theta}_{s_{1}}}\right\\}\right.$
$\displaystyle\quad\qquad\qquad+\left\\{(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1})+(P_{i1}-X_{ui})(\beta-1)\right.$
$\displaystyle\qquad\qquad\qquad\left.\left.\left(P^{\beta-2}_{i1}-\bar{P}^{\beta-2}_{i1}\right)\right\\}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{2}}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{1}}}\right]$

Applying CLT, when $k_{i}\to\infty$ for $i=1,\ldots,I,$
$\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}X_{ui}\xrightarrow{p}P_{i1}$.

$\displaystyle\noindent\text{Therefore,}\quad\frac{1}{k_{i}}\sum_{u_{i}=1}^{k_{i}}\frac{\partial^{2}H_{i\beta}(\bm{\theta})}{\partial\bm{\theta}_{s_{2}}\partial\bm{\theta}_{s_{1}}}\xrightarrow{p}(\beta+1)\left(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1}\right)\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{2}}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{1}}}$
$\displaystyle\text{and,}\quad\frac{\partial(D_{K,s_{1}}(\bm{\theta}))}{\partial\bm{\theta}_{s_{2}}}\xrightarrow{p}\sum_{i=1}^{I}\frac{k_{i}}{K}(\beta+1)\left(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1}\right)\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{2}}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s_{1}}}$

Consider $\bm{\theta}_{0}$ to be the true value of parameters, then apply the
Taylor series expansion $D_{k,s}(\bm{\theta})$ around $\bm{\theta}_{0}$

$\displaystyle D_{K,s}(\bm{\theta})$
$\displaystyle=D_{K,s}(\bm{\theta}_{0})+\sum_{m=1}^{2J}\left.\frac{\partial(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{m}}\right|_{\bm{\theta}=\bm{\theta}_{0}}(\hat{\bm{\theta}}_{m}-\bm{\theta}_{0m})$
$\displaystyle\quad+\frac{1}{2}\sum_{l=1}^{2J}\sum_{m=1}^{2J}\left.\frac{\partial^{2}(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{l}\partial\bm{\theta}_{s}}\right|_{\bm{\theta}=\bm{\theta}^{\ast}}(\hat{\bm{\theta}}_{l}-\bm{\theta}_{0{l}})(\hat{\bm{\theta}}_{m}-\bm{\theta}_{0m})$
(20)

Since, $D_{K,s}(\hat{\bm{\theta}}_{\beta})=0$ it can be written that,

$\displaystyle-\sqrt{K}D_{K,s}(\bm{\theta}_{0})$
$\displaystyle=\sqrt{K}\sum_{m=1}^{2J}\left(\left.\frac{\partial(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{m}}\right|_{\bm{\theta}=\bm{\theta}_{0}}\right.$
$\displaystyle\qquad+\left.\frac{1}{2}\sum_{l=1}^{2J}\left.\frac{\partial^{2}(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{l}\partial\bm{\theta}_{m}}\right|_{\bm{\theta}=\bm{\theta}^{*}}(\hat{\bm{\theta}}_{l}-\bm{\theta}_{0{l}})\right)(\hat{\bm{\theta}}_{m}-\bm{\theta}_{0m})$
(21)

Denoting,

$\displaystyle A_{s,m}$
$\displaystyle=\left.\frac{\partial(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{m}}\right|_{\bm{\theta}=\bm{\theta}_{0}}+\frac{1}{2}\sum_{l=1}^{2J}\left.\frac{\partial^{2}(D_{K,s}(\bm{\theta}))}{\partial\bm{\theta}_{l}\partial\bm{\theta}_{m}}\right|_{\bm{\theta}=\bm{\theta}^{*}}(\hat{\bm{\theta}}_{l}-\bm{\theta}_{0{l}})$

Then,

$\displaystyle A_{s,m}$
$\displaystyle\xrightarrow{p}\sum_{i=1}^{I}\frac{k_{i}}{K}(\beta+1)\left(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1}\right)\frac{\partial(P_{i1})}{\partial\bm{\theta}_{m}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s}}$
where $\displaystyle\;\frac{k_{i}}{K}\;\text{is finite
when}\;k_{i}\xrightarrow[]{}\infty.$ $\displaystyle\implies$ $\displaystyle
A_{\beta}\xrightarrow{p}(\beta+1)J_{\beta}(\bm{\theta}_{0})$ (22)

where,$\;A_{\beta}$ is the $2J\times 2J$ matrix with $(s,m)$th element as
$A_{s,m}$ and

$\displaystyle J_{\beta}(\bm{\theta}_{0})$
$\displaystyle=\left[\left(\sum_{i=1}^{I}\frac{k_{i}}{K}\left(P^{\beta-1}_{i1}+\bar{P}^{\beta-1}_{i1}\right)\frac{\partial(P_{i1})}{\partial\bm{\theta}_{m}}\frac{\partial(P_{i1})}{\partial\bm{\theta}_{s}}\right)_{s,m}\right]$

We can write,

$\displaystyle-\sqrt{K}D_{K,s}(\bm{\theta}_{0})$
$\displaystyle=\sqrt{K}\sum_{s=1}^{2J}(\hat{\bm{\theta}}_{s}-\bm{\theta}_{0s})A_{s,m}$

Therefore, it can be expressed as

$\displaystyle
T_{\beta}=\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})A_{\beta}$
$\displaystyle\implies\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})=A_{\beta}^{-1}T_{\beta}$
$\displaystyle\implies\sqrt{K}(\hat{\bm{\theta}}_{\beta}-\bm{\theta}_{0})=A_{\beta}^{-1}T_{\beta}\sim
N\left(0,J_{\beta}(\bm{\theta}_{0})^{-1}K_{\beta}(\bm{\theta}_{0})J_{\beta}(\bm{\theta}_{0})^{-1}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(\text{by
equation \eqref{a1}})$

where,

$\displaystyle K_{\beta}(\bm{\theta}_{0})$
$\displaystyle=\sum_{i=1}^{I}\frac{k_{i}}{K}K_{i\beta}(\bm{\theta})\;;\;\frac{k_{i}}{K}\;\text{is
finite, when}\;k_{i}\xrightarrow[]{}\infty.$

## References

  * [1] L. Wang, Inference of constant-stress accelerated life test for a truncated distribution under progressive censoring, Applied Mathematical Modelling 44 (2017) 743–757.  doi:https://doi.org/10.1016/j.apm.2017.02.011.
  * [2] A. Xu, J. Fu, Y. Tang, Q. Guan, Bayesian analysis of constant-stress accelerated life test for the weibull distribution using noninformative priors, Applied Mathematical Modelling 39 (20) (2015) 6183–6195.  doi:https://doi.org/10.1016/j.apm.2015.01.066.
  * [3] S. Roy, Bayesian accelerated life test plans for series systems with weibull component lifetimes, Applied Mathematical Modelling 62 (2018) 383–403.  doi:https://doi.org/10.1016/j.apm.2018.06.007.
  * [4] P. W. Srivastava, N. Jain, Optimum ramp-stress accelerated life test for m identical repairable systems, Applied Mathematical Modelling 35 (12) (2011) 5786–5793.  doi:https://doi.org/10.1016/j.apm.2011.05.003.
  * [5] N. Balakrishnan, M. H. Ling, Em algorithm for one-shot device testing under the exponential distribution, Computational Statistics & Data Analysis 56 (3) (2012) 502–509.  doi:https://doi.org/10.1016/j.csda.2011.09.010.
  * [6] 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.
  * [7] N. Balakrishnan, H. Y. So, M. H. Ling, A bayesian approach for one-shot device testing with exponential lifetimes under competing risks, IEEE Transactions on Reliability 65 (1) (2015) 469–485.  doi:10.1109/TR.2015.2440235.
  * [8] 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.
  * [9] N. Balakrishnan, E. Castilla, Em-based likelihood inference for one-shot device test data under log-normal lifetimes and the optimal design of a csalt plan, Quality and Reliability Engineering International 38 (2) (2022) 780–799.  doi:https://doi.org/10.1002/qre.3014.
  * [10] Y. Lan, L. M. Leemis, The logistic–exponential survival distribution, Naval Research Logistics (NRL) 55 (3) (2008) 252–264.  doi:https://doi.org/10.1002/nav.20279.
  * [11] S. Ali, S. Dey, M. H. Tahir, M. Mansoor, Two-parameter logistic-exponential distribution: Some new properties and estimation methods, American Journal of Mathematical and Management Sciences 39 (3) (2020) 270–298.  doi:https://doi.org/10.1080/01966324.2020.1728453.
  * [12] N. Balakrishnan, D. Kundu, Birnbaum-saunders distribution: A review of models, analysis, and applications, Applied Stochastic Models in Business and Industry 35 (1) (2019) 4–49.  doi:https://doi.org/10.1002/asmb.2348.
  * [13] 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.
  * [14] 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.
  * [15] N. Balakrishnan, E. Castilla, N. Martín, L. Pardo, Robust estimators and test statistics for one-shot device testing under the exponential distribution, IEEE transactions on Information Theory 65 (5) (2019) 3080–3096.  doi:10.1109/TIT.2019.2903244.
  * [16] 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.
  * [17] A. Basu, A. Mandal, N. Martin, L. Pardo, Testing statistical hypotheses based on the density power divergence, Annals of the Institute of Statistical Mathematics 65 (2013) 319–348.  doi:https://doi.org/10.1007/s10463-012-0372-y.
  * [18] 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.
  * [19] M. 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.
  * [20] Surveillance, epidemiology, and end results (seer) program (www.seer.cancer.gov) seer*stat database: Incidence - seer research data, 8 registries, nov 2021 sub (1975-2019) - linked to county attributes - time dependent (1990-2019) income/rurality, 1969-2020 counties, national cancer institute, dccps, surveillance research program, released april 2022, based on the november 2021 submission.
  * [21] J. Dik, M. De Gunst, The distribution of general quadratic forms in norma, Statistica Neerlandica 39 (1) (1985) 14–26.  doi:https://doi.org/10.1111/j.1467-9574.1985.tb01121.x.
  * [22] H. Jh, Adaptation in natural and artificial systems, Ann Arbor (1975).
  * [23] D. E. Golberg, Genetic algorithms in search, optimization, and machine learning, Addion wesley 1989 (102) (1989) 36.
  * [24] S.-L. Ou, L.-y. D. Liu, Y.-C. Ou, Using a genetic algorithm-based raroc model for the performance and persistence of the funds, Journal of Applied Statistics 41 (5) (2014) 929–943.  doi:https://doi.org/10.1080/02664763.2013.856870.
  * [25] R. L. Haupt, S. E. Haupt, Practical genetic algorithms, John Wiley & Sons, 2004\.
  * [26] M. Babaei, M. Sheidaii, Optimal design of double layer scallop domes using genetic algorithm, Applied Mathematical Modelling 37 (4) (2013) 2127–2138.  doi:https://doi.org/10.1016/j.apm.2012.04.053.
  * [27] C.-D. Liou, K.-H. Wang, M.-W. Liou, Genetic algorithm to the machine repair problem with two removable servers operating under the triadic (0, q, n, m) policy, Applied Mathematical Modelling 37 (18-19) (2013) 8419–8430.  doi:https://doi.org/10.1016/j.apm.2013.03.036.
  * [28] M. Soolaki, I. Mahdavi, N. Mahdavi-Amiri, R. Hassanzadeh, A. Aghajani, A new linear programming approach and genetic algorithm for solving airline boarding problem, Applied Mathematical Modelling 36 (9) (2012) 4060–4072.  doi:https://doi.org/10.1016/j.apm.2011.11.030.
  * [29] J. López-Campos, A. Segade, J. Fernández, E. Casarejos, J. Vilán, Behavior characterization of visco-hyperelastic models for rubber-like materials using genetic algorithms, Applied Mathematical Modelling 66 (2019) 241–255.  doi:https://doi.org/10.1016/j.apm.2018.08.031.
  * [30] B. Chakraborty, P. Chaudhuri, On the use of genetic algorithm with elitism in robust and nonparametric multivariate analysis, Austrian journal of statistics 32 (1&2) (2003) 13–27.  doi:https://doi.org/10.17713/ajs.v32i1\&2.447.
  * [31] N. Balakrishnan, M. H. Ling, Multiple-stress model for one-shot device testing data under exponential distribution, IEEE transactions on reliability 61 (3) (2012) 809–821.  doi:10.1109/TR.2012.2208301.
  * [32] 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.