Christophe Chesneau

Hassan S Bakouch

Pedro Luiz Ramos
email: pedrolramos@usp.br

A continuous probability

Keywords: Hazard rate, Two-parameter distributions, Reliability and statistical measures, Maximum Likelihood Estimation, Data applications

come   L'archive ouverte pluridisciplinaire

Introduction with motivations

In analysis of the lifetime data, monotone hazard rates are common. Such data can be modelled using the log-normal, Weibull and gamma distributions. The Weibull distribution is more popular than log-normal and gamma because the survival and hazard rate functions of the last two distributions have not a closed form and hence numerical integrations are required. Gupta and Kundu [START_REF] Gupta | Exponentiated exponential family: An alternative to gamma and weibull distributions[END_REF] introduced the exponentiated exponential (EE) distribution as an extension to the exponential distribution and also as an alternative to the gamma distribution. Further developments on the exponentiated exponential distribution can be seen in [START_REF] Gupta | Generalized exponential distribution: Existing results and some recent development[END_REF].

In many practical applications, continuous probability models that allow occurrence of zero values have vast importance, for example in forecast models when we observe the monthly rainfall precipitation, it is common in dry periods the non occurrence of precipitation, therefore the occurrence of zero values can be observed in different measures, such as the average, maximum and minimum. In survival analysis, we may observe data with instantaneous failure due construction problem. Another example of zero occurrence is the hydrologic data in arid and semiarid regions, like annual peak flow discharges. Moreover, zero occurrence can be met in many areas, e.g. manufacturing defects, medical consultations, hydrology, ecology and econometrics.

The above four distributions (log-normal, Weibull, gamma and exponentiated exponential) do not provide the characteristic of zero occurrence. Therefore, we introduce a new continuous probability model that allows occurrence of zero values and it shall be called Polynomial-exponential (PE) distribution. This model represents a powerful alternative to the mentioned distributions above and also an extension to exponential, linear exponential and other combinations of Weibull distributions as it shall be seen later. The model is specified in terms of the cumulative distribution function (cdf) as:

F (x) = 1 -e -λx x α -1
x-1 , x ∈ (0, +∞)\{1}, where α > 0, λ > 0 (and F (x) = 0 if x < 0) and Figure 1 presents plots of the cdf for different values of α and λ = 2. Note that, there is a discontinuity when x = 1, to avoid that we can considered the continuous extension for x = 1 given by

F (x) = 1 -e -λx x α -1 x-1 , x ∈ (0, +∞)\{1}, 1 -e -λα , x = 1. (1)
Moreover, a significant account of mathematical properties for the new distribution are provided and the hazard rate function has constant, increasing or decreasing shape, which making the PE distribution an alternative to the mentioned distributions above. Another attractive feature of the PE distribution is that it has closed form expressions for its cdf and hazard rate function, which is not the case, for the log-normal and gamma distributions. Further, the distribution has several particular sub-models. For α = 1, the PE distribution gives the exponential and when α = 2, it reduces to the one parameter linear exponential distribution. Also, when α is an integer, we can express F (x) as

F (x) = 1 -exp -λ α k=1
x k , x > 0, and hence the survival function

S(x) = 1 -F (x)
of the PE distribution is the product of the survival functions of Weibull distributions with parameters (λ, 1), (λ, 2), ..., and (λ, α) with respect to value of k, which means the distribution has an exponential general polynomial. Also, we can note F (x) given by (1) has the form:

F (x) = G(H(x)),
where G(x) is the cdf associated to the exponential distribution with parameter λ > 0 and

H(x) = x x α -1 x-1 (with H(1) = α) is a positive increasing function with H(0) = 0 and lim x→+∞ H(x) = +∞.
The inferential procedure for the parameters of PE distribution is presented using the maximum likelihood estimation. It is shown that one of the estimators can be obtained in closed-form and this allows us to obtain the estimates solving a simple one non-linear equation. The performance of the MLEs are compared using extensive numerical simulations.

The paper is organized as follows. In Section 2 we introduce the PE distribution. Section 3 is devoted to some of its mathematical properties. Estimations of the parameters via the maximum likelihood method are investigated in Section 4. A simulation analysis is given in Section 5. In Section 6 we apply our proposed model in two real data sets. Finally in Section 7 we conclude the paper.

Polynomial-exponential distribution

The associated probability density function (pdf) of the cdf given by the equation ( 1) is

f (x) = λ αx α+1 -(α + 1) x α + 1 (x -1) 2 e -λx x α -1 x-1 , x ∈ (0, +∞)\{1}, (2)
with the continuous extension for x = 1: f (1) = λα(α+1) 2 e -λα . Figure 2 presents some plot of the PE distribution for different values of α and λ, and showing various shapes of the density function with left skewness. As we know, many distributions such as log-normal, Weibull, gamma and exponentiated exponential, to list a few, do not allow occurrence of zero values. In this regard, the following remark shows that the PE distribution can be used as a model with occurrence of zero values.

Let us observe that f (0) = λ > 0 for all α > 0 and λ > 0. Further, it follows from equation (3) that f (x) ∼ λαx α-1 e -λx α as x → ∞. Therefore, the upper tail behavior of the pdf is a product of a polynomial power and an exponential polynomial power decay, both of them depend only on α. Obviously, larger values of α lead to faster decay of the upper tail, which interprets α as a shape parameter.

The hazard rate function (hrf) is given by

h(x) = λ αx α+1 -(α + 1) x α + 1 (x -1) 2 , x ∈ (0, +∞)\{1},
with the continuous extension for x = 1 given by h(1) = λα(α+1)

2

. The study of the behavior of the hazard function is not an easy task. Glaser's [START_REF] Glaser | Bathtub and related failure rate characterizations[END_REF] lemma is difficult to be implemented since η(t) = -d dt log(f (t)) does not has a simple form. However, from graphical analysis we observed that the hazard function presents an decreasing hazard rate for α > 1 and λ > 0, increasing hazard rate for α < 1 and λ > 0 and constant rate for α = 1, and for this purpose some plots of the hazard function with various values for the parameters α and λ are presented in Figure 3. Moreover, note that, h(x) ∼ λ as x → 0, and h(x) ∼ λαx α-1 as x → ∞. Hence, we conclude that, the lower tail of the hazard rate function is a constant, while its upper tail is a polynomial which allows for increasing, decreasing and constant hazard rate shapes.

Remark 1. Using the indicator function: 1 A (x) = 1 if x ∈ A and 0 elsewhere, we have following analytic expressions for F (x), f (x) and h(x):

F (x) = 1 -e -λx x α -1 x-1 1(0,+∞)\{1}(x) 1 -e -λα 1{1}(x) 1 (0,+∞) (x), f (x) = λ αx α+1 -(α + 1) x α + 1 (x -1) 2 e -λx x α -1 x-1 1(0,+∞)\{1}(x) × (3)
× λα(α + 1) 2 e -λα 1{1}(x)
1 (0,+∞) (x), and

h(x) = λ αx α+1 -(α + 1) x α + 1 (x -1) 2 1(0,+∞)\{1}(x) λα(α + 1) 2 1{1}(x)
1 (0,+∞) (x).

These analytic expressions will be useful in the next.

Mathematical properties

Some useful expansions

The result below presents a polynomial expansion of the cdf F (x) given by (1).

Proposition 3.1. We have the following expansion, for x ∈ (0, +∞)\{1},

F (x) = 1 - ∞ k=0 k =0 ∞ j=0 A k, ,j (-1) k x k+α +j 1 (0,1) (x) + x α -j 1 (1,+∞) (x) ,
where

A k, ,j = k -k j 1 k! (-1) +j λ k .
Proof of Proposition 3.1. First of all, let us now investigate an expansion for e -λx x α -1

x-1 by distinguishing the case x ∈ [0, 1) and the case x > 1.

• If x ∈ [0, 1), note that x α -1 x-1 = (1 -x α ) 1 1-x .
The exponential and binomial series give

e -λx x α -1 x-1 = ∞ k=0 1 k! -λx x α -1 x -1 k = ∞ k=0 1 k! (-λ) k x k (1 -x α ) k 1 (1 -x) k = ∞ k=0 k =0 ∞ j=0 k -k j 1 k! (-1) +j (-λ) k x k+α +j . • If x > 1, note that x x α -1 x-1 = -(1 -x α ) 1 1-1 x
. It follows from the exponential and binomial series that

e -λx x α -1 x-1 = ∞ k=0 1 k! -λx x α -1 x -1 k = ∞ k=0 1 k! λ k (1 -x α ) k 1 1 -1 x k = ∞ k=0 k =0 ∞ j=0 k -k j 1 k! (-1) +j λ k x α -j . Hence e -λx x α -1 x-1 = ∞ k=0 k =0 ∞ j=0 A k, ,j (-1) k x k+α +j 1 (0,1) (x) + x α -j 1 (1,+∞) (x) .
Therefore

F (x) = 1 -e -λx x α -1 x-1 = 1 - ∞ k=0 k =0 ∞ j=0 A k, ,j (-1) k x k+α +j 1 (0,1) (x) + x α -j 1 (1,+∞) (x) .
This complete the proof of Proposition 3.1.

The result below presents an expansion of the pdf f (x) given by ( 2) via polynomial and the exponential function e -λx , which will be important to ensure the permutation of sum and integral in several probabilistic quantities. Proposition 3.2. We have the following expansion, for x ∈ (0, +∞)\{1},

f (x) = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k R k, ,j (x)1 (0,1) (x) + S k, ,j (x)1 (1,+∞) (x) e -λx , where B k, ,j = k -k -2 j 1 k! (-1) +j λ k+1 , (4)
R k, ,j (x) = αx 2k+(α-1) +j+α+1 -(α + 1) x 2k+(α-1) +j+α + x 2k+(α-1) +j and S k, ,j (x) = αx k+(α-1) -j+α-1 -(α + 1) x k+(α-1) -j-2+α + x k+(α-1) -j-2 .

Proof of Proposition 3.2. Let us observe that

λx x α -1 x -1 = λx x α -1 x -1 -1 + λx = λx 2 x α-1 -1 x -1 + λx.
Therefore, we can express the pdf f (x) as

f (x) = λ(αx α+1 -(α + 1) x α + 1)e -λx × 1 (x -1) 2 e -λx 2 x α-1 -1 x-1 .
Let us now investigate an expansion for 1 (x-1) 2 e -λx 2 x α-1 -1

x-1 by distinguishing x ∈ [0, 1) and x > 1.

• If x ∈ [0, 1), it follows from the exponential and binomial series that

1 (x -1) 2 e -λx 2 x α-1 -1 x-1 = ∞ k=0 1 k! 1 (x -1) 2 -λx 2 x α-1 -1 x -1 k = ∞ k=0 1 k! (-λ) k x 2k (1 -x α-1 ) k 1 (1 -x) k+2 = ∞ k=0 k =0 ∞ j=0 k -k -2 j 1 k! (-1) +j (-λ) k x 2k+(α-1) +j .
• If x > 1, exponential and binomial series give

1 (x -1) 2 e -λx 2 x α-1 -1 x-1 = ∞ k=0 1 k! 1 (x -1) 2 -λx 2 x α-1 -1 x -1 k = ∞ k=0 1 k! λ k x k (1 -x α-1 ) k x -2 1 1 -1 x k+2 = ∞ k=0 k =0 ∞ j=0 k -k -2 j 1 k! (-1) +j λ k x k+(α-1) -j-2 . Hence 1 (x -1) 2 e -λx x α -1 x-1 = ∞ k=0 k =0 ∞ j=0 k -k -2 j 1 k! (-1) +j λ k × (-1) k x 2k+(α-1) +j 1 (0,1) (x) + x k+(α-1) -j-2 1 (1,+∞) (x) .
Owing to this equality, we obtain the desired expansion:

f (x) = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k x 2k+(α-1) +j 1 (0,1) (x) + x k+(α-1) -j-2 1 (1,+∞) (x) × αx α+1 -(α + 1) x α + 1 e -λx = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k R k, ,j (x)1 (0,1) (x) + S k, ,j (x)1 (1,+∞) (x) e -λx .
This ends the proof of Proposition 3.2.

Moments and moment generating function

Here and after, we consider a random variable X following the PE(α, λ) distribution with α > 0 and λ > 0.

We define the upper incomplete gamma function as Γ(s, x) = +∞ x t s-1 e -t dt and the lower incomplete gamma function as γ(s, x) = x 0 t s-1 e -t dt, s > 0, x ≥ 0.

Let r ≥ 0. Using the notations and the result of Proposition 3.2, the r-moments of X is given by

E(X r ) = ∞ -∞ x r f (x)dx = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k U k, ,j,r + V k, ,j,r ,
where U k, ,j,r = 1 0 x r R k, ,j (x)e -λx dx and V k, ,j,r = +∞ 1

x r S k, ,j (x)e -λx dx. We have

U k, ,j,r = α 1 0 x r+2k+(α-1) +j+α+1 e -λx dx -(α + 1) 1 0
x r+2k+(α-1) +j+α e -λx dx

+ 1 0 x r+2k+(α-1) +j e -λx dx = α γ(r + 2k + (α -1) + j + α + 2, λ) λ r+2k+(α-1) +j+α+2 -(α + 1) γ(r + 2k + (α -1) + j + α + 1, λ) λ r+2k+(α-1) +j+α+1 + γ(r + 2k + (α -1) + j + 1, λ) λ r+2k+(α-1) +j+1 .
On the other hand, we have

V k, ,j,r = α +∞ 1 x r+k+(α-1) -j+α-1 e -λx dx -(α + 1) +∞ 1
x r+k+(α-1) -j-2+α e -λx dx

+ +∞ 1 x r+k+(α-1) -j-2 e -λx dx = α Γ(r + k + (α -1) -j + α, λ) λ r+k+(α-1) -j+α -(α + 1) Γ(r + k + (α -1) -j -1 + α, λ) λ r+k+(α-1) -j-1+α + Γ(r + k + (α -1) -j -1, λ) λ r+k+(α-1) -j-1 .
The moment generating function of X is given by, for t < λ,

M (t) = E(e tX ) = ∞ -∞ e tx f (x)dx = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k U * k, ,j (t) + V * k, ,j (t) , where U * k, ,j (t) = 1 0 e tx R k, ,j (x)e -λx dx and V * k, ,j (t) = +∞ 1
e tx S k, ,j (x)e -λx dx. We have

U * k, ,j (t) = α 1 0 x 2k+(α-1) +j+α+1 e -(λ-t)x dx -(α + 1) 1 0 x 2k+(α-1) +j+α e -(λ-t)x dx + 1 0 x 2k+(α-1) +j e -(λ-t)x dx = α γ(2k + (α -1) + j + α + 2, λ -t) (λ -t) 2k+(α-1) +j+α+2 -(α + 1) γ(2k + (α -1) + j + α + 1, λ -t) (λ -t) 2k+(α-1) +j+α+1 + γ(2k + (α -1) + j + 1, λ -t) (λ -t) 2k+(α-1) +j+1 .
On the other hand, we have

V * k, ,j (t) = α +∞ 1 x k+(α-1) -j+α-1 e -(λ-t)x dx -(α + 1) +∞ 1
x k+(α-1) -j-2+α e -(t-λ)x dx

+ +∞ 1 x k+(α-1) -j-2 e -(λ-t)x dx = α Γ(k + (α -1) -j + α, λ -t) (λ -t) k+(α-1) -j+α -(α + 1) Γ(k + (α -1) -j -1 + α, λ -t) (λ -t) k+(α-1) -j-1+α + Γ(k + (α -1) -j -1, λ -t) (λ -t) k+(α-1) -j-1 .

On other means and moments

The following result proposes an expansion of the primitive t 0 x r f (x)dx, with t > 0. It will be useful in the next. Proposition 3.3. For any r ≥ 0 and t > 0, we have

t 0 x r f (x)dx = ∞ k=0 k =0 ∞ j=0 B k, ,j (-1) k U • k, ,j,r (t) + V • k, ,j,r (t)1 [1,+∞) (t) , (5)
where B k, ,j is defined by (4),

U • k, ,j,r (t) = α γ(r + 2k + (α -1) + j + α + 2, λ min(t, 1)) λ r+2k+(α-1) +j+α+2 -(α + 1) γ(r + 2k + (α -1) + j + α + 1, λ min(t, 1)) λ r+2k+(α-1) +j+α+1 + γ(r + 2k + (α -1) + j + 1, λ min(t, 1)) λ r+2k+(α-1) +j+1 and V • k, ,j,r (t) = α Γ(r + k + (α -1) -j + α, λt) λ r+k+(α-1) -j+α -(α + 1) Γ(r + k + (α -1) -j -1 + α, λt) λ r+k+(α-1) -j-1+α + Γ(r + k + (α -1) -j -1, λt) λ r+k+(α-1) -j-1 .
The proof of Proposition 3.3 follows from Proposition 3.2 with

U • k, ,j,r (t) = min(t,1) 0 x r R k, ,j (x)e -λx dx and V • k, ,j,r (t) = t 1 x r S k, ,j (x
)e -λx dx. The expressions of these integrals in terms of upper incomplete gamma function and the lower incomplete gamma function is obtained proceeding as Subsection 3.2.

Several crucial conditional moments use the integral t 0 x r f (x)dx for various values of r. The most useful of them are presented below. For any t > 0,

• The r-th conditional moments of X is given by,

E(X r | X > t) = 1 1 -F (t) +∞ t x r f (x)dx = 1 1 -F (t) E(X r ) - t 0 x r f (x)dx .
• The r-th reversed moments of X is given by

E(X r | X ≤ t) = 1 F (t) t 0 x r f (x)dx. Let µ = E(X).
• The mean deviations of X about µ is given by

δ = E(|X -µ|) = 2µF (µ) -2 µ 0 xf (x)dx
• The mean deviations of X about the median M is given by

η = E(|X -M |) = µ -2 M 0 xf (x)dx.
Residual life parameters can be also determined using E(X r ) and t 0 x r f (x)dx for several values of r. In particular,

• The mean residual life is defined as

K(t) = E(X -t | X > t) = 1 S(t) E(X) - t 0 xf (x)dx -t
and the variance residual life is given by

V (t) = V ar(X-t | X > t) = 1 S(t) E(X 2 ) - t 0 x 2 f (x)dx -t 2 -2tK(t)-[K(t)] 2 .
• The mean reversed residual life is defined as

L(t) = E(t -X | X ≤ t) = t - 1 F (t) t 0 xf (x)dx
and the variance reversed residual life is given by

W (t) = V ar(t -X | X ≤ t) = 2tL(t) -[L(t)] 2 -t 2 + 1 F (t) t 0 x 2 f (x)dx.

Stress-strength reliability

Let X be a random variable following the PE(α, λ 1 ) distribution with pdf denoted by f X (x) and Y be a random variable following the PE(α, λ 2 ) distribution with cdf denoted by F Y (x), with α > 0, λ 1 > 0 and λ 2 > 0. Then the stress-strength reliability is defined by R = P (X > Y ). Since the integral on (0, +∞) of the pdf (2) with parameters (α, λ 1 + λ 2 ) denoted by f * (x) is equal to one, we have

R = P (X > Y ) = +∞ 0 f X (x)F Y (x)dx = 1 - +∞ 0 λ 1 αx α+1 -(α + 1) x α + 1 (x -1) 2 e -(λ1+λ2)x x α -1 x-1 dx = 1 - λ 1 λ 1 + λ 2 +∞ 0 f * (x)dx = λ 2 λ 1 + λ 2 .
This result is of interest in a parametric estimation context; only λ 1 and λ 2 need to be estimated to have an estimation of R (the maximum likelihood estimators for λ 1 and λ 2 yield the maximum likelihood estimator for R by the plug-in method).

Order statistics distributions

We now introduce order statistics and present some of their properties in our mathematical framework (general results can be found, for instance, in [START_REF] David | Wiley series in probability and statistics[END_REF]). Let X 1 , X 2 , . . . , X n be n i.i.d. random variables following the PE(α, λ) distribution with α > 0 and λ > 0. Let us consider its order statistics as X 1:n , X 2:n , . . . , X n:n , i.e., for any i ∈ {1, . . . , n}, X i:n ∈ {X 1 , . . . , X n } with X 1:n ≤ . . . ≤ X n:n (so

X 1:n = X (1) = inf(X 1 , . . . , X n ) and X n:n = X (n) = sup(X 1 , . . . , X n )).
Let us now present some important distributions related to X 1:n , X 2:n , . . . , X n:n . Some important of them involving our distribution are presented below. The general expression of the cdf of X i:n is given by

F Xi:n (x) = n! (i -1)! (n -i)! n-i k=0 n -i k (-1) k i + k [F (x)] i+k , x ∈ R.
Hence, for any x ∈ (0, +∞)/{1}, we have

F Xi:n (x) = n! (i -1)! (n -i)! n-i k=0 n -i k (-1) k i + k 1 -e -λx x α -1
x-1 i+k .

For the case x = 1, we have

F Xi:n (1) = n! (i -1)! (n -i)! n-i k=0 n -i k (-1) k i + k 1 -e -λα i+k .
The general expression of the pdf of X i:n is given by

f Xi:n (x) = n! (i -1)! (n -i)! [F (x)] i-1 [1 -F (x)] n-i f (x), x ∈ R.
Thus, for any x ∈ (0, +∞)/{1}, we have

f Xi:n (x) = n! (i -1)! (n -i)! λ 1 -e -λx x α -1 x-1 i-1 αx α+1 -(α + 1) x α + 1 (x -1) 2 e -λ(n-i+1)x x α -1 x-1 .
For the case x = 1, we have

f Xi:n (1) = n! (i -1)! (n -i)!
1 -e -λα i-1 λα(α + 1) 2 e -λα(n-i+1) .

As in Remark 1, one can express f Xi:n (x) in a one form using 1 (0,+∞)\{1} (x) and 1 {1} (x). For i < j and x i < x j , the general expression of the joint pdf of (X i:n , X j:n ) is given by

f (Xi:n,Xj:n) (x i , x j ) = n! (i -1)! (n -j)! (j -i -1) [F (x i )] i-1 [F (x j ) -F (x i )] j-i-1 × ×[1 -F (x j )] n-j f (x i )f (x j ).
For the case (x i , x j ) ∈ (0, +∞) 2 /{(1, 1)}, we have

f (Xi:n,Xj:n) (x i , x j ) = n! (i -1)! (n -j)! (j -i -1) λ 2 1 -e -λxi x α i -1 x i -1 i-1 e -λxi x α i -1 x i -1 -e -λxj x α j -1 x j -1 j-i-1 × (αx α+1 i -(α + 1) x α i + 1)(αx α+1 j -(α + 1) x α j + 1) (x i -1) 2 (x j -1) 2 e -λxi
x α i -1

x i -1 -λ(n-j+1)xj

x α j -1

x j -1 .

The expression of f (Xi:n,Xj:n) (x i , x j ) for (x i , x j ) ∈ (0, +∞) 2 /{(1, 1)} can be set in a similar manner, using the values of F (1) and f (1).

In the following proposition, we provide the asymptotic distributions of the extreme values X 1:n and X n:n , and show that they are exponential and Gumbel distributions, respectively, which adapt the standards of the asymptotic distribution of extremes. Proposition 3.4. Let (X n ) n≥1 be a sequence of i.i.d. random variables following the P E(α, λ) distribution, then

• nX (1) n≥1 converges in distribution to a random variable X having the exponential distribution of parameter λ.

• X α (n) -log(n) λ n≥1
converges in distribution to a random variable X having the Gumbel distribution of parameters 0 and 1 λ . Proof of Proposition 3.4. Let us prove the two points in turn.

• Since X 1 , . . . , X n are i.i.d., using standard mathematical arguments, for x ∈ (0, +∞)/{n}, the cdf of nX (1) is given by

F nX(1) (x) = 1 -1 -F x n n = 1 -e -λx ( x n ) α -1 x n -1
.

So lim n→+∞ F nX(1) (x) = 1-e -λx = F X (x). This ends the proof of the first point. • Again, since X 1 , . . . , X n are i.i.d., using standard mathematical arguments, for

x ∈ -log(n) λ , +∞ / 1 -log(n) λ , the cdf of X α (n) -log(n) λ is given by F X α (n) -log(n) λ (x) = F x + log(n) λ 1 α n =   1 -e -λ(x+ log(n) λ ) 1 α x+ log(n) λ -1 ( x+ log(n) λ ) 1 α -1    n .
Therefore, when n → +∞, several equivalences give

F X α (n) -log(n) λ (x) ∼ e n log 1-e -λx n ∼ e -e -λx .
Hence lim n→+∞ F X α

(n) -log(n) λ (x) = e -e -λx = F X (x)
. The second point is proved.

Stochastic ordering

The ordering mechanism in life time distributions can be illustrate by the concept of stochastic ordering. See, for instance, [START_REF] Shaked | Stochastic orders and their applications[END_REF]. This subsection presents the basic of this concept, with a result using the proposed distribution. A random variable X is said to be stochastically smaller than a random variable Y in the

• stochastic order (X ≤ st Y ) if the associated cdfs satisfy: F X (x) ≥ F Y (x) for all x.
• hazard rate order (X ≤ hr Y ) if the associated hrfs satisfy: h X (x) ≥ h Y (x) for all x. • likelihood ratio order (X ≤ lr Y ) if the ratio of the associated pdfs given by fX (x) fY (x) decreases in x.

Important equivalences exist; when the supports of X and Y have a common finite left end-point, then we have: X ≤ lr Y ⇒ X ≤ hr Y ⇒ X ≤ st Y . Proposition 3.5. Let X be a random variable following the PE(α, λ 1 ) distribution with pdf denoted by f X (x) and Y be a random variable following the PE(α, λ 2 ) distribution with pdf denoted by f Y (x), with α > 0, λ 1 > 0 and λ 2 > 0. If λ 1 > λ 2 , then we have X ≤ lr Y .

Proof of Proposition 3.5. For any x ∈ (0, +∞)/{1}, we have

f X (x) f Y (x) = λ 1 λ 2 e (λ2-λ1)x x α -1
x-1

and fX (1) fY (1) = λ1 λ2 e (λ2-λ1)α . As mentioned in Introduction, the function

H(x) = x x α -1 x-1
is increasing function of x. Indeed, we have

H (x) = 1+αx α+1 -(1+α)x α (x-1) 2
and a study of

function shows that 1 + αx α+1 -(1 + α)x α ≥ 0 for any α > 0. Hence, if λ 1 > λ 2 , fX (x) fY (x)
decreases in x and X ≤ lr Y . Proposition 3.5 is proved.

Remark 2. One can prove that λ 1 > λ 2 implies that X ≤ hr Y,which follows immediately from the definition (2): for any x ∈ (0, +∞)/{1}, we have h

X (x) = λ 1 αx α+1 -(α+1)x α +1 (x-1) 2 ≥ λ 2 αx α+1 -(α+1)x α +1 (x-1) 2 = h Y (x), and for x = 1, we have h X (1) = λ1α(α+1) 2 ≥ λ2α(α+1) 2 = h Y (1).

Record values distributions

Let us now focus our attention on record values and present some of their properties in our mathematical context (general results can be found, for instance, in [START_REF] Ahsanullah | Record statistics[END_REF]). Let X 1 , X 2 , . . . , X n be n i.i.d. random variables following the PE(α, λ) distribution with α > 0 and λ > 0. We define a sequence of record times U (n) as follows:

U (1) = 1, U (n) = min{j; j > U (n -1), X j > X U (n-1) } for n ≥ 2.
We define the i-th upper record value by

R i = X U (i) , with R 1 = X 1 .
The general expression of the cdf of R i is given by

F Ri (x) = 1 -(1 -F (x)) i-1 k=0 [-log(1 -F (x))] k k! , x ∈ R.
Hence, for any x ∈ (0, +∞)/{1}, we have

F Ri (x) = 1 -e -λx x α -1 x-1 i-1 k=0 1 k! λx x α -1 x -1 k .
For the case x = 1, we have

F Ri (1) = 1 -e -λα i-1 k=0 (λα) k k! .
The general expression of the pdf of R i is given by

f Ri (x) = [-log(1 -F (x))] i-1 (i -1)! f (x), x ∈ R.
Hence, for any x ∈ (0, +∞)/{1}, we have

f Ri (x) = 1 (i -1)! λ i x i-1 (x α -1) i-1 (αx α+1 -(α + 1) x α + 1) (x -1) i+1 e -λx x α -1 x-1 .
Note that, for x = 1, we have

f Ri (1) = 1 (i -1)! (λα) i α + 1 2 e -λα .
The general expression of the joint pdf of (R 1 , . . . , R n ) is given by, for (x 1 , . . . ,

x n ) ∈ R n with x 1 < . . . < x n , f (R1,...,Rn) (x 1 , . . . , x n ) = f (x n ) n-1 k=1 h(x k ).
For the case (x 1 , . . . , x n ) ∈ (0, +∞) n /{(1, . . . , 1)}, we have 2 .

f (R1,...,Rn) (x 1 , . . . , x n ) = λ n+1 αx α+1 n -(α + 1) x α n + 1 (x n -1) 2 e -λxn x α n -1 xn -1 n-1 k=1 αx α+1 k -(α + 1) x α k + 1 (x k -1)
The expression of f (R1,...,Rn) (x 1 , . . . , x n ) for (x 1 , . . . , x n ) ∈ (0, +∞) n /{(1, . . . , 1)} can be set in a similar manner, using the values of f (1) and h(1). For i < j and x i < x j , the general expression of the joint pdf of (R i , R j ) is given by

f (Ri,Rj) (x i , x j ) = [-log(1 -F (x i ))] i-1 (i -1)! log 1-F (xi) 1-F (xj) j-i-1 (j -i -1)! h(x i )f (x j ).
For the case (x i , x j ) ∈ (0, +∞) 2 /{(1, 1)}, we have

f (Ri,Rj) (x i , x j ) = [-log(1 -F (x i ))] i-1 (i -1)! log 1-F (xi) 1-F (xj) j-i-1 (j -i -1)! h(x i )f (x j ) = 1 (i -1)! λ j x i x α i -1 x i -1 i-1 1 (j -i -1)! x j x α j -1 x j -1 -x i x α i -1 x i -1 j-i-1 × (αx α+1 i -(α + 1) x α i + 1)(αx α+1 j -(α + 1) x α j + 1) (x i -1) 2 (x j -1) 2 e -λxj
x α j -1

x j -1 .

The expression of f (Ri,Rj) (x i , x j ) for (x i , x j ) ∈ (0, +∞) 2 /{(1, 1)} can be set in a similar manner, using the values of F (1), f (1) and h(1).

Maximum likelihood estimation

Let X 1 , X 2 , . . . , X n be a random sample with common distribution the PE(α, λ) distribution with α > 0 and λ > 0. Let θ = (α, λ) be the parameter vector and x 1 , x 2 , . . . , x n be the observed values. Then the likelihood function associated to x 1 , . . . , x n is given by

L(θ) = n i=1 λ αx α+1 i -(α + 1) x α i + 1 (x i -1) 2 e -λxi x α i -1 x i -1 1(0,+∞)\{1}(xi) λα(α + 1) 2 e -λα 1{1} (xi)
.

For the set of simplicity, let us set

u i = 1 (0,+∞)\{1} (x i ), v i = 1 {1} (x i ),
and

n i=1 u i + n i=1 v i = n.
The log-likelihood function can be expressed as

(θ) = log(λ) n i=1 u i + n i=1 u i log αx α+1 i -(α + 1) x α i + 1 (x i -1) 2 -λ n i=1 u i x i x α i -1 x i -1 + n i=1 v i (log(λ) + log(α) + log(α + 1) -log(2)) .
The nonlinear log-likelihood equations ∂ (θ) ∂θ = 0 are given by

∂ (θ) ∂α = n i=1 u i x α i (x i -1)[1 + α log(x i )] -log(x i ) αx α+1 i -(α + 1) x α i + 1 -λ n i=1 u i x i x α i x i -1 log(x i ) + 1 α + 1 α + 1 n i=1 v i = 0, (6)
and

∂ (θ) ∂λ = n λ - n i=1 u i x i x α i -1 x i -1 = 0. (7)
Note that after some algebraic manipulations in [START_REF] Gupta | Generalized exponential distribution: Existing results and some recent development[END_REF] we have that

λ = n n i=1 u i x i x α i -1 xi-1 . (8)
Replacing ( 8) in ( 6) the maximum likelihood estimates of α and λ are determined by solving the one linear equation [START_REF] Gupta | Exponentiated exponential family: An alternative to gamma and weibull distributions[END_REF]. Since it does not admit any explicit solution, numerical procedures can be used. Under mild conditions the maximum likelihood estimators are asymptotically normal, with an asymptotic variance-covariance matrix depending on the Fisher information matrix. Crucial quantities to determine the entries of this matrix are the second partial derivatives of the log-likelihood function given by

∂ 2 (θ) ∂α 2 = n i=1 u i x α i log(x i )((x i -1)[2 + α log(x i )] -log(x i ))(αx α+1 i -(α + 1) x α i + 1) (αx α+1 i -(α + 1) x α i + 1) 2 - n i=1 u i x α i x α i [(x i -1)[1 + α log(x i )] -log(x i )] 2 (αx α+1 i -(α + 1) x α i + 1) 2 -λ n i=1 u i x i x α i x i -1 (log(x i )) 2 - 1 α 2 + 1 (α + 1) 2 n i=1 v i , and
∂ 2 (θ) ∂λ 2 = - n λ 2 , ∂ 2 (θ) ∂λ∂α = - n i=1 u i x i x α i x i -1 log(x i ).

Simulation analysis

In this section a simulation study is presented to compare the efficiency of the maximum likelihood method. The comparison is performed by computing the Bias and the mean square errors (MSE) given by

Bias(α i ) = 1 N N i=1 αi -α , MSE(α i ) = 1 N N j=1 (α i -α) 2 , Bias(λ i ) = 1 N N i=1 λi -λ , MSE(λ i ) = 1 N N j=1 ( λi -λ) 2 ,
where N is the number of estimates obtained through the MLE. The 95% coverage probability of the asymptotic confidence intervals are also evaluated. Here we expect that the most efficient estimation method returns both Bias and MSE closer to zero. Additionally, for a large number of experiments, using a 95% confidence level, the frequencies of intervals that covered the true values of α and λ should be closer to 95%. The programs can be obtained, upon request. The values of the PE were generated considering the following algorithm:

(1) Generate U i ∼ Uniform(0, 1), i = 1, . . . , n;

(2) Find x i from the solution of F (x i ) -u i = 0 , i = 1, . . . , n;

The simulation study is performed under the assumption (0.5, 2) and (4, 2), N = 100, 000 and n = (20, 35, . . ., 460). The chosen values allow us to obtain data with both increasing (α < 1) and decreasing (α > 1) hazard rate. It is important to point out that, the results of this simulation study were similar for different choices of α and λ. Figures 4 and5 present the Bias, the MSE and the coverage probability with a 95% confidence level of the estimates obtained through the MLE for different samples of size. From the obtained results, we can conclude that as there is an increase of n both Bias and MSE tend to zero, i.e., the estimator are asymptotic efficiency. Moreover, the coverage probability of the confidence levels tend to the nominal value assumed 0.95. Therefore, the MLE showed to be a good estimator for the parameters of the PE distribution.

Application to real data

In this section, we illustrate the flexibility of our proposed distribution by considering two real data sets. The results obtained from the PE distribution are compared with ones of the Weibull, Gamma, Lognormal and the EE distributions, and nonparametric survival function Here, different discrimination criterion are considered based on log likelihood function. Let k be the number of parameters to be fitted and θ the MLEs of θ, the discrimination criterion methods are respectively: The best model is the one which provides the minimum values of these criteria. The Kolmogorov-Smirnov (KS) test is also considered aiming to check the goodness of the fit for the models. This procedure is widely known and based on the KS statistic D n = sup x |F n (x) -F (x; θ)|, where sup x is the supremum of the set of distances, F n (x) is the empirical distribution function and F (x; θ) is the cdf of the fitted distribution. Under a significance level of 5% if the data comes from F (x; θ) (null hypothesis), the hypothesis is rejected if the P-value is smaller than 0.05.

The next subsections give a description of the used data and their analysis under the mentioned distributions above.

Air conditioning system data

The data have been presented by Proschan [START_REF] Proschan | Theoretical explanation of observed decreasing failure rate[END_REF] and further analyzed by Adamidis and Loukas [START_REF] Adamidis | A lifetime distribution with decreasing failure rate[END_REF]. Table 1 consists of the number of successive failures of the air conditioning system of each member of a fleet of 13 Boeing 720 jet airplanes.

Table 2 displays the MLEs, standard-error and 95% confidence intervals for α and λ. Table 3 presents the results of AIC, AICc, HQIC, CAIC criteria, for the compared In Figure 6, we have the TTT-plot, the survival function adjusted by the compared distributions and the non-parametric survival function.

Comparing the empirical survival function with the adjusted models we observe a goodness of the fit for the PE distribution, which is confirmed from different discrimination criterion methods as the PE distribution has the minimum value for all statistics and the largest for the P-value. Consequently, we conclude that the data related to the number of successive failures of the air conditioning system of each member of a fleet of 13 Boeing 720 jet airplanes can be described by the PE distribution. The adjusted models when compared to the empirical survival show a goodness of the fit for the PE distribution. Additionally, this result is corroborated by the P-value of the KS test. Therefore, our proposed distribution can be used to describe the data related to the total monthly rainfall during September at São Carlos.

Figure 1 .

 1 Figure 1. Cumulative density function shapes for PE distribution considering different values of α and fixed λ = 2.

Figure 2 .

 2 Figure 2. Density function shapes for PE distribution considering different values of α and λ.

Figure 3 .

 3 Figure 3. Hazard function shapes for PF distribution for λ = 2 and considering different values of α.

Figure 4 .

 4 Figure 4. Bias, MSEs related from the estimates of α = 3 and λ = 2 for N simulated samples under the MLE.

Figure 5 .

 5 Figure 5. Bias, MSEs related from the estimates of α = 0.9 and λ = 0.1 for N simulated samples under the MLE.

•

  Akaike information criterion AIC = -2l( θ; x) + 2k; • Corrected Akaike information criterion AICC = AIC +(2 k (k + 1))/(n -k -1); • Hannan-Quinn information criterion HQIC = -2 l( θ; x) + 2 k log (log(n)); • Consistent Akaike information criterion CAIC = -2 l( θ; x) + k (log(n) + 1).

Figure 7 .

 7 Figure 7. Survival function adjusted by the compared distributions and a non-parametric method considering the data set related to the total monthly rainfall during September at São Carlos.

Table 1 .

 1 Data set related to the number of successive failures of the air conditioning system of each member of a fleet of 13 Boeing 720 jet airplanes.

	194 15	41	29	33 181 413 14	58	37 100 65
	9	169 447 184 36 201 118 34	31	18	18	67
	57	62	7	22	34	90	10	60 186 61	49	14
	24	56	20	79	84	44	59	29 118 25 156 310
	76	26	44	23	62 130 208 70 101 208 74	57
	48	29 502 12	70	21	29 386 59	27 153 26
	326 55 320 56 104 220 239 47 246 176 182 33
	15 104 35	23 261 87	7	120 14	62	47 225
	71 246 21	42	20	5	12 120 11	3	14	71
	11	14	11	16	90	1	16	52	95	97	51	11
	4	141 18 142 68	77	80	1	16 106 206 82
	54	31 216 46 111 39	63	18 191 18 163 24
	50	44 102 72	22	39	3	15 197 188 79	88
	46	5	5	36	22 139 210 97	30	23	13	14
	359	9	12 270 603	3	104	2	438 50 254	5
	283 35	12 130 493 487 18 100	7	98	5	85
	91	43 230	3	130 102 209 14	57	54	32	67
	59 134 152 27	14 230 66	61	34
	distributions.

Table 2 .

 2 MLE, Standard-error and 95% confidence intervals for α and λ.

	θ	MLE S. error	CI 95% (θ)
	α	0.9010 0.01791 (0.7906; 1.0113)
	λ	0.0032 0.00003 (0.0076; 0.0282)

Table 3 .

 3 Results of AIC, AICc, HQIC, CAIC criteria and the p-value for the KS test for the compared distributions considering the number of successive failures of the air conditioning system of each member of a fleet of 13 Boeing 720 jet airplanes.

	Test	PE	Weibull Gamma Lognormal	EE
	AIC	2358.61 2359.17 2360.58	2361.76	2360.81
	AICc	2354.67 2355.23 2356.64	2357.82	2356.86
	CAIC 2367.34 2367.89 2369.30	2370.48	2369.52
	HQIC 2361.33 2361.89 2363.30	2364.470	2363.52
	P-value 0.63435 0.61336 0.37719	0.56520	0.33913

Table 6 .

 6 Results of KS test for the compared distributions considering the data set related to the total monthly rainfall during September at São Carlos.In the Figure7, the survival function adjusted by the compared distributions and the Kaplan-Meier estimator.

	Test	PE	Weibull Gamma Lognormal	EE
	P-value 0.66711 0.0000	0.0000	0.0000	0.0000

Figure 6. Survival function adjusted by the compared distributions and a non-parametric method considering the data sets related to the number of successive failures of the air conditioning system of each member of a fleet of 13 Boeing 720 jet airplanes.

Monthly rainfall data

In this subsection, we considered the data set firstly presented in Bakouch et al. [START_REF] Bakouch | Binomial-exponential 2 distribution: Different estimation methods with weather applications[END_REF].

The data set is related to the total monthly rainfall during September at São Carlos located in southeastern Brazil. Such city has an active industrial profile and high agricultural importance where the study of the behavior of dry and wet periods has proved to be strategic and economically significant its development. Table 4 presents the data related to the total monthly rainfall (mm) during September at São Carlos. Nadarajah and Haghighi [START_REF] Nadarajah | An extension of the exponential distribution[END_REF] observed that maximum likelihood estimate of the shape parameter is non-unique for the Gamma, Weibull and Generalized exponential distributions if data set consists of zeros and therefore none of these three distributions can fit this kind of data set. On the other hand the PE distribution is defined as x ≥ 0, which allow us to use the original values in the presence of zero. Table 5 displays the MLE, standard-error and 95% confidence intervals for α and λ. Table 6 presents the results of the P-value for the KS test for the compared distributions.


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

Carolina Meier-Hirmer

Michel Roussignol

Bivariate Gamma wear processes for track geometry modelling, with application to intervention scheduling

Keywords: Gamma wear process, L evy process, dependent stochastic processes, reliability, track maintenance

Bivariate Gamma wear processes for track geometry modelling, with application to intervention scheduling

Sophie Mercier, Carolina Meier-Hirmer, Michel Roussignol

Introduction

This paper is concerned with the intervention scheduling of a railway track, based on the observation of two indicators, which measure the deterioration of the geometry track. The railway track is considered as deteriorated when any of these two indicators is beyond a given threshold. The goal of the paper is the study of the intervention scheduling, which must ensure that, given some observations provided by inspection, the railway track will remain serviceable until the next intervention with a high probability.

The deterioration of the track geometry is characterized by the development of di erent representative parameters like, for example, the levelling of the track. Figure 1 shows the defects that are measured by two of these parameters: the longitudinal (NL) and transversal (NT) levelling indicators.

Longitudinal levelling (NL)

Transversal levelling (NT) Based on expert judgements, a Gamma process has been used in (Meier-Hirmer et al.

2009) both to model the evolution of the NL indicator and to plan maintenance actions.

As noted by J.M. van Noortwijk in his survey (van Noortwijk 2009), this stochastic process is widely used in reliability studies (see also [START_REF] Abdel-Hameed | A gamma wear process[END_REF][START_REF] Grall | Continuous { time predictive { maintenance scheduling for a deteriorating system[END_REF][START_REF] Zuckerman | Optimal replacement policy for the case where the damage process is a one-sided L evy process[END_REF]). Various domains of applications exist, such as civil engineering 3 [START_REF] Buijs | Time-dependent reliability analysis of ood defences using gamma processes[END_REF][START_REF] Cinlar | Stochastic process for extrapolating concrete creep[END_REF], highway engineering [START_REF] Nicolai | A comparison of models for measurable deterioration: An application to coatings on stell structures[END_REF] or railway engineering [START_REF] Meier-Hirmer | Maintenance optimisation for system with a Gamma deterioration process and intervention delay: application to track maintenance[END_REF]. Gamma processes are also used in other domains, such as nance [START_REF] Joshi | Intensity Gamma : a new approach to pricing portfolio credit derivatives[END_REF] or risk analysis [START_REF] Dufresne | Risk theory with the gamma process[END_REF]. All these papers use univariate Gamma processes.

In the present case, as the two indicators NL and NT are dependent, the use of a bivariate model is required. For this purpose, di erent processes might be used, such as

Bessel [START_REF] Newby | A bivariate process model for maintenance and inspection planning[END_REF] or L evy processes [START_REF] Kallsen | Characterization of dependence of multidimensional L evy processes using L evy copulas[END_REF]. In this paper, the approach of F.A. with parameters ( i ; 1) for i 2 f1; 2; 3g and from b 1 > 0, b 2 > 0, one de nes:

X (1) t = Y (1) t + Y (3) t =b 1 ; and X (2) t = Y (2) t + Y (3) t
=b 2 for all t 0:

The process (X t ) t 0 = X

t ; X

(2) t t 0 is then a homogeneous process in time with independent increments and it is a L evy process. The marginal processes of (X t ) t 0 are univariate Gamma processes with respective parameters (a i ; b i ), where a i = i + 3 for i = 1; 2.

For any bivariate L evy process, the correlation coe cient Xt of X

(1)

t and X

(2)

t is known to be independent of t. For a bivariate Gamma process, one obtains:

= Xt = 3 p a 1 a 2 and 1 = a 1 p a 1 a 2 ; 2 = a 2 p a 1 a 2 ; 3 = p a 1 a 2 .
This entails

0 max = min (a 1 ; a 2 ) p a 1 a 2 : (1)
See [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF]) section XI.3 for similar results on bivariate Gamma distributions.

This leads to two equivalent parameterizations of a bivariate Gamma process:

( 1 ; 2 ; 3 ; b 1 ; b 2 ) and (a 1 ; a 2 ; b 1 ; b 2 ; ).

November 29, 2010 10:17

Structure and Infrastructure Engineering swp0000

5

With the parameterization ( 1 ; 2 ; 3 ; b 1 ; b 2 ), the joint p.d.f. of X t is:

g t (x 1 ; x 2 ) = b 1 b 2 Z min(b1x1;b2x2) 0 f 1t;1 (b 1 x 1 x 3 ) f 2 t;1 (b 2 x 2 x 3 ) f 3t;1 (x 3 ) dx 3 ; = b 1 b 2 e b1x1 b2x2 ( 1 t) ( 2 t) ( 3 t) Z min(b1x1;b2x2) 0 (b 1 x 1 x 3 ) 1t 1 (b 2 x 2 x 3 ) 2t 1 x 3t 1 3 e x3 dx 3 (2) for x 1 0, x 2 0.
A L evy process is characterized both by its Laplace transform and by its L evy measure.

It is easy to compute their expressions for the bivariate Gamma process constructed above.

For x = (x 1 ; x 2 ) with x 1 0 and x 2 0, the Laplace transform of X(t) is equal to :

E e x X(t) = E e x1=b1Y (1) t (x1=b1+x2=b2)Y (3) t x2=b2Y (2) t = 1 1 + x 1 =b 1 1t 1 1 + x 1 =b 1 + x 2 =b 2 3t 1 1 + x 2 =b 2 2t = b 1 b 1 + x 1 (a1 p a1a2)t b 1 b 2 b 1 b 2 + x 1 b 2 + x 2 b 1 p a1a2t b 2 b 2 + x 2 (a2 p a1a2)t

:

Since the L evy measure X is caracterized by

E e x Xt = exp 0 B @t Z Z [0;+1) 2 e x y 1 X (dy) 1 C A for all x = (x 1 ; x 2 ) with x 1 0 and x 2 0, it is easy to obtain d X (y 1 ; y 2 ) = a 1 p a 1 a 2 y 1 e b1y1 dy 1 d 0 (y 2 ) + a 2 p a 1 a 2 y 2 e b2y2 d 0 (y 1 ) dy 2 + p a 1 a 2 y 1 e b1y1 dy 1 d b1y1=b2 (y 2 ) :
The tail integral function U X (x 1 ; x 2 ) = (]x 1 ; +1[ ]x 2 ; +1[) will also be used in the following. It is equal to:

U (x 1 ; x 2 ) = p a 1 a 2 Z +1 sup(b1x1;b2x2)
e y y dy:

(3)

6

A more general Gamma process may be constructed as a bivariate L evy process with univariate Gamma processes as marginal processes. For such a process, if

a 1 , b 1 , a 2 , b 2
are the parameters of the marginal Gamma processes, the linear correlation coe cient may be proved to satisfy:

0 Z Z [0;+1) 2 min r a 1 a 2 E(u 1 ); r a 2 a 1 E(u 2 ) du 1 du 2 where E(u) = R +1 u e y
y dy. If a 1 6 = a 2 , this upper bound is strictly greater than the upper bound (1) for the process constructed by trivariate reduction. Then the construction by trivariate reduction cannot provide all levels of dependence for Gamma processes since it does not cover the entire range of possible linear correlation coe cients. Nevertheless it is su cient for our applicative purpose here and in all the following a bivariate Gamma process stands for a process constructed by trivariate reduction.

Parameter estimation

The data used for the parameter estimation are values of the process increments for non overlapping time intervals on a single trajectory, and also on di erent independent trajectories. The data can then be represented as realizations of t j ; X

(1) j ; X

(2) j 1 j n where t j = t j s j stands for a time increment and

X (i) j = X (i) tj X (i)
sj for the associated i-th marginal increment (i = 1; 2). For di erent j, the random vectors X

(1) j ; X

(2) j are independent, but not identically distributed. The random variable X (i)

j (i = 1;
2) is Gamma distributed with parameters (a i t j ; b i ). The joint p.d.f. of the random vector

X (1) j ; X (2) j
is equal to g tj (:; :), with t j substituted to t in (2). In the same way as for the parameter estimation of an univariate Gamma process, both empirical and maximum likelihood methods are possible in the bivariate case.

Empirical estimators

Using E X (i) j

= ai bi t j and Var X (i) j

= ai b 2 i t j for i = 1; 2 and for all j, empirical estimators (â 1 ; b1 ; â2 ; b2 ) of (a 1 ; b 1 ; a 2 ; b 2 ) are given in [START_REF] Cinlar | Stochastic process for extrapolating concrete creep[END_REF] (4)

Using

Cov X

(1)

j ; X (2) j = p a 1 a 2 b 1 b 2 t j ;
a similar estimator ^ may be given for , with:

^ p â1 â2 b1 b2 = P n j=1 X (1) j â1 b1
t j X

(2) j â2 b2

t j P n j=1 t j 1 P n j=1 tj P n j=1 ( t j ) 2 :
(5)

These estimators satisfy:

E âi bi = a i b i ; E âi b2 i ! = a i b 2 i ; E ^ p â1 â2 b1 b2 = p a 1 a 2 b 1 b 2 :
If the time increments t j are equal, these estimators co• ncide with the usual empirical estimators in the case of i.i.d. random variables.

Maximum likelihood estimators

The parameter estimation of a univariate Gamma process is usually done by maximizing the likelihood function (see e.g. [START_REF] Meier-Hirmer | Maintenance optimisation for system with a Gamma deterioration process and intervention delay: application to track maintenance[END_REF]). With this method, estimators a i and b i (i = 1; 2) of the marginal parameters are computed by solving the equations:

a i b i = P n j=1 X (i) j P n j=1 t j and 0 @ n X j=1 t j 1 A ln a i P n j=1 t j P n j=1 X (i) j ! + n X j=1 t j ln X (i) j ( a i t j ) = 0;
where

(x) = d (x) dx (x) , (x) = Z 1 0 e u u x 1 du
for all x > 0 ( is the Digamma function).

In order to estimate all the parameters of the bivariate process ( 1 ; 2 ; 3 ; b 1 ; b 2 ) (which are here preferred to (a 1 ; b 1 ; a 2 ; b 2 ; )), the likelihood function associated with the data t j ; X

(1) j ; X

(2) j 1 j n can be written as

L( 1 ; 2 ; 3 ; b 1 ; b 2 ) = Q n j=1 g tj ( X (1)
j ; X

(2) j ) where g tj is provided by (2). However, because of the expression of the function g t (:; :), it seems complicated to optimize this likelihood function directly. An EM algorithm (see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF]) is then used, considering Y

(3) j = Y

(3) tj

Y

(3) sj 1 j n as hidden data. This procedure is still too complicated for estimating the ve parameters and does not work numerically. So, the procedure is restricted to the three parameters ( 1 ; 2 ; 3 ). For the parameters b 1 ; b 2 , the values b 1 ; b 2 computed using the maximum likelihood method for each univariate marginal process are taken.

For sake of simplicity, the realizations of

t j ; X (1) j ; X (2) j ; Y (3) j 1 j n are denoted by t j ; x (1) j ; x (2) j ; y
(3) j 1 j n in the following, the associated n-dimensional random vectors by X (1) ; X

(2) ; Y

(3) and the associated n-dimensional data vectors by x (1) ; x (2) ; y (3) .

The joint p.d.f. of the random vector X

(1)

t ; X (2) t ; Y (3) t is equal to: b 1 b 2 f 1t;1 (b 1 x 1 y 3 ) f 2 t;1 (b 2 x 2 y 3 ) f 3t;1 (y 3 ) = b 1 b 2 ( 1 t) ( 2 t) ( 3 t) e (b1x1+b2x2) (b 1 x 1 y 3 ) 1t 1 (b 2 x 2 y 3 ) 2t 1 y 3t 1 3 e y3 ;
with 0 y 3 min (b 1 x 1 ; b 2 x 2 ), x 1 > 0 and x 2 > 0.

Then, the log-likelihood function Q x (1) ; x (2) ; y (3) associated with the complete data

x (1) ; x (2) ; y (3) is derived:

Q x (1) ; x (2) ; y (3) = n (ln (b 1 ) + ln (b 2 )) n X j=1 (ln ( 1 t j ) + ln ( 2 t j ) + ln ( 3 t j )) b 1 n X j=1 x (1) j b 2 n X j=1 x (2) j + n X j=1 ( 1 t j 1) ln b 1 x (1) j y (3) j + ( 2 t j 1) ln b 2 x
(2) j y

(3) j + ( 3 t j 1) ln y

(3) j + y

(3) j

:

For the EM algorithm, the conditional log-likelihood of the complete data given the November 29, 2010 10:17

Structure and Infrastructure Engineering swp0000 9 observed data is needed:

E Q X (1) ; X (2) ; Y (3) j X (1) = x (1) ; X (2) = x (2) = n (ln (b 1 ) + ln (b 2 )) b 1 n X j=1 x (1) j b 2 n X j=1 x (2) j + n X j=1 ( 1 t j 1) E ln b 1 x
(1) j

Y

(3) j j X

(1)

j = x (1) j ; X (2) j = x (2) j + ( 2 t j 1) E ln b 2 x
(2) j

Y

(3) j j X

(1)

j = x (1) j ; X (2) j = x (2) j + ( 3 t j 1) E ln Y (3) j j X
(1)

j = x (1) j ; X (2) j = x (2) j +E Y (3) j j X (1) j = x (1) j ; X (2) j = x (2) j n X j=1 (ln ( 1 t j ) + ln ( 2 t j ) + ln ( 3 t j )) : (6)
Finally, the conditional probability density function of

Y (3) t given X (1) t = x 1 ; X (2) t = x 2
is equal to:

f 1t;1 (b 1 x 1 y 3 ) f 2t;1 (b 2 x 2 y 3 ) f 3 t;1 (y 3 ) R min(b1x1;b2x2) 0 f 1t;1 (b 1 x 1 x 3 ) f 2t;1 (b 2 x 2 x 3 ) f 3 t;1 (x 3 ) dx 3 = (b 1 x 1 y 3 ) 1t 1 (b 2 x 2 y 3 ) 2 t 1 y 3 t 1 3 e y3 R min(b1x1;b2x2) 0 (b 1 x 1 x 3 ) 1t 1 (b 2 x 2 x 3 ) 2t 1 x 3t 1 3 e x3 dx 3 ;
where 0 y 3 min (b 1 x 1 ; b 2 x 2 ), x 1 > 0 and x 2 > 0.

Step k of the EM algorithm consists of computing new parameter values

( (k+1) 1 ; (k+1) 2 ; (k+1) 3
) given the current values (

(k) 1 ; (k) 2 ; (k)
3 ) in two stages: stage 1: compute the conditional expectations in (6) using the current set

( (k) 1 ; (k) 2 ; (k)
3 ) of parameters, with:

f 1 j; k) 1 ; (k) 2 ; (k) 3 = E ln b 1 x (1) j Y (3) j j X (1) = x (1) j ; X (2) = x (2) j ; f 2 j; k) 1 ; (k) 2 ; (k) 3 = E ln b 2 x (2) j Y (3) j j X (1) = x (1) j ; X (2) = x (2) j ; f 3 j; k) 1 ; (k) 2 ; (k) 3 = E ln Y (3) j j X (1) = x (1) j ; X (2) = x (2) j ; h k) 1 ; (k) 2 ; (k) 3 = n X j=1 E Y (3) j j X (1) = x (1) j ; X (2) = x
(2) j : stage 2: take for (

(k+1) 1 ; (k+1) 2 ; (k+1) 3
) the values of ( 1 ; 2 ; 3 ) that maximize (6), 2). We notice that the parameters a (k) i stabilize more quickly than the parameter (k) . In other ways, the parameters

(k) i + (k)
3 for i = 1, 2 are quickly stable (about 5 iterations), but values of (k) i are much longer to stabilize (between 20 and 30 iterations). The conclusion of this section is that the estimation of the parameters (a i ; b i ) by empirical and maximum likelihood methods both give satisfactory results, with a slight preference to maximum likelihood estimation. Empirical estimators of have a good order of magnitude, but are sometimes not precise enough. Estimators of obtained by EM are always reasonable. The estimation of the three parameters ( 1 ; 2 ; 3 ) (column EM1) seems to give slightly better results than those obtained for 3 alone (column EM2). The results obtained by the EM algorithm for parameters a i (column EM1) are good, with a quality quite similar to those obtained by univariate maximum likelihood estimation. Finally, the EM algorithm does not seem very sensitive to initial values, at least if the initial value of 3 is not too small.

Intervention planning

A bivariate Gamma process

X t = X (1) t ; X (2) t
is now used to model the development of two deterioration indicators of a system. We assume that there exist thresholds s i (i = 1; 2) for each indicator, above which the system is considered to be deteriorated.

The system is not continuously monitored but only inspected at time intervals, with a perfect observation of the deterioration level. When one (or both) indicator(s) is (are) observed to be beyond its threshold, an intervention is undertaken. When both indicators are observed to be below their thresholds, no action is undertaken and a new inspection is planned. The time to next inspection ( ) must ensure with a high probability that neither X

(1) t nor X

(2) t go beyond their thresholds s i before the next inspection.

Let (x 1 ; x 2 ) 2 [0; s 1 [ [0; s 2 [ be the observed deterioration level at some inspection time, say at time t = 0 with no restriction. Also, let " 2]0; 1[ be some con dence level.

Using the bivariate Gamma process, the time to the next maintenance action B is chosen as the maximal value which ensures that the process stays below the thresholds with a probability greater than 1 " : B = max 0 such that P (x1;x2) X (1) < s 1 ; X (2) < s 2 1 " :

where P (x1;x2) stands for the conditional probability given X

(1) 0 ; X

(2) 0 = (x 1 ; x 2 ). dependence is ignored, the policy is not as safe as it should be. Also, the optimal time to the next preventive maintenance action is increasing with dependence ( B increases with ), which implies that the error made when considering separate models ( U ) is all the important that the components are less dependent. This also implies that the safest attitude, in case of an unknown correlation, is to consider both components as independent and chose = ? , where ? = max 0 such that P x1 X (1) < s 1 P x2 X (2) < s 2 1 " : Figure 4 shows that taking into account the single information x 1 = 0:4 as presently done at the SNCF may lead to too late maintenance actions. As an example, if x 2 = 0:4, one has B = 134:7 (and (2) = 152:9). The preventive maintenance action based only on NL is consequently scheduled 15 days too late. If x 2 = 0:5, one obtains B = 95:9 ( (2) = 97:5) and the maintenance action is undertaken 54 days too late. If x 2 = 0:6, one obtains B = 47:1 ( (2) = 47:2) and this is 103 days too late.

Application to track maintenance

Concluding this section, one can nally observe that if x 1 is not too close to x 2 , the value U = min (1) ; (2) seams reasonable for maintenance scheduling (see Figure 4), contrary to the currently used (1) , which may entail large delays in its planning (more than 100 days in our example). If x 1 is close to x 2 , the values of U and B have the same order of magnitude, with U > B however, so that the preventive maintenance action is again planned too lately (15 days in the example).

A bivariate Gamma process has been used to model the development of two deterioration indicators. Di erent estimation methods have been proposed for the parameters and tested on simulated data. Based on these tests, the best estimators seem provided by univariate likelihood maximization for the marginal parameters and by an EM algorithm for the correlation coe cient.

Preventive maintenance scheduling has then been studied for a system that deteriorates according to a bivariate Gamma process. In particular, it has been shown that, given an observed bivariate deterioration level, the optimal time to maintenance is increasing with dependence. It has been proven that the optimal time to maintenance is always shorter when taking into account the dependence between the two deterioration indicators than when considering them separately (or only considering one of them).

Finally, a bivariate Gamma process has been used to study a real track maintenance problem. The application shows that when both observed deterioration indicators are close to each other, the bivariate process gives safer results for maintenance scheduling than both univariate processes considered separately or one single univariate process, with the same order of magnitude in each case however. When the observed deterioration indicators are clearly di erent, considering one single univariate process as it is done in current track maintenance, may lead to clearly inadaquate results. The preventive maintenance action is scheduled too late and the track passes to a deteriorated state with a high probability. This application to real data of railway track deterioration hence shows the interest of a bivariate model for a correct de nition of a maintenance strategy.

Figure 1 .

 1 Figure 1. Levelling defects

Figures 2

 2 Figures 2(a), 2(b) and 2(c) show the evolution of a (k) i and (k) along the di erent steps of the EM algorithm (case EM1 in Table2). We notice that the parameters a

Figure 2 .

 2 Figure 2. Evolution of the parameters during the execution of the EM algorithm.

Figure 3 .

 3 Figure 3. B with respect to and U , for the four cases of Table 4

AFigure 4 .

 4 Figure 4. (1) , (2) and B with respect to x 2 with x 1 = 0:4

  Buijs, J.W. Hall, J.M. van Noortwijk and P.B. Sayers in[START_REF] Buijs | Time-dependent reliability analysis of ood defences using gamma processes[END_REF]) is used: a speci c L evy process called bivariate Gamma process is considered. This process is constructed from three independent univariate Gamma processes by trivariate reduction, and has univariate Gamma processes as marginal processes.

	for all x 2 [0; 1), and Laplace transform
				t
		E e s(Yt Y0) =	+ s
	for all s 0.
	Recall that E (Y t Y 0 ) = t , Var(Y t Y 0 ) = t 2 for all t 0 (see (van Noortwijk 2009)
	for more details).
	Following (Buijs et al. 2005), a bivariate Gamma process (X t ) t 0 = X t ; X (1) t (2)	t 0	is
	constructed by trivariate reduction: starting from three independent univariate Gamma
	(i) t processes Y	t 0
				Section 4 is devoted to the study
	of maintenance scheduling and to the comparison of the results based on the bivariate
	and on the univariate models. Finally, a bivariate Gamma process is tted to real data of
	railway track deterioration in Section 5 and it is shown that the maintenance scheduling
	based on the two available deterioration indicators is clearly safer than those based on a
	single one, or on both taken separately.
	2. The bivariate Gamma process
	Recall that an univariate (homogeneous) Gamma process (Y t ) t 0 with parameters > 0
	and > 0 is a process with homogeneous independent increments such that Y t Y 0 is
	Gamma distributed ( t; ) with probability density function (p.d.f.)
		f t; (x) =	t ( t)	x t 1 e x

It is the rst time that both NL and NT indicators are used conjointly to predict the optimal dates of interventions. The objective is to analyze the correlation between the two processes and to determine in what circumstances this bivariate process allows for a better prediction of the intervention times than the current univariate one, based only on the NL indicator.

The paper is organized in the following way: bivariate Gamma processes are introduced in Section 2. Empirical and maximum likelihood estimators for their parameters are provided in Section 3. An EM algorithm is proposed to carry out the maximum likelihood estimation. Both methods are tested on simulated data.

  and (van

	November 29, 2010 10:17 Structure and Infrastructure Engineering swp0000
										7
	Noortwijk and Pandey 2004), with:
	âi bi	=	P n j=1 X j (i) P n j=1 t j	and	âi b2 i	=	P n j=1 j=1 t j P n	X j (i) 1 P n j=1	tj	âi bi P n j=1 ( t j ) 2 : 2 t j

which here becomes:

( 1 t j 1) f 1 j;

(k) 1 ;

(k) 2 ;

(k) 3

+ ( 2 t j 1) f 2 j;

(k) 1 ;

(k) 2 ;

(k) 3

+ ( 3 t j 1) f 3 j;

(k) 1 ;

(k) 2 ;

(k) 3 n X j=1

(ln ( 1 t j ) + ln ( 2 t j ) + ln ( 3 t j )) + h (k)

1 ;

(k) 2 ;

(k) 3

:

The maximization procedure in stage 2 is done by solving the following equation with respect to i : @g 1 ; 2 ; 3 ;

(k) 2 ;

(k) 3

This EM algorithm provides estimates for 1 ; 2 ; 3 using the univariate ML estimates ( b 1 ; b 2 ) and consequently delivers estimates for a 1 ; a 2 and .

In the same way, it is possible to estimate only 3 by an EM algorithm using the univariate ML estimates ( a 1 ; b 1 ; a 2 ; b 2 ). In that case,

is the solution of the equation:

Tests on simulated data

We now test the previous methods on simulated data. 500 time increments (t j ) 1 j 500 are randomly chosen with similar magnitude as the data of track deterioration (the proposed methods will be used on these data in Section 5). Then, 500 values of a bivariate Gamma 1, 2 and3, each corresponding to independent series of data. In these tables, one can nd: the true values in column 2, the empirical estimators in column 3, the univariate maximum likelihood estimators of a 1 ; b 1 ; a 2 ; b 2 in column 4, the EM estimator of the three parameters a 1 ; a 2 ; in column 5, using the parameters b 1 ; b 2 previously estimated by the univariate maximum likelihood method (from column 4), and the second EM estimator of the parameter in column 6, using the estimated parameters a 1 ; b 1 ; a 2 ; b 2 from column 4.

The initial values for the EM algorithm are di erent for the three tables. For Table 1, the EM algorithm has been initiated with We have

where g t is the p.d.f. of X t (see Eq. ( 2)). This quantity is continuous and strictly decreasing in t. If we set Eq. ( 8) equal to 1 ", the time B may hence be numerically computed as the unique solution to this equation in t.

Without a bivariate model, a natural time to next maintenance action for the system is:

with

where P xi stands for the conditional probability given

where F ait;bi (x) is the cumulative distribution function of the distribution (a i t; b i ), the time (i) is computed by solving the equation F ai (i) ;bi (s i x i ) = 1 "; for i = 1; 2, and

Clearly we have :

The goal now is to study the di erence between U and B , and more generally, to understand the in uence of the dependence between both components of X t on B . We rst provide numerical examples. x 1 and x 2 , Table 4 gives the corresponding values for max (as provided by ( 1)) and the resulting (1) , (2) , U and B ( max ). and the resulting max , (1) , (2) , U and B ( max ).

The value of B is plotted with respect to in the Figures 3 for the di erent cases of Table 4, and the corresponding value of U is indicated.

In these gures, one can observe that both B = U and B < U are possible. One can also observe some monotony property of B with respect of the dependence between

X

(1) t and X

(2)

t . More speci cally, we have the following result: In conclusion to this section, one can see that using a bivariate model instead of two separate univariate models generally shortens the time to the next maintenance action ( B U ). If both processes are correlated, the frequency of inspections must then be increased in order to achieve the same safety level as in the uncorrelated case. So, if