Estimation and Reliability for a Special Type of Semi-Markov Process

Corresponding author: Andreas Makrides University o Rouen, France, University of Central Lancashire, Cyprus Email: amakrides@uclan.ac.uk Abstract: This work deals with multi-state systems modelled by means of semi-Markov processes. The sojourn times are seen to be independent not identically distributed random variables and assumed to belong to one of two different general classes of distributions so that the first class is closed under maxima and the second one is closed under minima. We obtain maximum likelihood estimators of the parameters of interest under several estimation schemes and investigate their asymptotic properties. Plug-in type estimators are furnished for various quantities related to the system under study.


Introduction
Let us consider a multi-state system with state space E = {1, 2,…, N}, defined on a probability space (Ω, , P) (cf. Natvig, 1982;El-Neveihi and Proschan, 1984). As it will be introduced in the next section the evolution of the system is assumed to follow a continuous time semi-Markov process.
The MSS could be investigated and analyzed with the use of Markov processes (cf. Limnios and Ouhbi, 2006;Limnios et al., 2005;Lisnianski et al., 2010;Lisnianski and Levitin, 2003). Markov processes are widely used for reliability analysis because the number of failures in arbitrary time intervals in many practical cases can be described as a Poisson process and time up to the failure and repair (or maintenance) is often Exponentially distributed. Such systems play a key role in reliability theory, engineering, economics and finance, geosciences and biomedicine with special interest lying in extreme or rare events like among others, natural disasters, total power supply failures, global economic crises, etc.
The main characteristic of this paper is that the sojourn times in a given state are assumed to belong to one of two different general classes of distributions, cf. Relations (3) and (30) respectively. The first class of distributions is closed under maxima and contains several distributions, like the Bernoulli distribution, the power function distribution and the extreme value Type I distribution.
The second class of distributions is closed under minima and includes the exponential, the Weibull, the Pareto, the Rayleigh and the Erlang truncated Exponential distribution (cf. Balasubramanian et al., 1991). Barbu et al. (2017) the class of distributions closed under minima has been considered. In this paper we deal with: (a) The class of distributions closed under maxima, where we consider several cases of no censoring and censoring at the beginning and/or at the end and appropriate estimators have been proposed; (b) the class of distributions closed under minima, where we consider censoring sojourn times at the beginning and/or at the end.
Our main objective of this paper is the proposal of parsimonious modeling for multistate systems, considering also a semi-Markov framework. Thus, we can have a useful and powerful tool with a reduced number of parameters, which can be of great importance from a practical point of view.
The outline of the paper is the following: In Section 2, a brief discussion about semi-Markov processes and multi-state systems regarding a family of distributions closed under maxima is introduced. Section 3 provides the likelihood function in addition with the maximum likelihood estimators of the parameters under investigation obtained under several statistical settings. In Section 4 a special case for the class of distributions closed under minima is established while the maximum likelihood estimators of the parameters of interest in this are also defined. Observe that: Let Tij be a potential time in state i before transit to the next state j. Let also Fij(t, ij) be the cumulative distribution function (cdf) and ij its m-dimensional parameter. The above distribution is assumed to be absolutely continuous with respect to the Lebesgue measure with density fij (t, ij).
The evolution of the model is governed by the following characteristic: The maximum value of Til identifies the next visited state, assuming that the process was in state i. Hence, under this setting the semi-Markov kernel takes the form:

Class of Distributions
As mentioned at the beginning, we consider two classes of distributions; the first is discussed below while the second one in Section 4. The first one is closed under maxima and the second one is closed under minima. More specifically, although the distribution functions Fij(; ij), have different parameters for different values of I, j = 1,…,N, they have the same form-pattern. Equivalently, the random variables involved in the above distributions are independent but not identically distributed (inid). Any member of the specific class of distribution functions with parameter a can be written in the following form (Balasubramanian et al., 1991): which is assumed to be absolutely continuous w.r.t. the Lebesgue measure with density f(t; a).
The above class of distributions includes the Bernoulli distribution, the power function distribution and the extreme value Type I distribution. A representative Vlad Stefan Barbu et al. / Journal of Mathematics and Statistics 2019, ■■ (■): ■■■.■■■ DOI: 10.3844/jmssp.2019.■■■.■■■ ■■ example comes from structural engineering where engineers are interested in stress and strain diagrams that graphically display the basic material characteristics when designing various types of constructions like bridges, highways or building. Strain is the elongation or contraction of a material per unit length of the material while stress is the applied load divided by the material area it is acting on. For measuring the extend of stress and strain, engineers often use (a) the tensile strength which is the amount of tensile stress that a material can resist before failing, (b) the compressive strength which is he amount of compressive stress that a material can resist before failing and (c) the ductile which is the ability of a material to be subjected to large strains before it ruptures or fails. In such cases the investigator is looking for the maximum time that represents the moment the material ruptures, fails or cracks (cf. Carreira and Chu, 1985;1986;Courtney, 2005).
The following Lemma, shows that the distribution function of the maximum order statistic falls also to class 3, in other words the above class of distributions is closed under maxima. This is stated in the next result.

Lemma 1
Let X1,…,XN be inid random variables such that Xi  F(x; ai) which belongs to class (3). Then the distribution function F (N) of the maximum order statistic X(N) belongs also to (3) (Balasubramanian et al. 1991).
Note that the dependence of semi-Markov kernel to aim is due to the fact that the parameter of the distribution function F(t) that appears in the semi-Markov kernel is im mE a   (Equation (4)).

Proposition 1
According the setting of the present section, the following quantities can be derived: and: By taking the limit: We also need the distribution of the maximum which is given by: Thus, by taking the derivative of the above equation, the pdf in (7) is easily derived.

Maximum Likelihood Estimation
We consider two statistical settings. In the first one we have only one sample path of the system, while in the second one we have several sample paths. On each situation we take into account several cases: One where all sojourn times are observed and one where the sojourn time in the last and/or in the first visited state can be right censored. In addition, the cases where right censoring appears only at the beginning or at the end are also considered. Right censoring at the end is associated with "the lost to follow up" while right censoring at the beginning is associated with the lack of information in the moment the action had began.

Maximum Likelihood Estimation for one Trajectory
Let denote the total observation time with M. For the scope of this section we introduce the processes Ni(t) which counts the number of visits to state i of the process  

The Case of no Censoring for One Trajectory
Firstly, we take into account one single sample path of a semi-Markov process {j0, x1, j1, x2,…, jN(M)} where the evolution of the system for 1 trajectory is given by the following likelihood: where,   1, k i x is the sojourn time in state i during the k th visit, k = 1,…,Ni(M). The superscript "1" indicates the use of a single trajectory; later we will have several trajectories.
Under the class of distributions (3), the likelihood above can be written as: The estimator of aij are obtained through derivation of the preceding log-likelihood w.r.t. aij, I, jE, that is: As for the estimators of Q, p, W and f can be derived using expressions (4) -(7). Let {j0,x1,j1,x2,…,jN(M), uM} a censored sample path, where uM := M-SN(M) is the last censored sojourn time in the last visited state. We follow the methodology proposed in Barbu et al. (2018) for parametric estimation of discrete-time semi-Markov processes. The contribution to the likelihood of a right censored time k, in state i is defined as: Censoring at the Beginning for One Trajectory Let now {x0,j0,x1,j1,x2,…,jN(M)}, where x1 is the first censored sojourn time in the first visited state. Note that the contribution to the likelihood is the same as before. Hence, the likelihood can be written as:

Censoring at the End for One Trajectory
Under the class of distributions in (3), by taking the derivative of the log-likelihood w.r.t. aij, i, jE, we get that the estimator   ij aM is given by solving the equation:

Maximum Likelihood Estimation for Several Trajectories
The following counting processes will be useful later:

The case of no Censoring for Several Trajectories
Let us consider L sample paths of a semi-Markov process, censoring. In this case the likelihood becomes: Note that for L = 1 the likelihood (17) reduces to the likelihood of the one trajectory case given in (11).
Under the class of distributions given in (3) the likelihood for L trajectories without censoring becomes: and we obtain the estimator of aij, namely   , ij a L M , which is given by: The MLE of the initial law is given by:

Censoring at the end for Several Trajectories
, l = 1,…,L, be the L censored sample paths of a semi-Markov process with censoring at time M, then the associated likelihood is: Note that, if the censoring time M in a certain trajectory l is a jump time, then for the corresponding observed censored time we have Consequently, the contribution to the likelihood of the associated term will be equal to 1. For this reason, if no censoring is involved, the uncensored likelihood given in (17) is just a particular case of (20). Note also that for L = 1 the likelihood given in (20) becomes the likelihood of one trajectory with censoring, given in (13).
The estimator   , ij a L M , under the class of distributions in (3), is established by solving the following equation: Regarding the estimator of the initial distribution, given by (19).

■■
k i x is the first censored sojourn time in state i as the first state, during the k th visit, k = 1,…, The estimator   , ij a L M in this case, is obtained by solving the following equation:

Censoring at the Beginning and at the End for Several Trajectories
We can combine the cases considered in Sections 3.2.2 and 3.2.3 so that all trajectories have censoring both at the beginning and at the end. Note through that, from practical point of view, it is important to properly write the case where only several trajectories have censoring at the beginning and/or at the end. This will be done in this section.
Consider the case of L trajectories where some of the sojourn times are censored either at the beginning or/and at the end. Given L censored sample paths of a semi-Markov process:  The likelihood in this case can be written as: Note that, in the above likelihood, due to the appearance of the counting processes , we do not take into account all the trajectories but only those that have censoring at the beginning and/or at the end. All others without censoring, are included in the term For the class of distributions given in (3), the logarithm of the likelihood given in (24) takes the form: log log log log log 1 log 1 .   For the uncensored case the above estimators have explicit forms by combining (27)

Proof
The proof follows along the lines of the proof of Theorem 1 (Barbu et al., 2017) by taking under consideration the convergence of the empirical estimators to the true transition probabilities of the embedded Markov chain, the law of large numbers and the fact that the MRC is regular.
The following Corollary follows immediately from the Proposition above.  (19), are strongly consistent, as L goes to infinity.
The above Corollary holds also for L = 1 with   ij M  given in (12) when M tends to infinity.

Special Case for the Class of Distribution Closed Under Minima
In this Section, a similar framework as in Barbu et al. (2017) is considered in order to investigate the case of censoring at the beginning for the following class of distribution functions which is closed under minima (30). More precisely, in the cited paper we have considered the case of no censoring and the case of censoring at time M for both situations of one sample path and several sample paths. However, here we deal with the general case where some of the sojourn times are censored either at the beginning and/or at the end for several trajectories.
Similarly, as in the case of class (3), the class of distributions given by: is absolutely continues w.r.t the Lebesgue measure and with density f(t; a). Members of the above class if distributions is the Geometric distribution, the Pareto distribution, the Weibull distribution, the Exponential, the Rayleigh and the Erlang truncated exponential.
The following Lemma shows that the class (30) is closed under minima (cf. Balasubramanian et al., 1991).

Lemma 2
Let X1,…,XN be inid random variables such that Xi  F(x; ai) which belongs to class (30). Then the distribution function F (1) of the minimum order statistic X(1) belongs also to (30).
The evolution of the model, in this case, is governed by the following characteristic: The minimum value of Til identifies the next visited state, assuming that the process was in state i. Thus, under this setting the semi-Markov kernel takes the form (cf. Barbu et al., 2017 Consider the case of L trajectories where some of the sojourn times are censored either at the beginning or/and at the end. Recall that   then the associated likelihood: For the class of distributions given in (30), the above likelihood takes the form:

Remark 1
In the case of no censoring, the estimator in (37) is reduced to the estimator given in Barbu et al. (2017) The proof is omitted due to the similarity with the proof of Theorem 1 (Barbu et al., 2017) and Proposition 2 at the end of Section 3.

Markov Renewal Function and Semi-Markov Transition Matrix
The Markov renewal function and the semi-Markov transition function are furnished and discussed in this section.
The Markov renewal function ij(t), i, jE, t≥0, is defined as (cf. Limnios and Oprișan, 2001 The semi-Markov transition matrix (function) is defined as: Limnios and Oprișan (2001) proved that the semi-Markov transition function satisfies the following Markov renewal equation (MRE): where, IN is the N × N identity matrix,  represents the convolution and "diag" is a diagonal matrix. Let where, (t) = (ij(t))i,jE and it is shown that (IN-Q) (-1) (t) = (t), where () (-1)

Reliability/Survival Analysis Indicators
For the purpose of this section we assume that the state space E can be splitted into two subsets, U (containing the functioning states) and D (containing the failure states), such as E = UD and E = UD = , where let say that U = {1,…, n} and D = {n +1,…,N}. Each matrix can be partitioned according to this state space partition.
This section is devoted to the construction of estimators for a variety of reliability indicators in the case of several censored paths at the beginning and/or at the end (Limnios and Ouhbi, 2003).

Reliability
The reliability at time t, R(t), is defined by: where, TD := inf{t|Zt  D} is the lifetime of the system. The following result presents the estimator of the reliability of a semi-Markov system in terms of estimators of basic quantities of a semi-Markov process.

Proposition 5
For a semi-Markov system, the estimator of the reliability at time t > 0 is given by (cf. Ouhbi and Limnios, 1996)

Lemma 3
The estimator of the failure rate of the system is given by:

Proposition 6
The estimator of the availability at time t > 0, for a semi-Markov system, is given by (Ouhbi and Limnios, 1996):

Maintainability
The maintainability of the system, M(t), is defined by: where, TU := inf{t|Zt  U} is the duration of repair.

Proposition 7
The estimator of the maintainability at time t > 0, for a semi-Markov system, is given by (cf. Limnios and Oprișan, 2001): is the mean sojourn time spent in state i. Note that for a regular and positive recurrent MRP, we have mi < , I  E (cf. Limnios and Oprișan, 2001). The estimator of mi can be obtained in two ways, firstly using an MLE considering the estimators obtained in Section 3 and secondly using the empirical estimator. We can estimate the mean sojourn time in state i in two different ways, namely using a plugin MLE (considering the estimators obtained in Section 3) or the empirical estimator, namely: where the estimators of aij are obtained in Section 3. For more details concerning the MTTF of a semi-Markov system, one can see (Limnios and Oprișan, 2001).

Proposition 8
Let a semi-Markov system, where the matrix     n UU I p M  is nonsingular. The estimator of MTTF can be obtained by one of the following two expressions: where, c is the scale parameter and is taken to be 100.
The initial distribution  is considered to follow the discrete Uniform distribution with parameters 1 and N.
In this work we study several cases, like:  Uncensored sample paths  Censored sample paths at the end (50% censored rate for the times associated with the end of the trajectory)  Censored sample paths at the beginning (50% censored rate for the times associated with the beginning of the trajectory)  Censored sample paths at the beginning and/or at the end (50% censored rate for the times associated with the beginning/end of the trajectory) The above cases can also be studied in the case of one trajectory. Variations of the above cases are also possible. Furthermore, one can consider other distributions for the sojourn times that belong to one of the two classes of distributions and investigate a variety of scenarios like the ones above. For the purpose of the simulation study the values are arbitrarily chosen in Table 1.

Several Uncensored Sample Paths
One way of illustrating the accuracy of the estimators is by providing the Squared Errors (S.E.) of the estimators for several values of L. The results are shown in Table 2. Observe that the accuracy for the transition probabilities of the embedded Markov chain is higher than that of the transition rates. As expected the squared errors are extremely good even for very small sample sizes.
Observe that the larger the number of trajectories, the smaller the squared errors of the estimators of aij and pij. On the other hand, according to Table 3, the smaller the value of t the most accurate the estimators of the transition probabilities of the semi-Markov process. In other words, we have more accurate predictions as we get closer to the present.

Several Censored Sample Paths at the Beginning and/or at the End
Let us first explain, the way we simulate the censoring at the beginning. Using a Bernoulli we randomly choose the trajectories that have censored sojourn time in the first visited state. For the trajectories that do have censoring at the beginning, we randomly cut the interval that have been computed as the first sojourn time in two parts using the Uniform distribution with parameters 0 and the first sojourn time that has been computed. The second part is considered to be the censored sojourn time in the first visited state. Of course, one can consider modifications of the above procedure.
As for the trajectories that have censored sojourn time in the last visited state, the same procedure is followed.
It should be noted that the censoring rate which is taken to be equal to 50% refers not to all sojourn times but only to the times associated with the beginning and/or the end of the trajectory. Consequently, the overall percentage of all sojourn times is much less than 50%.
The estimated values of the parameters for L = 1000 rejection with censoring at the beginning and/or at the end one pointed in the Table 4. Table 5 gives a comparison between the estimators of aij and pij obtained for different values of L. The estimators approach the true value as L increases. In addition, Table 6 shows that we obtain good estimators for values of t approaching the upper limit and even better estimators, in terms of the squared errors for small values of t. The estimators in the central part of the time interval appear to be less accurate as compared with the ones at the end points of the intervals. Table 7 provides the estimated values of aij, pij in the special case of one censored sample path at the beginning and at the end, where the observation time M is set to be 100000. For the specific example the number of visited states is 74031.

One Censored Sample Path at the Beginning and at the End
The accuracy of the estimators is studied using the square errors for several values of M. The results for M = 1000, 10000 and 100000 are shown in Tables 8 and 9.
The conclusions are similar to the ones in the case of several censored paths. According to Table 8 as M increases we obtain much better estimators. Table 9 gives better estimators for small values of t with a small reduction of the squared errors while we approach the upper limit of the time variable.

Influence of Taking into Account the Censoring
We would like now to examine the case of estimating the parameters where in the target case we have censoring (for example at the end). We take into account two different cases: (i) The case of estimating the parameters using the formula with censoring at the end; (ii) and estimating the parameters using (incorrectly) the formula without censoring, under the assumption that the first and the last observation times (0 and M) are jump times. This is done in order to investigate the influence of taking into account the censoring. For this procedure we consider both statistical settings with one trajectory and several trajectories.
The average number of visited states in L = 10 trajectories for the values of M in Table 10  Observe that, in the cases where the censored part represents approximately more than 1% of the average number of visited states, the target case gives better estimators than the case where we do not take into account the censored part. However, when the censored part represents less than 1% of the average number of visited states, the contribution of the censored part is not significant.
The number of visited states for the values of M in Table 11 is 7, 19, 147 and 1467 respectively. The length of the trajectory (the total non-censored sojourn time) is 421.7443, 999.8391, 9937.624 and 99988.39 respectively. Observe that, in the cases where the censored part represents the 12% of the number of visited states, the target case gives better estimators than the case where we do not take into account the censored part. However, when the censored part represents less than 12% of the number of visited states, the contribution of the censored part is not significant. Note that the huge difference between the two cases for M = 500 is predictable since the censored part is not to so small as compared to the total trajectory.

Simulations: Reliability
Let us consider the complete (full) case of several sample paths with censoring at the beginning and/or at the end in order to investigate the behavior of several reliability indices that have been proposed in Section 6.
According to Table 12 the Reliability estimator performs better for values of t that belong to the second and third quarter of its domain. However, as for the estimators of Availability and Maintainability, there is no significant difference for the entire spectrum of the time interval.

■■ Conclusion
As mentioned earlier the main objective of this work is the proposal of parsimonious modeling for multi-state systems, considering also a semi-Markov framework. This is a useful and powerful tool applicable to diverse scientific fields such as economics and finance, biomedicine, geosciences, engineering, etc. As future work we wish to implement the proposed methodology to geosciences and in particular to seismology using real seismic data. In such a case the seismic zones are considered as the states of a semi-Markov process and the proposed methodology is employed for estimating the transition probabilities of seismic events transit from a specific seismic zone to another one.