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

Authors: Qihong Duan; Zhiping Chen; Dengfu Zhao
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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 arem transient states in the system and that there are n failure time data. The devised algorithm only needs to compute the exponential of m×m upper triangular matrices for O(nm2) times in each iteration. Finally, the algorithm is applied to two real data sets, which indicates the practicality and efficiency of our algorithm.


---

## Body

## 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 CTMCX:={Xt:t≥0} with transient states 1,2,…,m, and one absorbing state m+1. The infinitesimal transition rates Q=(qij) of X satisfy (2.1)qm+1,m+1=0,qii<0,fori≠m.
The quantities qij,i≠m+1, are unknown. However, the characteristic matrix C=(cij), (2.2)cij={1,qij≠0,i≠j,0,otherwise,
which represents the structure of the system, is the known. Assume that the characteristic C is an upper triangular matrix, that is, cij=0 for i<j. We suppose that the initial probability function π=(π1,π2,…,πm+1) satisfies πm+1=0 and π1,π2,…,πm are unknown.Example 2.1 (The device of stages).
The system has a state transition diagram shown in Figure1. The characteristic matrix C=(cij) of this system is
(2.3)cij=1,forj=i+1,cij=0,forj≠i+1.
Hence C 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 Figure2. The characteristic matrix C=(cij) of this system is
(2.4)ci,j=1,forj=m+1,cij=0,forj≠m+1.
Hence C 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 nodem+1. The kth path, rk=(n1,n2,…,ns), is characterized by a set of nodes n1,n2,…,ns, representing the nodes visited before absorption, where cnini+1=1 for i=1,2,…,s-1 and cns,m+1=1. Since the characteristic matrix C and hence the infinitesimal transition rates Q, are upper triangular matrixes, any node will never appear again in a path. So the length of a path is no more than m+1, and hence the number of paths is finite.According to the above discussion, similarly to part b) of Theorem5.4 in [21], we can provide the following noncanonical construction of a Markov chain with the transition rates Q and the initial distribution π, respectively. The hidden variables, which comprise the noncanonical construction, will serve the EM algorithm in Section 4.Lemma 2.3.
Fork=1,2,…,n, let T1k,T2k,…,Tmk be a sequence of exponentially distributed random variables with parameters -q11,-q22,…,-qmm, respectively. For k=1,2,…,n, let J0k,J1k,J2k,…,Jmk be a sequence of random variables with discrete distributions described by probability functions
(2.5)qi(j)=P(Jik=j)={0,i=j,i≠0,qij-qii,i≠j,i≠0,πj,i=0.
Assume also that the sequences {J0k,Jjk,Tjk,j=1,2,…,m,k=1,2,…,n} comprise an independent family. Define
(2.6)Y0k=J0k,Ysk=JYs-1kk,fors≥1,Ys-1k≠m+1,Xtk=Y0k,0≤t<TY0kk,Xtk=Ys+1k,forTY0kk+TY1kk+⋯+TYskk≤t<TY0kk+TY1kk+⋯+TYs+1kk.
Then the constructed process Xk:={Xtk:t≥0} is Markovian with the initial distribution π and the infinitesimal transition rates Q. Moreover, Xk,k=1,2,…,n are independent.
For a sample ofJ:=(J0k,J1k,J2k,…,Jmk), Lemma 2.3 indicates that there is a unique path corresponding to the discrete-time embedded Markov chain (Y0k,Y1k,…). We denote this path as ℛ(J). Write j∈ℛ(J) if the node j appears in the path ℛ(J). Given a path ℛ(J), define the first node ℱJ(0)=J0k for this path, and the following node ℱJ(j)=Jjk for any node j∈ℛ(J).
A sample(τ1,τ2,…,τn) of the TTFs for a system is a family of independent identically distributed random variables. Precisely, the stochastic processes Xk,k=1,2,…,n, are independent, which have the same infinitesimal transition rates Q and the same initial probability function π. Then define(2.7)τk=inf{t≥0:Xtk=m+1}.
Using the above notations, we have
(2.8)τk=∑i∈ℛ(J)Tik,J=(J0k,J1k,J2k,…,Jmk).
Based on the observations τ1,τ2,…,τn and the known characteristic matrix C, our problem is to give estimators π̂ and Q̂ 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 superscriptk of the stochastic processes when the clarity is being held. Firstly, we recall the cumulative distribution function of the TTF τ, (3.1)Fτ(t)=P(τ≤t)=pm+1(t).
Here p(t)=(p1(t),p2(t),…,pm+1(t)) is the probability function of Xt, which satisfies the Kolmogorov's forward equation p'(t)=p(t)Q. From the solution to the equation in exponential form, we have
(3.2)pm+1(t)=πexp(tQ)1DT,
where the vector 1DT of dimensions m+1 is the transpose of 1D=(0,0,…,0,1). So the probability density function of the TTF τ is(3.3)fτ(t)=πexp(tQ)Q1DT.
The above methodology can be used to calculate other probabilities and expectations in this section.Lemma 3.1.
For the above noncanonical construction ofX in Lemma 2.3, we have for 1≤i≤m,
(3.4)P(i∈ℛ(J0,J1,…,Jm)∣τ=t)=1fτ(t)π(i)exp(tQi)Qi1PT,fort≥0.
Here Qi=(qrs′) is a matrix of dimensions (m+2)×(m+2) determined as follows:
(3.5)qrs′={qrs,r<i,s≤i,∑t=i+1m+1qrt,r<i,s=m+2,qrs,i≤r≤m+1,s≤m+1,0,otherwise.1P=(0,0,…,0,1,0) and π(i)=(π1(i),π2(i),…,πm+2(i)) are vectors of dimensions m+2,
(3.6)πs(i)={πs,s≤i,0,i<s<m+1,∑t=i+1m+1πt,s=m+2.Proof.
Define the events(3.7)At={i∈ℛ(J0,J1,…,Jm),τ≤t},t≥0.
Considering all possible realizations of J=(J0,J1,…,Jm), it follows from (2.8) that
(3.8)P(At)=∑ℛ(J)∋iP(J0=j0,Js=ℱ(J,s),s∈ℛ(J),∑s∈ℛ(J)Ts≤t).
Since the characteristic matrix C is an upper triangular matrix, the numbers of nodes along a path are monotonic increasing. Then if J=(j0,j1,…,jm) satisfies ℛ(J)∋i, we have j0≤i and
(3.9)ℱ(J,s)≤i,fors∈ℛ(J),s<i.
Considering another stochastic processX′ with the infinitesimal transition rates Qi and the initial probability function π(i), we can deduce a noncanonical construction consisting of J0′,J1′,…,Jm′,T1′,T2′,…,Tm′. It follows from (3.9) that(3.10)P(At)=∑ℛ(J')∋iP(J0′=j0,Js′=ℱ(J′,s),s∈ℛ(J′),∑s∈ℛ(J')Ts≤t).
For a path J′ of the stochastic process X', the property ℛ(J′)∋i means that the end of the path J′ is the absorbing state m+1. Similarly to (3.2), we have
(3.11)P(At)=P(Xt′=m)=π(i)exp(tQi)1PT.
The result follows from the above equation and
(3.12)P(At∣τ=t)=limϵ↓0P(At+ϵ)-P(At)P(t<τ≤t+ϵ).Similarly to the discussion in Lemma3.1, we have the following Lemma.Lemma 3.2.
For the noncanonical construction ofX in Lemma 2.3, 1≤i≤m, given i∈ℛ(J0,J1,…,Jm), the conditional density function of τ is
(3.13)fτ(t∣i∈ℛ(J))=1pi′π(i)exp(tQi)Qi1PT,fort≥0,
where J=(J0,J1,…,Jm) and pi′=P(i∈ℛ(J)).pi′ 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 Qi. However, since pi′ will be eliminated in the further propositions, it is not necessary to present the expression.Theorem 3.3.
For the noncanonical construction ofX in Lemma 2.3, 1≤i≤m, 1≤j≤m+1, and t≥0, we have
(3.14)P(j=ℱ(J,i),i∈ℛ(J0,J1,…,Jm)∣τ=t)={0,ci,j=0,1fτ(t)π(i)exp(tQij)Qij1PT,ci,j≠0.
Here 1P and π(i) are the same as those in Lemma 3.1. Qij=(q̃rs) is a matrix of dimensions (m+2)×(m+2) determined as follows:
(3.15)q̃rs={qrs,r≠i,qis,r=i,s=iorj,∑t=i+1m+1qit-qij,r=i,s=m+2,0,otherwise.Proof.
Ifci,j=0, it is obvious that P(i∈ℛ(J),j=ℱ(J,i),τ≤t)=0. Otherwise, by constructing an auxiliary stochastic process with the infinitesimal transition rates Qij and the initial probability function π(i), we can obtain, similarly to the proof of Lemma 3.1, that
(3.16)P(i∈ℛ(J),j=ℱ(J,i),τ≤t)=π(i)exp(tQij)1PT,
where J=(J0,J1,…,Jm). The required result can then be derived in the similar way as that in Lemma 3.1.Theorem 3.4.
For the noncanonical construction ofX in Lemma 2.3 and i≤m,t≥0, we have
(3.17)E[Ti∣i∈ℛ(J),τ=t]=-qiiexp(qiit)π(i)exp(tQi)Qi1PT×(tπm+1(i,T)+π(i,T)∫0t(t-u)exp(u(Mi-qiiI))duMi1PT),
where J=(J0,J1,…,Jm). I is an identity matrix and Mi=(Mrs) is a matrix of dimensions (m+2)×(m+2) determined as follows:
(3.18)Mrs={qrs,r<i,s<i,-qriqi,sqii,r<i,i<s≤m+1,∑t=i+1m+1qrt,r<i,s=m+2,qrs,i<r≤m+1,i<s≤m+1,0,otherwise.π(i,T)=(π1(i,T),π2(i,T),…,πm+2(i,T)) is vectors of dimensions m+2 determined as follows:
(3.19)πs(i,T)={πs,s<i,0,s=i,-πiqisqii,i<s≤m+1,∑t=i+1m+1πt,s=m+2.Proof.
It follows from (2.8) that given i∈ℛ(J),
(3.20)τ=Ti+∑j∈ℛ(J),j≠iTj,
where Ti and Ti=∑j∈ℛ(J),j≠iTj, are independent random variables. Given i∈ℛ(J), the conditional density function of Ti is
(3.21)fi(t)=-qiiexp(qiit),fort≥0.
Given i∈ℛ(J), suppose that the conditional density function of Ti is gi(t), where gi(t)=0 for t<0. Since the Jacobian of the transformation (x,y)→(x,x+y) is 1, the joint density of Ti,τ is given by f(t1,t2)=fi(t1)gi(t2-t1). Then for t≥0, we have
(3.22)E[Ti∣i∈ℛ(J),τ=t]=1fτ(t∣i∈ℛ(J))∫0tufi(u)gi(t-u)du=1fτ(t∣i∈ℛ(J))∫0t(t-u)fi(t-u)gi(u)du.
Here the conditional density function fτ(t∣i∈ℛ(J)), given i∈ℛ(J), is determined through (3.13) in Lemma 3.2.
As to the cumulative distribution function ofTi, given i∈ℛ(J), we consider another stochastic process X′ with the infinitesimal transition rates Mi and the initial probability function π(i,T). It follows from Lemma 2.3 that X′ has a noncanonical construction consisting of J0′,J1′,…,Jm′,T1′,T2′,…,Tm′. Similarly to the discussion in Lemma 3.1, we have (3.23)P(Ti≤t,i∈ℛ(J))=P(∑j∈ℛ(J'),j≠iTj≤t,theendofthepathJ′istheabsorbingstatem+1),
where J′=(J0′,J1′,…,Jm′). Hence
(3.24)P(Ti≤t,i∈ℛ(J))=π(i,T)exp(tMi)1PT.
Then it follows from (3.13) and (3.24) that
(3.25)P(Ti≤t∣i∈ℛ(J))=1pi′π(i,T)exp(tMi)1PT.
Since the probability in (3.25) may jump at t=0, the conditional density function corresponding to (3.24) is
(3.26)gi(t)=1pi′πm+1(i,T)δ(t)+1pi′π(i,T)exp(tMi)Mi1PT,fort≥0,
where δ(t) is the Dirac function. pi′=P(i∈ℛ(J)) is the same as that in Lemma 3.2. Then the integral in (3.22) becomes
(3.27)∫0t(t-u)fi(t-u)gi(u)du=-qiipi′texp(qiit)πm+1(i,T)-qiipi′exp(qiit)π(i,T)∫0t(t-u)exp(u(Mi-qiiI))duMi1PT.
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(m+2)×3(m+2) block triangular matrix:
(3.28)exp[t(OIOOOIOO(Mi-qii)I)]=(ItIH(t)OIG(t)OOF(t)).
Then the upper right block H(t) of this matrix gives the integral ∫0t(t-u)exp(u(Mi-qiiI))du. 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 i may never be reached. Since the estimation of Ti is nonsense under this case, we can assign any value to E[Ti∣i∈ℛ(J),τ=t]. We suggest -1/qii.The following corollary can be deduced from the above theorems and the tower property of conditional expectations.Corollary 3.7.
For the noncanonical construction ofX in Lemma 2.3 and 1≤i≤m,1≤j≤m+1,t≥0, we have
(3.29)E[Ti∣τ=t]=piE[Ti∣i∈ℛ(J),τ=t]-1-piqii,P(Ji=j∣τ=t)=P(j=ℱ(J,i),i∈ℛ(J)∣τ=t)-(1-pi)qijqii,P(J0=j∣τ=t)=1fτ(t)πieiexp(tQ)Q1DT.
Here ei is the unit basis vectors (0,0,…,1,0,…,0,0) with “1” in the ith position. J=(J0,J1,…,Jm) and pi=P(i∈ℛ(J0,J1,…,Jm)∣τ=t).Proof.
The first equation follows from Theorem3.4 and the tower property of conditional expectations. The second equation follows from
(3.30)P(Ji=j∣τ=t)=P(j=ℱ(J,i)∣i∈ℛ(J)τ=t)pi+P(Ji=j)(1-pi).
As to the third equation, we have
(3.31)P(J0=j,τ≤t)=P(J0=j)P(τ≤t∣J0=j).
From the solution to the Kolmogorov's forward equation in the exponential form, we obtain
(3.32)P(τ≤t∣J0=j)=πieiexp(tQ)Q1DT.
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 fθ for some hidden random variable Xi, where θ is an element of a parameter space Θ. Then based on a hidden complete data X1,…,Xn, a likelihood ℒ can be written as (4.1)ℒ(θ)=∏k=1nfθ(Xi).
Let a σ- field ℱ=σ{X1,…,Xn}. For a random variable Y∈ℱ, Eθ*Y denotes the expectation of Y with assumption that θ* is the parameter of Xis. Suppose another σ- field 𝒴⊂ℱ.Two iterative steps in the EM algorithm are as follows.(1) Expectation step: fix θ*=θt, then compute Q(·,θ*), where
(4.2)Q(θ,θ*)=Eθ*[log(ℒ(θ))∣𝒴].(2) Maximization step: Set θt+1=θ̂ as a best solution to the following optimization problem:(4.3)maxθ̂∈ΘQ(θ̂,θ*).Under this framework, we propose an EM algorithm which attempts to find a maximum likelihood estimate (MLE) of the parameterθ=(π,Q). Here 𝒴 is a σ- field generated from a sequence (τ1,τ2,…,τn) of TTFs. As for the hidden complete data, we choose T1k,T2k,…,Tmk and J0k,J1k,J2k,…,Jmk, where k=1,2,…,n. As showed in Lemma 2.3, no information is gained or lost by considering stochastic processes Xk,k=1,2,…,n, as hidden processes instead of T.. and J... 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(4.4)ℒ=∏k=1n(∏i=0mqi(Jik)∏i=1m(-qii)exp(qiiTik)).
The factors qi(Jik) in (4.4) represent the probability functions of the hidden random variables Jik. The factors -qiiexp(qiiTik) represent the density functions of Tik. Taking the logarithm of (4.4) yields(4.5)logℒ=∑k=1n(∑i=0m∑j=1mlog(qi(j))I(Jik=j)+∑i=1m(log(-qii)+qiiTik)),
where I(Jik=j)=1 if Jik=j and otherwise I(Jik=j)=0.A new parameter estimate, sayθ̂=(π̂,Q̂), is obtained as the parameter which maximizes the expected value of (4.5) given (τ1,τ2,…,τn) and a current parameter estimate θ*=(π,Q). The constrained maximization gives us the following new estimate:(4.6)q̂i(j)=1n∑k=1nPθ*(Ji=j∣τ=τk),q̂ii=-n∑k=1nEθ*(Ti∣τ=τk),
where the conditional probabilities and the conditional expectations are given in Corollary 3.7. Then it follows from (2.5) that(4.7)π̂j=q̂0(j),q̂ij=-q̂iiq̂i(j),fori≠j,i≠0.
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].
0.10.2111112367111218181818182132364045464750556063636767676772757982828384848485858585858686Using our EM algorithm withm=3, which means that there are transient states 1,2,3, and one absorbing state 4, we obtain (5.1)π̂=(0.000007,0.124512,0.875480,0.0),Q̂=(-1.0070041.0070040.00.00.0-1.0640190.0112951.0527240.00.0-0.0192700.0192700.00.00.00.0).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 (5.2)G(rn)=(∑i=1rTi:n+(n-r)Tr:n)∑i=1nTi:n
against r/n. Here r=1,…,n and Ti:n,i=1,…,n, 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 π̂,Q̂. Due to the sensitivity of TTT plot to Tn:n, 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].
Time intervalNumber of failuresTime intervalNumber of failures153229329436513625722816918108112212111313145155164171181Using our EM algorithm withm=3, we obtain (5.3)π̂=(0.998642,0.001358,0.0,0.0),Q̂=(-0.1729570.1655930.0073630.00.0-1.1645601.1596260.0049340.00.0-4.7670494.7670490.00.00.00.0).
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 withn failure time data, for a system with m transient states, the proposed EM algorithm only needs to compute the exponential of m×m upper triangular matrices for O(nm2) 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.



---

*Source: 242567-2010-09-21.xml*