An Expectation Maximization Algorithm to Model Failure Times by Continuous-Time Markov Chains

Abstract

In many applications, the failure rate function may present a bathtub shape curve. In this paper, an expectation maximization algorithm is proposed to construct a suitable continuous-time Markov chain which models the failure time data by the first time reaching the absorbing state. Assume that a system is described by methods of supplementary variables, the device of stage, and so on. Given a data set, the maximum likelihood estimators of the initial distribution and the infinitesimal transition rates of the Markov chain can be obtained by our novel algorithm. Suppose that there are 𝑚 transient states in the system and that there are 𝑛 failure time data. The devised algorithm only needs to compute the exponential of 𝑚×𝑚 upper triangular matrices for 𝑂(𝑛𝑚²) times in each iteration. Finally, the algorithm is applied to two real data sets, which indicates the practicality and efficiency of our algorithm.

1. Introduction

Among many quantitative analysis methods of system reliability, the state space representation has often been employed by reliability engineers, for example, [1–5]. The state space representation has been recognized by industrial IEC61511 standards [6]. This method assumes that the structure of the system is modelled with a continuous-time Markov chain (CTMC). In the state space representation, nodes represent states of the system and arcs represent the transitions between nodes. This method is well adapted to study the reliability of various systems and allows an exact analysis of their probability of failure.

In the state space representation method, the time interval between two consecutive transitions is a random variable of an exponential distribution. However, many studies have presented counterexamples. For some mechanical and electronic components, the failure rate function of time to failure (TTF) has a bathtub curve: a monotonically decreasing function initially, eventually becoming a constant, and finally changing to a monotonically increasing function after sufficient time elapses [7]. Because existing approaches do not provide a perfect solution to deal with this class of systems, several researchers have explored models with nonexponential distributions (see [8–10] and references therein). Among many approaches treating TTF with nonexponential distributions, the extended Markov-model [11] is recommendable. In the extended Markov-model, an operation state is divided into substates with different levels of failure rates, which result in a nonconstant failure rate of the operation state. That is, in these models, the bug of failure rate functions is fixed within the framework of the state space representation (see [12, 13] and references therein for some applications).

Methods of supplementary variables [14] and the device of stages [15] are two classical approaches of extended Markov-models. These two approaches allow one to model a wide range of experimental probability density functions through a set of additional stages in series or in parallel. For supplementary variables with given parameters, a formula for evaluating the occurrence of failures is derived in [16]. In fact, with suitable parameters, additional stages can approximate the behavior of many types of distributions in practice. Reference [17] is a recent application of these approaches.

Assume that the structure of a system is established. That is, the nodes and arcs of the system are given. For a experimental data set on the TTFs of the system, the bottleneck of the above two approaches is the estimation of parameters for the nodes and arcs. Since the system is modelled with a CTMC, the TTFs are phase-type distributed random variables. The problem is similar to phase-type fitting methods [18]. However, traditional phase-type fitting methods give the maximum likelihood estimators of the minimal representation [18], which may not correspond to the system with the expected structure. Moreover, from the standpoint of computational complexity, traditional methods are not suitable for the system with many paths.

Inspired by the theory of continuous-time Markov chains observed in continuous-time channels (see [19, 20]), we propose a new expectation maximization (EM) algorithm for data of TTFs in this paper. In Section 2 we specify the framework and hidden variables for the system. In Section 3 we derive several conditional probabilities and conditional expectations, which will be used in the EM algorithm. In Section 4 we present the new EM algorithm. In Section 5, we provide numerical results. Section 6 concludes the paper.

2. System Definition and Hidden Variables

Assume that a system is modelled by a CTMC 𝑋∶={𝑋_𝑡∶𝑡≥0} with transient states 1,2,…,𝑚, and one absorbing state 𝑚+1. The infinitesimal transition rates 𝑄=(𝑞_𝑖𝑗) of 𝑋 satisfy 𝑞_{𝑚+1,𝑚+1}=0,𝑞_𝑖𝑖<0,for𝑖≠𝑚. (2.1) The quantities 𝑞_𝑖𝑗,𝑖≠𝑚+1, are unknown. However, the characteristic matrix 𝐶=(𝑐_𝑖𝑗), 𝑐_𝑖𝑗=1,𝑞_𝑖𝑗≠0,𝑖≠𝑗,0,otherwise, (2.2) which represents the structure of the system, is the known. Assume that the characteristic 𝐶 is an upper triangular matrix, that is, 𝑐_𝑖𝑗=0 for 𝑖<𝑗. We suppose that the initial probability function 𝜋=(𝜋₁,𝜋₂,…,𝜋_𝑚+1) satisfies 𝜋_𝑚+1=0 and 𝜋₁,𝜋₂,…,𝜋_𝑚 are unknown.

Example 2.1 (The device of stages). The system has a state transition diagram shown in Figure 1. The characteristic matrix 𝐶=(𝑐_𝑖𝑗) of this system is 𝑐_𝑖𝑗=1,for𝑗=𝑖+1,𝑐_𝑖𝑗=0,for𝑗≠𝑖+1. (2.3) Hence 𝐶 is an upper triangular matrix.

Figure 1 

The state transition diagram of a device of stages.

Example 2.2 (The parallel system). The system has a state transition diagram shown in Figure 2. The characteristic matrix 𝐶=(𝑐_𝑖𝑗) of this system is 𝑐_𝑖,𝑗=1,for𝑗=𝑚+1,𝑐_𝑖𝑗=0,for𝑗≠𝑚+1. (2.4) Hence 𝐶 is also an upper triangular matrix.

Figure 2 

The state transition diagram of a parallel system.

Like the above two examples, the structures of many practical systems can be represented by upper triangular characteristic matrixes.

We define a path, the sequence of nodes visited from an initial node to the absorbing node 𝑚+1. The 𝑘th path, 𝑟_𝑘=(𝑛₁,𝑛₂,…,𝑛_𝑠), is characterized by a set of nodes 𝑛₁,𝑛₂,…,𝑛_𝑠, representing the nodes visited before absorption, where 𝑐_{𝑛𝑖𝑛𝑖+1}=1 for 𝑖=1,2,…,𝑠−1 and 𝑐_{𝑛𝑠,𝑚+1}=1. Since the characteristic matrix 𝐶 and hence the infinitesimal transition rates 𝑄, are upper triangular matrixes, any node will never appear again in a path. So the length of a path is no more than 𝑚+1, and hence the number of paths is finite.

According to the above discussion, similarly to part b) of Theorem 5.4 in [21], we can provide the following noncanonical construction of a Markov chain with the transition rates 𝑄 and the initial distribution 𝜋, respectively. The hidden variables, which comprise the noncanonical construction, will serve the EM algorithm in Section 4.

Lemma 2.3. For 𝑘=1,2,…,𝑛, let 𝑇^𝑘₁,𝑇^𝑘₂,…,𝑇^𝑘_𝑚 be a sequence of exponentially distributed random variables with parameters −𝑞₁₁,−𝑞₂₂,…,−𝑞_𝑚𝑚, respectively. For 𝑘=1,2,…,𝑛, let 𝐽^𝑘₀,𝐽^𝑘₁,𝐽^𝑘₂,…,𝐽^𝑘_𝑚 be a sequence of random variables with discrete distributions described by probability functions 𝑞_𝑖(𝑗)=𝑃(𝐽^𝑘_𝑖𝑞=𝑗)=0,𝑖=𝑗,𝑖≠0,_𝑖𝑗−𝑞_𝑖𝑖𝜋,𝑖≠𝑗,𝑖≠0,_𝑗,𝑖=0. (2.5) Assume also that the sequences {𝐽^𝑘₀,𝐽^𝑘_𝑗,𝑇^𝑘_𝑗,𝑗=1,2,…,𝑚,𝑘=1,2,…,𝑛} comprise an independent family. Define 𝑌^𝑘₀=𝐽^𝑘₀,𝑌^𝑘_𝑠=𝐽^𝑘_{𝑌𝑘𝑠−1},for𝑠≥1,𝑌^𝑘_𝑠−1𝑋≠𝑚+1,_𝑘𝑡=𝑌^𝑘₀,0≤𝑡<𝑇^𝑘_𝑌𝑘0,𝑋^𝑘_𝑡=𝑌^𝑘_𝑠+1,for𝑇^𝑘_𝑌𝑘0+𝑇^𝑘_𝑌𝑘1+⋯+𝑇^𝑘_𝑌𝑘𝑠≤𝑡<𝑇^𝑘_𝑌𝑘0+𝑇^𝑘_𝑌𝑘1+⋯+𝑇^𝑘_{𝑌𝑘𝑠+1}. (2.6) Then the constructed process 𝑋^𝑘∶={𝑋^𝑘_𝑡∶𝑡≥0} is Markovian with the initial distribution 𝜋 and the infinitesimal transition rates 𝑄. Moreover, 𝑋^𝑘,𝑘=1,2,…,𝑛 are independent.
For a sample of 𝐽∶=(𝐽^k₀,𝐽^𝑘₁,𝐽^𝑘₂,…,𝐽^𝑘_𝑚), Lemma 2.3 indicates that there is a unique path corresponding to the discrete-time embedded Markov chain (𝑌^𝑘₀,𝑌^𝑘₁,…). We denote this path as ℛ(𝐽). Write 𝑗∈ℛ(𝐽) if the node 𝑗 appears in the path ℛ(𝐽). Given a path ℛ(𝐽), define the first node ℱ_𝐽(0)=𝐽^𝑘₀ for this path, and the following node ℱ_𝐽(𝑗)=𝐽^𝑘_𝑗 for any node 𝑗∈ℛ(𝐽).
A sample (𝜏¹,𝜏²,…,𝜏^𝑛) of the TTFs for a system is a family of independent identically distributed random variables. Precisely, the stochastic processes 𝑋^𝑘,𝑘=1,2,…,𝑛, are independent, which have the same infinitesimal transition rates 𝑄 and the same initial probability function 𝜋. Then define𝜏^𝑘=inf{𝑡≥0∶𝑋^𝑘_𝑡=𝑚+1}. (2.7) Using the above notations, we have 𝜏^𝑘=_{𝑖∈ℛ(𝐽)}T^𝑘_𝑖,𝐽=(𝐽^𝑘₀,𝐽^𝑘₁,𝐽^𝑘₂,…,𝐽^𝑘_𝑚). (2.8) Based on the observations 𝜏¹,𝜏²,…,𝜏^𝑛 and the known characteristic matrix 𝐶, our problem is to give estimators 𝜋 and 𝑄 for the unknown parameters.

3. Conditional Probabilities and Conditional Expectations

To establish an EM algorithm for our problem, some conditional probabilities and conditional expectations will be calculated in this section. From now on, we will omit the superscript 𝑘 of the stochastic processes when the clarity is being held. Firstly, we recall the cumulative distribution function of the TTF 𝜏, 𝐹_𝜏(𝑡)=𝑃(𝜏≤𝑡)=𝑝_𝑚+1(𝑡). (3.1) Here 𝑝(𝑡)=(𝑝₁(𝑡),𝑝₂(𝑡),…,𝑝_𝑚+1(𝑡)) is the probability function of 𝑋_𝑡, which satisfies the Kolmogorov's forward equation 𝑝^′(𝑡)=𝑝(𝑡)𝑄. From the solution to the equation in exponential form, we have 𝑝_𝑚+1(𝑡)=𝜋exp(𝑡𝑄)1^𝑇_𝐷, (3.2) where the vector 1^𝑇_𝐷 of dimensions 𝑚+1 is the transpose of 1_𝐷=(0,0,…,0,1). So the probability density function of the TTF 𝜏 is𝑓_𝜏(𝑡)=𝜋exp(𝑡𝑄)𝑄1^𝑇_𝐷. (3.3) The above methodology can be used to calculate other probabilities and expectations in this section.

Lemma 3.1. For the above noncanonical construction of 𝑋 in Lemma 2.3, we have for 1≤𝑖≤𝑚, 𝑃(𝑖∈ℛ(𝐽₀,𝐽₁,…,𝐽_𝑚1)∣𝜏=𝑡)=𝑓_𝜏𝜋(𝑡)^(𝑖)exp(𝑡𝑄_𝑖)𝑄_𝑖1^𝑇_𝑃,for𝑡≥0. (3.4) Here 𝑄_𝑖=(𝑞^′_𝑟𝑠) is a matrix of dimensions (𝑚+2)×(𝑚+2) determined as follows: 𝑞^′_𝑟𝑠=𝑞_𝑟𝑠,𝑟<𝑖,𝑠≤𝑖,_{𝑚+1𝑡=𝑖+1}𝑞_𝑟𝑡𝑞,𝑟<𝑖,𝑠=𝑚+2,_𝑟𝑠,𝑖≤𝑟≤𝑚+1,𝑠≤𝑚+1,0,otherwise. (3.5)1_𝑃=(0,0,…,0,1,0) and 𝜋^(𝑖)=(𝜋₁^(𝑖),𝜋₂^(𝑖),…,𝜋^(𝑖)_𝑚+2) are vectors of dimensions 𝑚+2, 𝜋_𝑠^(𝑖)=𝜋_𝑠,𝑠≤𝑖,0,𝑖<𝑠<𝑚+1,_{𝑚+1𝑡=𝑖+1}𝜋_𝑡,𝑠=𝑚+2. (3.6)

Proof. Define the events 𝐴_𝑡={𝑖∈ℛ(𝐽₀,𝐽₁,…,𝐽_𝑚),𝜏≤𝑡},𝑡≥0. (3.7) Considering all possible realizations of 𝐽=(𝐽₀,𝐽₁,…,𝐽_𝑚), it follows from (2.8) that 𝑃(𝐴_𝑡)=_{ℛ(𝐽)∋𝑖}𝑃(𝐽₀=𝑗₀,𝐽_𝑠=ℱ(𝐽,𝑠),𝑠∈ℛ(𝐽),_{𝑠∈ℛ(𝐽)}𝑇_s≤𝑡). (3.8) Since the characteristic matrix 𝐶 is an upper triangular matrix, the numbers of nodes along a path are monotonic increasing. Then if 𝐽=(𝑗₀,𝑗₁,…,𝑗_𝑚) satisfies ℛ(𝐽)∋𝑖, we have 𝑗₀≤𝑖 and ℱ(𝐽,𝑠)≤𝑖,for𝑠∈ℛ(𝐽),𝑠<𝑖. (3.9)
Considering another stochastic process 𝑋^′ with the infinitesimal transition rates 𝑄_𝑖 and the initial probability function 𝜋^(𝑖), we can deduce a noncanonical construction consisting of 𝐽^′₀,𝐽^′₁,…,𝐽^′_𝑚,𝑇^′₁,𝑇^′₂,…,𝑇^′_𝑚. It follows from (3.9) that𝑃(𝐴_t)=_{ℛ(𝐽′)∋𝑖}𝑃(𝐽^′₀=𝑗₀,𝐽^′_s=ℱ(𝐽^′,𝑠),𝑠∈ℛ(𝐽^′),_{𝑠∈ℛ(𝐽′)}𝑇_𝑠≤𝑡). (3.10) For a path 𝐽^′ of the stochastic process 𝑋^′, the property ℛ(𝐽^′)∋𝑖 means that the end of the path 𝐽^′ is the absorbing state 𝑚+1. Similarly to (3.2), we have 𝑃(𝐴_𝑡)=𝑃(𝑋^′_𝑡=𝑚)=𝜋^(𝑖)exp(𝑡𝑄_𝑖)1^𝑇_𝑃. (3.11) The result follows from the above equation and 𝑃(𝐴_𝑡∣𝜏=𝑡)=lim_𝜖↓0𝑃(𝐴_𝑡+𝜖)−𝑃(𝐴_𝑡).𝑃(𝑡<𝜏≤𝑡+𝜖) (3.12)

Similarly to the discussion in Lemma 3.1, we have the following Lemma.

Lemma 3.2. For the noncanonical construction of 𝑋 in Lemma 2.3, 1≤𝑖≤𝑚, given 𝑖∈ℛ(𝐽₀,𝐽₁,…,𝐽_𝑚), the conditional density function of 𝜏 is 𝑓_𝜏1(𝑡∣𝑖∈ℛ(𝐽))=𝑝^′_𝑖𝜋^(𝑖)exp(𝑡𝑄_𝑖)𝑄_𝑖1^𝑇_𝑃,for𝑡≥0, (3.13) where 𝐽=(𝐽₀,𝐽₁,…,𝐽_𝑚) and 𝑝^′_𝑖=𝑃(𝑖∈ℛ(𝐽)).

𝑝^′_𝑖 in Lemma 3.2 can be expressed through the 1-step transition probability matrix of the embedded Markov chain of the continuous-time Markov chain with the infinitesimal transition rates 𝑄_𝑖. However, since 𝑝^′_𝑖 will be eliminated in the further propositions, it is not necessary to present the expression.

Theorem 3.3. For the noncanonical construction of 𝑋 in Lemma 2.3, 1≤𝑖≤𝑚, 1≤𝑗≤𝑚+1, and 𝑡≥0, we have 𝑃(𝑗=ℱ(𝐽,𝑖),𝑖∈ℛ(𝐽₀,𝐽₁,…,𝐽_𝑚)∣𝜏=𝑡)=0,𝑐_𝑖,𝑗1=0,𝑓_𝜏𝜋(𝑡)^(𝑖)exp(𝑡𝑄_𝑖𝑗)𝑄_𝑖𝑗1^𝑇_𝑃,𝑐_𝑖,𝑗≠0. (3.14) Here 1_𝑃 and 𝜋^(𝑖) are the same as those in Lemma 3.1. 𝑄_𝑖𝑗=(𝑞_𝑟𝑠) is a matrix of dimensions (𝑚+2)×(𝑚+2) determined as follows: 𝑞_𝑟𝑠=𝑞_𝑟𝑠𝑞,𝑟≠𝑖,_𝑖𝑠,𝑟=𝑖,𝑠=𝑖or𝑗,_{𝑚+1𝑡=𝑖+1}𝑞_𝑖𝑡−𝑞_𝑖𝑗,𝑟=𝑖,𝑠=𝑚+2,0,otherwise. (3.15)

Proof. If 𝑐_𝑖,𝑗=0, it is obvious that 𝑃(𝑖∈ℛ(𝐽),𝑗=ℱ(𝐽,𝑖),𝜏≤𝑡)=0. Otherwise, by constructing an auxiliary stochastic process with the infinitesimal transition rates 𝑄_𝑖𝑗 and the initial probability function 𝜋^(𝑖), we can obtain, similarly to the proof of Lemma 3.1, that 𝑃(𝑖∈ℛ(𝐽),𝑗=ℱ(𝐽,𝑖),𝜏≤𝑡)=𝜋^(𝑖)exp(𝑡𝑄_𝑖𝑗)1^𝑇_𝑃, (3.16) where 𝐽=(𝐽₀,𝐽₁,…,𝐽_𝑚). The required result can then be derived in the similar way as that in Lemma 3.1.

Theorem 3.4. For the noncanonical construction of 𝑋 in Lemma 2.3 and 𝑖≤𝑚,𝑡≥0, we have 𝐸[𝑇_𝑖𝑞∣𝑖∈ℛ(𝐽),𝜏=𝑡]=−^𝑖𝑖exp(𝑞_𝑖𝑖𝑡)𝜋^(𝑖)exp(𝑡𝑄_𝑖)𝑄_𝑖1^𝑇_𝑃×(𝑡𝜋^(𝑖,𝑇)_𝑚+1+𝜋^{(𝑖,𝑇)𝑡}₀(𝑡−𝑢)exp(𝑢(𝑀_𝑖−𝑞_𝑖𝑖𝐼))𝑑𝑢𝑀_𝑖1^𝑇_𝑃), (3.17) where 𝐽=(𝐽₀,𝐽₁,…,𝐽_𝑚). 𝐼 is an identity matrix and 𝑀_𝑖=(𝑀_𝑟𝑠) is a matrix of dimensions (𝑚+2)×(𝑚+2) determined as follows: 𝑀_𝑟𝑠=𝑞_𝑟𝑠−𝑞,𝑟<𝑖,𝑠<𝑖,_𝑟𝑖𝑞_𝑖,𝑠𝑞_𝑖𝑖,𝑟<𝑖,𝑖<𝑠≤𝑚+1,_{𝑚+1𝑡=𝑖+1}𝑞_𝑟𝑡𝑞,𝑟<𝑖,𝑠=𝑚+2,_𝑟𝑠,𝑖<𝑟≤𝑚+1,𝑖<𝑠≤𝑚+1,0,otherwise. (3.18)𝜋^(𝑖,𝑇)=(𝜋₁^(𝑖,𝑇),𝜋₂^(𝑖,𝑇),…,𝜋^(𝑖,𝑇)_𝑚+2) is vectors of dimensions 𝑚+2 determined as follows: 𝜋_𝑠^(𝑖,𝑇)=𝜋_𝑠−𝜋,𝑠<𝑖,0,𝑠=𝑖,_𝑖𝑞_𝑖𝑠𝑞_𝑖𝑖,𝑖<𝑠≤𝑚+1,_{𝑚+1𝑡=𝑖+1}𝜋_𝑡,𝑠=𝑚+2. (3.19)

Proof. It follows from (2.8) that given 𝑖∈ℛ(𝐽), 𝜏=𝑇_𝑖+_{𝑗∈ℛ(𝐽),𝑗≠𝑖}𝑇_𝑗, (3.20) where 𝑇_𝑖 and 𝑇^𝑖=∑_{𝑗∈ℛ(𝐽),𝑗≠𝑖}𝑇_𝑗, are independent random variables. Given 𝑖∈ℛ(𝐽), the conditional density function of 𝑇_𝑖 is 𝑓_𝑖(𝑡)=−𝑞_𝑖𝑖exp(𝑞_𝑖𝑖𝑡),for𝑡≥0. (3.21) Given 𝑖∈ℛ(𝐽), suppose that the conditional density function of 𝑇^𝑖 is 𝑔_𝑖(𝑡), where 𝑔_𝑖(𝑡)=0 for 𝑡<0. Since the Jacobian of the transformation (𝑥,𝑦)→(𝑥,𝑥+𝑦) is 1, the joint density of 𝑇_𝑖,𝜏 is given by 𝑓(𝑡₁,𝑡₂)=𝑓_𝑖(𝑡₁)𝑔_𝑖(𝑡₂−𝑡₁). Then for 𝑡≥0, we have 𝐸[𝑇_𝑖1∣𝑖∈ℛ(𝐽),𝜏=𝑡]=𝑓_𝜏(𝑡∣𝑖∈ℛ(𝐽))^𝑡₀𝑢𝑓_𝑖(𝑢)𝑔_𝑖=1(𝑡−𝑢)𝑑𝑢𝑓_𝜏(𝑡∣𝑖∈ℛ(𝐽))^𝑡₀(𝑡−𝑢)𝑓_𝑖(𝑡−𝑢)𝑔_𝑖(𝑢)𝑑𝑢. (3.22) Here the conditional density function 𝑓_𝜏(𝑡∣𝑖∈ℛ(𝐽)), given 𝑖∈ℛ(𝐽), is determined through (3.13) in Lemma 3.2.
As to the cumulative distribution function of 𝑇^𝑖, given 𝑖∈ℛ(𝐽), we consider another stochastic process 𝑋^′ with the infinitesimal transition rates 𝑀_𝑖 and the initial probability function 𝜋^(𝑖,𝑇). It follows from Lemma 2.3 that 𝑋^′ has a noncanonical construction consisting of 𝐽^′₀,𝐽^′₁,…,𝐽^′_𝑚,𝑇^′₁,𝑇^′₂,…,𝑇^′_𝑚. Similarly to the discussion in Lemma 3.1, we have 𝑃(𝑇^𝑖≤𝑡,𝑖∈ℛ(𝐽))=𝑃(_{𝑗∈ℛ(𝐽′),𝑗≠𝑖}𝑇_𝑗≤𝑡,theendofthepath𝐽^′istheabsorbingstate𝑚+1), (3.23) where 𝐽^′=(𝐽^′₀,𝐽^′₁,…,𝐽^′_𝑚). Hence 𝑃(𝑇^𝑖≤𝑡,𝑖∈ℛ(𝐽))=𝜋^(𝑖,𝑇)exp(𝑡𝑀_𝑖)1^𝑇_𝑃. (3.24) Then it follows from (3.13) and (3.24) that 𝑃(𝑇^𝑖1≤𝑡∣𝑖∈ℛ(𝐽))=𝑝^′_𝑖𝜋^(𝑖,𝑇)exp(𝑡𝑀_𝑖)1^𝑇_𝑃. (3.25) Since the probability in (3.25) may jump at 𝑡=0, the conditional density function corresponding to (3.24) is 𝑔_𝑖1(𝑡)=𝑝^′_𝑖𝜋^(𝑖,𝑇)_𝑚+11𝛿(𝑡)+𝑝^′_𝑖𝜋^(𝑖,𝑇)exp(𝑡𝑀_𝑖)𝑀_𝑖1^𝑇_𝑃,for𝑡≥0, (3.26) where 𝛿(𝑡) is the Dirac function. 𝑝^′_𝑖=𝑃(𝑖∈ℛ(𝐽)) is the same as that in Lemma 3.2. Then the integral in (3.22) becomes ^𝑡₀(𝑡−𝑢)𝑓_𝑖(𝑡−𝑢)𝑔_𝑖𝑞(𝑢)𝑑𝑢=−^𝑖𝑖𝑝^′_𝑖𝑡exp(𝑞_𝑖𝑖𝑡)𝜋^(𝑖,𝑇)_𝑚+1−𝑞_𝑖𝑖𝑝^′_𝑖exp(𝑞_𝑖𝑖𝑡)𝜋^{(𝑖,𝑇)𝑡}₀(𝑡−𝑢)exp(𝑢(𝑀_𝑖−𝑞_𝑖𝑖𝐼))𝑑𝑢𝑀_𝑖1^𝑇_𝑃. (3.27) The desired result follows from the above equation and Lemma 3.2.

Remark 3.5. The integral in (3.17) can be efficiently computed by using an approach developed by Van Loan [22]. To apply that approach here, we first compute the exponential of the following 3(𝑚+2)×3(𝑚+2) block triangular matrix: exp[𝑡𝑂𝐼𝑂𝑂𝑂𝐼𝑂𝑂(𝑀_𝑖−𝑞_𝑖𝑖)𝐼]=𝐼𝑡𝐼𝐻(𝑡)𝑂𝐼𝐺(𝑡)𝑂𝑂𝐹(𝑡). (3.28) Then the upper right block 𝐻(𝑡) of this matrix gives the integral ∫^𝑡₀(𝑡−𝑢)exp(𝑢(𝑀_𝑖−𝑞_𝑖𝑖𝐼))𝑑𝑢. Moreover, the matrix exponentials can be efficiently performed by using the diagonal Padé approximation with repeated squaring, as recommended in [22, 23].

Remark 3.6. Sometimes, the denominator of (3.17) is 0. For example, under some characteristic matrices, a node 𝑖 may never be reached. Since the estimation of 𝑇_𝑖 is nonsense under this case, we can assign any value to 𝐸[𝑇_𝑖∣𝑖∈ℛ(𝐽),𝜏=𝑡]. We suggest −1/𝑞_𝑖𝑖.

The following corollary can be deduced from the above theorems and the tower property of conditional expectations.

Corollary 3.7. For the noncanonical construction of 𝑋 in Lemma 2.3 and 1≤𝑖≤𝑚,1≤𝑗≤𝑚+1,𝑡≥0, we have 𝐸[𝑇_𝑖∣𝜏=𝑡]=𝑝_𝑖𝐸[𝑇_𝑖∣𝑖∈ℛ(𝐽),𝜏=𝑡]−1−𝑝_𝑖𝑞_𝑖𝑖,𝑃(𝐽_𝑖=𝑗∣𝜏=𝑡)=𝑃(𝑗=ℱ(𝐽,𝑖),𝑖∈ℛ(𝐽)∣𝜏=𝑡)−(1−𝑝_𝑖)𝑞_𝑖𝑗𝑞_𝑖𝑖,𝑃(𝐽₀1=𝑗∣𝜏=𝑡)=𝑓_𝜏𝜋(𝑡)^𝑖𝑒_𝑖exp(𝑡𝑄)𝑄1^𝑇_𝐷. (3.29) Here 𝑒_𝑖 is the unit basis vectors (0,0,…,1,0,…,0,0) with “1” in the 𝑖th position. 𝐽=(𝐽₀,𝐽₁,…,𝐽_𝑚) and 𝑝_𝑖=𝑃(𝑖∈ℛ(𝐽₀,𝐽₁,…,𝐽_𝑚)∣𝜏=𝑡).

Proof. The first equation follows from Theorem 3.4 and the tower property of conditional expectations. The second equation follows from 𝑃(𝐽_𝑖=𝑗∣𝜏=𝑡)=𝑃(𝑗=ℱ(𝐽,𝑖)∣𝑖∈ℛ(𝐽)𝜏=𝑡)𝑝_𝑖+𝑃(𝐽_𝑖=𝑗)(1−𝑝_𝑖). (3.30) As to the third equation, we have 𝑃(𝐽₀=𝑗,𝜏≤𝑡)=𝑃(𝐽₀=𝑗)𝑃(𝜏≤𝑡∣𝐽₀=𝑗). (3.31) From the solution to the Kolmogorov's forward equation in the exponential form, we obtain 𝑃(𝜏≤𝑡∣𝐽₀=𝑗)=𝜋_𝑖𝑒_𝑖exp(𝑡𝑄)𝑄1^𝑇_𝐷. (3.32) Similarly to (3.12), we can see that the third equation follows from the above two equations and (3.2).

4. EM Algorithm

The EM algorithm is a two-step iterative process for computing the maximum likelihood estimate. This process is terminated when some imposed measure of convergence for the sequence of estimators is satisfied. A filter-based form of the EM algorithm for a Markov chain was presented in [24] and was developed in [25]. Here we review firstly the filter-based EM framework.

Let 𝜃 index a given family of probability density functions 𝑓_𝜃 for some hidden random variable 𝑋_𝑖, where 𝜃 is an element of a parameter space Θ. Then based on a hidden complete data 𝑋₁,…,𝑋_𝑛, a likelihood ℒ can be written as ℒ(𝜃)=^𝑛_𝑘=1𝑓_𝜃(𝑋_𝑖). (4.1) Let a 𝜎- field ℱ=𝜎{𝑋₁,…,𝑋_𝑛}. For a random variable 𝑌∈ℱ, 𝐸^𝜃∗𝑌 denotes the expectation of 𝑌 with assumption that 𝜃^∗ is the parameter of 𝑋_𝑖s. Suppose another 𝜎- field 𝒴⊂ℱ.

Two iterative steps in the EM algorithm are as follows.

(1) Expectation step: fix 𝜃^∗=𝜃_𝑡, then compute 𝑄(⋅,𝜃^∗), where 𝑄(𝜃,𝜃^∗)=𝐸^𝜃∗[log(ℒ(𝜃))∣𝒴]. (4.2)

(2) Maximization step: Set 𝜃_𝑡+1=𝜃 as a best solution to the following optimization problem:𝜃max_∈Θ𝑄(𝜃,𝜃^∗). (4.3)

Under this framework, we propose an EM algorithm which attempts to find a maximum likelihood estimate (MLE) of the parameter 𝜃=(𝜋,𝑄). Here 𝒴 is a 𝜎- field generated from a sequence (𝜏¹,𝜏²,…,𝜏^𝑛) of TTFs. As for the hidden complete data, we choose 𝑇^𝑘₁,𝑇^𝑘₂,…,𝑇^𝑘_𝑚 and 𝐽^𝑘₀,𝐽^𝑘₁,𝐽^𝑘₂,…,𝐽^𝑘_𝑚, where 𝑘=1,2,…,𝑛. As showed in Lemma 2.3, no information is gained or lost by considering stochastic processes 𝑋^𝑘,𝑘=1,2,…,𝑛, as hidden processes instead of 𝑇^._. and 𝐽^._.. However, our choice of the hidden random variables results in a simple E-step, and it leads to an explicit M-step for estimating the parameters.

In this case, the complete data likelihood ℒ can be written asℒ=^𝑛_𝑘=1(^𝑚_𝑖=0𝑞_𝑖(𝐽^𝑘_𝑖)^𝑚_𝑖=1(−𝑞_𝑖𝑖)exp(𝑞_𝑖𝑖𝑇^𝑘_𝑖)). (4.4) The factors 𝑞_𝑖(𝐽^𝑘_𝑖) in (4.4) represent the probability functions of the hidden random variables 𝐽^𝑘_𝑖. The factors −𝑞_𝑖𝑖exp(𝑞_𝑖𝑖𝑇^𝑘_𝑖) represent the density functions of 𝑇^𝑘_𝑖. Taking the logarithm of (4.4) yieldslogℒ=^𝑛_𝑘=1(^𝑚𝑚_𝑖=0𝑗=1log(𝑞_𝑖(𝑗))𝐼(𝐽^𝑘_𝑖=𝑗)+^𝑚_𝑖=1(log(−𝑞_𝑖𝑖)+𝑞_𝑖𝑖𝑇^𝑘_𝑖)), (4.5) where 𝐼(𝐽^𝑘_𝑖=𝑗)=1 if 𝐽^𝑘_𝑖=𝑗 and otherwise 𝐼(𝐽^𝑘_𝑖=𝑗)=0.

A new parameter estimate, say 𝜃=(𝜋,𝑄), is obtained as the parameter which maximizes the expected value of (4.5) given (𝜏¹,𝜏²,…,𝜏^𝑛) and a current parameter estimate 𝜃^∗=(𝜋,𝑄). The constrained maximization gives us the following new estimate:𝑞_𝑖1(𝑗)=𝑛^𝑛_𝑘=1𝑃^𝜃∗(𝐽_𝑖=𝑗∣𝜏=𝜏^𝑘𝑞),_𝑖𝑖𝑛=−∑^𝑛_𝑘=1𝐸^𝜃∗(𝑇_𝑖∣𝜏=𝜏^𝑘), (4.6) where the conditional probabilities and the conditional expectations are given in Corollary 3.7. Then it follows from (2.5) that𝜋_𝑗=𝑞₀(𝑗),𝑞_𝑖𝑗=−𝑞_𝑖𝑖𝑞_𝑖(𝑗),for𝑖≠𝑗,𝑖≠0. (4.7) These expressions constitute the M-step of our EM algorithm.

5. Two Applications

In this section, we present two examples to illustrate the results derived in the previous sections.

The first data set is a widely used sample of 50 components taken from [26], which possess a bath-shaped failure rate property. These data consist of times to failure of 50 devices put on life test at time 0. The data is shown in Table 1.

Table 1 

Failure time data from [26].

Using our EM algorithm with 𝑚=3, which means that there are transient states 1,2,3, and one absorbing state 4, we obtain .𝜋=(0.000007,0.124512,0.875480,0.0),𝑄=−1.0070041.0070040.00.00.0−1.0640190.0112951.0527240.00.0−0.0192700.0192700.00.00.00.0 (5.1)

In many applications of the failure rate function, the total time on test (TTT) plot [26] was used to show the efficiency of the selected model. The TTT plot is obtained by plotting 𝑟𝐺(𝑛)=(∑^𝑟_𝑖=1𝑇_𝑖∶𝑛+(𝑛−𝑟)𝑇_𝑟∶𝑛)∑^𝑛_𝑖=1𝑇_𝑖∶𝑛 (5.2) against 𝑟/𝑛. Here 𝑟=1,…,𝑛 and 𝑇_𝑖∶𝑛,𝑖=1,…,𝑛, are the order statistics of the sample. To check the efficiency of the result, Monte Carlo simulations have also been carried out. We generated 500 independent TTFs, with the estimated parameters 𝜋,𝑄. Due to the sensitivity of TTT plot to 𝑇_𝑛∶𝑛, the simulated TTFs are truncated by the maximum of the data set. The TTT plots for the first data set and the truncated simulated TTFs are shown in Figure 3.

Figure 3 

The TTT plots for the first data set and the simulated TTFs. The circles represent the real data and the solid points represent the simulated results with the estimated parameters.

Another real set of lifetime data set, taken from [27], is the failure data collected during unit testing. During a unit testing phase, the number of failures in a certain time interval (of equal length) is recorded. The total number of units tested is equal to 311. The data is given in Table 2.

Table 2 

Failure time data from [27].

Using our EM algorithm with 𝑚=3, we obtain ,.𝜋=(0.998642,0.001358,0.0,0.0)𝑄=−0.1729570.1655930.0073630.00.0−1.1645601.1596260.0049340.00.0−4.7670494.7670490.00.00.00.0 (5.3) The TTT plots for this data set and the Monte Carlo simulated TTFs are given in Figure 4.

Figure 4 

The TTT plots for the second data set and the simulated TTFs. The dashed boxes represent the real data and the solid points represent the simulated results with the estimated parameters.

6. Conclusions

One great advantage with the extended Markov-models is that these models can tackle operation states of nonexponential distributions under the framework of the state space representation. Based on the theories of continuous-time Markov chains and EM algorithms, this paper deals with the problem of parameters estimation for extended Markov-models. Given a data set with 𝑛 failure time data, for a system with 𝑚 transient states, the proposed EM algorithm only needs to compute the exponential of 𝑚×𝑚 upper triangular matrices for 𝑂(𝑛𝑚²) times in each iteration. The effectiveness of the proposed EM algorithm is well demonstrated through two real world examples.

Two real applications show the practicality and efficiency of our new EM algorithm. Nevertheless, one theoretical limitation about EM algorithms is that its iterations sequence converges to a stationary point of the likelihood. Only in the case that the likelihood function is convex, the global convergence from any initial value can be ensured. However, our experience with many examples is against this theoretical drawback. In few practical cases, the stationary points are saddle points. As a simple perturbation can lead the algorithm to diverge from the saddle point, this phenomenon is thus not a problem in practice. However, to avoid that the limit point is only a local maximum of the likelihood, the initial value must be carefully selected in some cases. We suggest to use moments of related random variables to construct the initial estimates.

Of course, it will be interesting to investigate conditions ensuring the convergence to the global maximum of the likelihood. This is left for future research.

Acknowledgments

Dengfu Zhao and Qihong Duan acknowledge the funding of the National Natural Science Foundation of China (Grant number 50977073) and the Fundamental Research Funds for the Central Universities. Zhiping Chen and Qihong Duan acknowledge the funding of the National Natural Science Foundation of China (Grant number 70971109).