Laurent Bordes
email: laurent.bordes@univ-pau.fr

Sophie Mercier
email: sophie.mercier@univ-pau.fr

Emmanuel Remy
email: emmanuel.remy@edf.fr

Emilie Dautrême
email: emilie.dautreme@edf.fr

Partially observed

Keywords: Competing risks, Composite likelihood, Gamma process, Poisson process, Reliability theory, Residual life

published or not. The documents may come

Introduction and motivating example 1.Introduction

Competing risks are largely used in survival or reliability analysis when several causes of death or failure are in competition. In a reliability study, the competing risk model can be seen as a series system, which fails whenever one of its component fails. To be more specic, let X 1 , . . . , X k be the respective lifetimes of a k-component system with competing risks. The lifetime of the system then is X = min{X 1 , . . . , X k } and the index C such that X = X C is known as the failure cause. The failure cause is generally observed together with the system lifetime, or it may be partially observed, with C ∈ C ⊂ {1, . . . , k}. The later case is known as a masked cause of failure. In order to prevent identiability issues that may occur when lifetimes X j 's are dependent, the lifetimes often are considered as independent. Also, the number of competing causes is generally xed, see Crowder [START_REF] Crowder | Classical Competing Risks[END_REF] for an overview on classical competing risks models and their application. In presence of dependency (see [START_REF] Somboonsavatdee | Parametric inference for multiple repairable systems under dependent competing risks[END_REF] for an example where the dependence is due to the sharing of a common frailty random variable by the competing risks), it is generally the cause specic hazard rates that are estimated. The number of competing causes may also sometimes be random, see [START_REF] Balakrishnan | COM-Poisson Cure Rate Models and Associated Likelihood-based Inference with Exponential and Weibull Lifetimes[END_REF] for an example in biostatistics. It seems however that the case of a random number of competing causes with non identically distributed component lifetimes has not been much addressed in the previous survival or reliability analysis literature.

In this paper, we are interested in a system subject to several competing deterioration sources (defects), which appear one at a time and next independently propagate. These defects could be for instance pits due to corrosion [START_REF] Kuniewski | Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection[END_REF], fatigue cracks [START_REF] Guida | A gamma process model for the analysis of fatigue crack growth data[END_REF], or any measurable indications of deterioration which initiate over time and next propagate. In the present applicative context, these defects correspond to aw indications, which are not further specied, due to condentiality issues. However we use this vocabulary in the remaining of the paper, and we will speak of "aw indication" or "indication" for short. The competing deterioration sources (the aw indications) initiate at random times T 1 < T 2 < • • • , which are the points of a counting process N t = ∞ i=1 1 {T i ≤t} t≥0 . Once initiated at time T j , the jth indication propagates according to a non decreasing process Z (j) = (Z (j) t ) t≥0 , which is set equal to zero on [0, T j ). The system is considered to be out of order as soon as the value of one among all the initiated aw indications has reached a critical threshold, known in advance. This leads to a competition between the N t indications present at time t. Then we call competing degradation process the bivariate process (N t , Z t ) t≥0 where Z t = max{Z [START_REF] Abdel-Hammed | A gamma wear process[END_REF] t , . . . , Z (Nt) t }. This denomination comes from the fact that the crossing time of the degradation level by (Z t ) t≥0 is nothing but the minimum of the crossing times of processes Z (1) , . . . , Z (Nt) for the same level. This connection between multiple degradation processes and competing risks has been used to derive new lifetime models. In [START_REF] Singpurwalla | On competing risk and degradation processes[END_REF] some examples are provided by using Brownian Maximum Processes or gamma processes for the degradation process, whereas in [START_REF] Haghighi | Parametric degradation model with multiple competing risks[END_REF] the authors use degradation processes of type Z (j) t = g(t; A j ) where g is a known function and A i a random vector.

Though the propagation is assumed to be similar for each of the aw indications, their initiation times are dierent, leading to non identically distributed Z (j) 's and consequently, to a competition between a random number of non identical degradation processes. A similar approach has been proposed in Kuniewski et al. [START_REF] Kuniewski | Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection[END_REF], in order to evaluate damages caused by corrosion of industrial systems in the oil industry. Their model combines a non homogeneous spatial Poisson process for the initiation process (N t ) t≥0 with a homogeneous gamma process for the propagation of the defect. The parameters of both Poisson and gamma processes are assumed to be known, up to a proportionality constant for the Poisson process. The paper develops a Bayesian method for estimating the proportionality constant and predictive indicators, based on partial information data including imperfect defects detection. In [START_REF] Caballé | A conditionbased maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes[END_REF], Caballé et al. propose a condition-based maintenance of a system subject to fatal shocks and to deterioration due to pitting corrosion, with a similar model as [START_REF] Kuniewski | Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection[END_REF] for the latter. In the same way, Cha and Castro envision sudden failures in competition with soft failures due to a similar deterioration model in [START_REF] Cha | A stochastic failure model with dependent competing risks and its applications to condition-based maintenance[END_REF], where they study some positive dependence properties between the competing failure times and some condition-based maintenance policy. A very close deterioration model is also used in [START_REF] Castro | Age-based preventive maintenance for passive components submitted to stress corrosion cracking[END_REF][START_REF] Velázquez | Statistical modelling of pitting corrosion: Extrapolation of the maximum pit depth-growth[END_REF] (among others).

The paper is organized as follows. First, we introduce in the next sub-section the industrial example which motivated the study. The probabilistic model is described in Section 2, where the conditional probability density function of Z t given N t is calculated, together with additional prognostic indicators of interest for the application. All the calculations are relegated to Appendices AC. In Section 3, a parametric estimation procedure is proposed, which accounts for the incomplete information data coming from the eld, as described in the next sub-section. Section 4 is devoted to numerical developments, including a validation of the estimation method through a Monte Carlo study as well as the application to the real data set. Some concluding remarks are given in the last section.

1.2

Motivating industrial application EDF, the world's biggest electricity generator, performs in-service inspections of the passive components within its electric power plants in order to ensure that their degradation is lower than a critical level and guarantee the safety and the availability of the installations. These examinations are regular but non periodic, and their times dier from one component to the other. They allow to collect successive degradation measurements for each component, from which the point is to predict the degradation propagation and estimate the residual operation time upon which the critical degradation level (denoted by ) is reached.

By an examination, the measurement process does not always give a perfect image of the degradation. Indeed, because of technical limitations, too small aw indications are detected but their size cannot be measured. Moreover, if several competing aw indications have initiated on one component, the measurement process can count the number of existing aw indications but only the size of the largest one can be measured. For a better understanding of the available information, let us now be more specic on the examination procedure, which is made through two non-destructive ultrasonic testing processes with dierent objectives and performances:

• The rst process aims at detecting aw indications: it gives a binary response "presence"

or "absence" of aw indications, together with the number of initiated indications.

• The second one is conceived to measure aw indications. The measurement rst requires the detection of an indication by this second process, which is less sensitive than the rst one and only detects aw indications whose measurement exceeds a known xed threshold. Also, due to the use of two dierent technologies for this second process, there are several (always known) thresholds under which the aw indication cannot be measured. Finally, when several aw indications are observable, only the largest one is measured by the second testing process.

We assume that all aw indications are detected by the rst testing process and that the second testing process gives the exact measurement of the largest indication when observable.

The available data concern 228 components that are supposed to be independent and identical.

The commissioning date is known for each component as well as its inspection times. Each previously described two-step examination provides one of the possible following information about the degradation level:

• No aw indication is detected on the component: no indication is initiated at the inspection time.

• At least one aw indication is detected by the rst testing process at the inspection time but no aw indication is detected by the second testing process: the number of initiated indications is known and the measurement of the largest aw is lower than a deterministic known value (left censoring).

• The second testing process detects and measures aw indications: the number of initiated indications is known, as well as the measurement of the largest indication on the component (which is therefore higher than a known value).

In order to make accurate degradation predictions, it is essential to take into account this partial information coming from the eld. It is important to mention that no other source of information is available to predict the evolution of the indications with time. In particular, the underlying physical phenomena are complex and no representative mechanical model exists: a structural reliability approach [START_REF] Lemaire | Structural Reliability[END_REF] is thus not feasible. The measurements from the inspections are the only available material and a statistical analysis of these data is the unique way to assess the residual lifetime of the components.

The object of the present study is to compute predictive indicators such as the yearly expected number of initiated aw indications, the expected annual growth of one aw indication, quantiles of the residual operating time given some observations at a given time. Another applicative object is the study of the inuence of the number of initiated aw indications on the residual operating time, in order to know whether it is worth or not to gather this information from the eld. A last point of interest from an industrial point of view is to know whether it is worth taking into account censoring. To be more specic, a rst possibility for estimating the model parameters naturally is to take into account the censoring nature of the data. A second possibility is to use the so-called imputation method deployed within survival analysis studies. In this method, the censoring level of any censored data is considered to be the data itself. Here, considering a largest aw indication which is known to be smaller than a censoring level C, the imputation method consists in considering that the size of the largest aw indication is exactly equal to C. (And this is done for all the censored data). Of course, this imputation method leads to an over-estimation of the real sizes of the (largest) aw indications and hence introduces some bias in the estimated parameters. However, it may be of interest in an industrial context because it leads to simpler estimation procedures. More important, the method has the good property of being conservative, in the sense that the results over-estimate the degradation propagation and hence provide pessimistic prognostic indicators.

(Such a method is particularly adapted as soon as safety is concerned). The point here is to know whether the imputation method provides very dierent results from an estimation procedure taking into account the censoring nature of the data and hence whether it is worth taking into account censoring. The answer to this question requires the development of the two dierent estimation procedures, for comparison purpose.

Z t = max 1≤j≤Nt Z (j) t if N t ≥ 1, 0 elsewhere. If 0 ≤ t 1 < • • • < t m
are the inspection times of the component (dierent from one component to another and non informative), data for one component are (n j , z j ) 1≤j≤m with n j = N t j (ω) and z j = Z t j (ω) for an experiment ω. These data correspond to partial observation N t j , Z t j 1≤j≤m of the competing degradation process (N t , Z t ) t≥0 .

In order to account several possible levels of censoring at time t, we set C t to be the known deterministic censoring level at time t. We also introduce a censoring indicator:

D t = 0 if Z t ≤ C t , 1 otherwise.
The data now are of the form (n j , u j , d j ) 1≤j≤m with n j = N t j (ω), u j = U t j (ω), d j = D t j (ω), and correspond to one observation of N t j , U t j , D t j 1≤j≤m with

U t =    Z t if Z t > C t , C t if Z t ≤ C t and N t > 0, 0 otherwise.
An example of censored trajectory is provided in Figure 1 for the competing degradation process. The censoring level is assumed to be constant for sake of clarity. Four aw indications initiate at time T 1 (blue), T 2 (red), T 3 (green) and T 4 (black), respectively. At the beginning, the largest indication (namely Z t ) is the blue one, next the red one and nally the green one. This is depicted through black/white dots, where black dots means that Z t is below the censoring level (so that U t = C t ) and white dots refers to an observable Z t beyond the censoring level (so that U t = Z t ). The observable counting process (N t ) t≥0 is simultaneously increased by one at each initiation time T i , i = 1, . . . , 4.

Assumptions. We assume that the degradation processes are identical for all the components and that they are independent between components. The independence is justied by the fact that the components are in dierent production units throughout France so that they do not interfere one with the other. We also assume that once initiated, the processes dening the propagation of all aw indications present on one component are independent and identically distributed (i.i.d.). Indeed there are only few aw indications at the same time on a component (otherwise the component is changed) and the size of a component is large (thus limiting some line

Z (4) t t Z (1) t Z (2) t Z (3) t N t Z t T 1 T 2 T 3 T 4 × × × × censoring level observed Z t censored Z t 1 2 3 4 Figure 1:
An example of trajectory for the competing degradation process potential mechanical interactions between aw indications), so that it seems mostly natural to assume that the development of the dierent aw indications are similar and independent one from the other. Note that this i.i.d. assumption is pretty classical when dealing with similar problems, see for instance [START_REF] Caballé | A conditionbased maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes[END_REF], [START_REF] Castro | Age-based preventive maintenance for passive components submitted to stress corrosion cracking[END_REF], [START_REF] Castro | A condition-based maintenance for a system subject to multiple degradation processes and external shocks[END_REF], [START_REF] Kuniewski | Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection[END_REF], [START_REF] Scarf | A stochastic model of crack growth under periodic inspections[END_REF] or [START_REF] Velázquez | Statistical modelling of pitting corrosion: Extrapolation of the maximum pit depth-growth[END_REF] among lots of others. For each component, the number of aw indications is assumed to follow a Poisson process (N t ) t≥0 , where N t is Poisson distributed with parameter Λ(t) and distribution denoted by P(Λ(t)), see [START_REF] Ross | Stochastic processes[END_REF] for more details on Poisson processes. In the sequel, we either suppose that this process is homogeneous (HP) with Λ(t) = λt, or non homogeneous (NHP) with Λ(t) = αt β (powerlaw process), where λ, α, β > 0. The rate of this Poisson process is denoted by λ(t), with λ(t) = λ for the HP process and λ(t) = Λ (t) = αβt β-1 for the NHP process. When a aw indication appears, we assume that its measurement increases according to a gamma process with parameters (A(t), b), e.g. see [START_REF] Abdel-Hammed | A gamma wear process[END_REF][START_REF] Kahle | Degradation Processes in Reliability[END_REF][START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF][START_REF] Ye | Stochastic modelling and analysis of degradation for highly reliable products[END_REF] for more details on the use of gamma processes for modeling deterioration. Here again, we envision the homogeneous case (HG) with A(t) = at and the non homogeneous case (NHG) with A(t) = ηt γ , where a, η, γ > 0.

To specically dene the measurements of the various indications present on a component at time t, we introduce a sequence X (j) j∈N * of independent gamma processes, all with param- eters (A(t), b). For each j ∈ N * and each t ≥ 0, the random variable X where

(j) t is gamma distributed with probability density function (p.d.f.) f A(t),b (x) = b A(t) x A(t)-1 e -bx 1 R + (x) /γ (A(t)), cumu- lative distribution function (c.d.f.) denoted by F A(t)
FA(t),b = 1 -F A(t),b . With the previous notations, we have E X (j) t = A(t)/b and var X (j) t = A(t)/b 2 .
For each 1 ≤ j ≤ N t , the dierence t -T j corresponds to the elapsed time at time t since the initiation of the j-th indication. We then set Z

(j) t = X (j) t-T j
to be the measurement of the j-th aw indication at time t. If j > N t , the j-th indication has not yet been initiated at time t and we put Z (j) t = 0. Setting (t -T j ) + = max (t -T j , 0) and remembering that X

(j) 0 = 0, this can be summed up into Z (j) t = X (j) (t-T j ) + for all j ∈ N * .
As only the largest aw indication can be measured by an examination, without censoring, a measure at time t hence is an observation of

Z t = max j∈N * Z (j) t = max 1≤j≤Nt X (j) t-T j if N t ≥ 1, 0 otherwise.
We now derive some probabilistic results, to be used later on for inference purpose and computation of the predictive indicators.

Probabilistic results

Let us rst remark that the maximum among the Z (j) t 's will generally not correspond to the same trajectory over time. A direct consequence is that the distribution of an increment Z t -Z t 0 (with t 0 < t) is not an easy thing to write down. Also, the process (Z t ) t≥0 does not have independent increments. The obtaining of probabilistic results hence requires some care, especially for the derivation of the conditional distribution of (Z t , N t ) given (Z t 0 , N t 0 ) given in Proposition 2.

We rst provide the conditional distribution of Z t when the number of aw indications is known at time t.

Proposition 1 Given that N t = 0, then Z t is conditionally almost surely null (namely P(Z t = 0|N t = 0) = 1). Otherwise, for each n ≥ 1, the conditional distribution of Z t given that N t = n has the following cumulative distribution function (c.d.f.)

F Zt|Nt (z|n) = t 0 F A(y),b (z) λ (t -y) dy Λ (t) n , (1)
and its probability distribution function (p.d.f.) with respect to the Lebesgue measure is given by

f Zt|Nt (z|n) = n (Λ (t)) n t 0 f A(y),b (z) λ (t -y) dy t 0 F A(y),b (z) λ (t -y) dy n-1 (2)
for all z ≥ 0.

Proof. See Appendix A.

The joint distribution of (Z t , N t ) is easy to derive as the product of the conditional distribution of Z t given N t cross the distribution of N t (see also Equation (B.3) in Appendix B). It admits the following density function

f (Nt,Zt) (n, z) = e -Λ(t) if n = 0, e -Λ(t) (n-1)! t 0 f A(y),b (z) λ (t -y) dy t 0 F A(y),b (z) λ (t -y) dy n-1 if n > 0 (3)
with respect to

δ 0 (dn) × δ 0 (dz) + +∞ k=1 δ k (dn) × dz.
We now provide conditional distribution results of the competing degradation process at a time t > 0 given its value at a previous time t 0 ∈ [0, t).

Proposition 2 For t 0 ∈ [0, t) and z ≥ 0 we have

P (Z t ≤ z|Z t 0 = 0) = e -(Λ(t)-Λ(t 0 ))+ t t 0 F A(t-v),b (z)λ(v)dv = e -t t 0 FA(t-v),b (z)λ(v)dv . ( 4
)
For t 0 ∈ (0, t), 0 < z 0 ≤ z and n ≥ n 0 ≥ 1, we have

P (Z t ≤ z, N t = n|Z t 0 = z 0 , N t 0 = n 0 ) = e -(Λ(t)-Λ(t 0 )) × t 0 0 F A(t-u)-A(t 0 -u),b (z -z 0 ) f A(t 0 -u),b (z 0 ) λ (u) du t 0 0 f A(t 0 -y),b (z 0 ) λ (y) dy × t 0 0 z 0 0 F A(t-u)-A(t 0 -u),b (z -x) f A(t 0 -u),b (x) dx λ (u) du t 0 0 F A(t 0 -y),b (z 0 ) λ (y) dy n 0 -1 × t t 0 F A(t-v),b(z) λ (v) dv n-n 0 (n -n 0 )! . (5)
Proof. See Appendix B.

Note that a similar expression to (4) is provided by Equation ( 18) in [START_REF] Kuniewski | Sampling inspection for the evaluation of time-dependent reliability of deteriorating systems under imperfect defect detection[END_REF] in the specic case of a homogeneous gamma process with t 0 = 0.

Writing P (Z t ≤ z|Z t 0 = z 0 , N t 0 = n 0 ) = +∞ n=n 0 P (Z t ≤ z, N t = n|Z t 0 = z 0 , N t 0 = n 0 ) ,
one immediately gets the following corollary.

Corollary 1 For 0 < t 0 < t, 0 < z 0 ≤ z and n 0 ≥ 1, we have

P (Z t ≤ z|Z t 0 = z 0 , N t 0 = n 0 ) = e -t t 0 FA(t-v),b (z) λ (v) dv × t 0 0 F A(t-u)-A(t 0 -u),b (z -z 0 ) f A(t 0 -u),b (z 0 ) λ (u) du t 0 0 f A(t 0 -y),b (z 0 ) λ (y) dy × t 0 0 z 0 0 F A(t-u)-A(t 0 -u),b (z -x) f A(t 0 -u),b (x) dx λ (u) du t 0 0 F A(t 0 -y),b (z 0 ) λ (y) dy n 0 -1
.

If necessary, one could also derive from the previous results the conditional distribution of Z t given Z t 0 (and next the distribution of an increment Z t -Z t 0 ) by writing

P (Z t ≤ z|Z t 0 = z 0 ) = ∞ n 0 =1 P(Z t ≤ z|Z t 0 = z 0 , N t 0 = n 0 ) P(Z t 0 = z 0 , N t 0 = n 0 ) P(Z t 0 = z 0 )
for all z 0 > 0 and using the joint distribution of (Z t 0 , N t 0 ) provided by Equation (B.3) in Appendix B. This is not used in the following and consequently not made explicit here.

Prognostic indicators

The initiation time U of the rst aw indication veries P(

U > t) = P(N (t) = 0) = e -Λ(t) = e -αt β
. It is therefore distributed according to a Weibull distribution with scale parameter α -1/β and shape parameter β. Mean, variance and quantiles of the distribution of U are given in Table 1.

ε µ c = α -1/β γ (1 + 1/β) σ c 2 = α -2/β γ (1 + 2/β) -µ c 2 q c,ε = α -1/β (-log (1 -ε)) 1/β
In our model, the propagation of a aw indication follows a gamma process with parameters (A(t), b) and the expected propagation between times s and t (for 0 < s < t < ∞) is (A(t)-A(s))/b. For a homogeneous gamma process with parameters (a, b), the yearly expected propagation is 365 a/b since the time unit is the day.

Remembering that stands for the critical degradation level of the system, we set τ to be the crossing time of level , with

τ = inf{t ≥ 0 : Z t > }.
Given t 0 ≥ 0 and an observation (Z t 0 , N t 0 ) = (z 0 , n 0 ) at time t 0 with z 0 ∈ [0, ) (so that τ > t 0 ), the residual life of a component at time t 0 is identically distributed as the conditional distribution of τ -t 0 given (Z t 0 , N t 0 ) = (z 0 , n 0 ). Given ε in (0, 1), a prognostic indicator of major interest for the industrial application is the quantile t ( ,t 0 ,z 0 ,n 0 ) ε of order ε of the residual life at time t 0 given that Z t 0 = z 0 and N t 0 = n 0 . For z 0 = 0, we necessarily have n 0 = 0 so that we write t

( ,t 0 ,0) ε instead of t ( ,t 0 ,0,0) ε .
The quantile t

( ,t 0 ,z 0 ,n 0 ) ε veries P τ -t 0 < t ( ,t 0 ,z 0 ,n 0 ) ε |Z t 0 = z 0 , N t 0 = n 0 = ε, or equivalently P Z t ( ,t 0 ,z 0 ,n 0 ) ε +t 0 ≤ |Z t 0 = z 0 , N t 0 = n 0 = 1 -ε,
using the fact that {τ ≥ t} = {Z t ≤ }. Based on Equation (4) or Corollary 1 according to whether z 0 = 0 or z 0 > 0, respectively, it is now easy to write down the equation fullled by the quantile t ( ,t 0 ,z 0 ,n 0 ) ε . As will be seen later on, the model that we nally retain for the application is the combination of the NHP process with the HG process (denoted by NHPHG model) so that the equation fullled by the quantile t ( ,t 0 ,z 0 ,n 0 ) ε is given only under this assumption.

Proposition 3 Assume that (N t , Z t ) t≥0 follows the NHPHG model. For z 0 = 0 the quantile t

( ,t 0 ,0) ε is obtained by solving t ( ,t 0 ,0) ε 0 Fay,b ( ) t (t 0 ,0) ε + t 0 -y β-1 dy = - log (1 -ε) αβ . (6)
For z 0 ∈ (0, ) the quantile t

( ,t 0 ,z 0 ,n 0 ) ε is obtained by solving log F a t ( ,t 0 ,z 0 ,n 0 ) ε ,b ( -z 0 ) -αβ t ( ,t 0 ,z 0 ,n 0 ) ε 0 (u + t 0 ) β-1 Fa t ( ,t 0 ,z 0 ,n 0 ) ε -u ,b ( ) du + (n 0 -1) × log t 0 0 u β-1 z 0 0 F a t ( ,t 0 ,z 0 ,n 0 ) ε ,b ( -x) f a(t 0 -u),b (x) dx du = log (1 -ε) + (n 0 -1) log t 0 0 F a(t 0 -y),b (z 0 ) y β-1 dy . (7)
Proof. See Appendix C.

3 Estimation procedure

Estimation principle

Remember that our data are independent observations of random vectors N t j , Z t j 1≤j≤m , pos- sibly censored, where the number of inspections m and the inspection times t j are dierent from one component to another. As the increments of (N t , Z t ) t≥0 are not independent, it is a real challenge to obtain the joint distribution of N t j , Z t j 1≤j≤m . It is however possible to write down a recursive formula, which unfortunately is of no real help for the computation of the joint distribution, due to its complexity. The usual maximum likelihood method cannot hence be used. Using the joint density function of (N t j , Z t j ) provided by [START_REF] Caballé | A conditionbased maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes[END_REF], it is possible to write down a composite likelihood function (whose formal denition has been introduced by Lindsay [START_REF] Lindsay | Composite likelihood methods[END_REF]) based on the product of the likelihood of each observation with

L θ| (n j , z j ) 1≤j≤m = 1≤j≤m f (Nt j ,Zt j ) (n j , z j ; θ) ,
where θ stands for the set of parameters to be estimated. In that case, the optimization is made with respect to all parameters in one single step (with ve parameters for the NHP-NHG model for instance). This method has been tested on simulated data. However, the results

were not very stable and the following two-step procedure has been preferred for its numerical robustness.

Let us recall that the parameters to be estimated are:

• Parameter(s) of the Poisson process: θ P = λ (HP) or θ P = (α, β) (NHP),

• Parameters of the gamma process:

θ G = (a, b) (HG) or θ G = (η, γ, b) (NHG).
As a rst step, we classically estimate the parameter θ P of the Poisson process by the usual maximum likelihood method, based on the uncensored observation of N t j 1≤j≤m and on the joint distribution of N t j 1≤j≤m . The estimator is denoted by θP .

As a second step, the parameters of the gamma processes are next estimated, plugging the estimator of the Poisson parameter into some composite likelihood function (see [START_REF] Cox | A note on pseudolikelihood constructed from marginal densities[END_REF] and [START_REF] Varin | An overview of composite likelihood methods[END_REF]), based on the conditional distribution of the Z t j 's given N t j 's. More specically, we consider the following composite likelihood:

L θ G | (z j ) 1≤j≤m , (n j ) 1≤j≤m ; θP = 1≤j≤m s.t. n j ≥1 f Zt j |Nt j z j |n j ; θ G , θP ,
where f Zt j |Nt j z j |n j ; θ G , θP is the conditional p.d.f. of Z t j given that N t j = n j with respect to Lebesgue measure for n j ≥ 1 (see Proposition 1), and where θ P is replaced by θP .

Observations for which n j = 0 do not contain any information on θ G because N t j = 0 implies that Z t j = 0. That is why only data such that n j ≥ 1 are involved in the likelihood function.

In practice, N independent components are observed. Adding exponent (i) to both processes and observations related to the i-th component, the log-composite-likelihood can be written as

θ G |z, n, t; θP = N i=1 1≤j≤m (i) s.t. n (i) j ≥1 log f Z t (i) j |N t (i) j z (i) j |n (i) j ; θ G , θP where z = z (i) j 1≤j≤m (i) 1≤i≤N , n = n (i) j 1≤j≤m (i) 1≤i≤N and t = t (i) j 1≤j≤m (i) 1≤i≤N
.

If in addition we take into account censoring, we get:

θ G |u, n, t, d; θP (8) = N i=1 1≤j≤m (i) s.t. n (i) j ≥1 d (i) j log f Z t (i) j |N t (i) j u (i) j |n (i) j ; θ G , θP + 1 -d (i) j log F Z t (i) j |N t (i) j u (i) j |n (i) j ; θ G , θP with u (i) j = z (i) j if d (i) j = 1 and u (i) j = c (i)
j otherwise, and d = d

(i) j 1≤j≤m (i) 1≤i≤N
.

Fitting the Poisson process

Parameter θ P of the Poisson process can be estimated from data (n, t) by maximizing the log-likelihood function:

(θ P |t, n) = N i=1 log P ∩ m (i) j=1 N t (i) j = n (i) j
.

Using the independent increments of (N t ) t≥0 , we easily get:

log P ∩ m j=1 N t j = n j ∝ -Λ (t m ) + m-1 j=0 (n j+1 -n j ) log (Λ (t j+1 ) -Λ (t j )) for all 0 = t 0 < t 1 < • • • < t m and 0 = n 0 ≤ n 1 ≤ • • • ≤ n m ,
where ∝ means "is equal to" up to an additive constant independent of the parameter of interest (here θ P ).

The log-likelihood can then be written as

(θ P |t, n) ∝ - N i=1 Λ t (i) m (i) + N i=1 m (i) -1 j=0 n (i) j+1 -n (i) j log Λ t (i) j+1 -Λ t (i) j
, where we set t

(i) 0 = n (i) 0 = 0, for all 1 ≤ i ≤ N .
In case of a non homogeneous Poisson process, we have Λ (t) = αt β and θ P = (α, β). Fixing β and solving ∂L ∂α (α, β|t, n) = 0 leads to

α(β) = N i=1 n (i) m (i) N i=1 t (i) m β . (9)
We replace α by α(β) in the log-likelihood function, and then we look for the maximizer β of

(β|t, n) ≡ (α(β), β|t, n) (10)
∝ -

N i=1 n (i) m (i) log N i=1 t (i) m (i) β + N i=1 m (i) -1 j=0 n (i) j+1 -n (i) j log t (i) j+1 β -t (i) j β .
Using a numerical optimization method, we rst obtain β; then we set α = α( β).

Fitting the gamma process

The parameter θ G of the gamma process is estimated from (u, n, t, d) by maximizing the logcomposite-likelihood function given by [START_REF] Crowder | Classical Competing Risks[END_REF]. Setting λ(t) = λ(t; θP ) and plugging formulae (1) and (2) into (8), we obtain

θ G |u, n, t, d; θP ∝ N i=1 1≤j≤m (i) s.t. n (i) j ≥1 d (i) j log t (i) j 0 f A(y),b u (i) j λ t (i) j -y dy + n (i) j -d (i) j log t (i) j 0 F A(y),b u (i) j λ t (i) j -y dy . (11)
In case of a homogeneous gamma process, A(y) is replaced by ay in [START_REF] Haghighi | Parametric degradation model with multiple competing risks[END_REF]. Parameters (a, b) are estimated by (â, b) which maximize a, b|u, n, t, d; θP . In case of a non homogeneous gamma process, A(y) is replaced by ηy γ in [START_REF] Haghighi | Parametric degradation model with multiple competing risks[END_REF]. Parameters (η, γ, b) are estimated by η, γ, b which maximizes η, γ, b|u, n, t, d; θP .

Estimating the prognostic indicators

Prognostic indicators of Section 2.3 are estimated by substituting the unknown parameters of both the initiation and degradation process by their estimated values. For instance, the estimator t( ,t 0 ,z 0 ,n 0 )

ε (resp. t( ,t 0 ,0) ε ) of t ( ,t 0 ,z 0 ,n 0 ) ε (resp. t ( ,t 0 ,0) ε
) is obtained by solving (7) (resp. ( 6))

replacing parameters (a, b, α, β) by (â, b, α, β).

4 Numerical study

Monte Carlo simulations

This section is devoted to a simulation study which beyond checking the validity of our estimation method aims at quantifying the inuence of the censoring process on the estimators behaviour. The next section deals with the application to the EDF data. As already mentioned, EDF data are well tted by combining the NHP process with the HG process. That is why the present section is limited to this specic NHPHG model although the estimation method has also been validated by a similar Monte Carlo study for the three other models.

In the EDF data, we have at our disposal 228 time data sets, which each corresponds to the This procedure provides us with a simulated censored data set based on the observation of N components, which is next used to derive an estimation of each parameter. This procedure is independently repeated 500 times, and condence intervals (CI) are nally constructed from empirical quantiles based on these 500 estimations. All the computations of this section have been carried out with MATLAB [16] on a dual core laptop computer. The optimization is made with the function fmincon through an interior-point algorithm and the BFGS method (please see MATLAB documentation on the fmincon function for more details and references).

As announced in the introduction, two dierent estimation methods are tested, rst through the imputation method (namely considering censoring levels as observed measurements) and next taking into account the censoring.

A rst set of parameters is considered for both methods: parameters of the NHP process Estimation without taking into account censoring (imputation method). We use the log-composite-likelihood function [START_REF] Haghighi | Parametric degradation model with multiple competing risks[END_REF] with all d (i) j 's equal to 1. Results are provided in Table 2. As for the parameters (α, β) of the NHP process, we observe that they are well estimated, both from a bias and standard deviation point of view. Whatever the sample size is, the condence intervals always contain the true values of the parameters. This good behaviour is consistent with the fact that the imputation method does not aect the number of indications at observation times and consequently does not aect the estimation procedure of the Poisson process parameters.

The results concerning the HG process are not convincing. Indeed, even if empirical means are not too far from their true values, we can see that there is some bias on the estimates of parameters a and b. As a matter of fact, condence intervals are either very large (for small N ) or do not even contain the true values of the parameters (for large N ). Though the bias can be seen to be less important on the mean rates per unit time of both mean (a/b) and variance (a/b

2 ) than on parameters a and b, it is clear that the imputation method leads to an overestimation of the sizes of the aw indications. It consequently does not seem to be appropriate for the present study and it is not considered any more in the following. Note however that for all parameters, standard deviations are divided by two when the sample size is multiplied by four, which is an indicator of a "rootofN " convergence rate with an asymptotic bias. 4 and5 As a conclusion to this simulation study, it is clearly preferable to take censoring into consideration and doing so, the estimation procedure seems to provide reliable results (which can be more or less long to obtain according to the parameters).

Application to EDF data

EDF data concern N = 228 components. For condentiality reasons the unit of the aw indication measure is not specied. Each component has been inspected several times. The total number of inspection times is 1,695. Our approach requires both numerical integration and numerical optimization for the computation and optimization of the log-composite-likelihood function provided by [START_REF] Haghighi | Parametric degradation model with multiple competing risks[END_REF]. To validate the estimation results, the MATLAB program used in the previous simulation study is used together with another one written in R [START_REF] Prabhu | R: A language and environment for statistical computing[END_REF]. The two programs provide very similar results. In addition, condence intervals are computed via standard bootstrap method, using 1,000 sets of 228 trajectories uniformly drawn (with replacement) from the 228 trajectories of our data set. For each bootstrap sample, estimates are provided for the parameters. Empirical mean, standard deviation and both 90% and 95% condence intervals are next derived for each parameter, based on the 1,000 estimation results. The 90% (resp. 95%) condence interval corresponds to [q 0.05 , q 0.95 ] (resp. [q 0.025 , q 0.975 ]), where q α is the empirical α-quantile based on the 1,000 estimation results. As often, this method is quite time consuming. Note that, based on the strong evidence provided by the previous Monte Carlo study for taking into account censoring, all the results of this section are computed under this basis.

Choosing between HP or NHP for the initiation process. Estimation results concerning the Poisson process (HP and NHP) are given in Table 6 (top) together with nonparametric bootstrap results. Note that whatever the model is, the estimates are within the bootstrapped condence interval. For each model, Table 6 (bottom) next gives the yearly estimated average number of new indications. For the HP process, it is equal to λ × 365 whereas for the NHP process it is equal to α × (365 k) β -(365 (k -1)) β for the k-th year (remember that the time unit is the day). In Table 6 (top), we can see that the 90% and 95% bootstrap condences intervals for β do not contain the value 1, advocating for the selection of the NHP for the initiation process. Also, EDF experts consider that the rate of initiation of new indications should be increasing over time, leading here again to the choice of the NHP based on Table 6 (bottom). Choosing between HG and NHG for the propagation process. Estimation results are given in Table 7 (left) for the parameters of both HG and NHG degradation processes, considering both HP and NHP initiation processes. Table 7 (right) provides the annual mean growth of one aw indication considering the four models. As we can see, whatever the initiation process is (HP or NHP), in case of a NHG propagation process (with A(t) = ηt γ ), the parameter γ is very small (γ < 10 -8 ) resulting in a strongly non-linear growth for one aw indication, which essentially expands during the rst year and next remains mostly constant. Based on EDF specialists considerations, we hence suggest to retain the homogeneous version of the gamma process (combined with an NHP process). However, this point might require further investigation. For example, a similar bootstrap method as for the initiation process model (HP vs. NHP) might be used, providing a bootstrap condence interval for the unknown parameter γ. This would allow to test the null hypothesis γ = 1 (meaning that the gamma process is homogeneous). Unfortunately, we have faced the problem of huge calculation times and numerical instabilities. Thus, this approach could not be achieved. Estimation results for the selected NHPHG model. Estimation results concerning the selected NHPHG model are available in Table 8. The parameters are reestimated using nonparametric bootstrap, which allows to obtain condence intervals as well as standard deviations for the dierent parameters. Estimation results for the parameter b of the homogeneous gamma process are rather stable as well as for the parameter β of the initiation process. The results are less stable for the parameters α (initiation process) and a (homogeneous gamma process). (90, 0, 0, 0) (90, 20, 0, 0) (90, 25, 30, 3) ε of the residual life τ 90 -t 0 (in years) with respect of ε ∈ (0, 1) for three dierent values of (t 0 , z 0 , n 0 ) 3: Estimated εquantiles of the residual life τ 90 -t 0 (in years) given z 0 = n 0 = 0 at time t 0 > 0 with respect of t 0 for ε ∈ {0.75, 0.9, 0.95}

Concluding remarks

We proposed a stochastic model for the competing degradation process (N t , Z t ) t≥0 describing initiation and propagation of degradation on a passive component of an electric power plant. This model associates a Poisson process together with a gamma process, where both processes may be homogeneous or not. This leads to four dierent models, which have been tted to EDF data. The nal retained model is based on a bootstrap condence interval for the shape parameter of the non homogeneous Poisson process and on EDF experts considerations for the gamma process. It is a combination of a NHP process for the initiation of the aw indications and a process for describing the growth of the indications.

The Monte Carlo simulation study shows that the sample size of the EDF data is sucient to guarantee the quality of estimation results. The study also shows that taking into account censoring is necessary to avoid some bias on the estimates of the gamma process parameters. This bias disappears when accounting for censoring.

The 95% condence intervals obtained on the EDF data may appear as rather large for some parameters such as a or α. However, looking at the condence intervals for the indicators of interest (mean and standard deviation of the initiation time of the rst aw indication, yearly mean growth of one aw indication, εquantiles of the residual life) provided in Tables 9, 10 and 11, we can see that they are not so large. The estimated indicators hence appear as accurate enough for helping the decision process (e.g. deciding whether the component should be replaced or not). The inuence of the number of initiated aw indications on a component residual life is also made clear by the study. It is hence of a major importance to go on collecting these data from the eld, and to take this degradation indicator into account for prognostic purpose.

From a statistical point of view, studying the asymptotic behavior of our parametric esti-mators and selecting one model are challenging objectives. The later point requires to use some information criteria like AIC (see [START_REF] Varin | A note on composite likelihood inference and model selection[END_REF]) or BIC (see [START_REF] Gao | Composite likelihood bayesian information criteria for model selection in high-dimensional data[END_REF]), previously developed in the composite likelihood approach as dened by Lindsay [START_REF] Lindsay | Composite likelihood methods[END_REF]. The use of these criteria would require the preliminary adaptation of our estimation method, which is based on a two-step procedure plus a composite likelihood method (see [START_REF] Gao | Composite likelihood bayesian information criteria for model selection in high-dimensional data[END_REF] for some recent results concerning BIC selection model criteria for composite likelihood approach).

Another challenging statistical issue would be to integrate measurement errors in our model, with two possible levels of errors. The rst level deals with the number of initiated aw indications, where typically, too small indications may be unrevealed by the testing processes.

The second level is linked to the sizing process, which may be perturbed by measurement errors.

Finally, the possibility to include a time inuence on the degradation dynamic could be considered by including a time eect on the degradation rate. For instance one could consider a similar model where the ith degradation process has a shape parameter depending on its initiation time T i : a = exp(θT i ) where θ ∈ (0, +∞). Note however that this model requires further investigation since the T i 's are generally interval censored. whatever t < t 0 or t ≥ t 0 . Now, let us set P(•) = P(•|N t 0 = 0). Under P, it is easy to check that the counting process (N t ) t≥0 still has independent increments. Based on the previous remark, (N t ) t≥0 hence is a Poisson process under P, with rate function λ(•). For t ≥ t 0 and n ≥ 0, we consequently have P(Z t ≤ z|N t 0 = 0, N t = n) = P(Z t ≤ z|N t = n) Second we consider the case z 0 > 0 (or equivalently: n 0 > 0). Let (t, n, z) such that t 0 < t, n 0 ≤ n and z 0 ≤ z. We have:

P (Z t 0 ≤ z 0 , N t 0 = n 0 , Z t ≤ z, N t = n) = P n 0 i=1 X (i) t 0 -T i ≤ z 0 ∩ n i=1 X (i) t-T i ≤ z ∩ {N t 0 = n 0 } ∩ {N t = n} = P n 0 i=1 X (i) t 0 -T i ≤ z 0 , X (i) t-T i ≤ z ∩ n i=n 0 +1 X (i)
t-T i ≤ z ∩ {N t 0 = n 0 } ∩ {N t = n} (using the convention that if n 0 = n, the corresponding empty intersection is equal to Ω).

Conditioning by σ (N ) t≥0 , X (i) for 1 ≤ i ≤ n 0 , we get

P (Z t 0 ≤ z 0 , N t 0 = n 0 , Z t ≤ z, N t = n) = E n 0 i=1 1 X (i) t 0 -T i ≤z 0 ,X (i) t-T i ≤z n i=n 0 +1 F A(t-T i ),b (z) 1 {Nt 0 =n 0} 1 {Nt=n} . Writing X (i) t-T i = X (i) t 0 -T i + X (i) t-T i -X (i) t 0 -T i
and using the independent increments of X (i) , we obtain through conditioning by σ (N ) t≥0 , X Given that {N t 0 = n 0 , N t = n}, the random vector (T 1 , . . . , T n 0 ) is conditionally identically distributed as the order statistic Y (1) , . . . , Y (n 0 ) of the same i.i.d. Y i 's as in (A.1) whereas the random vector (T n 0 +1 , . . . , T n ) can be seen to be conditionally identically distributed as the order statistic (V (n 0 +1) , . . . , V (n) ), where V n 0 +1 ,. . . , V n are i.i.d. copies of a random variable V with p.d.f. provided by

(i) t 0 -T i for 1 ≤ i ≤ n 0 P (Z t 0 ≤ z 0 , N t 0 = n 0 , Z t ≤ z, N t = n) = E n 0 i=1 1 X (i) t 0 -T i ≤z 0 F A(t-T i )-A(t 0 -T i ),b z -X (i) t 0 -T i n i=n 0 +1
f V (v) = λ (v) Λ (t) -Λ (t 0 ) 1 [t 0 ,t] (v) . (B.1)
Also, the vectors (T 1 , . . . , T n 0 ) and (T n 0 +1 , . . . , T n ) are conditionally independent. This provides

P (Z t 0 ≤ z 0 , Z t ≤ z|N t 0 = n 0 , N t = n) = E n 0 i=1 ϕ 1 (Y i ) n i=n 0 +1 ϕ 2 (V i ) = (E (ϕ 1 (Y ))) n 0 (E (ϕ 2 (V ))) n-n 0 = t 0 0 ϕ 1 (y) λ (y) Λ (t 0 ) dy n 0 × t t 0 ϕ 2 (v) λ (v) Λ (t) -Λ (t 0 ) dv n-n 0
and based on the independent increments of the Poisson process (N t ) t≥0 , we get P (Z t 0 ≤ z 0 , N t 0 = n 0 , Z t ≤ z, N t = n) = P (Z t 0 ≤ z 0 , Z t ≤ z|N t 0 = n 0 , N t = n) P (N t 0 = n 0 ) P (N t -N t 0 = n -n 0 ) = t 0 0 ϕ 1 (y) λ (y) dy

Appendix C Proof of Proposition 3

First case: z 0 = 0. Based on Equation (4), the quantile t ( ,t 0 ,0) ε veries t 0 +t ( ,t 0 ,0) ε t 0 Fa t ( ,t 0 ,0) ε +t 0 -v ,b ( ) λ (v) dv = log (1 -ε)

which may easily be written as [START_REF] Cha | A stochastic failure model with dependent competing risks and its applications to condition-based maintenance[END_REF].

Second case: z 0 > 0. We write t ε = t ( ,t 0 ,z 0 ,n 0 ) ε for sake of simplicity. Based on Corollary 1 and in case of a homogeneous gamma process, the quantile t ε is the solution of which may easily be written as [START_REF] Cox | A note on pseudolikelihood constructed from marginal densities[END_REF].

  ,b and survival function denoted by FA(t),b

  sequence of observation times for one component. In all the simulation study, the observation times of one component are chosen among these 228 time data sets up to a xed multiplicative constant. This rescaling of the time data through a multiplicative constant allows us to prevent from numerical instability and overow. The same rescaling is used for the application to EDF data in the next sub-section. Several possibilities are considered for the number N of observed components, with the following choice of time data sets: for N = 228, we simply take all the available 228 (rescaled) time data sets; for N = 4 × 228 = 912, each of the 228 time data sets is used four times; for N = 25, 50 and 100, the N time data sets are chosen randomly among the 228 available data sets, by sampling without replacement. Once the observation times aregiven for each component, the degradation samples are generated according to the NHPHG model. In order to be as close as possible from the EDF data, the simulated data are next censored with a xed censoring level c (which means that any indication measurement with size below c is left censored by c). In each of the following examples, the censoring level c is chosen such that the proportion of censored data is around 9%, which is similar to the EDF data.

  with cumulative intensity Λ(t) = αt β are equal to (α, β) = (1, 1.5). Parameters of the HG process are equal to (a, b) = (1, 2). The censoring level is c = 0.2.

  for N ∈ {25, 50, 100, 228}. The results mostly share the same good quality as for the rst parameter set. The computation times are a little longer however: considering the estimation of one parameter set for N = 228 components, it takes about 70 central processing unit (c.p.u.) time for the rst parameter set, 95 c.p.u. for the second one and 220 c.p.u. for the third one. This is due to the optimization step for the HG parameters which takes more time, apparently because of a atter objective function to maximize.

Figure 2 :

 2 Figure 2: Estimated εquantile t(90,t 0 ,z 0 ,n 0 )

= t 0 Ft 0 FF

 00 A(y),b (z) λ(t -y)dy Λ(t) n where the second equality holds by Proposition 1. Then we deriveP(Z t ≤ z|Z t 0 = 0) = n≥0 P(Z t ≤ z|N t 0 = 0, N t = n)P(N t = n) = n≥0 A(y),b (z) λ(t -y)dy Λ(t) n e -Λ(t) Λ (t) n n! = e -(Λ(t)-Λ(t 0 ))+ t-t 0 0 A(y),b (z)λ(t-y)dy = e -tt 0 FA(t-v),b (z)λ(v)dv .

FF

  A(t-T i ),b (z) 1 {Nt 0 =n 0} 1 {Nt=n} = E n 0 i=1 ϕ 1 (T i ) n i=n 0 +1 ϕ 2 (T i ) 1 {Nt 0 =n 0} 1 {Nt=n} where ϕ 1 (u) = E 1 X (i) t 0 -u ≤z 0 F A(t-u)-A(t 0 -u),b z -X A(t-u)-A(t 0 -u),b (z -x) f A(t 0 -u),b (x) dx, ϕ 2 (u) = F A(t-u),b (z) .

ϕ 2 FF 0 0FF

 20 (v) λ (v) dv n-n 0 e -Λ(t) n 0 ! (n -n 0 )! after simplication.Substituting ϕ 1 and ϕ 2 by their expressions, we now haveP (Z t 0 ≤ z 0 , N t 0 = n 0 , Z t ≤ z, N t = n) A(t-u)-A(t 0 -u),b (z -x) f A(t 0 -u),b (x) dx λ (u) du A(t-v),b (z) λ (v) dv n-n 0 e -Λ(t) n 0 ! (n -n 0 )! and P (Z t 0 ∈ dz 0 , N t 0 = n 0 , Z t ≤ z, N t = n) = t A(t-u)-A(t 0 -u),b (z -z 0 ) f A(t 0 -u),b (z 0 ) λ (u) du × A(t-u)-A(t 0 -u),b (z -x) f A(t 0 -u),b (x) dx λ (u) du

F 0 F

 0 A(t-v),b (z) λ (v) dv n-n 0 e -Λ(t) (n 0 -1)! (n -n 0 )! dz 0 . (B.2)Using Appendix A or (B.2), one easily getsP (Z t 0 ∈ dz 0 , N t 0 = n 0 ) = e -Λ(t 0 ) (n 0 -1)! t 0 A(t 0 -y),b (z 0 ) λ (y) dy

  t 0 -y),b (z 0 ) λ (y) dy dz 0 (B.3) from where we nally derive Equation (5) by dividing (B.2) by (B.3).

logF 0 0F

 0 {F a tε,b ( -z 0 )} -tε 0 λ (t 0 + u) Fa(tε-u),b ( ) du + (n 0a tε,b ( -x) f a(t 0 -u),b (x) dx du -(n 0 -1) log t a(t 0 -y),b (z 0 ) λ (y) dy = log (1 -ε) ,

  The points of the counting process (N t ) t≥0 stand for the initiating times of the aw indications on the component and are denoted by T 1 , T 2 , . . . , with 0 < T 1 < T 2 < • • • (almost surely). In case N t ≥ 1, we set Z with 1 ≤ j ≤ N t . The measurement of the largest aw indication at time t hence is:

	2 Competing degradation model and prognostic indicators
	2.1	Competing stochastic degradation processes
	We here specify the stochastic model for the initiation and development of aw indications on
	one single component. We rst set N t to be the random number of indications already initiated
	at time t. (j) t	to be the (random) measurement of the j-th aw indication
	at time t,

Table 1 :

 1 Main characteristics of the initiation time of the rst aw indication

	Mean	Variance	Quantile of order

Table 2 :

 2 Monte Carlo estimations for the NHP-HG model based on 500 simulated censored samples of size N (censoring level c = 0.2), without considering censoring (imputation method).Estimation taking into account censoring. We rst provide results in Table3for the same parameter sets as for the previous imputation method. The quality of the estimation results is now good for both processes (NHP and HG). The previous bias on the HG process parameters almost disappears and, contrary to the results of Table2, it decreases as the sample size increases. Condence intervals are well centered on the true parameters values. Standard deviations are still divided by two as the sample size is multiplied by four, which means that the "rootofN " convergence rate remains valid. Two other parameter sets are next considered in Tables

		N	a	b	a/b	a/b 2	α	β
	value		1	2	0.5	0.25	1	1.5
		25	1.37	2.57	0.53	0.22	1.00	1.50
		50	1.37	2.58	0.53	0.21	1.00	1.50
	mean	100	1.34	2.54	0.53	0.21	1.00	1.50
		228	1.31	2.49	0.52	0.21	1.00	1.50
		912	1.28	2.44	0.52	0.22	1.00	1.50
		25	0.39	0.62	0.03	0.05	0.12	0.05
		50	0.30	0.47	0.02	0.03	0.09	0.04
	stand. err. 100	0.21	0.35	0.02	0.03	0.07	0.03
		228	0.14	0.23	0.01	0.02	0.05	0.02
		912	0.11	0.18	0.01	0.03	0.02	0.01
		25	[0.81, 2.13] [1.67, 3.84] [0.47, 0.58] [0.15, 0.30] [0.81, 1.18] [1.43, 1.59]
		50	[0.94, 1.93] [1.86, 3.45] [0.49, 0.57] [0.16, 0.27] [0.86, 1.15] [1.44, 1.57]
	90% CI	100 [1.03, 1.70] [2.03, 3.15] [0.50, 0.55] [0.17, 0.25] [0.90, 1.11] [1.46, 1.55]
		228 [1.09, 1.55] [2.12, 2.89] [0.50, 0.54] [0.18, 0.24] [0.93, 1.08] [1.47, 1.53]
		912 [1.17, 1.40] [2.28, 2.65] [0.51, 0.53] [0.20, 0.23] [0.97, 1.04] [1.48, 1.51]
		25	[0.76, 2.22] [1.53, 3.91] [0.47, 0.59] [0.14, 0.33] [0.78, 1.26] [1.41, 1.60]
		50	[0.89, 2.07] [1.81, 3.69] [0.48, 0.58] [0.15, 0.28] [0.83, 1.17] [1.43, 1.58]
	95% CI	100 [0.99, 1.76] [1.94, 3.27] [0.49, 0.56] [0.16, 0.26] [0.87, 1.14] [1.45, 1.56]
		228 [1.06, 1.60] [2.08, 2.97] [0.50, 0.55] [0.18, 0.25] [0.91, 1.09] [1.46, 1.54]
		912 [1.15, 1.42] [2.23, 2.67] [0.51, 0.54] [0.20, 0.23] [0.96, 1.04] [1.48, 1.52]

Table 3 :

 3 Monte Carlo estimations for the NHP-HG model based on 500 simulated censored samples of size N considering censoring (censoring level c = 0.2).

		N	a	b	a/b	a/b 2	α	β
	value		1	2	0.5	0.25	1	1.5
		25	1.14	2.22	0.51	0.24	1.00	1.50
		50	1.08	2.13	0.51	0.24	1.00	1.50
	mean	100	1.03	2.04	0.50	0.25	1.01	1.50
		228	1.03	2.04	0.50	0.25	1.00	1.50
		912	0.99	1.99	0.50	0.25	1.00	1.50
		25	0.37	0.59	0.04	0.05	0.13	0.06
		50	0.23	0.38	0.02	0.04	0.09	0.04
	stand. err. 100	0.15	0.25	0.02	0.03	0.07	0.03
		228	0.12	0.20	0.01	0.04	0.04	0.02
		912	0.09	0.15	0.01	0.03	0.02	0.01
		25	[0.65, 1.87] [1.44, 3.39] [0.44, 0.57] [0.16, 0.33] [0.80, 1.22] [1.41, 1.59]
		50	[0.76, 1.49] [1.59, 2.83] [0.47, 0.55] [0.19, 0.30] [0.87, 1.16] [1.44, 1.57]
	90% CI	100 [0.79, 1.29] [1.66, 2.50] [0.47, 0.53] [0.21, 0.29] [0.91, 1.11] [1.45, 1.54]
		228 [0.88, 1.20] [1.80, 2.33] [0.48, 0.52] [0.22, 0.28] [0.93, 1.08] [1.47, 1.53]
		912 [0.92, 1.07] [1.86, 2.14] [0.49, 0.51] [0.24, 0.27] [0.97, 1.04] [1.48, 1.51]
		25	[0.57, 2.05] [1.27, 3.61] [0.43, 0.58] [0.15, 0.37] [0.76, 1.27] [1.40, 1.61]
		50	[0.70, 1.61] [1.52, 2.93] [0.46, 0.55] [0.18, 0.32] [0.82, 1.19] [1.43, 1.58]
	95% CI	100 [0.76, 1.37] [1.61, 2.64] [0.47, 0.54] [0.20, 0.30] [0.88, 1.14] [1.45, 1.55]
		228 [0.85, 1.25] [1.75, 2.41] [0.48, 0.52] [0.22, 0.28] [0.92, 1.09] [1.47, 1.54]
		912 [0.90, 1.10] [1.83, 2.17] [0.49, 0.51] [0.23, 0.27] [0.96, 1.05] [1.48, 1.52]

Table 4 :

 4 Monte Carlo estimations for the NHP-HG model based on 500 simulated censored samples of size N considering censoring (censoring level c = 6). .40, 3.20] [1.49, 2.97] [0.94, 1.07] [0.35, 0.64] [1.79, 2.26] [0.96, 1.05] 228 [1.55, 2.64] [1.60, 2.55] [0.96, 1.04] [0.40, 0.61] [1.86, 2.14] [0.97, 1.03]

		N	a	b	a/b	a/b 2	α	β
	value		2	2	1	0.5	2	1
		25	2.35	2.31	1.01	0.48	2.01	1.00
	mean	50	2.22	2.19	1.01	0.49	2.00	1.00
		100	2.13	2.11	1.00	0.49	2.01	1.00
		228	2.03	2.02	1.00	0.50	2.00	1.00
		25	0.86	0.75	0.06	0.14	0.24	0.05
	stand. err.	50	0.66	0.58	0.05	0.14	0.18	0.04
		100	0.47	0.41	0.04	0.13	0.12	0.02
		228	0.29	0.25	0.03	0.09	0.08	0.02
		25	[1.13, 3.97] [1.21, 3.67] [0.90, 1.10] [0.29, 0.76] [1.62, 2.41] [0.92, 1.08]
	90% CI	50	[1.33, 3.55] [1.42, 3.37] [0.94, 1.07] [0.31, 0.67] [1.72, 2.29] [0.95, 1.06]
		100 [1.48, 2.93] [1.56, 2.83] [0.95, 1.06] [0.37, 0.62] [1.83, 2.22] [0.96, 1.04]
		228 [1.60, 2.51] [1.65, 2.45] [0.96, 1.04] [0.42, 0.59] [1.88, 2.12] [0.97, 1.03]
		25	[1.05, 4.00] [1.16, 3.75] [0.88, 1.13] [0.28, 0.80] [1.56, 2.49] [0.90, 1.11]
	95% CI	50	[1.29, 3.79] [1.34, 3.62] [0.92, 1.09] [0.29, 0.71] [1.67, 2.33] [0.94, 1.07]
		100 [1

Table 5 :

 5 Monte Carlo estimations for the NHP-HG model based on 500 simulated censored samples of size N considering censoring (censoring level c = 7.5).

		N	a	b	a/b	a/b 2	α	β
	value		2	1	2	2	3	0.75
		25	2.30	1.13	2.02	1.93	3.00	0.75
	mean	50	2.15	1.07	2.01	1.95	3.01	0.75
		100	2.08	1.03	2.01	1.97	3.00	0.75
		228	2.03	1.01	2.00	2.01	3.01	0.75
		25	0.75	0.3	0.10	0.51	0.30	0.04
	stand. err.	50	0.49	0.21	0.07	0.35	0.22	0.03
		100	0.35	0.15	0.05	0.25	0.16	0.02
		228	0.24	0.11	0.06	0.38	0.10	0.01
		25	[1.25, 3.89] [0.66, 1.76] [1.86, 2.20] [1.21, 2.94] [2.50, 3.55] [0.69, 0.82]
	90% CI	50	[1.43, 3.09] [0.75, 1.46] [1.89, 2.12] [1.43, 2.60] [2.63, 3.35] [0.70, 0.80]
		100 [1.58, 2.69] [0.82, 1.29] [1.92, 2.10] [1.60, 2.39] [2.74, 3.26] [0.71, 0.78]
		228 [1.70, 2.44] [0.86, 1.19] [1.95, 2.06] [1.69, 2.28] [2.83, 3.19] [0.73, 0.77]
		25	[1.16, 3.96] [0.62, 1.82] [1.81, 2.23] [1.17, 3.04] [2.41, 3.61] [0.68, 0.83]
	95% CI	50	[1.32, 3.26] [0.69, 1.55] [1.86, 2.15] [1.35, 2.74] [2.59, 3.43] [0.70, 0.81]
		100 [1.49, 2.82] [0.77, 1.36] [1.91, 2.11] [1.53, 2.51] [2.71, 3.30] [0.71, 0.79]
		228 [1.63, 2.52] [0.84, 1.23] [1.94, 2.07] [1.67, 2.33] [2.80, 3.22] [0.72, 0.78]

Table 6 :

 6 Estimation and bootstrap results based on 1,000 bootstrapped samples of size 228 for the parameters of the Poisson process (top: rows 3 to 7) and estimated yearly average number of new indications for one component (bottom: rows 8 to 12)

	Model	HP		NHP
	Parameter	λ (×10 -4 )	α (×10 -9 )	β
	Estimation	1.56	1.22	2.28
	Mean	1.57	2.24	2.28
	St. dev.	0.17	2.60	0.13
	90% conf. int.	[1.31, 1.85]	[0.18, 7.56] [2.08, 2.49]
	95% conf. int.	[1.28, 1.91]	[0.12, 9.13] [2.05, 2.54]
			Year 1 8.47 × 10 -4	Year 6	1.71 × 10 -2
	Yearly average		Year 2 3.27 × 10 -3	Year 8	2.55 × 10 -2
	number of new 5.71 × 10 -2	Year 3 6.26 × 10 -3	Year 10	3.45 × 10 -2
	indications		Year 4 9.62 × 10 -3	Year 12	4.44 × 10 -2
			Year 5 1.33 × 10 -2	Year 15	5.92 × 10 -2

Table 7 :

 7 Estimation results for the propagation process on EDF data and the dierent models

	(left: columns 2 to 4), and the corresponding estimated expected annual growth for one aw
	indication (right: columns 5 and 6)
	Model	Estimates of the parameters of the gamma process	Expected annual growth of one aw indication
	HPHG		â	b		1.12
			2.44 × 10 -4	0.08
	NHPHG		â	b		1.97
			8.63 × 10 -4	0.16
	HPNHG	η	γ	b	Year 1	11.61
		0.47 9.22 × 10 -9	0.04	Year 2	7.43 × 10 -8
	NHPNHG	η	γ	b	Year 1	10.93
		0.36 2.24 × 10 -10	0.03	Year 2	1.70 × 10 -9

Table 8 :

 8 Estimation (second line) and bootstrap (third line) results on EDF data based on

	1, 000 bootstrapped samples of size 228 for the NHPHG model parameters
		α (×10 -9 )	β	â (×10 -4 )	b	â/ b (×10 -3 ) â/ b2 (×10 -2 )
	Estimates	1.22	2.28	8.63	0.16	5.48	3.48
	Mean	2.24	2.28	11.28	0.20	5.49	3.37
	St. dev.	2.60	0.13	12.36	0.20	0.42	1.05
	[q 0.05 , q 0.95 ]	[0.18, 7.56] [2.08, 2.49] [5.15, 19.55] [0.10, 0.33] [4.83, 6.20]	[1.73, 4.99]
	[q 0.025 , q 0.975 ] [0.12, 9.13] [2.05, 2.54] [4.80, 26.72] [0.10, 0.44] [4.69, 6.31]	[1.35, 5.48]

Table 11 :

 11 Estimation of t, the quantiles of order ε of the hitting time of level = 90 for z 0 = 0 at t 0 = 0 and inference results based on 10 3 bootstrapped samples of size 228 [q 0.05 , q 0.95 ] Mean [St. dev.] [q 0.05 , q 0.95 ]

		(90,0,0)
		ε
	ε	0.9		0.95
	Data	75.65	79.87
		Mean [St. dev.]
	Bootstrap	75.89 [4.57]	[69.66, 83.66]	80.09 [5.25]	[73.41, 88.58]

Table 12 :

 12 Estimated εquantiles of the residual life τ 90 -t 0 given various values of (z 0 , n 0 ) at time t 0 = 25 and ε ∈ {0.75, 0.9, 0.95} t 0 z 0 n 0 t

				( ,t 0 ,z 0 ,n 0 ) 0.75	t ( ,t 0 ,z 0 ,n 0 ) 0.90	t 0.95 ( ,t 0 ,z 0 ,n 0 )	t 0 z 0 n 0 t ( ,t 0 ,z 0 ,n 0 ) 0.75	t	( ,t 0 ,z 0 ,n 0 ) 0.90	( ,t 0 ,z 0 ,n 0 ) t 0.95
	25	20	2	36.681	41.480	44.283	25	40	2	29.611		34.494	37.473
	25	20	4	32.264	36.660	39.225	25	40	4	26.562		30.974	33.580
	25	20	6	29.381	33.580	35.981	25	40	6	24.380		28.716	31.008
	25	20	8	27.052	31.402	33.652	25	40	8	22.631		26.859	29.214
	25	20	10	25.119	29.683	31.904	25	40	10	21.046		25.395	27.647
	25	30	2	33.404	38.362	41.142	25	50	2	25.141		29.981	32.880
	25	30	4	29.622	34.023	36.558	25	50	4	22.954		27.470	30.054
	25	30	6	27.070	31.311	33.637	25	50	6	21.252		25.560	28.051
	25	30	8	25.093	29.258	31.571	25	50	8	19.749		24.076	26.345
	25	30	10	23.341	27.655	29.876	25	50	10	18.374		22.781	25.125

Acknowledgements The authors thank the Associate Editor and the two Referees for their very careful reading of the paper and for their deep questions, which have helped us to highly improve the paper. This work was supported by the French National Research Agency (ANR), AMMSI project, ref. ANR-2011-BS01-021.

Additional prognostic indicator evaluation. All the computations of this paragraph have been carried out with MATLAB. The estimation results concerning the initiation time of the rst aw indication are given in Table 9. As we saw in Section 2.3, the initiation time is Weibull distributed. The estimated mean is close to 19 years with a standard deviation estimated to 9 years. From the bootstrap study, the results seem to be reliable. [q 0.05, q 0.95 ] [q 0.025, q 0.975 ] We nish this section with estimation results for the quantile of the residual life τ 90 -t 0 given (Z t 0 , N t 0 ) = (z 0 , n 0 ) for the critical level = 90. First, Figure 2 depicts the estimated quantile curve ε → t(90,t 0 ,z 0 ,n 0 ) ε of the residual life as a function of ε (the time scale is given in years) for three dierent values of ( , t 0 , z 0 , n 0 ). We can see that, for each value of (t 0 , z 0 , n 0 ), the εquantile increases with ε. Looking at the boundary values for ε (0 + and 1 -), we can also see that the εquantile quickly increases (towards ∞) when ε approaches to 1 -and quickly decreases (towards 0) when ε goes to 0 + . These observations are consistent with what might have been expected. Second, we consider components for which no aw indication has appeared at the end of the observation period (denoted by t 0 ), that is components for which z 0 = 0, and thus n 0 = 0, at time t 0 > 0. Corresponding quantiles t (90,0,0) ε are given for several values of ε (ε ∈ {0.9, 0.95}) in Table 11. The mean and standard deviation values are calculated using the non-parametric bootstrap method. Figure 3 plots similar quantiles with respect to t 0 for ε equal to 0.75, 0.9 and 0.95. Table 12 nally provides estimated εquantiles of the residual life τ 90 -t 0 given various values of (z 0 , n 0 ) at time t 0 = 25 and ε ∈ {0.75, 0.9, 0.95}. Looking at Figure 3 and Table 12, we can see that, as expected and for a given ε, the εquantile is decreasing with respect of t 0 , z 0 and n 0 . Also, we can conclude from these results that the residual life of a component is clearly inuenced by the observed number of initiated aw indications, which consequently should continue to be collected.

Appendix A Conditional distribution of Z t given N t = n

The case n = 0 is clear so that only the case n ≥ 1 is considered. For n ≥ 1 and z ≥ 0, we have

Given that N t = n, the random vector (T 1 , . . . , T n ) is known to be identically distributed as the order statistics Y (1) Based on the independence between the X (i) 's and N t , we obtain

t-Y (i) ≤ z does not depend on the order of the Y (i) 's for the last equality. Conditioning by the Y i 's and remembering that F A(u),b is the common c.d.f. of the X (i) u 's (and using the independence between the X (i) 's and the Y i 's), we now have

n based on the i.i.d. property of the Y i 's for the second line and on the fact that the common p.d.f. of the Y i 's is provided by (A.1) for the last line.

The derivation of the p.d.f. with respect to Lebesgue measure of the conditional distribution of Z t given N t = n is now clear, which achieves this proof.

Appendix B Proof of Proposition 2

First we consider the case z 0 = 0 (or equivalently: n 0 = 0). Let us rst observe that given N t 0 = 0, the random variable N t is conditionally Poisson distributed with parameter Λ(t) = (Λ(t) -Λ(t 0 )) + = t t 0 λ(u) du,