Uniqueness of Maximum Likelihood Estimators for a Backup System in a Condition-Based Maintenance

Abstract

A parameter estimation problem for a backup system in a condition-based maintenance is considered. We model a backup system by a hidden, three-state continuous time Markov process. Data are obtained through condition monitoring at discrete time points. Maximum likelihood estimates of the model parameters are obtained using the EM algorithm. We establish conditions under which there is no more than one limitation in the parameter space for any sequence derived by the EM algorithm.

1. Introduction

Suppose a backup system is represented by a continuous time homogeneous Markov chain 𝑋={𝑋_𝑡∶𝑡≥0} with a state space 𝑆={0,1,2}. States 0, 1, and 2 are the healthy state, unhealthy state, and the failure state, respectively. Assume that the system is in a healthy state at time 0, and the transition rate matrix is given by Λ=−1/𝜃₁,1/𝜃₁,00,−1/𝜃₂,1/𝜃₂0,0,0, (1) where 𝜃_𝑖∈(𝛼,+∞) for 𝑖=1,2 are unknown. Here 𝛼>0 is a known extreme edge of the parameters. Suppose the system is observed at time points 0,Δ,2Δ,…, where [𝑘Δ,(𝑘+1)Δ] is a benchmark interval. While the system is failed at an inspection, a new system replaces it. Let two processes 𝑌={𝑌(𝑘)∶𝑘=1,…,𝑛}, and let 𝑅={𝑅(𝑘)∶𝑘=1,…,𝑛} be a record of the system. Here 𝑅(𝑘)=1 if the system is failed during (𝑘−1,𝑘), and 𝑅(𝑘)=0 otherwise. And 𝑌(𝑘)=𝑋_𝑘Δ. As a path of 𝑋 is a stepped right-continuous function, 𝑌(𝑘)=0 when a replacement occurs at time 𝑘Δ. Moreover, we set 𝑌(0)=0 for convenience. The process 𝑅 represents the replacement of the system, and the process 𝑌 represents observable information of the system collected through condition monitoring.

The maximum likelihood estimates (MLE) of the model parameters for such models have been studied by [1, 2]. As the stochastic processes, 𝑌,𝑅 are not Markov processes and the sample path of 𝑋 is not observable, the likelihood function of incomplete data (𝑌,𝑅) is complex. Hence, it is difficult to obtain directly the MLE 𝜃=argmax_𝜃∈Θln𝐿(𝜃∣𝑌,𝑅), where 𝜃=(𝜃₁,𝜃₂) and Θ={(𝑥,𝑦)∶𝑥,𝑦>𝛽}. Here 𝛽>0 is a prearranged constant. Both [1, 2] suggest the EM algorithm (see e.g., [3, 4]). Let 𝜃(0) be initial values of the unknown parameters. The EM algorithm works as follows.

The 𝐸 step. For 𝑛>0, compute the following pseudo-likelihood function: 𝑄𝜃^′∣𝜃(𝑛)=𝐸_𝜃(𝑛)ln𝐿𝜃^′∣𝐷∣𝑌,𝑅. (2) Here 𝐷 is the complete data set of the process 𝑋. The forms of the complete data set may be different for different purpose. For example, the forms are different in [1, 2].

The 𝑀 step. Choose 𝜃(𝑛+1)=𝑀(𝜃(𝑛)). Here 𝑀(𝜃)=argsup_𝜃′∈Θ𝑄𝜃^′∣𝜃. (3) The 𝐸 and 𝑀 steps are repeated. According to the theory of EM algorithms (Theorem  1 in [3]), 𝐿(𝜃(𝑛+1)∣𝑌,𝑅)≥𝐿(𝜃(𝑛)∣𝑌,𝑅) for any given initial value 𝜃(0) and 𝑛=0,1,…. It is clear that if an MLE 𝜃 in Θ is one of these fixed-points when it exists.

In this paper, we consider the uniqueness of the MLE 𝜃. As the likelihood function ln𝐿(𝜃∣𝑌,𝑅) of incomplete data is complex, we do not follow the classical method by which the uniqueness of a MLE is demonstrated by establishing the global concavity of the log-likelihood function (see e.g., [5–7].) Alternatively, we investigate conditions under which the operator 𝑀 is a contraction. The conditions ensure that the MLE is unique if it exists. Moreover, the conditions implies that there is not more than one limitations in Θ for different sequences derived by the EM algorithm. For the complete data set we present in the next section, we have the following main theorem of this paper.

Theorem 1. There is not more than one fixed-point of the operator 𝑀 in Θ provided that ▵<(√2/2)𝛼 and the record {𝑌,𝑅} has at least one replacement.

2. Complete Data Set

To establish the expression of operators 𝑄, we present our construction of the path of 𝑋 and the complete data set 𝐷 of the process 𝑋. Suppose that random variables 𝑇^𝑖_𝑘, 𝑖=1,2, 𝑘=0,1,… are independent, 𝑇^𝑖_𝑘 has an exponential distribution, and 𝐸𝑇^𝑖_𝑘=𝜃_𝑖/Δ. As every state of 𝑋 has an exponential duration distribution, we may construct a path of 𝑋 through the approach introduced by Theorem  5.4 in [8]. The path of 𝑋 restricted on 𝑡∈[0,Δ) has the following form: 𝑋_𝑡1_0≤𝑡<Δ=𝑗forΔ^𝑗_𝑖=1𝑇^𝑖₁≤𝑡<Δ^𝑗+1_𝑖=1𝑇^𝑖₁∧1,𝑗=0,1, (4) and 𝑋_𝑡1_{0≤𝑡≤Δ}=2 for Δ(𝑇¹₁+𝑇²₁)≤𝑡<Δ. If 𝑚=lim_{𝑡↑𝑘Δ}𝑋_𝑡≠2 for 𝑘≥1, we can construct the following path of 𝑋 restricted on 𝑡∈[𝑘Δ,(𝑘+1). Consider 𝑋_𝑡1_{𝑘≤𝑡<𝑘+1}=𝑗,Δ∑^𝑗_𝑖=𝑚𝑇^𝑖_𝑘≤𝑡−𝑘Δ<Δ∑^𝑗+1_𝑖=𝑚𝑇^𝑖_𝑘∧1,𝑗≠2,2,Δ∑²_𝑖=𝑚𝑇^𝑖_𝑘≤𝑡−𝑘Δ<Δ. (5) If lim_{𝑡↑𝑘Δ}𝑋_𝑡=2, that is, the system is failed during (𝑘Δ,(𝑘+1)Δ), a new system replaces it at time (𝑘+1)Δ and a path of 𝑋 restricted on 𝑡∈[𝑘Δ,(𝑘+1)Δ) is presented by setting 𝑚=1 in (5). In this paper, the complete data set 𝐷={𝑇^𝑖_𝑘∶𝑖=1,2,𝑘=0,1,…,𝑛−1}.

Our choice of the complete data set 𝐷 ensures a simple form of the operator 𝑄. The log-likelihood function has the following form: ln𝐿𝜃^′∣𝐷=2𝑛lnΔ−𝑛ln𝜃^′₁+ln𝜃^′₂−Δ^𝑛−1_𝑘=0𝑇¹_𝑘𝜃^′₁+𝑇²_𝑘𝜃^′₂. (6) And the Markovian property of the process 𝑋 implies that 𝐸_𝜃𝑇^𝑖_𝑘∣𝑌=𝐸_𝜃𝑇^𝑖_𝑘∣𝑌(𝑘−1),𝑌(𝑘),𝑅(𝑘), (7) where 𝑌(0)=0.

In Table 1, we denote by 𝑛_𝑖 the number of different values of the triple (𝑌(𝑘−1),𝑌(𝑘),𝑅(𝑘)). For example, 𝑛₁ is the number of triple (𝑌(𝑘−1),𝑌(𝑘),𝑅(𝑘))=(0,0,0). It follows from (6) that 𝑄𝜃^′∣𝜃=2𝑛lnΔ−𝑛ln𝜃^′₁+ln𝜃^′₂−Δ⁵_𝑖=1𝑛_𝑖𝑀¹_𝑖𝜃^′₁+𝑀²_𝑖𝜃^′₂. (8)

Table 1 

The number for values of (𝑌(𝑘−1),𝑌(𝑘),𝑅(𝑘)).

Here the following forms of 𝑀^𝑚_𝑖 for 𝑚=1,2 follows from (7) and Table 1. Consider 𝑀^𝑚₁=𝐸_𝜃𝑇^𝑚₁∣𝑇¹₁𝑀≥1,_𝑚2=𝐸_𝜃𝑇^𝑚₁∣𝑇²₁𝑀≥1,_𝑚3=𝐸_𝜃𝑇^𝑚₁∣𝑇²₁𝑀<1,_𝑚4=𝐸_𝜃𝑇^𝑚₁∣𝑇¹₁<1,𝑇¹₁+𝑇²₁𝑀≥1,_𝑚5=𝐸_𝜃𝑇^𝑚₁∣𝑇¹₁+𝑇²₁<1. (9)(10)(11)(12)(13)

For any given 𝜃, it is obvious that the function 𝑄(𝜆₁,𝜆₂)=𝑄(1/𝜆₁,1/𝜆₂∣𝜃) is a concave function for 0<𝜆₁,𝜆₂<∞. Hence there is a unique vector 𝜃^′=𝑀(𝜃) satisfying 𝜕𝑄𝜃^′∣𝜃𝜕𝜃^′_𝑖=0,𝑖=1,2. (14) Moreover, it follows from the definition (3) that if 𝑀(𝜃)∈Θ, 𝑀(𝜃)=𝑀(𝜃). Therefore, every fixed-point of 𝑀 in Θ is a fixed-point of 𝑀 and vice versa. In this paper, we will prove Theorem 1 through studying the number of fixed-points of 𝑀 in Θ. Here we present the from of 𝑀(𝜃)=(𝑀₁(𝜃),𝑀₂(𝜃)) derived by (14). Consider 𝑀_𝑚(𝜃)=Δ𝑛⁵_𝑖=1𝑛_𝑖𝑀^𝑚_𝑖,𝑚=1,2. (15)

3. Two Lemmas

Lemma 2. Consider the following: 𝜕𝑀^𝑚₄𝜕𝜃_𝑖≤Δ𝜃²_𝑖,𝑖=1,2,𝑚=1,2. (16)

Proof. Writing 𝜎={(𝑥₁,𝑥₂)∣𝑥₁,𝑥₂≥0,𝑥₁+𝑥₂≤1} and 𝑉=_𝜎Δ²𝑥_𝑚𝜃₁𝜃₂𝑒^{(−Δ/𝜃1)𝑥1−(Δ/𝜃2)𝑥2}𝑑𝑥₁𝑑𝑥₂,𝐷=_𝜎Δ²𝜃₁𝜃₂𝑒^{(−Δ2/𝜃1)𝑥1−(Δ/𝜃2)𝑥2}𝑑𝑥₁𝑑𝑥₂, (17) we have 𝑀^𝑚₄=𝑉/𝐷. Hence, 𝜕𝑀^𝑚₄𝜕𝜃_𝑖=𝑀^𝑚₄𝑉⁻¹𝜕𝑉𝜕𝜃_𝑖−𝐷⁻¹𝜕𝐷𝜕𝜃_𝑖. (18)
As 0≤𝑥_𝑚≤1 on 𝜎, we have that 0≤𝑉≤𝐷 and 0≤𝑀^𝑚₄≤1. Therefore, 𝜕𝑀^𝑚₄𝜕𝜃_𝑖≤𝑉⁻¹𝜕𝑉𝜕𝜃_𝑖−𝐷⁻¹𝜕𝐷𝜕𝜃_𝑖. (19) Moreover, we have 𝜕𝑉𝜕𝜃_𝑖𝑉=−𝜃_𝑖+Δ³𝜃²_𝑖^𝜎𝑥_𝑚𝑥_𝑖𝜃₁𝜃₂𝑒^{−(Δ/𝜃1)𝑥1−(Δ/𝜃2)𝑥2}𝑑𝑥₁𝑑𝑥₂. (20) As 0≤𝑥_𝑖≤1 on the region 𝜎, it follows from the definition to 𝑉 that there is 0<𝜂₁<1 such that 𝜕𝑉𝜕𝜃_𝑖𝑉=−𝜃_𝑖+𝜂₁Δ𝑉𝜃²_𝑖. (21) Similarly, there is 𝜂₂∈(0,1) such that 𝜕𝐷𝜕𝜃_𝑖𝐷=−𝜃_𝑖+𝜂₂Δ𝐷𝜃²_𝑖. (22) The result follows from (19), (20), and (22).

Lemma 3. Let a function 𝐴∶𝑅→𝑅 be 𝐴(0)=1/2 and 𝐴(𝑥)=1/(1−𝑒^𝑥)+1/𝑥 for 𝑥≠0. We have 0≤𝐴^′(𝑥)≤−1/12 for any 𝑥∈𝑅.

Proof. By a routine analysis, we may obtain that 𝐴^′(0)=−1/12 and for 𝑥≠0, 𝐴^′(𝑥)=𝑒^−𝑥/2−𝑒^𝑥/2−2−𝑥⁻². (23) Moreover, we may obtain that for 𝑥>0, 𝐴^′′(𝑥)=𝐴₂(𝑥)𝑥⁻³𝑒^𝑥−1⁻³, (24) where 𝐴₂(𝑥)=(𝑒^𝑥−1)³−𝑥³𝑒^𝑥(1+𝑒^𝑥). We have 𝐴^′₂(𝑥)=−𝑒^𝑥𝐴₃(𝑥), (25) where 𝐴₃(𝑥)=3𝑥²(1+𝑒^𝑥)+𝑥³(1+2𝑒^𝑥)−6(𝑒^𝑥−1)². A routine analysis may confirm that 𝐴₃⁽⁴⁾(𝑥)=𝑒^𝑥𝐴₄(𝑥), where 𝐴₄(𝑥)=96+96𝑥+27𝑥²+2𝑥³−96𝑒^𝑥. Then it follows from 𝑒^𝑥≤1+𝑥+𝑥²/2+𝑥³/6 that 𝐴₄(𝑥)≤0. Moreover, as 𝐴₃(0)=𝐴^′₃(0)=𝐴₃^′′(0)=𝐴₃⁽³⁾(0)=0, we have that for any 𝑥>0, there is 𝜉∈[0,𝑥] such that 𝐴₃(𝑥)=1𝐴24³⁽⁴⁾(𝜉)𝑥⁴=1𝐴24⁴(𝜉)𝑒^𝜉𝑥⁴≤0. (26) Then it follows from 𝐴₂(0)=0 and (25) that 𝐴₂(𝑥)≥0. Hence, it follows from (24) that 𝐴^′′(𝑥)≥0 for any 𝑥>0. Therefore, for any 𝑥>0−112=𝐴^′(0)≤𝐴^′(𝑥)≤lim_𝑥→+∞𝐴^′(𝑥)=0. (27) The result follows from the above formula, and the fact that 𝐴^′(𝑥) is an even function.

4. Proof of the Main Theorem

We may derive that 𝑀(𝜃) is a contraction by investigating the Jacobian matrix 𝜕(𝑀₁,𝑀₂)/𝜕(𝜃₁,𝜃₂).

By a routine analysis, it follows from (12) that 𝑀¹₂Δ=1+𝐴𝜃₁−Δ𝜃₂,𝑀²₂𝜃=1+²ΔΔ+𝐴𝜃₁−Δ𝜃₂. (28) Hence, it follows from Lemma 3 that 𝜕𝑀¹₂𝜕𝜃₂=Δ𝜃²₂𝐴^′Δ𝜃₁−Δ𝜃₂≤Δ12𝜃²₂. (29) Similarly, we have 𝜕𝑀¹₂𝜕𝜃₁≤Δ12𝜃²₁. (30) As Δ<𝛼<𝜃₂, it follows from (28) and Lemma 3 that 𝜕𝑀²₂𝜕𝜃₁+𝜕𝑀²₂𝜕𝜃₂=1Δ. (31)

It follows from (11) that 𝑀²₅=𝐴(𝜃₂/Δ), and 𝑀¹₅=𝜃₁/Δ. Therefore, 𝜕𝑀²₅𝜕𝜃₁=0,𝜕𝑀¹₅𝜕𝜃₁=0,𝜕𝑀¹₅𝜕𝜃₂=1Δ. (32) From the definition (11), we have that 𝑀²₅=𝐴(Δ/𝜃₂). Then, it follows from Lemma 3 that 𝜕𝑀²₅𝜕𝜃₂≤Δ12𝜃²₂. (33)

It follows from (9) and (10) that 𝑀¹₁=1+𝜃₁/Δ, 𝑀²₁=𝜃₂/Δ, and hence for 𝑚=1,2, 𝜕𝑀^𝑚₁𝜕𝜃₁+𝜕𝑀^𝑚₁𝜕𝜃₂≤1Δ,𝜕𝑀^𝑚₃𝜕𝜃₁+𝜕𝑀^𝑚₃𝜕𝜃₂≤1Δ. (34)

Write 𝑆=max{𝑆₁,𝑆₂}<1, where for 𝑖=1,2, 𝑆_𝑖=sup_𝜃∈Θ𝜕𝑀_𝑖𝜕𝜃₁+𝜕𝑀_𝑖𝜕𝜃₂. (35) As 𝜕𝑀_𝑖/𝜕𝜃_𝑗 are continuous on the convex set Θ, it follows from Theorem  5.19 in [9] that for 𝑥=(𝑥₁,𝑥₂),𝑦=(𝑦₁,𝑦₂)∈Θ, there exist 𝑧∈Θ such that 𝑀(𝑥)𝑦=𝜕𝑀−𝑀¹,𝑀₂𝜕𝜃₁,𝜃_2𝑧𝑥−𝑦. (36) Therefore, we have that ‖𝑀(𝑥)−𝑀(𝑦)‖≤𝑆‖𝑥−𝑦‖, where ‖𝑥‖=max{|𝑥₁|,|𝑥₂|}. Hence, there is not more than one fixed-point of 𝑀 in Θ when 𝑆<1. As every fixed-point of 𝑀 in Θ is a fixed-point of 𝑀, we will prove the theorem by indicating 𝑆<1.

It follows from (29) to (35), (15), and Lemma 2 that 𝑆₁≤𝑛₁𝑛+𝑛₃𝑛+𝑛₅𝑛+𝑛₄𝑛2Δ²𝛼²+𝑛₂𝑛Δ²6𝛼²,𝑆₂≤𝑛₁𝑛+𝑛₃𝑛+𝑛₂𝑛+𝑛₄𝑛2Δ²𝛼²+𝑛₅𝑛Δ²12𝛼². (37)(38) Then it follows from Δ<√2𝛼/2 that 𝑆₁≤1 and 𝑆₂≤1.

The record has at least one replacement. That is, there is 𝑘≥1 such that 𝑅(𝑘)=1. As a new system replaces the old failed system at time 𝑘Δ, we have that 𝑌(𝑘)=0. Now the theorem will be accomplished in two cases.

Case 1. The case 𝑌(𝑘−1)=0. We have that 𝑛₄≠0. It follows from (37) and (38) that 𝑛𝑆≤1−⁴𝑛1−2Δ²𝛼²<1. (39)

Case 2. The case 𝑌(𝑘−1)≠0. According to the condition-based maintenance policy, a failed system is replaced at an inspection. Hence 𝑌(𝑘−1)=1 follows from 𝑌(𝑘−1)≠0. As 𝑌(𝑘−1)=1, 𝑌(𝑘)=0, and 𝑅(𝑘)=1, we have that 𝑛₅≠0. Assume that the last replacement occurs at 𝑠Δ if any replacement occurs before (𝑘−1)Δ. And we write 𝑠=0 when there is no replacement before (𝑘−1)Δ. Now we study the sequence 𝑌(𝑠+1),…,𝑌(𝑘−1), which consists of digits 0 and 1, corresponding healthy and unhealthy states of a system without replacement. It is obvious that there is 𝑠+1≤𝑡≤𝑘−1 such that 𝑌(𝑡)=1 and 𝑌(𝑡−1)=0. Then, 𝑛₂≠0 follows from 𝑌(𝑡−1)=0, 𝑌(𝑡)=1 and 𝑅(𝑡)=0. It follows from (37) and 𝑛₂≠0 that 𝑆₁𝑛≤1−²𝑛Δ1−²6𝛼²<1. (40) Moreover, it follows from (38) and 𝑛₅≠0 that 𝑆₂𝑛≤1−⁵𝑛Δ1−²12𝛼²<1. (41)

5. Example and Discussion

We will apply the EM algorithm to a simulation dataset. Based on this example, we will show the efficiency and accuracy of the EM algorithm. Moreover, by this example, we will show some limitations and shortcomings of Theorem 1. In this example, we make ensembles of 10³ consecutive inspection of a simulating backup system defined by (1). The true parameters are 𝜃₁=10 and 𝜃₂=10/3, which is adopted from [2].

We describe first the process of iterations described by the EM algorithm (2) and (3). For a given couple of initial value 𝜃(0) in a given parameter space Θ, we may derive 𝑀(𝜃(0)) from (15). If 𝑀(𝜃(0))∈Θ, it is the unique solution to (3), and we have that 𝜃(1)=𝑀(𝜃(0)). If we are fortunate and can repeat the operation again and again, then we obtain a sequence 𝜃(𝑛)∈Θ,𝑛=1,2,….

It follows from the expression (8) of 𝑄(𝜃^′∣𝜃) that 𝑄 is continuous with respect to both 𝜃^′ and 𝜃. Similar to the discussion of the Theorem  1 in [4], we can prove that if the limitation of the sequence 𝜃(𝑛)∈Θ,𝑛=1,2,… exists and is also in Θ, then the limitation is a fixed-point of the operator 𝑀. Theorem 1 shows that the limitation is unique for all such sequences.

In the first experiment, we run the EM algorithm for different initial values. In this experiment, we set 𝛽=1 and the parameter space Θ={(𝑥,𝑦)∶𝑥,𝑦>𝛽}. We run the algorithm for 64 couples of initial values which are chosen randomly from 2 to 12. For each couple of initial value, we run the algorithm for 200 iterations. Figure 1 draws the final estimations of the parameters for initial values. We can see that the algorithm converges to the same result for a great range of initial values. As Theorem 1 points out that the number of fixed points is not more than 1. So we can conclude that there is a unique fixed points, and hence it is the MLE of the model parameters on the parameter space Θ={(𝑥,𝑦)∶𝑥,𝑦>1}.

Figure 1 

Estimated parameters for different initial values.

Sometimes, the above procedure of 𝜃(𝑛),𝑛=1,2,… must stop without the output of the estimated parameters. In general, for a 𝜃(𝑛)∈Θ, if 𝑀(𝜃(𝑛)) derived from (15) is not an element of Θ, then 𝑀(𝜃(𝑛))≠𝑀(𝜃(𝑛)). For this case, the solution 𝑀(𝜃(𝑛)) to (3) is on the boundary of Θ. As we do not obtain the explicit expression 𝑀(𝜃(𝑛)) for this case, the procedure is aborted.

In the following second experiment, for the same dataset of the first experiment, we run the EM algorithm for another parameter space Θ={(𝑥,𝑦)∶𝑥,𝑦>𝛽^′} with 𝛽^′=12. We run the algorithm for 32 couples of initial values which are chosen randomly from 20 to 80. As we predict, the procedure is aborted for every couple of initial values. For these initial values, Figure 2 draws the maximal iteration numbers before the procedure is aborted.

Figure 2 

Maximal iteration number for different initial values.

As we know, MLE 𝜃 in Θ, when it exists, is one of fixed-points of the operator 𝑀. However, there may be other fixed-points of 𝑀, such as stationary points of 𝑄. Theorem 1 provides a sufficient condition under which such fixed-point does not exist. Our first experiment and Figure 1 illustrate this fact. However, in some cases, there is not a fixed point of 𝑀 on the parameter space Θ. Such a phenomenon occurs when we set a wrong parameter space as in the second experiment.

6. Conclusion

A parameter estimation problem for a backup system has been considered. We established an EM algorithm, which can be used to iteratively determine the maximum likelihood estimators given observations of the system at discrete time points. It has been found that for any initial values, the sequence derived by the EM algorithm converges to a unique point when the limitation belongs to the specified parameter space.

Acknowledgment

This research was supported by the National Natural Science Foundation of China Grant Nos. 50977073 and 70971109.