Stage Life Testing with Missing Stage Information : an EM-Algorithm Approach

Vol. 20, No. 01, pp 123-152 DOI:10.52547/jirss.20.1.123

Stage Life Testing with Missing Stage Information - an EM- Algorithm Approach

Erhard Cramer¹, and Benjamin Laumen¹

1Institute of Statistics, RWTH Aachen University, D-52062 Aachen, GERMANY.

Received: 23/11/2020, Revision received: 07/02/2021, Published online: 03/04/2021 Abstract. We consider a stage life testing model and assume that the information at which levels the failures occurred is not available. In order to find estimates for the lifetime distribution parameters, we propose an EM-algorithm approach which interprets the lack of knowledge about the stages as missing information. Furthermore, we illustrate the implementation difficulties caused by an increasing number of stages.

The study is supplemented by a data example as well as simulations.

Keywords.EM-Algorithm, Exponential Distribution, Missing Information, Progressive Censoring, Stage Life Testing, Weibull Distribution

MSC:62N05, 62F10.

1 Introduction

The notion of stage life testing (SLT) has been proposed in Laumen (2017) and Laumen and Cramer (2019b, 2021a) as an extension of progressive Type-I censoring (for a version with random stage changing times, see Laumen and Cramer (2021b)). The approach provides models that allow to incorporate additional life time information

Corresponding Author: Erhard Cramer (erhard.cramer@rwth-aachen.de) Benjamin Laumen (benjamin.laumen@rwth-aachen.de)

of progressively censored objects by performing additional testing of the removed items (for comments in this direction, cf., e.g., Balakrishnan and Aggarwala 2000, p. 3, Balakrishnan et al. 2011, p. 336, Balakrishnan and Cramer (2014, 2021), Cramer (2017)).

In fact, it is assumed that the progressively censored objects are further tested but under different conditions (called stages) whereas the remaining items are continued to be monitored under the initial conditions. An illustration of this concept is depicted in Figure 1 fork−1 stage changing timesτ1 <· · ·< τk−1and an effectively applied stage changing plan (r^?₁, . . . ,r^?_k−1). The experimental design of such a life test requires that r^?_j objects are randomly withdrawn from the life test at time τj, 1 ≤ j ≤k−1. Notice that the stage changing plan may be specified in different ways. The life span test is terminated when either the last remaining object in the life test fails or the last item is removed from the test. Notice that we do not assume that the experiment is stopped at timeτk−1as is commonly done in progressive Type-I censoring (cf. progressive Type-I censoring with fixed censoring times discussed in Laumen and Cramer 2019a).

τ1 τi τk−1

s₀ s₁ s_i s_k−1

R^?₁ =r^?₁

R^?_i =r^?_i

R^?_k−1=r^?_k−1

D1 D2 Di D_i+1 Dk−1 Dk

time

stage

Figure 1: Illustration of k-step SLT with stage-changing timesτ1 < · · · < τk−1 and an effectively applied stage changing plan (r^?₁, . . . ,r^?_k−1).

As can be seen from Figure 1, the sample is split at change-timeτ1in these items

which are tested under the initial conditions (stages₀), and those which are tested on stages₁with a possibly different load (which, of course, may be higher or lower). This process is continued for the items remaining on stages0at the following stage changing timesτ2, . . . , τk−1. The adaption of the load is modelled by the cumulative exposure model approach (see, e.g., Kundu and Ganguly (2017)). Furthermore, it should be mentioned that SLT can also be interpreted as a modified simple step-stress model (see Balakrishnan (2009), Kundu and Ganguly (2017)) where only a proportion of the items tested under initial conditions is selected for testing under other stress conditions. In this regard, SLT can also be interpreted as a model of accelerated life testing. For details, we refer to Laumen and Cramer (2019b, 2021a).

This paper is organized as follows. In Section 2, we introduce briefly the SLT model and recall some results presented in Laumen and Cramer (2019b, 2021a). In Section 3, we address maximum likelihood estimation in the SLT model under missing stage information and present the likelihood fork-step SLT. Afterwards, we illustrate the EM-algorithm approach for the 2- and 3-step SLT. In particular, we start with an exponential distribution on both stages in Section 3.2.1. In Section 3.2.2, we consider the combination of a Weibull and an exponential distribution for the stages of the SLT model. In Section 4.1, we provide an illustrative example. Finally, we present the results of a simulation study in Section 4.2.

2 SLT Model

Our discussion is based on the notation ofk-step SLT order statistics as introduced in Laumen and Cramer (2021a). Assume that n identical objects with iid lifetimes X₁, . . . ,Xnare placed on a life test at the initial stages0. At thejth prefixed stage-change timeτj,R^?_j ≥0 of the surviving items are randomly withdrawn (if possible) from the sample and further tested on the changed stages_j, the remaining objects are further tested under the conditions of stages₀, 1≤ j≤k−1. The life test is terminated when allnobjects have failed.

The (random) numbers of failures observed on the initial stages₀in the intervals (−∞, τ1], (τ1, τ2], . . . , (τk−2, τk−1], (τk−1,∞),

are denoted byD₁,D₂, . . . ,D_k−1,D_k, whereτ0=−∞,τk=∞(see Figure 1).

LetM=D•k=Pk

j=1D_jbe the total number of observations failed on levels₀where D•0 =0. Arranging the data according to the stage levelss₀,s₁, . . . ,s_k−1, the observed

failure times on these levels are denoted by Y_h,Dh =

YD•h−1+1:M:n, . . . ,YD•h:M:n

denote the ordered failure times observed in the interval (τh−1, τh] on stages₀,h= 1, . . . ,k. Notice thatY_1,D1, . . . ,Y_k,Dk forms a progressively Type-I censored sample with fixed censoring times as discussed in Laumen and Cramer (2019a). This connection is reflected by the notation used;

Z_j,R?

j =

Z_j,1:R?

j, . . . ,Z_j,R? j:R^?_j

withY_D•j:M:n≤τj<Z_j,1:R?

j denote the ordered failure times observed on stages_j, j=1, . . . ,k−1.

The order statistics on stage s₀ and the order statistics on the stages s₁, . . . ,s_k−1 are represented by the vectorsY =

Y_1,D1, . . . ,Y_k,Dk

andZ= Z_1,R?

1, . . . ,Z_k−1,R^?_k−1

, respectively.

The complete sample is given by (Y,Z). Notice that the partitioning of the sample induces the assignment of failures to stages. But, in the following, we assume that the information about this assignment is not available. In order to describe the situation, we introduce random indicators Σ1, . . . ,Σn which provide the information about the stage, that is,

Σi =



















0, objectihas failed on stages₀ 1, objectihas failed on stages₁

...

k−1, objectihas failed on stages_k−1

, i=1, . . . ,n;

indicates whether an object has failed on stages₀, . . . ,s_k−1. Thus, each observationX_iis accompanied by an indicatorΣiso that we observe a pair (X_i,Σi) whereΣiprovides the information about the stage. By considering the order statisticsX^∗=(X_1:n, . . . ,X_n:n) of the sampleX, the stage indicators can be interpreted as a concomitant (see David and Nagaraja (1998), Bairamov and Eryılmaz (2006), Izadi and Khaledi (2007), Balakrishnan and Cramer (2014)), that is, we get the (bivariate) ’ordered’ sample

(X1:n,Σ[1:n]), . . . ,(Xn:n,Σ[n:n]).

Due to the construction of the sample, we know that X_h:n = Y_h:M:n and Σ[h:n] = 0, h=1, . . . ,D₁. For brevity, we subsequently writeX^∗_i =X_i:n,Σ^∗_i = Σ[i:n], 1≤i≤n. In the present discussion, we use the following notation and assumptions:

Throughout the manuscript, we use the notationw_j =(w₁, . . . ,w_j) for the vector of jcomponentsw₁, . . . ,w_jas well asw•j =Pj

i=1w_ifor their partial sum.

The random censoring numberR^?_j is generated by a (deterministic) function%jof the failures observed beforeτj, that is,R^?_j =%j(D_j), 1≤ j≤k−1. In the following, these functions may be chosen according to the needs of the experimenter. In Laumen and Cramer (2019b, 2021a) two options to generate the withdrawal numberR^?_j, j=1, . . . ,k−1, have been proposed:

%j(d_j)=









 jπj·

n−d•j−Pj−1

i=1%i(d_i)k

, Type-P,

minn

R⁰_j,max{n−d•j−Pj−1

i=1%i(d_i),0}o

, Type-M, (2.1) where d_j = (d1, . . . ,dj), 1 ≤ j ≤ k, and btc is defined as the largest integer not exceeding t ∈ R. The proportionsπj ∈ (0,1) as well as the numbers R⁰_j ∈ N, 1 ≤ j ≤ k− 1, are specified in advance, respectively. These choices can be interpreted as follows:

◦ Type-P: Atτj, a (prefixed) proportionπjof the surviving objects is selected for testing on stages_j, j=1, . . . ,k−1.

◦ Type-M: The second way to generateR^?_j, j = 1, . . . ,k−1, is similar to the censoring procedure of progressive censoring with fixed failure times (see Laumen and Cramer (2019a)). Given a prefixed number R⁰_j, it is intended to select atτjas many items as possible (at mostR⁰_j) for testing on stages_j,

j=1, . . . ,k−1.

By construction,D_kis a (deterministic) function ofD₁, . . . ,D_k−1

and (R^?₁, . . . ,R^?_k−1) , that is,

D_k=n−D•k−1−R^?_•k−1=n−D•k−1−

k−1

X

i=1

%i(D_i).

The support of (D₁, . . . ,D_k) is represented by the set D_(k)=n

a_k∈N^k₀

ai ≤max{n−a^•i−1−r^?•i−1(a_i−1),0}, i=1, . . . ,k−1,

a_k =max{n−a•k−1−r^?_•k−1(a_k−1),0}o

; P^∗=(π1, . . . , πk−1) denotes the proportional stage-changing plan (see Type-P);

R∗⁰ = (R⁰₁, . . . ,R⁰_k−1) denotes the initially planned stage-changing plan (see Type- M);

F_jdenotes the absolutely continuous cumulative distribution function with density function f_jon stages_j, j∈ {0, . . . ,k−1};

The cumulative exposure model is supposed to hold, that is, the distribution function F_j on stage j is connected to the baseline distribution function F₀ as follows: For 1≤ j≤k−1, we have valuesv₁, . . . ,v_k−1such thatv_jis the solution of the equation

F₀(τj)=F_j(v_j).

Hence, given the stage-changing time τj, the cumulative distribution function and the corresponding probability density functions of a test unit on stages_jare given by

F_0,j(t)=











F₀(t), t≤τj

Fj(t+vj−τj), τj <t, f_0,j(t)=











f₀(t), t≤τj

fj(t+vj−τj), τj <t. (2.2) Details on the cumulative exposure model can be found in Kundu and Ganguly (2017, Chapter 2).

Laumen and Cramer (2021a) have obtained the joint density function ofk-step SLT order statistics as given in Theorem 2.1.

Theorem 2.1. LetY_1,D₁, . . . ,Y_k,Dk andZ_1,R^?

1, . . . ,Z_k−1,R^?_k−1 be k-step SLTOSs and let Fi be an absolutely continuous cumulative distribution function with density function f_i, i∈ {0, . . . ,k− 1}. Further, let−∞=τ0< τ1 <· · ·< τk−1< τk =∞.

Then, the joint density function f_1...n^Y,Z,DkofY=

Y_1,D1, . . . ,Y_k,Dk ,Z=

Z_1,R?

1, . . . ,Z_k−1,R^?_k−1

andD_k = (D₁, . . . ,D_k) (w.r.t. the product of the n dimensional Lebesgue measure and the k dimensional counting measure) is given by

f^Y,Z,Dk(y,z,d_k)= Yk

j=1

n−d•j−1−r^?_•j−1 d_j

! d_j!

( ^d•j Y

i=d^•j−1+1

f₀(y_i:m:n)11_(τj−1,τj](y_i:m:n) )

× (k−1

Y

a=1

r^?_a!

r^?_a

Y

b=1

f_a(z_a,b:r?

a +v_a−τa)11_(τa,∞)(z_a,b:r? a)

)

, (2.3)

ford_k ∈D_(k), y=

y_1,d1, . . . ,y_k,dk

with y_j,dj =

y_d•j−1+1:m:n, . . . ,y_d•j:m:n

, 1≤ j≤k, z=

z_1,r?

1, . . . ,z_k−1,r^?_k−1

with z_a,r?

a =

z_a,1:r^?_a, . . . ,z_a,r^?_a:r^?_a

, 1≤a≤k−1.

Notice that (Y,Z,D_k) are determined by (X,Σ) and vice versa. Therefore, using the above notation, the density function in (2.3) can equivalently be written as

f^X,Σ(x,σ)=cY

i:σi=0

f₀(x_i)

k−1

Y

j=1

Y

i:σi=j

f_j(x_i+v_j−τj), (2.4)

with d_j = P

i:σi=011_(τj−1,τj](x_i), j = 1, . . . ,k, and normalizing constant c = c(d_k). The corresponding density function of the ordered data is given by

f^X∗,Σ∗(x^∗,σ^∗)=c^∗ Y

i:σ^∗_i=0

f₀(x^∗_i)

k−1

Y

j=1

Y

i:σ^∗_i=j

f_j(x^∗_i +v_j−τj), (2.5)

withd_j =P

i:σ^∗_i=011_(τj−1,τj](x^∗_i),j=1, . . . ,k, and normalizing constantc^∗=c^∗(d_k).

3 Maximum Likelihood Estimation in SLT under Missing Stage Information

3.1 k-step SLT

Given a statistical model with parameter vectorθ = (ϑj)_j=0,...,k−1 ∈ Θ =

×

likelihood function is obtained from (2.5) as L(ϑ0, . . . ,ϑk−1|x^∗,σ^∗)= f_ϑ^X∗,Σ∗

0,...,ϑk−1(x^∗,σ^∗)∝ Y

i:σ^∗_i=0

f_0,ϑ0(x^∗_i)

k−1

Y

j=1

Y

i:σ^∗_i=j

f_j,ϑj(x^∗_i +v_j−τj), (3.1)

with density functions f_j,ϑj, j=0, . . . ,k−1. Assumingσ^∗as known, that is, we know which failure occurred on which stage, likelihood inference has been discussed for exponential and Weibull distributions in Laumen and Cramer (2019b, 2021a). However, if the stage informationΣ^∗ = σ is not available, the respective likelihood is obtained from (3.1) by summing over all possible values ofσso that the corresponding likelihood is obtained as marginal density function ofX^∗. It reads

L_MI(ϑ0, . . . ,ϑk−1|x^∗)=X

σ

f_ϑ^X∗,Σ∗

0,...,ϑk−1(x^∗,σ)∝X

σ

Y

i:σi=0

f_0,ϑ0(x^∗_i)

k−1

Y

j=1

Y

i:σi=j

f_j,ϑj(x^∗_i +v_j−τj). (3.2)

Notice that the normalizing constant does not depend on the failure assignmentσ. Due to the sum representation, direct optimization of the likelihood (3.2) may be quite hard.

However, treating the values ofΣ^∗=σas missing information, allows to address the maximization problem by an EM-algorithm type approach. In case of an exponential distributionExp(ϑ) with meanϑ > 0, that is, the probability density function and the cumulative distribution function are given by

f(x)= 1

ϑe^−x/ϑ11_(0,∞)(x), F(x)=

1−e^−x/ϑ

11_(0,∞)(x), ϑ >0,x∈R.

this approach is particularly useful since the resulting MLEs (under complete information) are available in a closed form representation. Suppose the lifetime distributions on stage s_jare exponential with meanϑj, j=0, . . . ,k−1. Then, one getsv_j =τ1ϑj

ϑ0. Furthermore, from (3.1), the likelihood function is given by

L(ϑ0, . . . , ϑk−1|x^∗,σ)

=c^∗ϑ⁻₀^m

k−1

Y

a=1

ϑ⁻a^r?j exp (

−

k−1

X

a=1

r^?_aτa

1 ϑ0

− 1 ϑa

− 1 ϑ0

Xn j=1

x^∗_j11{0}(σj)−

k−1

X

a=1

1 ϑa

Xn j=1

x^∗_j11{a}(σj) )

. (3.3) Using results of Laumen and Cramer (2021a), the corresponding MLEs are given by

bϑ0= 1 M

Xⁿ

i=1

X^∗_i11{0}(Σi)+

k−1

X

j=1

R^?_j τj

, (3.4)

and, providedR^?_h >0,

bϑh= 1 R^?_h

Xn i=1

(X^∗_i −τh)11{h}(Σi), h∈ {1, . . . ,k−1}. (3.5)

3.2 EM-Algorithm for Two-Step SLT 3.2.1 Exponential-Exponential Case

Let the sample (X^∗₁,Σ^∗₁), . . . ,(X^∗_n,Σ^∗_n) be incomplete in the sense that the information on which stage the observed failures occurred is missing, that is, the values of the stage

indicatorsΣ^∗’s are not available. The design of the life test is still known, i.e., we know n, τ1, π1, andR⁰₁, respectively. Furthermore, the ordered failure timesx^∗₁ ≤ . . . ≤ x^∗_n are observed. This situation is depicted in Figure 2. The question marks indicate that the assignment of the failure to the stage is not known. Notice that the failure timesx^∗_d

1+1, . . . ,x^∗_ncould have been observed on stagess₀ors₁. In order to estimate the unknown parametersϑ0andϑ1, we utilize an EM-algorithm (see, e.g., Dempster et al.

1977).

x^∗₁ x^∗_d

1 τ1 x^∗_d

1+1 x^∗_d

1+2 x^∗_d

1+3 x^∗_n

s₀

s₁ ^r1

d2=n−d1−r1

d1

? ? ? ?

time

stage

Figure 2: 2-step SLTOSs with missing information. The observationsx^∗_d

1+1, . . . ,x^∗_ncould have been observed on stagess0ors1.

Using the observationsx^∗₁, . . . ,x^∗_nand the design of the life test, the quantitiesd₁,d₂, andr^?₁ are obtained as

d₁= Xn

i=1

11₍−∞,τ1](x^∗_i), r^?₁ =%(d1)=











bπ1·(n−d₁)c, Type-P minn

n−d₁,R⁰₁o

, Type-M, d₂=n−d₁−r^?₁. Remark1. (i) Whend₁=n, we know that all observed failures occurred on stages₀.

Hence,bϑ0can be determined from equation (3.4) and the MLEϑb1does not exist.

(ii) When d₁ ≥ 0, d₂ > 0, and r^?₁ = 0, we know that all observed failures after τ1

occurred on stages₀. Thus, the MLEϑb1does not exist andbϑ0can be determined from equation (3.4).

(iii) Whend₁ > 0, d₂ = 0, and r^?₁ > 0, we know that all observed failures after τ1

occurred on stages₁. Therefore, bϑ0 andϑb1 can be obtained from (3.4) and (3.5), respectively.

In order to define the EM-algorithm in the present situation, we consider first the scenario of available stage information (X^∗,Σ^∗). Letθ = (ϑ0, ϑ1) andθ^(t) =

ϑ^(t)₀ , ϑ^(t)₁ , t ∈ N0. Using equation (3.3), the likelihood function for the complete data (X^∗,Σ^∗) = (x^∗,σ^∗) can be written in the form

L(θ|x^∗,σ^∗)=c^∗

d1

Y

i=1

f0(x^∗_i) Yn j=d1+1

h11{0}(σ^∗_j) f0(x^∗_j)+11{1}(σ^∗_j)f1(x^∗_j+v1−τ1)i

= c^∗ ϑ^d₀¹ exp

(

− 1 ϑ0

d₁

X

i=1

x^∗_i )

exp ( n

X

j=d1+1

11{0}(σ^∗_j)

− x^∗_j ϑ0

−log(ϑ0) )

×exp ( n

X

j=d1+1

11{1}(σ^∗_j)

− 1 ϑ1

x^∗_j+τ1 ϑ1

ϑ0

−τ1

−log(ϑ1) )

,

where Pn

j=d₁+1σ^∗_j = r^?₁. Notice that this simplification is possible since σ^∗_j ∈ {0,1} for k = 2. This yields the log-likelihood function for the complete data (X^∗,Σ^∗) = (x^∗,σ^∗) given by

`(θ|x^∗,σ^∗)= log(c^∗)−d₁log(ϑ0)− 1 ϑ0

d1

X

i=1

x^∗_i + Xn j=d1+1

11{0}(σ^∗_j)

− x^∗_j ϑ0

−log(ϑ0)

+ Xn j=d1+1

11{1}(σ^∗_j)

− 1 ϑ1

x^∗_j+τ1ϑ1

ϑ0

−τ1

−log(ϑ1)

. (3.6)

In order to perform theE-stepof the EM-algorithm, we have to calculate the expectation of`(θ|X^∗,Σ^∗) w.r.t.P^{Σ∗|X∗=x∗}

θ^(t) for the current estimateθ^(t). Thus, we get with (3.6) Q

θ;θ^(t)

=E_θ(t)

h`

θ|X^∗,Σ^∗

X^∗=x^∗i

= log(c^∗)−d₁log(ϑ0)− 1 ϑ0

d1

X

i=1

x^∗_i + Xn j=d1+1

P^(t)_0,j

− x^∗_j ϑ0

−log(ϑ0)

+ Xn

j=d₁+1

P^(t)_1,j

− 1 ϑ1

x^∗_j+τ1ϑ1

ϑ0

−τ1

−log(ϑ1) ,

where

P^(t)_s,j=E_θ(t)

h11{s}(Σ^∗_j)

X^∗=x^∗i

=P_θ(t)

Σ^∗_j=s

X^∗=x^∗

= f^Σ

∗ j,X^∗ θ^(t) (s,x^∗)

f^X∗

θ^(t)(x^∗) , s∈ {0,1}, j∈ {d₁+1, . . . ,n}, does not depend onθ. To proceed further, we introduce the sets

S_n= (

σ^∗∈ {0,1}ⁿ

σ^∗₁ =0, . . . , σ^∗_d1 =0, Xn j=d1+1

σ^∗_j =r^?₁ )

,

S_n^j,0=S_n∩n σ^∗_j=0o

and S_n^j,1 =S_n∩n σ^∗_j =1o

, j∈ {d₁+1, . . . ,n}. Hence, we have

f^Σ

∗ j,X^∗

θ^(t) (0,x^∗)= 1 ϑ^(t)₀ exp

− x^∗_j ϑ^(t)₀

Y^d1

h=1

1 ϑ^(t)₀ exp

− x^∗_h ϑ^(t)₀

× X

σn^∗∈S_n^j,0

Yn

i=d1+1 i,j

f^0,1

θ^(t)(x^∗_i, σ^∗_i)

!

, (3.7a)

f^Σ

∗ j,X^∗

θ^(t) (1,x^∗)= 1 ϑ^(t)₁ exp

− 1 ϑ^(t)₁

x^∗_j+τ1

ϑ^(t)₁ ϑ^(t)₀ −τ1

Y^d1

h=1

1 ϑ^(t)₀ exp

− x^∗_h ϑ^(t)₀

× X

σn^∗∈S_n^j,1

Yn

i=d₁+1 i,j

f_θ^0,1_(t)(x^∗_i, σ^∗_i)

!

, (3.7b)

j∈ {d₁+1, . . . ,n}, and f_θ^X_(t)^∗(x^∗)=

d1

Y

h=1

1 ϑ^(t)₀ exp

− x^∗_h ϑ^(t)₀

X

σn^∗∈S_n

Yn

i=d₁+1

f_θ^0,1_(t)(x^∗_i, σ^∗_i)

!

, (3.8)

with f^0,1

θ^(t)(x^∗_i, σ^∗_i)=11^{₀^}(σ^∗_i) 1 ϑ^(t)₀ exp

− x^∗_i ϑ^(t)₀

+11^{₁^}(σ^∗_i) 1 ϑ^(t)₁ exp

− 1 ϑ^(t)₁

x^∗_i +τ1

ϑ^(t)₁ ϑ^(t)₀

−τ1

.

For theM-step, we have to maximizeQ θ;θ^(t)

w.r.t.θ. First, we have Q

θ;θ^(t)

=Qe₀(ϑ0)+Qe₁(ϑ1), where

Qe0(ϑ0)= log(c^∗)−

d1+d^(t)₂

log(ϑ0)− 1 ϑ0

d1

X

i=1

x^∗_i − 1 ϑ0

Xn j=d1+1

P^(t)_0,jx^∗_j−r^?₁^(t)τ1

ϑ0 ,

Qe₁(ϑ1)= −r^?₁^(t)log(ϑ1)− 1 ϑ1

Xn

j=d₁+1

P^(t)_1,jx^∗_j+r^?₁^(t)τ1

ϑ1 , with

d^(t)₂ = Xn j=d1+1

P^(t)_0,j=E_θ(t)

h Xⁿ

j=d1+1

11{0}

Σ^∗_j

X^∗=x^∗i

=n−d₁−r^?₁ and r^?₁^(t)= Xn j=d1+1

P^(t)_1,j =r^?₁.

Thus, the updated estimatesϑ^(t+1)₀ andϑ^(t+1)₁ in the (t+1)th iteration step are given by ϑ^(t₀⁺¹⁾= 1

n−r^?₁

d1

X

i=1

x^∗_i + Xn j=d1+1

P^(t)_0,jx^∗_j+r^?₁τ1

!

and ϑ^(t₁⁺¹⁾= 1 r^?₁

Xn j=d1+1

P^(t)_1,jx^∗_j−r^?₁τ1

! ,

respectively. Possible initial values for the EM-algorithm are given by ϑ⁽⁰⁾₀ = 1

d₁

d1

X

i=1

x^∗_i +(n−d₁)τ1

!

and ϑ⁽⁰⁾₁ = 1 n−d₁

Xn

j=d₁+1

x^∗_j−(n−d₁)τ1

! ,

whend₁>0, and by ϑ⁽⁰⁾₀ =ϑ⁽⁰⁾₁ = 1

n Xn

j=1

x^∗_j−nτ1

! ,

whend₁ =0.

Alternatively, the MLEs under missing stage information (IMLE) can be computed by direct maximization of f^X

∗ n

θ^(t) w.r.t.θ^(t)(see equations (3.2) and (3.8)). Notice that this function is the marginal density function ofX^∗. In Section 4.2, we compare the results of both approaches showing that they lead almost to the same estimates.

3.2.2 Exponential-Weibull Case

We assume the same situation of missing information as in Section 3.2.1 but with Weibull lifetimes on stage s0. The probability density function and the cumulative distribution function of the Weibull distributionWei(ϑ, β) are given by

f(x)= β

ϑx^β−1e^−xβ/ϑ11_(0,∞)(x), F(x)=

1−e^−xβ/ϑ

11_(0,∞)(x), ϑ >0, β >0, x∈R.

We assume that the lifetimes on stages₀ are Wei(ϑ0, β)-distributes whereas they are Exp(ϑ1)-distributed on stages₁. Therefore,v₁=τ^β₁^ϑ_ϑ¹₀ and, withθ=(ϑ0, β, ϑ1) andθ^(t)= ϑ^(t)₀ , β^(t), ϑ^(t)₁

,t∈N0, the log-likelihood function for the complete data (X^∗,Σ^∗)=(x^∗,σ^∗) is given by

`(θ|x^∗,σ^∗)=log(c^∗)+d₁log(β)−d₁log(ϑ0)− 1 ϑ0

d1

X

i=1

(x^∗_i)^β+(β−1)

d1

X

i=1

log(x^∗_i)

+ Xn j=d1+1

11{0}(σ^∗_j)

− (x^∗_j)^β

ϑ0 +log(β)−log(ϑ0)+(β−1) log(x^∗_j)

+ Xn

j=d₁+1

11{1}(σ^∗_j)

− 1 ϑ1

x^∗_j+τ^β₁ϑ1

ϑ0

−τ1

−log(ϑ1) .

Hence, we get

Q θ;θ^(t)

= log(c^∗)+d₁log(β)−d₁log(ϑ0)− 1 ϑ0

d1

X

i=1

(x^∗_i)^β+(β−1)

d1

X

i=1

log(x^∗_i)

+ Xn j=d1+1

P^(t)_β,0,j

− (x^∗_j)^β

ϑ0 +log(β)−log(ϑ0)+(β−1) log(x^∗_j)

+ Xn

j=d₁+1

P^(t)_β,1,j

− 1 ϑ1

x^∗_j+τ^β₁ ϑ1

ϑ0

−τ1

−log(ϑ1) ,

where the weights

P^(t)_β,s,j= f^Σ

∗ j,X^∗ θ^(t) (s,x^∗)

f_θ^X_(t)^∗(x^∗) , s∈ {0,1}, j∈ {d₁+1, . . . ,n}, do not depend onθ. Further, with obvious changes, the functions f^Σ

∗ j,X^∗

θ^(t) and f^X∗

θ^(t)are as defined in (3.7) and (3.8), respectively. Moreover, we have to maximizeQ

θ;θ^(t) w.r.t.

θ. First, we have

Q θ;θ^(t)

=Qe₀(ϑ0, β)+Qe₁(ϑ1), where

Qe0(ϑ0, β)= log(c1)+(n−r^?₁) log(β) +(β−1)

d1

X

i=1

log(x^∗_i)+(β−1) Xn j=d1+1

P^(t)_β,0,j log(x^∗_j)

−(n−r^?₁) log(ϑ0)− 1 ϑ0

d₁

X

i=1

(x^∗_i)^β− 1 ϑ0

Xn j=d1+1

P^(t)_β,0,j(x^∗_j)^β−r^?₁τ^β₁ ϑ0 ,

Qe₁(ϑ1)= −r^?₁ log(ϑ1)− 1 ϑ1

Xn j=d1+1

P^(t)_β,1,jx^∗_j+ r^?₁τ1

ϑ1 .

The updated estimates ofϑ0andϑ1in the (t+1)th iteration step are given by ϑ^(t₀⁺¹⁾

β^(t+1)

= 1 n−r^?₁

d₁

X

i=1

(x^∗_i)^β(t+1) + Xn j=d1+1

P^(t)_β,0,j(x^∗_j)^β(t+1) +r^?₁τ^β₁^(t+1)

| {z }

=A_β(t+1),say,

!

, (3.9)

and

ϑ^(t+1)₁ = 1 r^?₁

Xn j=d1+1

P^(t)_β,1,jx^∗_j−r^?₁τ1

! ,

respectively. The partial derivative ofQ θ;θ^(t)

w.r.t.βis given by

∂Q

θ;θ^(t)

∂β = d₁

β − 1

ϑ0 d1

X

i=1

(x^∗_i)^βlog(x^∗_i)

| {z }

=Bβ,say,

+

d1

X

i=1

log(x^∗_i)− 1 ϑ0

r^?₁τ^β₁log(τ1)

− 1 ϑ0

Xn

j=d₁+1

P^(t)_β,0,jh

(x^∗_j)^βlog(x^∗_j)i

| {z }

=Cβ,say,

+ Xn

j=d₁+1

P^(t)_β,0,j 1

β +log(x^∗_j) =0.

Equation (3.9) plugged into the above equation yields 1

β^(t+1) −B_β(t+1) +C_β(t+1) +r^?₁τ^β₁^(t+1)log(τ1)

A_β(t+1) + 1

n−r^?₁

d1

X

i=1

log(x^∗_i)+ Xn j=d1+1

P^(t)_β,0,j log(x^∗_j)

!

=0.

(3.10) Therefore, we get the update variableβ^(t+1)ofβfor the (t+1)th iteration step by solving equation (3.10) numerically forβ^(t+1). As shown in Section A.1, equation (3.10) has a unique solution forβ^(t+1)>0. Possible initial values for the EM-algorithm are given by

ϑ⁽⁰⁾₀ β⁽⁰⁾

= 1 d₁

d1

X

i=1

(x^∗_i)^β(0) +(n−d1)τ^β₁⁽⁰⁾

!

and ϑ⁽⁰⁾₁ = 1 n−d₁

Xn j=d1+1

x^∗_j−(n−d1)τ1

! ,

whend₁ >0, whereβ⁽⁰⁾is the (unique) numerical solution of 1

β⁽⁰⁾ − (n−d₁)τ^β₁⁽⁰⁾log(τ1)+Pd1

i=1(x^∗_i)^β(0)log(x^∗_i) (n−d₁)τ^β₁⁽⁰⁾ +Pd1

i=1(x^∗_i)^β(0)

+ 1 d₁

d₁

X

i=1

log(x^∗_i)=0.

Ford1 =0, one may choose the initial values β⁽⁰⁾ =1 and ϑ⁽⁰⁾₀ =ϑ⁽⁰⁾₁ = 1

n Xn

j=1

x^∗_j−nτ1

! .

3.3 EM-Algorithm for Three-Step SLT

Fork>2, the situation is more involved. In order to illustrate the additional difficulties caused by more stages, we discuss the case of two fixed stage-change times, i.e., the case k = 3. An extension to more stages may be developed in the same manner. In particular, the EM-algorithm proposed in Section 3.2.1 is extended to a second stage- change time. As above, we assume that the design of the life test is still known, that is, n,τ1,τ2,P^∗=(π1, π2), andR∗⁰=(R⁰₁,R⁰₂) are given. The sample is given by the ordered failure timesx^∗₁ ≤ · · · ≤ x^∗_n only. Since the estimators forϑ1 andϑ2 do not exist if no failures have been observed, our calculations are conditional onR^?₁ > 0 and R^?₂ > 0.

ForR^?₁ =0 and/orR^?₂ =0, the following EM-algorithm can be used with the necessary adjustments.

Using the samplex^∗₁, . . . ,x^∗_nand the design of the life test, we have

d₁= Xn

i=1

11₍−∞,τ1](x^∗_i) and r^?₁ =%1(d₁)=











bπ1·(n−d₁)c, Type-P, minn

n−d₁,R⁰₁o

, Type-M.

Further, we introduce the following counters:

Number of observations in the interval (τ1, τ2]:b₁=Pn

i=111_(τ1,τ2](x^∗_i);

Number of observations in the interval (τ2,∞):b₂ =Pn

i=111_(τ2,∞)(x^∗_i);

Number of observations on stages₁in the interval (τ2,∞):r₁ =r^?₁ +d₂−b₁. Note that we can generally not determined2from the available information. Therefore, we have to consider all possible valuesd₂∈ {max{b₁−r^?₁,0}, . . . ,b₁}in our calculations so that

r^?₂ =%2(d₂)=











bπ2·(n−d1−d2−r^?₁)c, Type-P minn

n−d₁−d₂−r^?₁,R⁰₂o

, Type-M











>0.

The situation is illustrated in Figure 2. Note that the failure timesx^∗_d

1+1, . . . ,x^∗_d

1+b₁ could have been observed on stagess0 ors1, whereas the failure timesx^∗_d

1+b₁+1, . . . ,x^∗_n could have been observed on stagess0,s₁, ors2. Thus, in the intervals (τ1, τ2] and (τ2,∞), we have two and three options for each observation to allocate the observed data, respectively.