Julien Chiquet
email: julien.chiquet@utc.fr

Nikolaos Limnios
email: nikolaos.limnios@utc.fr

A Method to Compute the Transition Function of a Piecewise Deterministic Markov Process with Application to Reliability

Keywords: Piecewise deterministic Markov Process, Markov Renewal Process, Transition function, Reliability, Survival Analysis

We study the time evolution of an increasing stochastic process governed by a first order stochastic differential system. This defines a particular piecewise deterministic Markov process (PDMP). We consider a Markov renewal process (MRP) associated to the PDMP and its Markov renewal equation (MRE) which is solved in order to obtain a closed-form solution of the transition function of the PDMP. It is then applied in the framework of survival analysis to evaluate the reliability function of a given system. We give a numerical illustration and we compare this analytical solution with the Monte-Carlo estimator.

Introduction

In the literature of stochastic differential systems, a large part is dedicated to those involving diffusion type processes, that lead to the theory of stochastic differential equations (see, e.g., Øksendal, 2003). Another class of stochastic models has arisen to describe the random evolution of processes that do not involve diffusion type motion, but rather the mixture of deterministic motions and random jumps. Such an idea has been suggested by various authors, yet Davis provided a major contribution to this approach : he gave the underlying theory for the class of stochastic models called piecewise deterministic Markov processes (PDMP) in [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF], farther developed in [START_REF] Davis | Markov models and Optimization[END_REF] and [START_REF] Jacobsen | Point Process Theory and Applications[END_REF]. In [START_REF] Koroliuk | Stochastic Systems in Merging Phase Space[END_REF] the limit theory is studied, in a functional setting. [START_REF] Lapeyre | Introduction to Monte-Carlo Methods for Transport and Diffusion Equations[END_REF] worked with transport processes, which are a special case of PDMP, in order to give a stochastic interpretation of the transport equations used in physics for the modeling of the motion of particles. This family of processes are nowadays much used in several applications, e.g., reliability analysis [START_REF] Devooght | Probabilistic dynamics as a tool for dynamic psa[END_REF] or insurance [START_REF] Embrechts | Stochastic Processes: Theory and Methods[END_REF]. In [START_REF] Chiquet | Estimating stochastic dynamical systems driven by a continuous-time jump markov process[END_REF] The main contribution of this paper is to provide a closed-form solution for the transition function of a PDMP defined through a first order differential equation.

For this purpose, we associate to this PDMP a Markov renewal process (MRP), thus we will refer to results from this theory (see e.g. [START_REF] Limnios | Semi-Markov processes and Reliability[END_REF].

The outline is the following: in section 2, we give the model settings and standard notation for MRP. We then build and solve a Markov renewal equation (MRE) for the transition function of the PDMP. In section 3, an application for survival analysis is studied, where the numerical computation of the reliability is detailed.

We compare the result obtained through the direct resolution of the MRE with the empirical estimator obtained with the Monte-Carlo method.

Let (Z t , t ∈ R + ) be a real-valued stochastic process starting almost surely from z > 0 and governed by the following first order differential system

Żt = C(Z t , X t ), Z 0 = z, (1)
with the following assumptions:

A. 1 The process (X t , t ∈ R + ) is an irreducible Markov process with a countable state space E, an initial distribution α i = P(X 0 = i), and a matrix generator A = (a ij ) i,j∈E such that a ij ≥ 0, for all i = j, and a ii = -a i = -k∈E,i =k a ik .

A

: E -→ R, for x, y ∈ R + and i ∈ E, such that |C(x, i) -C(y, i)| ≤ f (i) |x -y|.
Each path Z t (ω) is built in a piecewise manner according to the function C and a path X t (ω) of the jump Markov process. That is, for t ∈ R + , the Cauchy problem

Żt (ω) = C(Z t (ω), i), Z 0 (ω) = z, i ∈ E, (2)
has a unique solution built on the successive intervals [S n (ω), S n+1 (ω)), where (S n , n ∈ N) is a random sequence describing the jump times of X t . For any t < S 1 (ω), we denote by ϕ z,i (t) the solution of (2), where X 0 (ω) = i, i.e., ϕ z,i (t) is the solution before the first jump time of X t , conditionally to the starting value of (Z 0 , X 0 ) = (z, i). We also assume that Z 0 and X 0 are independent. In the sequel, we focus on the transition function P , defined by

P ij (z, B, t) := P z,i (Z t ∈ B, X t = j), i, j ∈ E, B ∈ B, (3)
where B is a subset of B, the Borel σ-field of R + and

P z,i (•) := P(•|Z 0 = z, X 0 = i).
A Markov process is a special MRP, thus we may associate to (Z t , X t ) the extended

MRP (ζ n , Y n , S n , n ∈ N) such as ζ n = Z Sn , Y n = X Sn , n ∈ N. The process (Y n , S n ) is a standard MRP, while (ζ n , Y n , S n
) is an extended one. In the homogeneous case the associated semi-Markov kernel Q is defined, for t > 0, by

Q ij (z, B, t) := P z,i (ζ 1 ∈ B, Y 1 = j, S 1 -S 0 ≤ t). (4)
The Stieltjes-convolution is denoted by " * ", hence, the convolution of Q with a measurable function φ on the space

R + × E is (Q * φ) ij (z, B, t) = k∈E R + t 0 Q ik (z, dy, ds)φ kj (y, B, t -s),
for i, j ∈ E, z > 0 and B ∈ B. In the same way, the successive n-fold convolutions of the semi-Markov kernel are defined recursively. For n = 0, 1,

Q (0) ij (z, B, t) = 1 {i=j} 1 B (z)1 R + (t), Q (1) ij (z, B, t) = Q ij (z, B, t),
where 1 B (x) is the indicator function, i.e., 1

B (x) = 1 if x ∈ B, 0 otherwise. For n ≥ 2, the n-fold convolution turns to Q (n) ij (z, B, t) := (Q * Q (n-1) ) ij (z, B, t) = k∈E R + t 0 Q ik (z, dy, ds)Q (n-1) kj (y, B, t -s). (5)
The Markov renewal function Ψ, which plays a central role, is defined by

Ψ ij (z, B, t) = n≥0 Q (n) ij (z, B, t). (6)
In the case at hand, we have

(ζ n , Y n , S n ) a normal MRP, that is, Ψ ij (z, B, t) < ∞
for any fixed t > 0, z > 0, B ∈ B and i, j ∈ E.

A MRE has the following form

Θ ij (z, B, t) = g ij (z, B, t) + (Q * Θ) ij (z, B, t), ( 7
)
where g is a known function defined on R + ×E ×R + and Θ is the unknown function.

The solution (see e.g [START_REF] Koroliuk | Stochastic Systems in Merging Phase Space[END_REF] is given by

Θ ij (z, B, t) = (Ψ * g) ij (z, B, t). (8)
Let us now take advantage of Markov renewal theory to build a solvable MRE for P .

For this purpose, we first need to calculate Q. This is done in the following Lemma.

Lemma 1 The semi-Markov kernel Q of the MRP (ζ n , Y n , S n ) verifies, for i = j, Q ij (z, B, dt) = a ij e -a i t δ ϕ z,i (t) (B)dt, (9)
where δ x (B) is the Dirac distribution, equal to 1 if x ∈ B, 0 otherwise.

PROOF.

Assuming that S 0 = 0, and conditioning on definition (4), we get

Q ij (z, B, dt) = P z,i (Y 1 = j, S 1 ∈ dt)P z,i (ζ 1 ∈ B|Y 1 = j, S 1 = t).
First, since S 1 and Y 1 are independent, we have

P z,i (Y 1 = j, S 1 ∈ dt) = a ij e -a i t dt
from usual results of Markov theory. Second, Z t is fully characterized by ϕ z,i (t) before the first jump time S 1 , thus

P z,i (ζ 1 ∈ B|X 1 = j, S 1 = t) = P z,i (Z t ∈ B) = δ ϕ z,i (t) (B)
. Indeed, the probability P z,i (Z t ∈ B) is zero everywhere, excepted for the time points where B is reached. We hence get the expected result.

Proposition 2 The transition function P of (Z t , X t ) is governed by the MRE

P ij (z, B, t) = g ij (z, B, t) + (Q * P ) ij (z, B, t), whose solution is P ij (z, B, t) = (Ψ * g) ij (z, B, t
) and where

g ij (z, B, t) = e -a i t 1 {i=j} 1 B (ϕ z,i (t)). ( 10
)
PROOF. It is convenient to make appear S 1 in (3). Hence,

P ij (z, B, t) = P z,i (Z t ∈ B, X t = j, S 1 > t) P 1 + P z,i (Z t ∈ B, X t = j, S 1 ≤ t) P 2 .
Before the first jump of X t , i = j and Z t evolves according to ϕ z,i (t). Thus, P 1 = e -a(i)t 1 B (ϕ t (z, i))1 {i=j} . From Total Probability Theorem, it holds for P 2 that

P 2 = k∈E k =i t 0 P z,i (Z t ∈ B, X t = j|Y 1 = k, S 1 = s)P z,i (Y 1 = k, S 1 ∈ ds).
As long as P z,i (Y 1 = k, S 1 ∈ ds) = a ik e -a i s ds, and noticing that P z,i (Z t ∈ B, X t = j|Y 1 = k, S 1 = s) = P kj (ϕ z,i (s), B, t -s), then P 2 is fully characterized. Finally,

P ij (z, B, t) = e -a i t 1 {i=j} 1 B (ϕ z,i (t)) + k∈E k =i a ik t 0
e -a(i)s P kj (ϕ z,i (s), B, t -s)ds, which may be written, with Q given by ( 9) and g given by ( 10), as

P ij (z, B, t) = g ij (z, B, t) + k∈E R + t 0 Q ik (z, dy, ds)P kj (y, B, t -s).
This last equation is of the general form (7), whose solution is given by (8).

Application to reliability

As an application, let us now study system (1) under assumptions A.1, A.2 and E a finite state space. The system remains reliable until Z t reaches a critical threshold ∆ > z > 0 at random time τ , which is the failure time, with a cumulative distribution function F τ . A schematic view of the parallel evolution of both components Z t and X t can be seen on Figure 1, where three typical paths of the PDMP are plotted. Of course, since, C > 0, Z t is an increasing process. Looking towards reliability

State Space R + of Z t Time Z t (ω 1 ) Z t (ω 2 ) Z t (ω 3 ) failure boundary {∆} State space E of X t Time X t (ω 1 ) X t (ω 2 ) X t (ω 3 ) τ(ω 1 ) τ(ω 2 ) τ(ω 3 )
Fig. 1. Schematic view of the PDMP -three paths of (Z t , X t ) analysis, we put U = [z, ∆) and D = [∆, ∞) the respective sets of good states and failure states for Z t . Due to the continuous, increasing evolution of Z t , it necessarily passes through ∆. The reliability function R is defined as follows:

R(t) = P((Z t , X t ) ∈ U × E) = i,j∈E α i P ij (z, U, t) = 1 -F τ (t).
Through Proposition 2, P is known. Hence R (as well as F τ ) is fully characterized:

R(t) = i,j∈E α i × (Ψ * g) ij (z, U, t). ( 11
)
Note that the kernel Q can be calculated at a given time point t > 0 for the Borel subset U , by integrating ( 9). As a matter of fact, it holds that

Q ij (z, U, t) = a ij t 0 e -a i s 1 U (ϕ z,i (s))ds = a ij a i 1 -e -a i min(t z,i (∆),t) ,
where t z,i (∆) = inf {t ≥ 0 : ϕ z,i (t) ≥ ∆} .

Numerical implementation

The numerical calculation of R successively requires the kernel Q, the n-fold convolutions Q (n) for each n ≥ 0, the Markov renewal function Ψ built upon the Q (n) and the transition function P , by a convolution between g and Ψ. Since convolution products are time consuming, any simplification would mean a great time saving.

By (9), the n-fold convolution (5) of Q turns to

Q (n) ij (z, B, t) = k∈E k =i a ik t 0 e -a i s Q (n-1) kj (ϕ z,i (s), B, t -s)ds, (12)
hence removing the integral on R + , thanks to the Dirac distribution. Since our main interest is the reliability, we compute P just for the subset B ≡ U , that is,

P ij (z, U, t) = U t 0 Ψ ij (z, dy, ds)e -a j (t-s) 1 U (ϕ y,j (t -s)). (13)
Indeed, the sum on E has been removed thanks to the structure of g. Furthermore, the integration on y ∈ R + is limited on U since 1 U (ϕ y,j (t -s)) is zero elsewhere.

Theses functions have to be properly discretized to achieve the numerical computation. In the following, a function with an upper-index "#" means its discretized version. This discretization must be operated on both intervals U = [z, ∆) and [0, t], thus we set two numerical partitions {z = y 0 < y

1 < • • • < y < • • • < y L = ∆ -} and {0 = t 0 < t 1 < • • • < t m < • • • < t M = t}.
Both L and M , being the respective numbers of discretization steps for [z, ∆) and [0, t], have to be sufficiently large. When L, M → ∞ each numerical function tends to the associated "true" one. For instance, when L, M → ∞, then Q # → Q uniformly w.r.t a given matrix norm, by example, ||Q|| = max i,j Q ij (z, y, t) with t, z, B fixed. Hence, the numerical version of ( 13) is

P # ij (z, U, t) = y ∈[z,∆) tm∈(0,t] ∆ yt Ψ # ij (z, y , t m )e -a j (t-tm) 1 ϕ y ,j (t-tm) (U ),
where ∆ yt Ψ # ij (z, y , t m ) is the only unknown, which stands for the numerical evaluation of Ψ(z, dy, ds) in ( 13). By the definition (6), it can be evaluated through

∆ yt Ψ # ij (z, y , t m ) = n≥0 ∆ yt Q #(n) (z, y , t m ). ( 14
)
The derivative ∆ yt Q #(n) is estimated by finite differences on y and t: 15) is obtained by the numerical version of (12):

∆ yt Q #(n) (z, y , t m ) = [Q #(n) (z, y , t m ) -Q #(n) (z, y -1 , t m )] -[Q #(n) (z, y , t m-1 ) -Q #(n) (z, y -1 , t m-1 )]. (15) Each element Q #(n) in (
Q #(n) ij (z, y , t m ) = k∈E k =i a ik tm∈(0,t] e -a i tm Q #(n-1) kj (ϕ z,i (t m ), y , t -t m )∆t m ,
with ∆t m = t m -t m-1 , the time-step discretization. Finally, we point out that the sum on the n-fold convolutions of the kernel in ( 14) is truncated from the rank n * , provided that ||Q #(n * ) || < ε. We put ε a small real number, chosen closed to the machine precision. Note that the integer n * is finite since Q

#(n) ij (z, y, t) ---→ n→∞ 0 for a normal MRP with fixed values of i, j ∈ E, t > 0, z > 0 and y ∈ [z, ∆].

Numerical results

To check the validity of our results, let us study the process Z t governed by Żt

= aZ t × X t , Z 0 = z, (16)
with a = 0.0075, z = 10, ∆ = 50. The process X t is a three-state jump Markov process with E = {1, 2, 3} , a matrix generator and an initial law given by

A =           -0.02 0.02 0 0.027 -0.03 0.003 0.01 0 -0.01           , α = 2 3 1 3 0 .
The coupled process (Z t , X t ) defines a PDMP whose transition function and reliability function can be computed through Proposition 2 and (11), based upon the numerical implementation described above for which we put M = L = 100 points of discretization. Another way to compute the reliability is the Monte-Carlo method, which consists in simulating a large number of paths of Z t and counting when the state {∆} is reached or not. We use the empirical estimator applied on K = 10000 paths (Z k t ) k=1,...,K simulated through Monte-Carlo techniques, that is

R(t) = 1 K K k=1 1 {Z k t <∆}
. This estimator is compared with the direct calculus of R through the Markov renewal argument developed in this paper. Results can be found on Figure 2 where the Monte-Carlo estimator is used as a reference to check the validity both of the theoretical results and the numerical implementation. One can see that the reliability curves obtained with the two methods are quite close.

  ; Chiquet et al. (to appear), we took advantage of a PDMP for the modeling of degradation mechanisms arising in a structure subject to random environmental effects.

Fig. 2 .

 2 Fig. 2. Reliability estimation with two different methods

  . 2 The function C : R + × E -→ R + is measurable, strictly positive and Lipschitz w.r.t. the first argument, i.e, there is a function f

Because models based on system (1) are encountered in the literature for many different fields of applications, we think that the reliability results of this paper may be useful for engineering studies, such as structural mechanics (see, e.g., Chiquet et al., to appear). Yet, to be more suitable to real phenomena, the case of PDMP which are not only monotone should be considered, as further developments.


Julien Chiquet
email: julien.chiquet@utc.fr

Nikolaos Limnios

Mohamed Eid

Piecewise deterministic Markov processes applied to fatigue crack growth modelling

Keywords:

Markov renewal process Averaging principle Maximum likelihood estimation Fatigue crack growth

In this paper, we use a particular piecewise deterministic Markov process (PDMP) to model the evolution of a degradation mechanism that may arise in various structural components, namely, the fatigue crack growth. We first derive some probability results on the stochastic dynamics with the help of Markov renewal theory: a closed-form solution for the transition function of the PDMP is given. Then, we investigate some methods to estimate the parameters of the dynamical system, involving Bogolyubov's averaging principle and maximum likelihood estimation for the infinitesimal generator of the underlying jump Markov process. Numerical applications on a real crack data set are given.

Introduction

Let us consider a structure subject to a main degradation process (Z t , t ∈ R + ) which increases randomly on R * + = (0, ∞) until it reaches a critical value > 0, meaning the failure of the structure. It is assumed that Z t is an observable process whose sample paths are obtained from experimental feedback. The reliability of such a system may be expressed by

R(t) = P(Z t < ), (1)
that is the probability that the process Z t does not reach the failure boundary on the whole observation interval [0, t].

The dynamical evolution of the continuous-time, real-valued process Z t is modelled by a first order stochastic dynamical system of the following form

dZ t dt = C(Z t , X t ), ( 2
)
where C is a function with the appropriate existence and uniqueness properties. The general form of C is known from physical assumptions on the real phenomenon. The process Z t describes the time evolution of this phenomenon whose modelling through classical differential equations is not possible due to, for instance, some random environmental effects or some inhomogeneity in material of the structure. The stochastic process X t with state space E is thus introduced to handle with the randomness of the system.

In our study, X t is assumed to be a jump Markov process. The use of jump processes in the framework of stochastic dynamical systems modelling seems interesting from a physical point of view. As a matter of fact, it fits well with the mathematical description of a degradation mechanism, i.e., a process observed on a continuous state space, evolving continuously in time, but whose growth rate is likely to change through some "chocs", even very small, occurring at discrete instants of time. Both intensities and time occurrences of the chocs are random, as described by the jump Markov process.

The coupled process (Z t , X t ) has the Markov property and no further assumption is required for Z t . The process (Z t , X t ) belongs to the class of stochastic models sometimes called piecewise deterministic Markov processes (PDMP). PDMP are non-diffusion type stochastic models, combined with the mixture of deterministic and jump motions. [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of non-diffusion stochastic models[END_REF][START_REF] Davis | Markov Models and Optimization[END_REF] gave a full description of these stochastic processes as well as the basic associated stochastic calculus. Eq. ( 2) has also been referred in the literature as a stochastic version of the transport equations, used to model some particle systems (see, e.g., [START_REF] Papanicolaou | Asymptotic analysis of transport processes[END_REF][START_REF] Lapeyre | M éthodes de Monte Carlo pour les Équations de Transport et de Diffusion[END_REF], thus (Z t , X t ) is also known as a transport process. [START_REF] Ezhov | Markov processes homogeneous with respect to the second component I[END_REF] also used similar processes under the name of Switching processes. [START_REF] Koroliuk | Stochastic Systems in Merging Phase Space[END_REF] gave some asymptotic results corresponding to a large class of processes including PDMP. Note also that [START_REF] Jacobsen | Point Process Theory and Applications[END_REF] recently studied non-homogeneous PDMP and built the associated likelihood. We want to give in this paper some practical estimation methods for a particular PDMP dedicated to degradation modelling.

As a standard example of a degradation mechanism, the fatigue crack growth problem is handled in this paper to validate the models. Fig. 1 represents the extensive set of crack growth data obtained by [START_REF] Virkler | The statistical nature of fatigue crack propagation[END_REF] during rigorously controlled laboratory measurements. Roughly speaking, the data set consists in 68 specimens of the same material containing initially the same artificially made 9 mm-crack, that is Z 0 = 9 mm. These samples were exposed to the same cyclic loading. The sample paths represent the evolution in time of the crack size for the 68 specimens, whose measures were performed until they reach the critical value = 49. 8 mm. One can see the large variability of the crack size despite a well controlled experiment. Thus, stochastic dynamical systems were found well adapted to the modelling of fatigue crack growth (see, e.g., [START_REF] Lin | A stochastic theory of fatigue crack propagation[END_REF][START_REF] Bogdanoff | Probabilistic Models of Cumulative Damage[END_REF]Sobczyk, 1993;[START_REF] Tanaka | Importance sampling simulation for a stochastic fatigue crack growth model[END_REF], for an introduction).

The paper is organized as follows: in the next section, we present the family of stochastic dynamical systems involving jump Markov processes of the general form (2) and the settings associated to the modelling of degradation mechanisms. We use Markov renewal theory to obtain a solution of the transition function of the system, with which an expression of the reliability can be given. The third section is dedicated to statistical inferences for such a system, more particularly concerning the estimation of the generator of the jump Markov process X t , whose paths cannot be directly observed. The only available data concern the process Z t , whose paths are defined on random length time intervals. Finally, numerical applications for the crack growth problem are given.

The model settings

Let (Z t , t ∈ R + ) be an increasing stochastic process on R * + , with initial condition Z 0 = z. Looking toward reliability analysis, we define for the process Z t a set of working states U = [z, ), with 0 < z < , and a set of down states D = [ , ∞). Its evolution is continuous in time with values in the real positive axis, thus Z t necessarily passes through the point while reaching the set of down states.

The evolution of Z t is modelled through the following stochastic dynamical system

dZ t dt = C(Z t , X t ), Z 0 = z, (3)
with the following assumptions:

A. 1. The process (X t , t ∈ R + ) is an irreducible Markov process with a finite state space E, initial law i = P(X 0 = i), for all i ∈ E, stationary distribution = ( i ) i∈E and a matrix generator A = (a ij ) i,j∈E such as a ij 0 for all i = j, and

a ii = -a i = - k∈E,i =k a ik .
A. 2. The function C : R * + × E -→ R * + is measurable, strictly positive and of the Lipschitz type, that is, there is a function

f : E -→ R such that, for x, y ∈ R * + and i ∈ E, |C(x, i) -C(y, i)| f (i)|x -y|. A. 3. We have P(Z 0 = z) = 1, with 0 < z < .
Because X t is Markov, it can be easily seen that the coupled process (Z t , X t ) is a Markov homogeneous process on the state space R * + × E, provided that Z 0 and X 0 are independent. Obviously, the process Z t evolves in a deterministic way between the jumps of the process X t , occurring at random discrete instants of time

S 0 = 0 < S 1 < • • • < S n < • • •.
Hence, between the jumps, the evolution of a path of Z t is fully determined by a classical deterministic differential system. For instance, say that X t = i for t < S 1 where S 1 stands for the first jump time of X t . Then the Cauchy problem

dZ t dt = C(Z t , i), Z 0 = z,
has a unique solution denoted by z,i (t), that is (Z t = z,i (t)) 0 t<S 1 . We can thus build in a piecewise manner a unique solution Z t for system (3) on the successive jump time intervals [S n , S n+1 ) of X t .

There is much to say about the Markov process (Z t , X t ). First, its generator B (see, e.g., [START_REF] Lapeyre | M éthodes de Monte Carlo pour les Équations de Transport et de Diffusion[END_REF] is given by

Bf (z, i) = C(z, i) j jz f (z, i) + j∈E a ij f (z, j),
with i ∈ E, z ∈ R * + and f a measurable and differentiable function on the first argument. Unfortunately, the infinitesimal generator cannot always be exploited easily numerically. Thus, it can be more suitable to focus on an other important function of (Z t , X t ), namely the transition probability function P(t):

P ij (z, B, t) := P(Z t ∈ B, X t = j|Z 0 = z, X 0 = i).
(4)

It is defined for all z ∈ R * + , i, j ∈ E and B ∈ B, with B being the Borel -field of R * + . The transition function permits to give a simple expression of the probability distribution (t) of (Z t , X t ), using a conditional probability development, thus characterizing well the PDMP:

j (B, t) = P(Z t ∈ B, X t = j) = i∈E R * + P(Z t ∈ B, X t = j|Z 0 = y, X 0 = i) z (dy) i = i∈E i P ij (z, Bt),
where z is the Dirac distribution at point z.

Expression (1) of the reliability can be advantageously defined as a function of the coupled process (Z t , X t ) rather than as an expression of Z t only. Concerning (Z t , X t ), U × E is the set of working states while D × E is the set of down states. Moreover, we consider here that the system is not a reparable one. Indeed, we have

R(t) = P((Z t , X t ) ∈ U × E) = j∈E j (U, t) = i,j∈E i P ij (z, U, t).
(5)

Please cite this article as:

JSPI : 3675 ARTICLE IN PRESS

We also introduce the hitting time to D for the process Z t , that is the failure time, defined by

= inf{t 0 : Z t ∈ D} = inf{t 0 : (Z t , X t ) ∈ D × E}. (6)
The distribution function F of is directly linked to the reliability through

F (t) = 1 -R(t).
A Markov process is a special case of a Markov renewal process, thus we can associate to (Z t , X t , t ∈ R + ) the extended Markov renewal process (Z n , J n , S n , n ∈ N), where

Z n = Z Sn , J n = X Sn , n ∈ N,
with (S n ) n 0 being the jump times of X t . The Markov renewal function (t), which represents at time t the mean number of times a set of R × E is visited by the process (Z t , X t ), starting from a given point (z, i), is given by

ij (z, B, t) = E z,i [#transitions in B, j] = n 0 Q (n) ij (z, B, t),
where Q is the semi-Markov kernel defined by

Q ij (z, B, t) = P(Z n+1 ∈ B, J n+1 = j, S n+1 -S n t|Z n = z, J n = i).
The successive n-fold convolution Q (n) for t 0 are defined in a Stieltjes sense, that is

Q (0) ij (z, B, t) = 1 {i=j} 1 B ( z,i (t)), Q (1) ij (z, B, t) = Q ij (z, B, t), t > 0,
and, for n 2,

Q (n) ij (z, B, t) = k∈E R + t 0 Q ik (z, dy, ds)Q (n-1) kj (y, B, t -s).
It follows from results of Markov renewal theory (e.g., [START_REF] Limnios | Semi-Markov Processes and Reliability[END_REF]) the following proposition.

Proposition 1. The transition function P(t) is given by the following Markov renewal equation

P ij (z, B, t) = g ij (t; z, B) + k∈E R + t 0 Q ik (z, dy, ds)P kj (y, B, t -s), ( 7
)
where

Q ij (z, B, dt) = a ij e -ait z,i (t) (B) dt, and
g ij (z, B, t) = 1 {i=j} 1 B ( z,i (t))e -ait ,
where 1 B (x) is the indicator function of the subset B, i.e., 1 B (x) = 1 if x ∈ B and 0 otherwise. Solution of (7) is obtained through a Markov renewal argument, and

P ij (z, B, t) = R+ t 0 ij (z, dy, ds)g jj (y, B, t -s). ( 8
)
Expression ( 8) can be computed numerically, thus with (5), the reliability R(t) of the system as well as the cumulative distribution function of the failure time F can also be computed. An example of a numerical computation of R(t), based upon expression (8), is detailed in [START_REF] Chiquet | Modelling the reliability of degradation processes through Markov renewal theory[END_REF] for a given system with general form (3).

As statistical inference for the dynamical system (3) is the main concern of this paper, we shall consider K independent sample paths of Z t obtained from lab experiments. The paths are denoted by (Z k t ) k=1,...,K

, which are supposed to be the only available data for estimating the whole dynamical system. The observation procedure is the following: each path Z k t is observed on some random time intervals [0, k ∧ T obs ], with the k 's being K independent copies of the failure time defined in (6). It means that each path of Z t is observed until it reaches the failure threshold , provided that < T obs . Here, T obs is a real positive constant representing the maximal time of observation, beyond which the lab experiment become too expensive. If we always can afford to observe the system until it reaches , then we just put T obs = ∞, what is the case for the Virkler data set seen on Fig. 1. We point out that T obs is the same for any path. We also denote by t k the successive realizations of the random variables k ∧ T obs . Fig. 2 shows a schematic representation of the observation procedure for three paths of Z t and the corresponding paths of X t , which are obviously defined on the same random time intervals as the corresponding Z k t 's.

Please cite this article as:

Estimation

In this section we develop the scheme used for the estimation of the system, where we just have some K paths of the observable degradation process Z t . The estimation scheme of (3) splits in two parts:

(1) First, we study the asymptotic behavior of the system using the Bogolyubov's averaging principle (see [START_REF] Bogolyubov | Asymptotics Methods in the Theory of Nonlinear Oscillations[END_REF]. With this technique, the random component "disappears" and we obtain a deterministic system that corresponds to the mean behavior of (3). A regression analysis is performed on the asymptotic, deterministic system so that we estimate the fixed parameters in function C, whose general form is assumed to be known. (2) Second, we turn back to the remaining estimation of the random part of the system, that is, of the jump Markov process X t .

To this purpose, we estimate the paths of X t which are not directly observed, as well as its state space. We then build the associated likelihood function, keeping in mind that these paths are defined on random time intervals. The generator of X t is estimated by maximizing an approached likelihood function.

Bogolyubov's averaging principle

An approximation of Z t is obtained by analyzing system (3) in a series scheme (see, e.g., [START_REF] Koroliuk | Stochastic Systems in Merging Phase Space[END_REF], that is by studying the weak convergence, when → 0, of

dZ t dt = C(Z t , X t/ ), Z 0 = z, ( 9
)
where X t is assumed to be ergodic. In fact, the change of scale t → t/ is performed for X t in order to see the behavior of the dynamical system when the random component X t just adds the information it would add after a very long time of observation of (9), since t/ → ∞ when → 0. This so-called averaging approximation was first introduced by Bogolyubov (see [START_REF] Bogolyubov | Asymptotics Methods in the Theory of Nonlinear Oscillations[END_REF] who showed that (9) converges weakly when → 0 to the following deterministic system

d z t dt = C( z t ), z 0 = z, ( 10
)
with z t the limit deterministic process and C a mean function defined by

C(z) = lim T→∞ 1 T T 0 C(z,

JSPI : 3675 ARTICLE IN PRESS

In the particular case where X t is an ergodic jump Markov process with a stationary law , we have

C(z) = i∈E C(z, i) i .
Through this averaging technique, we have a limit deterministic system (10) associated to stochastic differential system (3). The fixed parameters appearing in the function C are the same as the ones appearing in C but in (10) the random part was "eliminated ": with the K sample paths Z k t , we can perform a classical regression analysis on (10) to estimate the fixed parameters of C. By this mean, we obtain the mean deterministic process, that is z t , associated to the dynamical system. The only unknowns that remain are the initial law , the generator A and the state space E of the random process X t . The stationary law will be obtained from A through

A = 0, i∈E i = 1.

Estimating the paths and the state space of the jump Markov process

Prior to any estimation of A or , some representations of the paths (X k t ) k=1,...,K are needed. To this purpose, it is assumed that there exists a function G from expression (3) giving X t as a function of Z t and its derivative denoted by Żt = dZ t /dt. Such a function does not always exist. Indeed, when the stochastic process X t is a linear additive or a multiplicative term in the function C, we may easily find the corresponding function G. This is the case for the so-called Paris-Erdogan model, which will be used to describe the crack growth mechanism of the Virkler data. Hence, we may obtain a first estimation for the X k t 's with the following relationship:

X k t = G(Z k t , Żk t ), (11)
where the derivative of Z t can be estimated through

Żk t = Z k t+ t -Z k t t , (12)
with t being the time discretization step of the data set. Due to the way the X k t are computed, that is through expression (11), we obtain some noisy paths taking their values in R, which values may be quite nearby. Our model requires a finite state space for the underlying jump Markov process X t with some piecewise shape paths, thus it is appropriate to "regroup" the values which are very close from each other to an unique state. By these means, we build a state space which is an approximation of the real state space. Two cases may occur: first, if the state space E of the process X t is assumed to be known, each point of the X k t 's is rounded to the closest element of E, according to a given distance (e.g., euclidean). Second, if E is not known, which is the most common case, a non-hierarchical clustering method can be used to fit as well as possible with the X k t 's. In Chiquet et al. (2006), we applied the classical K-means algorithm to estimate the state space of X t . The application of one of this clustering methods leads to K approximated paths X k t , still defined on the random time interval [0, k ∧ T obs ]. The X k t are then used for further estimations linked to the jump Markov process X t , namely for estimating and A.

Maximum likelihood method

Concerning the initial law of the jump Markov process X t , it is assumed that X 0 is independent of the generator and of the censoring time . We thus have the classical empirical estimator ( see, e.g., [START_REF] Sadek | Nonparametric estimation of reliability and survival function for continuous time finite Markov processes[END_REF]

) i = 1 K K k=1 1 {X k 0 =i} . ( 13
)
It turns out that the estimation problem of the infinitesimal generator A requires further attention. In order to build the likelihood function for A on the X k t 's, we need to give an expression of the density function of a path of X t defined on a random time interval [0, ]. To this purpose, we first remind fundamental results of statistical inference for jump Markov processes, obtained by [START_REF] Billingsley | Statistical Inference for Markov Processes[END_REF] and [START_REF] Albert | Estimating the infinitesimal generator of a continuous time, finite state Markov process[END_REF], who built the density function on the space of sample paths defined on the same fixed time interval [0, ] with a real positive constant. The following notation is required in the sequel: • (S n , n ∈ N), the jump times of the process X t , with S 0 = 0; • W n = S n+1 -S n , for n ∈ N, the sojourn times in the successively visited states; • (J n , n ∈ N), the so-called embedded Markov chain of the Markov process X t , that is J n = X Sn , for n ∈ N;

• N( ), the number of jumps in [0, ], that is the largest integer for which S N( ) ;

Please cite this article as: Chiquet, J., et al., Piecewise deterministic Markov processes applied to fatigue crack growth modelling. J. Statist. Plann. Inference ( 2008), doi: 10.1016/j.jspi.2008.05.034

• N ij ( ) the number of transitions from state i to state j observed in [0, ], that is,

N ij ( ) = N( ) k=1 1 {J k-1 =i,J k =j} ;
• V i ( ) the total length of time spent in state i on [0, ], that is,

V i ( ) = 0 1 {Xs=i} ds.
Hence, the history H of the process X t on the fixed time interval [0, ] can be represented as an ordered sequence:

H = ((J 0 , W 0 ), ... , (W N( )-1 , W N( )-1 ), (J N( ) , U )),
where U = -S N( ) . It means that the path starts at state J 0 at time zero, remains in J 0 for W 0 units of time, makes a jump to J 1 , and so on, until it makes the final jump to J N( ) where it stays for U units of time, that is until the censoring time .

The probability density function f A associated to a sample path defined on a fixed time interval [0, ], was given by [START_REF] Albert | Estimating the infinitesimal generator of a continuous time, finite state Markov process[END_REF]. We put h = H ( ) a realization of the history of X t associated to a given sample path. It can be shown that

f A (h ) = (x 0 ) i,j∈E i =j a n ij ( ) ij e -a ij vi( ) .
The lower case letters n ij ( ), v i ( ) have been used to write down the realizations associated to the corresponding random variables. These results can be slightly generalized to get an expression of the density for a path defined on a random time interval [0, ], with the random censoring variable. Let us put H the history of X t when observed on the random interval [0, ]. A realization h = H ( ) ( ) of the history associated to a sample path of X t censored by has the same structure as in the fixed time case. Yet the density f A, of a randomly censored path h is written as a function of the conditional density f A| and of the density function f of through

f A, (h ) = f A| (h | = )f ( ). ( 14
)
Let us now build the likelihood function based upon the density f A, , yet under the observation scheme of Fig. 2, that is when the paths are not necessarily observed until the random failure time , but until ∧ T obs , that is still random. We have K paths of X k t estimated from Section 3.2, defined on [0, t k ], with t k being a realization of k ∧ T obs . We associate K ordered sequences h k t k to the estimated paths X k t 's. Let us also introduce the following function, associated to the kth path:

I k (t k ) = f (t k ) i f k T obs , P( > T obs ) if k > T obs .
By expression (14) of f A, and using the function I , it turns out that the corresponding log-likelihood L A for the generator A is

L A (K) = K k=1 log f A| (h k t k | k = t k ) + K k=1 log I k (t k ). ( 15
)
The second term in the right hand of (15) can be written

K k=1 log I k (t k ) = k:t k <T obs log f (t k ) + k:t k =T obs log P( > t k ),
and if we denote by L the number of paths censored by T obs , with L K, we have for these L paths P( > t k ) = 1 -F (T obs ). Hence, the log-likelihood (15) becomes

L A (K) = K k=1 log f A| (h k t k | k = t k ) + k:t k <T obs log f (t k ) + L log(1 -F (T obs )). (16)
The above expression ( 16) is quite general. Let us here treat the case where the failure time and the process X t are independent, which is obviously not true when dealing with the dynamical system (3), but can be considered as a first approximation. If so, the

JSPI : 3675 ARTICLE IN PRESS

structure of the likelihood is almost the same as in [START_REF] Albert | Estimating the infinitesimal generator of a continuous time, finite state Markov process[END_REF]. The conditional density f A| turns to f A and the log-likelihood L 0 A to be maximized is just

L 0 A (K) = K k=1 log f A (h k t k ) = K k=1 log i,j∈E i =j a n k ij (t k ) ij e -a ij v k i (t k ) ,
where n k ij and v k i are, respectively, the number of transitions from i to j and the time spent in i for the kth observed history h k . If we denote n ij (K) and v i (K) the corresponding variables for all of the K paths, that is

n ij (K) = K k=1 n k ij (t k ), v i (K) = K k=1 v k i (t k ),
then we have

L 0 A (K) = i,j∈E i =j n ij (K) log a ij -a ij v i (K).
For an element a ij with i = j, we get

a ij (K) = n ij (K) v i (K) . ( 17
)
This is a straightforward generalization of the maximum likelihood estimator obtained in [START_REF] Billingsley | Statistical Inference for Markov Processes[END_REF] and [START_REF] Albert | Estimating the infinitesimal generator of a continuous time, finite state Markov process[END_REF], computed here on several sample paths with random time lengths. In Chiquet and Limnios (2006), estimator ( 17) was proved to be consistent with an asymptotic normality property, even when defined on randomly censored sample paths, in the independent case.

In the following numerical application, we use the estimator for the generator of the underlying jump Markov process of the dynamical system, in the framework of fatigue crack growth analysis. This is obviously an approximation because we consider the failure time and the process X t independent, which is not the case as seen in expression (6) of the failure time. Although good numerical results are obtained with this approximation, it could be significantly improved by introducing the dependency in the estimation procedure.

Application to crack growth modelling

To handle with the fatigue crack growth problem, some empirical laws have been suggested in fracture mechanics to describe the mean behavior of a crack. One of the most versatile one is the law given by [START_REF] Paris | A critical analysis of crack propagation laws[END_REF], describing deterministically the evolution of the crack size function z t by

dz t dt = C( K(z t )) n , (18)
where C and n are some material constants and K(z) = Y(z) S √ z is the so-called stress intensity factor range. S is the stress range of the load applied on the specimen, while Y(z) is a structure geometry function that is approximately equal to 1 when considering large structure. With this assumption, Eq. ( 18) turns to the much simpler form

dz t dt = a(z t ) b , ( 19
)
where a = C( S) n n/2 and b = n/2. We suggest to randomized ( 19) with a right-hand side multiplicative factor given by a jump Markov process X t , thus we model the phenomenon by the following stochastic dynamical system

dZ t dt = a(Z t ) b × X t , Z 0 = z, ( 20
)
with z being the initial crack size.

Applying the Bogolyubov's averaging principle described in Section 3.1, we have

dZ t dt = a(Z t ) b × X t/ , Z 0 = z,
where > 0, a series parameter. The asymptotic behavior, when → 0, gives where a 0 = a i∈E i i. The following stochastic dynamical system is equivalent to (20)

d z t dt = a 0 ( z t ) b , z 0 = z, (21)
dZ t dt = a 0 (Z t ) b × v(X t ), Z 0 = z, ( 22
)
where v(X t ) is a jump Markov process with the same generator as X t , but a state space E that is linked to the state space E of X t through the linear mapping v defined by v(x) = x/ i∈E i i, x ∈ E. We consider system (22) rather than (20) so that we "normalize" the process X t to v(X t ) by dividing by the state mean i i i. The term a 0 absorbs the compensatory term for the normalization. Hence, practically, the stationary law does not need to be computed. Moreover, the system is more suitable to regression analysis. As a matter of fact, the parameters a 0 and b can be estimated by a least squares method on the data set, by taking the logarithm on both sides of (21). Hence, ln Żt = ln a 0 + b ln Z t + ln v(X t ).

The available data, that is the Z k t 's, can be represented as a N-sample composed with all of the data points from the crack growth curves:

{(t i , Z ti ) : i = 1, ... , N}.
This N-sample can be transformed in {(ln Z ti , ln Żti )}, with Żti being obtained through (12). Introducing the notation x i = ln Z ti , y i = ln Żti and i = ln v(X ti ), we have the following classical regression problem:

y i = ln a 0 + bx i + i , i = 1, ... , N. ( 23
)
When → 0, that is in the series scheme of the Bogolyubov's principle described in Section 3.1, system (22) tends to (21). It also means that, on average, v(X t ) is 1. The term "average" must be understood here in the series scheme sense of the Bogolyubov's principle. Thus, the ln v(X ti )'s are "close to 0" in the regression equation ( 23) and correspond to the so-called residuals i . Minimizing 2 i , we get the estimators associated to the well-known least squares method:

a 0 = exp(y -bx), b = x i y i -Nxy x 2 i -Nx 2 , with x = 1 N x i , y = 1 N y i .
Some paths ( v(X t )) k , k = 1, ... , K, are then obtained explicitly by "reversing" the dynamical system, that is through the function G given by (11). In the case at hand, we have

( v(X t )) k = 1 a 0 (Z k t ) -b × Żk t .
We then apply the K-means clustering method, initially suggested by MacQueen (1967) ( see, e.g. Jain and Dubes, 1988, for implementation) to reduce the state space of the ( v(X t )) k , thus obtaining some ( v(X t )) k 's suitable for any further estimation of the jump Markov process. For the Virkler data, the histogram of a 18-state space obtained with the K-means method is represented on Fig. 3. The values of the estimated state space of E are also put in the figure . One can see that it has a kind of discrete log-normal dispersion profile, with values distributed around one, what fits well the physical point of view: when v(X t )=1, we are in the mean behavior of the asymptotic Paris-Erdogan law (21), corresponding to the assumption made during the least squares regression. Finally, the initial law and the infinitesimal generator A are estimated with ( 13) and ( 17), adapted to the jump Markov process v(X t ). Some crack growth curves can thus be simulated with the dynamical system (22) as long as all of the model parameters have been estimated. With Monte Carlo techniques, some important functions can then be computed. For instance, the cumulative distribution function of Z t is interesting for industrial care, especially for optimizing a maintenance strategy. It represents the probability distribution of the crack size in time. We estimate it for a fixed instant of time t through In the same way, we may estimate the reliability defined here by (1) through

F t (x) = P(Z t x) = 1 M M m=1 1 { m t x} , x ∈ [z, ], (24)
R(t) = 1 M M m=1 1 { m t < } . ( 25
)
Results are compared with the data points on Fig. 5.

Concluding remarks and perspectives

Concerning the generator A of the jump Markov process X t , the estimator used in the previous numerical application is an approximation, since we have used the maximum likelihood estimator for K sample paths censored by a random variable independent of the process. We pointed out that it was not the case in the framework of dynamical system developed in the present paper, as long as the random censoring is also the failure time of the system. In Section 3, we gave the general form of the log-likelihood to be maximized. One can see that the density function of appears within (16), and an important issue is to give an formal expression of f if we want to include it in the maximization process. This could be achieved using Proposition

Fig. 1 .

 1 Fig. 1. The Virkler data set for the fatigue crack growth problem.

Fig. 3 .Fig. 4 .

 34 Fig. 3. 18-State space estimated for v(Xt ) through the K-means method.

  1 and consecutive results, what the authors are working on. Please cite this article as: Chiquet, J., et al., Piecewise deterministic Markov processes applied to fatigue crack growth modelling. J. Statist. Plann. Inference (2008), doi: 10.1016/j.jspi.2008.05.034

  Observation scheme of the dynamical system.

	point {Δ}	Z t 2	Z t 3	Z t 1
	State space of Z t
	z				Time
	State space E of X t	2 X t	X t 3	X t 1
		t 2	t 3	t 1 = T obs	Time
	Fig. 2.

Chiquet, J., et al., Piecewise deterministic Markov processes applied to fatigue crack growth modelling. J. Statist. Plann. Inference (2008), doi: 10.1016/j.jspi.2008.05.034

Please cite this article as: Chiquet, J., et al., Piecewise deterministic Markov processes applied to fatigue crack growth modelling. J. Statist. Plann. Inference (2008), doi: 10.1016/j.jspi.2008.05.034