Julien Worms
email: julien.worms@uvsq.fr

Rym Worms
email: rym.worms@u-pec.fr

New estimators of the extreme value index under random right censoring, for heavy-tailed distributions

Keywords: Extreme value index, Tail inference, Random censoring, Kaplan-Meier integration Mathematics Subject Classification (2000) 62G32 (Extreme value statistics)

This paper presents new approaches for the estimation of the extreme value index in the framework of randomly censored (from the right) samples, based on the ideas of Kaplan-Meier integration and the synthetic data approach of S. Leurgans (1987). These ideas are developed here in the heavy tail case and for the adaptation of the Hill estimator, for which the consistency is proved under first order conditions. Simulations show good performance of the two approaches, with respect to the only existing adaptation of the Hill estimator in this context.

Introduction

Estimating the extreme value index is an important problem in extreme value statistics. A distribution function (d.f.) F is said to be in the maximum domain of attraction of H γ (noted F ∈ D(H γ )) with H γ (x) := exp -(1 + γx) -1/γ for γ = 0 and 1 + γx > 0 exp(-exp(-x)) for γ = 0 and x ∈ R , if there exist two normalizing sequences (a n ) ⊂ R + and (b n ) ⊂ R such that (for every x ∈ R)

F n (a n x + b n ) n→∞ -→ H γ (x).
If we observe a sample (X i ) i≤n with common distribution function F ∈ D(H γ X ), with γ X > 0, a classical estimator of the extreme value index γ X is the so-called Hill estimator γX,Hill :=

1 k n kn i=1 log X n-i+1,n X n-kn,n
where X 1,n ≤ • • • ≤ X n,n are the ascending order statistics associated to the X-sample and k n the sample fraction to keep from this sample. However, in a certain number of applications, such as survival analysis, reliability theory, insurance . . ., the variable of interest X is not necessarily completely available. This is the case in the presence of random right censoring. Examples of censored data with apparent heavy tails can be found in [Gomes and Neves (2011)] and [Einmahl et. al. (2008)].

The usual way to model this situation is to introduce a random variable C, independent of X, such that only Z = X ∧ C and δ = I X≤C are observed. The observed variable δ determines whether X has been censored or not. It is common sense that any classical estimator of the extreme value index (such as the Hill or the Moment estimator) is not consistent for estimating γ X if it is naively computed from the Z-sample (indeed, it estimates the extreme value index associated to the Z-sample, denoted by γ in the sequel).

Recently, [Beirlant et. al. (2007)] and [Einmahl et. al. (2008)] proposed an adaptation of classical extreme value index estimators in the case of right random censoring, therefore providing (to the best of our knowledge) the first methodological papers on this subject; their method will be presented in subsection 2.1.

In this paper, we propose two different approaches to deal with the estimation of γ X , relying on more natural heuristics in this randomly censored sample framework : one of these amounts to consider Kaplan-Meier integrals, and the other on ideas coming from censored regression. Given the combination of difficulties coming from extreme values and censoring, we will restrict ourselves to the application of these approaches to the adaptation of the Hill estimator, in the heavy tailed case, and to the consistency of these adapted estimators. It is however more than likely that our ideas adapt to other, more efficient, estimators of the extreme value index, and to other domains of attraction.

In Subsection 2.1, we define the framework more precisely and recall the existing methodology cited above. In Subsections 2.2 and 2.3, we present our methodology and state the consistency results for the adaptation of the Hill estimator. Section 3 is devoted to a small simulation study, and Section 5 to the proofs. A conclusion is provided in Section 4.

Notations : in the whole paper, the sign := denotes an equality defining the quantity on the left side, f (t -) denotes lim s↑t f (s), and f ← the general inverse of the function f . I A (x) is equal to 1 if x ∈ A and 0 if x ∈ A. The definition of a regularly varying function f of order α (noted f ∈ RV α ) is recalled in the Appendix.

Methodology

The framework and a general existing methodology

We consider in this paper two independent i.i.d. non-negative samples (X i ) i≤n and (C i ) i≤n with respective continuous distribution functions F and G (with end-points τ F and τ G , where τ F := sup{x, F (x) < 1}). In the context of randomly right-censored observations, one actually only observes, for i ≤ n,

Z i = X i ∧ C i and δ i = I Xi≤Ci .
We denote by H the distribution function of the Z-sample, satisfying

1 -H = (1 -F )(1 -G) and by Z 1,n ≤ • • • ≤ Z n,
n the associated order statistics. In the whole paper, δ 1,n , . . . , δ n,n denote the δ's corresponding to Z 1,n , . . . , Z n,n , respectively ([Stute (1995)] call them "concomitant" to the order statistics).

If F and G are assumed to be in the maximum domains of attraction D(H γ X ) and D(H γ C ) respectively, where γ X and γ C are real numbers, then this implies that H ∈ D(H γ ), for some γ ∈ R. [Einmahl et. al. (2008)] considered the following three most interesting cases : case 1:

γ X > 0 , γ C > 0 in this case γ = γ X γ C γ X + γ C case 2: γ X < 0 , γ C < 0 , τ F = τ G in this case γ = γ X γ C γ X + γ C case 3: γ X = γ C = 0 , τ F = τ G = +∞ in this case γ = 0.
The general existing method, appeared first in [Beirlant et. al. (2007)] and developed in [Einmahl et. al. (2008)], is to consider any consistent estimator γZ of the extremal index γ applied to the Z-sample and divide it by the proportion p of non-censored observations in the tail (i.e. in the k n largest observations of the Z-sample). That is, an adaptation of an extreme value index estimator in the presence of random right censoring is :

γc X := γZ p , where p := 1 k n kn j=1 δ n-j+1,n . (1)
It is proved in [Einmahl et. al. (2008)] that p consistently estimates p := γ C γ X +γ C and therefore γc X consistently estimates γ/p = γ X (obtaining as well the asymptotic normality, if it holds for the sequence γZ ). This method provides flexibility as to the choice of the estimator of γ Z . Up to now, the only alternative to it, in this context, is the adaptation of the ML estimator based on the excesses over a threshold, developed in [Beirlant et. al. (2010)], which we shall however not detail here.

We now present an alternative path for estimating the extreme value index γ X , based on ideas which are well-known in the survival analysis literature.

First approach

The starting point of the first new approach is the following well known result, which is the basis of censored regression methods (for instance, an early reference is [Koul et. al. (1981)]) : if φ is some nonnegative real function,

E δ 1 -G(Z) φ(Z) = E[φ(X)]. (2)
It is readily proved : since Z = X when δ = 1,

E δ 1 -G(Z) φ(Z) = I x≤c φ(x) 1 -G(x) dF (x) dG(c) = φ(x)(1 -G(x)) -1 +∞ x dG(c) dF (x) = φ(x) dF (x) = E[φ(X)].
In the context of extreme value statistics, the idea is to take advantage of this property and of the fact that some tail parameters of the distribution of X can be approached by the expectation of some function of X, therefore opening the way to their estimation. In this paper we will illustrate it in the context of heavy-tailed distributions, and for the estimation of the extreme value index, assuming that we are in the first of the three situations presented in paragraph 2.1 :

F ∈ D(H γ X ) , G ∈ D(H γ C ) with γ X > 0 and γ C > 0 (3)
which, as noted earlier, implies that H ∈ D(H γ ) with

γ = γ X γ C γ X + γ C .
In this case, it is well known that (see Remark 1.2.3 in [Haan and Ferreira (2006)] for instance)

lim t→+∞ E [log(X/t) | X > t] = γ X . (4)
If (k n ) is a sequence of integers satisfying, as n tends to +∞,

k n → +∞ and k n = o(n) (5)
then we can define a random version of φ(x) = (1 -F (t)) -1 log(x/t)I x>t , with random threshold t = Z n-kn,n φn (x) := 1

1 -Fn (Z n-kn,n ) log x Z n-kn,n I x>Z n-kn,n . (6)
Consequently, by combining (2) and ( 4) with this function φn , our first adaptation of the Hill estimator comes, valid in situation (3), γKM X,Hill :=

1 n(1 -Fn (Z n-kn,n )) kn i=1 δ n-i+1,n 1 -Ĝn (Z - n-i+1,n ) log Z n-i+1,n Z n-kn,n , (7)
where Fn and Ĝn (the Kaplan-Meier estimators of F and G, respectively) are defined as follows : for t < Z n,n ,

1-Fn (t) = Zi,n≤t n -i n -i + 1 δi,n and 1-Ĝn (t) = Zi,n≤t n -i n -i + 1 1-δi,n
.

Note that we take Ĝn (Z - n-i+1,n ) instead of Ĝn (Z n-i+1,n ), in the definition of γKM X,Hill , because 1 -Ĝn (Z n,n ) can be zero or undefined.

The following theorem provides the consistency of this estimator. For this purpose, we need two additional assumptions on the behavior of the function p • H ← : if p(z) := P(δ = 1|Z = z),

1 k n kn i=1 p H ← 1 - i n -p P -→ c ∈ R (8) sup (s,t)∈Cn |p(H ← (t)) -p(H ← (s))| → 0, for all C > 0 (9)
where

C n = {(s, t) such that s < 1 , 1 -k n /n ≤ t < 1 , |t -s| ≤ C √ k n /n }.
Theorem 1 Under assumptions (3), ( 5), ( 8), ( 9), if we additionally assume that

-log(k n /n) k n = O(n -δ ), ( 10
)
for some δ > 0 if γ X < γ C and for some δ ≥ γ X -γ C γ X + γ C if γ X ≥ γ C , then, as n tends to +∞, γKM X,Hill P -→ γ X .
Remark 1 Inspection of the proof (see the treatment of P n at the end of subsection 5.1.2) reveals that any sequence (k n ) such that k n ≥ cn a (for some a ∈]0, 1[, some constant c > 0, and n large) is suitable for the theorem to hold (whether γ X < γ C or not, and without assuming (10)).

Remark 2 Assumptions (8) and (9) have similarities with (but are weaker than) those used in [Einmahl et. al. (2008)]. The latter have been proved, in [Brahimi et. al. (2013)], to be unnecessary for obtaining the asymptotic normality of the adapted Hill estimator γc X,Hill := γZ,Hill /p (see (1)). The empirical process techniques used in [Brahimi et. al. (2013)] do not seem to be applicable in our setting, since our estimator is a nonconstantly-weighted version of the Hill estimator

γKM Z,Hill = 1 k n kn i=1 1 -Ĝn (Z n-kn,n ) 1 -Ĝn (Z n-i+1,n ) δ n-i+1,n log Z n-i+1,n Z n-kn,n whereas γc Z,Hill = 1 k n kn i=1 1 p log Z n-i+1,n Z n-kn,n
Remark 3 We have made the choice of a random threshold and fixed number of relative excesses, because it seemed closer to what is done in practice and, secondly, the other path went with its own difficulties. This other choice was to consider a deterministic threshold t n , and then (up to the estimation of F at t n , which causes no problem) write γKM

X as a proper Kaplan-Meier integral φ n (x) d Fn (x), with deterministic φ n (x) = (1-F (t n )) -1 log(x/t n )I x>tn (
with our choice, the function φn is random). However, this function φ n is intrinsically unbounded, has a "sliding towards infinity" support, and is dependent on n : we found no way to deal with this, using the Kaplan-Meier integration tools known in the literature.

Second approach

Our second approach, alternative to the Kaplan-Meier integral approach presented in the previous paragraph, is based on ideas of [Leurgans (1987)], who developed a "synthetic data" strategy in censored regression problems (see [START_REF] Delecroix | [END_REF]] for a more recent reference to this method). The starting point of this second approach is the following result :

if φ and ψ are two nonnegative R + → R functions such that x 0 ψ(t) dt = φ(x), then E Z 0 ψ(t) 1 -G(t) dt = E[φ(X)]. (11)
Indeed,

E Z 0 ψ(t) 1 -G(t) dt = +∞ 0 +∞ t dH(z) ψ(t) 1 -G(t) dt = +∞ 0 (1 -F (t))ψ(t) dt = +∞ 0 +∞ t dF (x) ψ(t) dt = +∞ 0 x 0 ψ(t) dt dF (x) = E[φ(X)]. Consequently, if φ(x) = x 0 ψ(z) dz, then 1 n n i=1 Zi 0 ψ(z) 1 -Ĝn (z -) dz should correctly estimate E[φ(x)]
. This estimator can be rewritten using the special form of function ψ and the piecewise constant form of the Kaplan-Meier estimator Ĝn : noting Z 0,n = 0, and rk(Z i ) the (ascending order) rank of the observation Z i in the Z-sample, we have indeed

Zi 0 ψ(z) 1 -Ĝn (z -) dz = rk(Zi) j=1 ]Zj-1,n,Zj,n] ψ(z) 1 -Ĝn (z -) dz = rk(Zi) j=1 1 1 -Ĝn (Z j-1,n ) Zj,n Zj-1,n ψ(z) dz = rk(Zi) j=1 φ(Z j,n ) -φ(Z j-1,n ) 1 -Ĝn (Z j-1,n ) .
Considering, once again, the function φn introduced in (6), we can now define our second new adaptation of the Hill estimator, valid in situation (3) γLeurg X,Hill :=

1 n n i=1 i j=1 φn (Z j,n ) -φn (Z j-1,n ) 1 -Ĝn (Z j-1,n )
which turns out to be, after some simplifications,

γLeurg X,Hill = 1 n(1 -Fn (Z n-kn,n )) kn i=1 1 1 -Ĝn (Z - n-i+1,n ) i log Z n-i+1,n Z n-i,n .
(12) We note that, while γKM X,Hill appeared as a weighted version of the classical form of the Hill estimator (mean of the log relative excesses log(Z n-i+1,n /Z n-kn,n )), our second candidate γLeurg X,Hill is a weighted version (but with weights which are always non null) of the mean of the so-called log spacings i log(Z n-i+1,n /Z n-i,n ), i.e. the other form of the Hill estimator.

The following theorem provides the consistency of this estimator, under less restrictive conditions than Theorem 1.

Theorem 2 Under assumptions (3), ( 5) and (10), then, as n tends to +∞, γLeurg X,Hill P -→ γ X .

Finite sample behaviour

In this section, we present some graphs (issued from an extensive study in the heavy tail framework) corresponding to the comparison, in terms of observed bias and mean squared error (MSE) of our new estimators γKM X,Hill and γLeurg X,Hill (defined by ( 7) and ( 12)) with other adapted estimators for γ X (see (1)) : the adapted Hill estimator γc X,Hill and the adapted version of the biascorrected Hill estimator (introduced and studied in [START_REF] Caeiro | [END_REF]] and [Gomes et. al. (2008)]) γc X,M VRB , where MVRB stands for minimum-variance reduced-bias.

We consider two classes of heavy-tailed distributions :

-Burr(β, τ, λ) with d.f. 1 -( β β+x τ ) λ , which extreme value index is 1 λτ . -Frechet(γ) with d.f. exp(-x -1/γ ), which extreme value index is γ.

For each considered distribution, 2000 random samples of length n = 500 were generated ; median bias and MSE of the four above-mentioned estimators are plotted against different values of k n , the number of excesses used.

We considered three cases : a Burr distribution censored by another Burr distribution (Fig. 1), a Fréchet distribution censored by another Fréchet distribution (Fig. 2) and a Fréchet distribution censored by a Burr distribution (Fig. 3). In each case, we considered a situation with γ X < γ C (subfigure (a)), which corresponds to a weak censoring in the tail, and the reverse situation with γ X > γ C (subfigure (b)), which corresponds to a strong censoring.

It seemed more natural to present separately, below, the comparison of our estimators γKM X,Hill and γLeurg X,Hill with γc X,Hill on one hand, and with γc X,M VRB on the other hand.

From the three situations presented above, it seems that our new estimators perform better than the former adapted Hill estimator, in the weak censoring case, both in term of bias and MSE. It is not surprising that, in the strong censoring case, results become worse for all estimators but (clearly) to a lesser extent for ours. Moreover, in the strong censoring situation, γLeurg X,Hill seems to have systematically the best behavior. Other simulations not presented here confirm this phenomenon.

Turning now to the comparison of our estimators with the adaptation of the MVRB estimator, first note that intuitively this should end with a better performance for the latter, since in the uncensored situation it possesses better theoretical and empirical properties than the Hill estimator, on which our estimators are based. It is indeed the case when the bias-reduction is successful (as in Figure 2(a)), but surprisingly, our adaptations of the Hill estimator

Conclusion

In this paper, we have introduced two new approaches for the estimation of the extreme value index in the case of randomly censored observations, based on ideas coming from the censored regression literature. The estimation problem of the e.v.i. in the censoring framework had been addressed in very few papers before, and our methodology, though we have applied it to the adaptation of the Hill estimator in the heavy tail case only, has some potential for more applications, either for other estimators (moment based estimators on first place) and other maximum domains of attraction, or maybe for the estimation of other tail parameters. This work thus forms a basis for future research in this recently studied area of censored extremes, which can prove much useful in applications, as was showed in the review paper [Gomes and Neves (2011)]. For the moment, technical problems prevent us from obtaining asymptotic normality results and rigorous evaluation of the variance, but simulations show that our two versions of the Hill estimator perform quite well with respect to the existing version, in terms of MSE, even in the apparently less favorable case of heavy censoring in the tail (γ X > γ C ).

Proofs

First, note that in several occasions in the next pages, reference will be made to Proposition 1 : it is stated in the Appendix.

Proof of Theorem 1

Since (1 -Fn (t))(1 -Ĝn (t)) = k n /n for t = Z n-kn,n , by introducing

A n := 1 -Ĝn (Z n-kn,n ) 1 -G(Z n-kn,n ) and C in := 1 -G(Z n-i+1,n ) 1 -Ĝn (Z - n-i+1,n ) we can write γKM X,Hill = A n 1 k n kn i=1 C in W in ,
where

W in := δ n-i+1,n log Z n-i+1,n Z n-kn,n 1 -G(Z n-kn,n ) 1 -G(Z n-i+1,n ) .
Therefore, we have the decomposition γ KM X,Hill = A n (W n + R n ) where

W n := 1 k n kn i=1 W in and R n := 1 k n kn i=1 (C in -1)W in .
First of all, relying on Theorem 2 in [Csörgő (1996)], continuity of G entails that A n P -→ 1, as n → +∞ (please note that, unfortunately, this theorem will not be sufficient for controlling the quantities C in , because they involve really extreme observations ; in section 5.1.2 we show how this difficulty is circumvented).

Consequently, it remains to prove that W n P -→ γ X and R n P -→ 0 : this is the purpose of the next two subsections.

Proof of W

n P -→ γ X
Let us first introduce the following notation, used throughout the rest of the proof, Zi,n := Z n-i+1,n Z n-kn,n .

Under assumption (3) (which implies that 1 -G ∈ RV -1/γ C and 1 -H ∈ RV -1/γ ), setting > 0, we can apply Potter bounds (29) stated in Proposition

1, to the function 1/(1 -G) ∈ RV γ -1 C
, and to t = Z n-kn,n P -→ +∞ and

x = Zi,n ≥ 1. We thus obtain for n sufficiently large, the following bounds for W in :

(

1-)δ n-i+1,n log( Zi,n )( Zi,n ) γ -1 C -≤ W in ≤ (1+ )δ n-i+1,n log( Zi,n )( Zi,n ) γ -1 C + . (13)
Therefore, it remains to prove that both the mean of the lower bound, and that of the upper bound, converges in probability to a quantity arbitrary close to γ X when is taken close to 0. We consider the case of the upper bound only, the lower bound being similar.

Recall that

p := γ C γ X + γ C
which is (see [Einmahl et. al. (2008)]) the limit of p(z) = P(δ = 1|Z = z) when z → ∞. We intend to rely on the closeness of the δ n-i+1,n to i.i.d. Bernoulli(p) random variables, independent of the log-spacings Zi,n .

Mimicking what was proposed in [Einmahl et. al. (2008)], we use the fact that the original (Z i , δ i ) i≤n are identically distributed as (Z i , δ i ) i≤n , where δ i = I Ui≤p(Z i ) and (U i ) i≤n denotes an i.i.d. sequence of standard uniform variables, independent of a given sequence (Z i ) i≤n of i.i.d. variables having H as their c.d.f. We thus carry on the proof by considering now that δ i is related to Z i by δ i = I Ui≤p(Zi) .

We then define (where the U i below is the same as the one in the above definition of δ i ) δi = I Ui≤p which are Bernoulli(p) distributed and independent of the sequences (Z i ) and ( Zi,n ). Note that we define ( δ1,n , . . . , δn,n ) as the rearrangements of the δi corresponding to the order induced by the order statistics (Z 1,n , . . . , Z n,n ) : these are however still independent of the sequences (Z i ) or ( Zi,n ).

According to (13), it comes

(1 + ) -1 W n ≤ 1 k n kn i=1 δn-i+1,n log( Zi,n )( Zi,n ) γ -1 C + + 1 k n kn i=1 (δ n-i+1,n -δn-i+1,n ) log( Zi,n )( Zi,n ) γ -1 C + =: J 1 n + J 2 n ( 14
) (i) Let us first prove that J 1 n P -→ g( ), with lim ↓0 g( ) = γ X .
Let Y 1,n , . . . , Y n,n be the ascending order statistics of n i.i.d standard Pareto random variables Y 1 , . . . , Y n with distribution function 1 -1/x, for x > 1. By independence of the δj,n and the Zi,n , if U is the quantile function associated to

H (i.e. U (t) = H ← (1 -1/t)) then Zi,n log( Zi,n ) d = U (Y n-i+1,n ) U (Y n-kn,n ) log U (Y n-i+1,n ) U (Y n-kn,n ) . ( 15
)
Under assumption (3), U is regularly varying with index γ and Proposition 1 can be applied to U . Taking for each 1 29) and their logarithm yield, for some given > 0 and n sufficiently large,

≤ i ≤ k n , t = Y n-kn,n and x = Y n-i+1,n /Y n-kn,n , bounds (
(1-) γ -1 C + H in ≤ U (Y n-i+1,n ) U (Y n-kn,n ) γ -1 C + log U (Y n-i+1,n ) U (Y n-kn,n ) ≤ (1+ ) γ -1 C + K in (16
) where , setting α -:= (γ -)(γ -1 C + ) and α + := (γ + )(γ -1 C + ),

H in = Y n-i+1,n Y n-kn,n α - log(1 -) + (γ -) log Y n-i+1,n Y n-kn,n K in = Y n-i+1,n Y n-kn,n α + log(1 + ) + (γ + ) log Y n-i+1,n Y n-kn,n
However, it is known that

(Y n-i+1,n /Y n-kn,n ) 1≤i≤kn d = ( Ỹi,kn ) 1≤i≤kn , (17)
where Ỹ1,kn , . . . , Ỹkn,kn are the ascending order statistics of k n i.i.d random variables Ỹ1 , . . . , Ỹkn with standard Pareto distribution. Thanks to the independence of the δj,n and the Y i , it follows that J 1 n is bounded above by some variable which equals in distribution

1 k n kn i=1 I Ui≤p Ỹ α + i (log(1 + ) + (γ + ) log( Ỹi )). ( 18
)
Independence of the (U i ) and ( Ỹi ) and the law of large numbers then yields that (since 0 < α + < 1) lim sup J 1 n is bounded above, in probability, by

p × log(1 + ) 1 1 -α + + (γ + ) 1 (1 -α + ) 2 .
Dealing now with the lower bound of ( 16), one can handle H in similarly and obtain a lower bound in probability for lim inf J 1 n : straightforward computations show that both bounds converge to γ X as and tend to 0. This thus concludes part (i).

(ii) It remains to prove that J

2 n P -→ 0. Let us put T n-i+1,n := log( Zi,n )( Zi,n ) γ -1
C + and recall that

J 2 n = 1 k n kn i=1 (δ n-i+1,n -δn-i+1,n )T n-i+1,n .
Let p > 1 and q > 1 such that 1 p + 1 q = 1. Hölder's Inequality gives :

|J 2 n | ≤ 1 k n kn i=1 δ n-i+1,n -δn-i+1,n q 1/q 1 k n kn i=1 (T n-i+1,n ) p 1/p = 1 k n kn i=1 δ n-i+1,n -δn-i+1,n 1/q 1 k n kn i=1 (T n-i+1,n ) p 1/p
So it remains to prove that 1 kn kn i=1 |δ n-i+1,n -δn-i+1,n | = o P (1) and that

1 kn kn i=1 (T n-i+1,n ) p = O P (1)
, for an appropriate p > 1. On one hand, according to (15), ( 16) and ( 17) (as in part (i)), we have to prove that

1 k n kn i=1 Ỹ α + i (log(1 + ) + (γ + ) log( Ỹi )) p = O P (1)
This is the case, using the law of large numbers, as soon as we take p

< 1 + γ C γ X (in this case α + p < 1, so E(( Ỹ α + 1 log( Ỹ1 )) p ) is finite).
On the other hand, recall that δ i,n = I Ui≤p(Zi) and δi,n = I Ui≤p . Then,

1 k n kn i=1 δ n-i+1,n -δn-i+1,n ≤ 1 k n kn i=1 I Ui≤p(Zn-i+1,n) -I Ui≤p(H ← (1-i n )) + 1 k n kn i=1 I Ui≤p(H ← (1-i n )) -I Ui≤p =: T 1,k + T 2,k .
Following the same lines as in [Einmahl et. al. (2008)] (p218) for the treatment of their terms T 1,k and T 2,k (which are different from ours in the rates but nevertheless quite similar), we show using a result of [Chow and Teicher (1997)] (p356), that both T 1,k and T 2,k tend to 0 in probability, thanks to assumptions (8) and ( 9). This concludes the proof of (ii) and therefore of Wn

P -→ γ X . 5.1.2 Proof of R n P -→ 0
As mentioned at the beginning of the proof, difficulties arise as to the control of the Kaplan-Meier estimate of G in the tail (here it takes the form of the variables C i,n , not to be mistaken with the censoring variables C i ) : we will circumvent them via a device known in the survival analysis literature. Let us define, for some > 0,

C(t) := t 0 dG(x) (1 -G(x)) 2 (1 -F (x))
and

h in := ( C(Z n-i+1,n )) -1 2 -.
We readily have

|R n | ≤ T 1 n T 2 n ,
where

T 1 n := sup 1≤i≤kn √ n|h in (C in -1)| and T 2 n := 1 k n kn i=1 W in h -1 in n -1 2
and we are going to prove that T 1 n = O P (1) and

T 2 n = o P (1). First remind that C i,n is the value of the function t → (1 -G(t))/(1 - Ĝn (t -)) at t = Z n-i+1
,n , and consequently, by continuity of G and C [Gill (1983)] applies and therefore the process

T 1 n ≤ sup t≤Zn,n √ n( C(t)) -1 2 -Ĝn (t -) -G(t) 1 -Ĝn (t -) ≤ sup t<Zn,n √ n( C(t)) -1 2 -Ĝn (t) -G(t) 1 -Ĝn (t) . Since +∞ 0 h 2 (t)d C(t) < ∞ for the function h(t) = ( C(t)) -1 2 -, Theorem 2.1 of
√ n( C(t)) -1 2 -Ĝ(t) -G(t) 1 -Ĝ(t) t<Zn,n converges in distribution. As a consequence, T 1 n = O P (1). It remains to prove that T 2 n = o P (1). First, from the definition of h in and C, since (1 -H) = (1 -F )(1 -G) we clearly have h -1 in < -log(1 -G(Z n-i+1,n )) 1 -H(Z n-i+1,n ) 1 2 + .
Moreover, according to (13), for some given > 0, we have T 2 n ≤ (1 + )P n Q n (n large) where

P n := n -1 2 -log(1 -G(Z n-kn,n )) 1 -H(Z n-kn,n ) 1 2 + Q n := 1 kn kn i=1 1-H(Z n-kn ,n ) 1-H(Zn-i+1,n) 1 2 + log(1-G(Zn-i+1,n)) log(1-G(Z n-kn,n )) 1 2 + log( Zi,n ) Zi,n γ -1 C +
Under assumption (3), 1/(1 -H) is regularly varying with index 1/γ and -log(1 -G) is slowly varying : therefore, for some given > 0 and n sufficiently large, the application of Proposition 1 to these functions yields the following upper bounds

1 -H(Z n-kn,n ) 1 -H(Z n-i+1,n ) ≤ (1 + ) Zγ -1 + i,n log(1 -G(Z n-i+1,n )) log(1 -G(Z n-kn,n )) ≤ (1 + ) Z i,n
and consequently

Q n ≤ (1 + ) 1+2 1 k n kn i=1 Zβ i,n log( Zi,n ), (19)
where β = (2γ) -1 + γ -1 C + , for some > 0 (arbitrarily small when and are closer to 0).

Concerning P n , we have (n/k n )(1 -H(Z n-kn,n )) P -→ 1 and therefore 1 -H(Z n-kn,n ) = k n n U n , and -log(1 -H(Z n-kn,n )) = -log k n n V n
for some sequences U n and V n tending to 1, in probability, as n → +∞. Using the inequality 1 -H ≤ 1 -G, it follows that, for some W n P -→ 1,

P n ≤ n -1 2 {-log(k n /n)/(k n /n)} 1 2 + W n = n {-log(k n /n)/k n } 1 2 + W n .
(20) In order to prove that P n Q n = o P (1), we have to distinguish the case γ X < γ C from the case γ X ≥ γ C (respectively weak and strong censoring in the tail).

(i) Case γ X < γ C Proceeding as in subsection 5.1.1 by using ( 15), ( 16) and ( 17), we obtain that the upper bound of Q n in ( 19) is O P (1) via the law of large numbers as soon as E( Ỹ1 β(γ+ ) log( Ỹ1 )) is finite, i.e. β(γ + ) < 1 (which is the case since γ X < γ C ). This proves Q n = O P (1). Finally, under assumption (10) on k n , (20) implies that P n P -→ 0 as soon as < δ/(4 -2δ). This concludes the proof of

P n Q n = o P (1) in this case. (ii) Case γ X ≥ γ C
We use, as in (i), equations ( 15), ( 16) and ( 17) to treat Q n . But, since β(γ + ) is not < 1 but ≥ 1, the upper bound of Q n in ( 19) is no longer O P (1). In this case, we rely on the Marcinkiewicz-Zygmund law of large numbers to obtain

1 k q n kn i=1 Ỹi β(γ+ ) log( Ỹi ) = o P (1),
where q := β(γ + ) + δ, for some δ > 0. This proves that Q n = o P (k q-1 n ). It remains to prove that k q-1 n P n = O P (1). From (20) and condition (10) on k n , we have

k q-1 n P n ≤ k n n q-1 n +q-1-δ 2 -δ .
So, and δ being arbitrary close to 0, k q-1 n P n = o P (1) as soon as δ ≥ γ X -γ C γ X +γ C , thanks to assumption (5). Note that if k n is of the order of n a , for some a ∈]0, 1[, the control of P n does not require condition (10).

Proof of Theorem 2

Similarly to γKM X,Hill , the estimator γLeurg X,Hill can be written as follows :

γLeurg X,Hill = A n (W n + R n ),
where A n ,C in and R n are defined as before (see Subsection 5.1) but now

W in := ξ in 1 -G(Z n-kn,n ) 1 -G(Z n-i+1,n ) and ξ in := i log Z n-i+1,n Z n-i,n .
Recall that A n P -→ 1 (see beginning of the proof of Theorem 1). In Section 5.2.2, we prove that R n P -→ 0. Let us first deal with W n .

Proof of W

n P -→ γ X
Recall first the notation used in Subsection 5.1 Zi,n := Z n-i+1,n Z n-kn,n .

Let η > 0 and > 0. Using Potter bounds (29) for 1/(1-G), which is regularly varying of order γ -1 C , from the definition of the W in we first obtain for small enough

P(W n -γ X > η) ≤ P(k -1 n kn i=1 Zγ -1 C + i,n ξ in -γ X > η 2 ), P(γ X -W n > η) ≤ P(γ X -k -1 n kn i=1 Zγ -1 C - i,n ξ in > η
2 ). Let us now consider constants c > 1 and c < 1, both close to 1, and α > 0 and α > 0 both close to γ/γ C , which values will be specified in the proof of Lemma 1 below. Using positivity of ξ in , it comes

P(W n -γ X > η) ≤ P   max i≤kn Zγ -1 C + i,n cu -α i > 1   +P c k -1 n kn i=1 u -α i ξ in -γ X > η 2 (21) P(γ X -W n > η) ≤ P   min i≤kn Zγ -1 C - i,n c u -α i < 1   +P γ X -c k -1 n kn i=1 u -α i ξ in > η 2 (22)
where

u i = u in := i/(k n + 1), for 1 ≤ i ≤ k n .
Convergence in probability of W n to γ X thus comes from the combination of ( 21), ( 22) and the following two Lemmas, by letting go to 0 in the end : Lemma 1 is applied with θ = γ -1 C + and θ = γ -1 C -, whereas Lemma 2 is applied twice with a = α and a = α (both close to γ/γ C which is < 1), noticing that γ/(1 -γ/γ C ) equals γ X .

Lemma 1 Let θ and θ > 0. There exist constants c > 1, c < 1 both arbitrarily close to 1, and α > 0, α > 0, arbitrarily close to γθ and to γθ respectively, such that

lim n→∞ P max i≤kn Zθ i,n cu -α i > 1 = lim n→∞ P min i≤kn Zθ i,n c u -α i < 1 = 0. Lemma 2 If 0 < a < 1, then 1 k n kn i=1 u -a i ξ in P -→ γ 1 -a . If a > 1 , then for any δ > 0, 1 k a+δ n kn i=1 u -a i ξ in P -→ 0.
Lemmas 1 and 2 are proved one after the other below. Note that Lemma 2 is proved by using what could be called a first order version of the techniques used in [Beirlant et. al. (2007)] (where a second order approximation of the log-spacings ξ i,n is obtained).

Proof of Lemma 1

Once again (as in Subsection 5.1.1), we introduce Y 1,n , . . . , Y n,n the ascending order statistics of n i.i.d standard Pareto random variables, in order to have Zi,n

d = U (Y n-i+1,n ) U (Y n-kn,n ) .
Applying Potter bounds (29) to U ∈ RV γ , it comes, for some given > 0, and n large enough,

(U (Y n-i+1,n )/U (Y n-kn,n )) θ ≤ (1 + ) θ Y n-i+1,n Y n-kn,n α (U (Y n-i+1,n )/U (Y n-kn,n )) θ ≥ (1 -) θ Y n-i+1,n Y n-kn,n α
where α = (γ + )θ and α = (γ -)θ .

Using (17) and introducing

V 1,kn ≤ • • • ≤ V kn,
Zθ i,n cu -α i > 1 ≤ P (1 + ) θ max i≤kn (V i,kn /u i ) -α > c , P min i≤kn Zθ i,n c u -α i < 1 ≤ P (1 -) θ min i≤kn (V i,kn /u i ) -α < c .
Relying on max 1≤i≤kn |V i,kn -u i | P -→ 0 (uniform consistency of the uniform empirical quantile process), we readily have, for any given β > 1,

lim n→∞ P(E n ) = 1 where E n := { ∀ 1 ≤ i ≤ k n , β -1 u i ≤ V i,kn ≤ βu i } (23)
Therefore, for given constants > 0 and β > 1, if we set c = β α (1 + ) θ and c = β -α (1 -) θ , then both probabilities appearing in Lemma 1 are bounded above by 1-P(E n ), which achieves the proof in view of (23). As announced, by choosing appropriate values of and β, the constants c, c , α, α , are respectively arbitrarily close to 1, 1, γθ and γθ .

Proof of Lemma 2

We proceed very similarly as in [Beirlant et. al. (2002)], therefore some details will be ommited. Let E 1 , . . . , E n be i.i.d. Exp(1) random variables.

Then Zi,n d = U (exp(E n-i+1,n ))/U (exp(E n-i,n )). The first order Potter-type bounds for U ∈ RV γ stated in Proposition 1 thus yield : for some given > 0, n large enough and 1 ≤ i ≤ k n ,

ξ i,n = i log Z n-i+1,n Z n-i,n d = γ i(E n-i+1,n -E n-i,n ) + i(B k,n (i) -B k,n (i + 1)) d =: γ ξ i + β i,n (24)
where log(1-)-(E n-i+1,n -E n-kn,n ) ≤ B k,n (i) ≤ log(1+ )+ (E n-i+1,n -E n-kn,n )

and the Rényi representation was used to derive (24), with ξ 1 , . . . , ξ kn denoting independent Exp(1) variables.

Using the fact that (E n-i+1,n -E n-kn,n ) 1≤i≤kn d = (-log V i,kn ) 1≤i≤kn , where V 1 ,. . . ,V kn are independent standard uniform variables, and using relation (23), we obtain (this is in fact the first order version of Theorem 2.1 in [Beirlant et. al. (2002)]), uniformly in i ∈ {1, ..., k n }, kn j=i βj,n j

= |B k,n (i)| = o P (log + (1/u i )) (25)
(with log + (x) = max(1, log x) and recalling that u i stands for i/(k + 1)). Regular application of the law of large numbers for triangular arrays of independent random variables yields (cf [Chow and Teicher (1997)] ; details are ommited), when 0 < a < 1 Suppose now that a > 1 and let us prove similarly (27) : we need to be more cautious since u 0 s -a ds is no longer defined. We have n is bounded by (1 + )P n Q n , for some > 0, where P n is defined as in 5.1.2 but, in the definition of Q n , the factor log Zi,n needs to be replaced by ξ in . The same arguments as before allow us to write, for some given , > 0 and n sufficiently large,

Q n ≤ (1 + ) (1+2 ) 1 k n kn i=1 Zβ i,n ξ in ,
where, as in Subsection 5.1.2, β = (2γ) -1 +γ -1 C + , for some > 0 (as small as needed when and are set close to 0). We now proceed as in Subsection 5.2.1 to control T 2 n and end the proof of R n = o P (1).

  Fig. 1 Comparison of bias and MSE for γKM X,Hill (solid), γLeurg X,Hill (thick), γc X,Hill (dotted) and γc X,M VRB (dashed) for a Burr distribution censored by another Burr distribution : (a) γ X = 1/2 and γ C = 1 (weak censoring), (b) γ X = 1 and γ C = 1/2 (strong censoring)

  Fig. 2 Comparison of bias and MSE for γKM X,Hill (solid), γLeurg X,Hill (thick) γc X,Hill (dotted) and γc X,M VRB (dashed) for a Fréchet distribution censored by another Fréchet distribution : (a) γ X = 1/2 and γ C = 1 (weak censoring), (b) γ X = 1 and γ C = 1/2 (strong censoring)

  0 < a < 1. The trick to show the negligibility (26) is to write (where u 0 = 0)

  s -a ds + |1 -a|(k n + 1) proof is identical to the case of the first theorem. As in Subsection 5.1.2, we have |R n | ≤ T 1 n T 2 n where T 1 n is left unchanged and is O P (1). The factor T 2

Let η > 0 and consider constants c > 1 close to 1 and α > 0 close to 1 2 + γ γ C .

The fact that T 2 n = o P (1) , then, comes from the combination of ( 28) and Lemmas 1 and 2, with delails given below.

First, Lemma 1 is applied with θ = β, which implies that α = (γ + )β and thus the first term of the right hand-side of (28) tends to 0.

Next, Lemma 2 is applied with a = α and needs to be combined with some rate for the factor P n : as in subsection 5.1.2, we need here to distinguish the case γ

In this case, α < 1 and therefore Lemma 2 implies that 1 kn kn i=1 u -α i ξ in = O P (1). Moreover, we have already proved in subsection 5.1.2 that P n = o P (1), and consequently the second term of the right hand-side of (28) tends to 0. (ii) case γ X ≥ γ C

In this case, α > 1 and therefore Lemma 2 implies that, for any given δ > 0, 1

. Therefore, it remains to prove that k α+δ -1 n P n = O P (1) : this has already been done in subsection 5.1.2 (where, there, 1/p denoted β(γ + ) + δ where δ was arbitrary close to 0, as is δ here).

Appendix

Definition 1 An ultimately positive function f : R + → R is regularly varying (at infinity) of order α ∈ R, if

This is noted f ∈ RV α . If α = 0, f is said to be slowly varying.

Proposition 1 (See [Haan and Ferreira (2006)] Proposition B.1.9) Suppose f ∈ RV α . If x > 0 and δ 1 , δ 2 > 0 are given, then there exists t 0 = t 0 (δ 1 , δ 2 ) such that for t ≥ t 0 and tx ≥ t 0 ,

).

If x ≥ 1 and > 0, then there exists t 0 = t 0 ( ) such that for every t ≥ t 0 ,


Noureddine Kettaf
email: nour-eddine.kettaf@uha.fr

Hafid Abouaissa
email: a.abouaissa@uha.fr

Thang Vu

Pascal Lorenz
email: p.lorenz@uha.fr

An Efficient QoS Routing Algorithm for Solving MCP in Ad hoc Networks

Providing guaranteed quality of service (QoS) in wireless networks is a key issue for deploying multimedia applications. To support such a QoS, an arduous problem concerning how to find a feasible end to end path to satisfy multiple QoS constraints should be studied. In general, multi-constrained path selection, with or without optimization, is an NP-complete problem that cannot be exactly solved in polynomial time. Approximation algorithms and heuristics with polynomial and pseudopolynomial time complexities are often used to deal with this problem. However, existing solutions suffer either from excessive computational complexities that cannot be used for multimedia applications in ad hoc networks characterized by mobility and performance constraints (e.g., limited energy, wireless medium, etc.). Recently a promising heuristic algorithm H_MCOP using a non linear Lagrange relaxation path

functions has demonstrated an improvement in its success rate and in finding feasible paths. However, the H_MCOP is not suitable for ad hoc networks and has not exploited the full capability that a Lagrange relaxation could offer. In this paper, we propose an efficient multi-constrained path heuristic called E_MCP, which exploits efficiently the Lagrange relaxation and enhances the path search process to be adequate to mobile ad hoc networks. Using extensive simulations on random mobile network with correlated and uncorrelated link weights, we show that the same level of computational complexity, E_MCP can achieve a higher success ratio of finding feasible paths.

Keywords: Multiple constraints, path selection, QoS routing, ad hoc.

1-Introduction

Routing in ad hoc networks has been an active research area for many years. Much of the original work was motivated by mobile applications, such as on-line live movies and video conferencing. However, most routing protocols proposed to date purely function on best effort basis with no attempt to provide any quality of service (QoS).

In essence, the purpose of QoS routing is to find a feasible path, which has sufficient unused resources to guarantee a certain service level agreement (SLA) between the service provider and the user applications. For example, the delay sensitive applications such as real time voice and video require the data flow to be received at the destination within a certain bounded time. In order to provide such QoS it is imperative that mobile ad hoc networks ensure QoS support in terms of bandwidth, delay, jitter, etc.

In general, routing consists of two basic tasks: distributing the state information of the network and searching this information for a feasible path [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF]. In this paper, we focus on the second task and assume that the network information has been disseminated throughout the network using a QoS-based routing protocol (e.g., ACOR [START_REF] Kettaf | Admission Control enables Routing Protocol -ACOR[END_REF]). Specifically, each link in the network is associated with multiple parameters which can roughly classified into additive [START_REF] Alles | ATM internetworking[END_REF], [START_REF] Wang | On the complexity of quality of service routing[END_REF]. For additive parameters (e.g., delay, jitter), the cost of an end to end path is given by the sum of the individual link values along that path [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF]. However, the cost of the path with respect to (w.r.t) a non-additive parameter, such as bandwidth, is determined by the value of that constraint at the bottleneck link. It is known that constraints associated with non-additive parameters can be easily dealt with a pre-processing step by pruning all links that do not satisfy these constraints [START_REF] Wang | Bandwidth-delay based routing algorithms[END_REF].

Recently, a diversity of QoS routing algorithms incorporating a variety of constraints have been proposed [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF] [6] [START_REF] Szviatovszki | Lagrange relaxation based method for the qos routing problem[END_REF]. For unicast routing, the Multi-Constrained Optimal Path (MCOP) and Multi-Constrained Path (MCP) are the most notorious ones for their NP-complete property [START_REF] Reeves | A distributed algorithm for delay-constrained unicast routing[END_REF]. A MCOP problem needs to find a path with the minimal cost subject to one or more path constraints, while a MCP problem needs to find a path subject to two or more path constraints without necessarily finding an "optimal" solution. Regarding the two classes of problems in wired networks, much work has been done, in particular when the number of constraints is small [START_REF] Gang | An efficient approximate algorithm for delay-costconstrained qos routing[END_REF], [START_REF] Cheng | A new qos routing framework for solving MCP[END_REF].

In ad hoc networks, where the topology changes dynamically, QoS routing is even more challenging. Hence, in this paper we propose an enhanced algorithm based on H_MCOP heuristic (Heuristic Multi-Constrained Optimal Path) [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF], namely E_MCP (Efficient Multi-Constrained Path) to solve the MCP problem in mobile ad hoc networks. To the best of our knowledge, none of the previous published papers deal with the MCP or the MCOP problem in ad hoc networks. Our heuristic aims to enhance the H_MCOP and explore the entire capability of the Lagrange relaxation technique to achieve a high success ratio of finding feasible paths in mobile ad hoc networks. The rest of the paper is organized as fallows. In section 2, we describe the H_MCOP heuristic. Then the proposed heuristic is presented in section 3. The formal specification and justification in section 4, and extensive simulation results are evaluated in section 5. In section 6, we briefly review some related works in wired networks, and finally our main conclusions are drawn in section 7.

2-The H_MCOP heuristic algorithm

In this section, we describe the heuristic algorithm H_MCOP proposed in [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF] Assume a network is modeled by a directed graph G = (V,E), where V is the set of nodes and E is the set of links. Each link is associated with

E j i ∈ ) , (
is associated with a primary cost parameter ) , ( j i c and K additive QoS parameters

) , ( j i w k , ; , 2 , 1 K k K = all parameters are non-negative. Given K constraints k c , ; , 2 , 1 K k K =
the problem is to find a path p from a source node s to a destination node d such that:

(i) k p j i k k c j i w p w ≤ = ∑ ∈ ) , ( ) ( ) , ( for ; , 2 , 1 K k K = and
(ii)

∑ ∈ = p j i j i c p c ) , ( ) , ( ) (
is minimized over all feasible paths satisfying (i).

The search of a feasible path with H_MCOP is done by approximating the non-linear cost function [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF].

λ λ λ λ ) ) ( ( ... ) ) ( ( ) ) ( ( ) ( 2 2 1 1 m m C p w C p w C p w p f + + + = (1)
Where 1 ≥ λ .

As mentioned above, the H_MCOP runs two slightly modified versions of Dijkstra's algorithm. First, it computes in the backward direction the shortest paths from every node to d w.r.t. a linear combination of all weights. Second, in the forward direction, H_MCOP starts from s and discovers every node u based on a path p, where p is a heuristically determined complete s-d path that is obtained by concatenating the already traveled sub-path, from s to u and the estimated remaining sub-path from u to d. Since the algorithm considers complete paths before reaching d, it can foresee some feasible paths during the search. If paths seem feasible, the algorithm can switch to explore these feasible paths for the one with minimum length. A pseudocode of the H_MCOP algorithm is shown in A1.

Algorithm A1:

The heuristic algorithm H_MCOP for the MCOP problem:

H_MCOP (G = (V,E), s, d, C m , m = 1, 2, …,K) Reverse_Dijkstra (G = (V,E),d) if r[u] > K then return failure /*no feasible paths*/ endif Look_Ahead_Dijkstra (G = (V,E),s) if G m [d] ≤ C m , m = 1, 2, …,K then
return the path /*a feasible path is found*/ endif return failure Among the existing work in the context of QoS routing algorithms, it is believed that the H_MCOP has the best performance in terms of finding feasible paths and low computational complexities for wired networks. However, in wireless mobile networks such as ad hoc networks characterized by node mobility and resource constraints have limited the application of the H_MCOP. Hence, in this paper, we investigate a new heuristic which adapts the H_MCOP algorithm and enhances its Lagrange relaxation procedure to provide a solution for the MCP problem in mobile ad hoc networks. By general, the proposed solution is applicable to any number of constraints, irrespective of their nature and interdependence [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF].

3-The efficient multi-constrained path heuristic (E_MCP)

We now present our heuristic algorithm E_MCP, which attempts to the MCP problem in ad hoc networks. E_MCP adopts the basic idea of the H_MCOP but by using an enhanced Lagrange relaxation procedure to be adequate to the mobile environment of ad hoc networks. Specifically, H_MCOP determines the sub-path from s to u by Look_Ahead_Dijkstra w.r.t. a complete path s-d. At this stage, a predecessor of a node v is determined as follows. If there is a neighbouring node u which satisfies the path feasibility from s to d, the intermediate node u is the predecessor, otherwise, the predecessor will be a node u that minimizes (.)

λ f
. The sub-path from u to d is determined using the linear cost function (.) 1 f .

3.1-Nonlinear cost function for MCP

We consider the non linear cost function (1) [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF] for any path p from the source to the destination. For a given 1 ≥ λ suppose there is an algorithm χ that returns a path p by minimizing the cost function [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF]. Then, the following bounds on the performance of algorithm χ can be established.

Theorem 1: Consider the MCP problem and assume that there is at least one feasible path p* in the network. Let p be a path that minimizes the cost function λ g for a given

1 ≥ λ . Then, (i) k k c p w ≤ ) (
for at least one k, and

λ k k Kc p w ≤ ) ( (ii)
for all other k's.

Proof: if the returned path p is feasible, then from (1) the above bounds are correct. Assume that p is not feasible. Since the algorithm returns the path p (and not the feasible path p*), it must be true that

*) ( ) ( p g p g λ λ ≤ In addition, since k k c p w ≤ ) (
for all k's, we have

K p w k ≤ *) ( Thus, K p w k ≤ ) ( (2)
If

k k c p w > ) ( for all k's, then K p g > ) ( λ . Since this contradicts (2), we must have k k c p w ≤ ) (
for at least one k, and the bound in part (i) is correct. Note that if

K p g > ) ( λ , then it is guaranteed that there is no feasible path p* in G because of at least one k, k k c q w > ) (
for any path q.

To prove the part (ii), assume to the contrary that for at least one constraint j c we have

λ j j Kc p w ≤ ) ( , so that K c p w j j > λ ) ) ( ( . It readily follows that K c p w p g j j > > λ λ ) ) ( ( ) (
, which contradicts (2). Hence, part (ii) is proved.

Corollary 1: as λ increases, the likelihood of finding a feasible path also increases.

Proof: follows immediately from theorem 1. Therefore, to increase the possibility of finding a feasible path, it makes sense to set λ to its largest value, i.e.,

.

∞ → λ

In order to provide a practical computational model for ∞ → λ , we can replace the cost function

) ( lim ) ( * p g p g λ λ ∞ → =
by another cost function that does not explicitly involve λ but that achieves the same ordering of candidate as * g .

3.2-Proposed heuristic for MCP in ad hoc networks

E_MCP adopts the basic idea of H_MCOP, but by using an enhanced Lagrange relaxation procedure to adapt it to the mobile environment of ad hoc networks.

Lagrange relaxation has the property of providing bounds on the value of the optimal objective function and, frequently, of quickly generate good, though not necessarily optimal, solutions with associated performance guarantees [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF]. The major modifications are related to the use of the same cost function in both directions and the evaluation of sub-path instead of the complete path.

K constraints k c , . , , 2 , 1 K k K =
For each node u, the algorithm maintains the following labels: λ g For that it uses Reverse_Dijkstra [START_REF] Patek | Simple alternate routing for differentiated services networks[END_REF] with some modifications to the relaxation procedure as shown in algorithm A3. If there is a possibility that the network contains feasible paths, a heuristic search procedure called Look_Ahead_Dijkstra is executed in the forward direction. Then, the algorithm checks

], [u r ], [u R k ], [u r π ], [u g ], [u G k ] [u g π and ] [u c , . , , 2 , 1 K k K = Label ] [u r represents the cost of the shortest path from u to d w.r.t. the cost function λ g . Labels ] [u R k , K k , , 2 , 1 K
E_MCP (G = (V,E), s, d, C m , m = 1, 2, …,K) Reverse_Dijkstra (G = (V,E),d, m w , m = 1, 2, …,K) if , ] [ m m C s R ≤ K m , , 2
[d R k , , , , 2 , 1 K k
K s r > ] [
to determine the possibility of finding feasible paths (please see Theorem 1).

If there is a possibility that the network contains feasible paths, Look_Ahead_Dijkstra is executed in the forward direction to identify whether there is other paths q, which improve the performance over the returned path p using the information provided by the above Reverse_Dijkstra. To implement Look_Ahead_Dijkstra, we need a slight modification to the relaxation procedure of Dijkstra's algorithm [START_REF] Lee | Routing subject to quality of service constraints in integrated communication networks[END_REF]. The preference rule takes as input two nodes and their labels. Then, it selects a path of these nodes such that minimizes the primary cost function if foreseen s-d paths passing through these nodes are feasible; otherwise, it selects the path that minimizes the objective function .

Reverse_Dijkstra_Relax (u,v) if ∑ = + > K m m m m C v u w v R u r 1 ) ) , ( ] [ ( ] [ λ then ∑ = + = K m m m m C v u w v R u r 1 ) ) , ( ] [ ( : ] [ λ ) , ( ] [ : ] [ v u w v R u R m m m + = for K m , , 2 , 1 K = u v r = : ] [ π endif
λ g Look_Ahead_Dijkstra_Relax (u,v)
Let tmp be a temporary node

λ ∑ = + = K m m j m C v u w u G tmp f 1 ) ) , ( ] [ ( : ] [ ) , ( ] [ : ] [ v u w u G tmp G m m m + = for K m , , 2 , 1 K = ] [ : ] [ v R tmp R j j = for K m , , 2 , 1 K = if ( ) ) , ( _ _ Pr tmp v tmp best the efer = then ] [ : ] [ tmp f v f = ] [ : ] [ tmp G v G m m = u v f = : ] [ π endif

4-Justification of E_MCP heuristic algorithm

We aim with E_MCP to exploit the full capability that a Lagrange relaxation technique could offer to solve the MCP problem in a mobile environment such as ad hoc networks.

Prefer_the_best (a,b) if ), ,..., 2 , 1 ( K m = ∀ m m m C a R a G ≤ + ] [ ] [ then return ( a ) if ), ,..., 2 , 1 ( K m = ∀ m m m C b R b G ≤ + ] [ ] [ then return ( b ) if ] [ ] [ b f a f < then return ( b )
Respectively, it is necessary to justify the proposed Lagrange process and the method that the E_MCP determines a path according to the nonlinear cost function

(.). λ f
Consider a special case where the following condition is satisfied:

∑ = Δ = k m m p x 1 ) ( (3)
Where Δ is a fixed value,

m m m C p w p x ) ( ) ( =
and p is any path between two nodes, say from u to v.

Under this condition, we need to determine the path from u to v that minimizes (.).

λ f
The problem can be described as the following constrained optimization problem:

Minimize ∑ = = k m m p x p f 1 )) ( ( ) ( λ λ
Subject to the condition in (3).

It can be further converted to the following unconstrained optimization problem:

Minimize ) ) ( ( )) ( ( ) ( 1 1 ∑ ∑ = = Δ - + = K m m K m m p x p x p g σ λ λ
Where σ is a Lagrange multiplier.

To solve the unconstrained optimization problem, we let

, 0 ) ( ) ( = ∂ ∂ p x p h m , ,..., 2 , 1 K m = , 0 ) ( = ∂ ∂ σ p h Which lead to , ) ( K p x m Δ = . ,..., 2 , 1 K m =
Remark: the derivation does not apply the case of . 1 = λ We demonstrate by the derivation above that under the condition in Eq.3, the optimal path would be the path that has the equal ratio of each weight to the corresponding upper bound. In real situation, Eq.3 does not hold normally, but we can infer that the optimal solution would be a path that has the following two characteristics.

C1: it has a relatively small summation of all ratios, i. Similarly, for the sub-path v-d expect the weight is accumulated from v to d. Hence, using the sub-paths search in this fashion should achieve the highest probability of finding a feasible MCP solution.

We believe that this method could achieve good performance when all the constraints are quite loose, but when 0 C becomes relatively tight, H_MCOP is likely to make mistakes. To enlighten this observation, consider the case where all constraints are quite loose except that 0 C is just slightly greater than ) ( 0 0 p w , where 0 p is the shortest path between s and d w.r.t. weight 0 w , which implies the availability of a very limited number of feasible paths.

5-Performance evaluation

In this section, we simulate our algorithm in ad hoc network environment. The results are for a network of 50 nodes randomly on a square 670m flat space. To preserve connectivity, every node has 12 neighbours on average. Nodes move following the Random Trip Model [START_REF] Jean-Yves | Perfect simulation and stationarity of a class of mobility models[END_REF]. We associate two randomly generated weights with each link ) , ( j i . As shown in Table 1, these weights are selected from uniform distributions under several types of correlation between them.

If there is positive correlation, we assume that both weights are selected from uniform distributions with either small mean or large mean. If there is negative correlation, we assume that one of weights is selected from a uniform distribution with small mean while the other is selected from another uniform distribution with large mean. If there is no specific correlation, we assume both weights are independently selected from uniform distributions. The primary cost of a link ) , ( j i is taken as

) , ( j i c ~]. 50 , 1 [ uniform
The source and destination of a request are randomly generated such that the minimum hop-count between them is at least three.

Positive correlation

No correlation Negative correlation We contrast the performance of E_MCP and MH_MCOP, which is a modified H_MCOP to be MCP based, using the success ratio SR, which refers to the fraction of connection requests for which feasible paths are found by the given approximation algorithm or heuristic in several scenarios.

The results reported here are averaged over several runs. In each run, random graphs are generated according to the mobility model. For each instance of a random graph, ten independent realizations of link weights are generated using different random seeds. Finally for each instance of a random graph with given link weights, about 5~20 connection requests are generated for graphs with 10, 25, and 50 nodes, respectively. When link weights are positively correlated, the path weights also become positively correlated [START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF], and thus a nonlinear approximation heuristic has a better success ratio of finding feasible path. However, if the link weights are negatively correlated, there will be more paths for which ) ( ) ( We now contrast our study on the performance comparison using the k-shortest paths for E_MCP and MH_MCOP. Comparing the equal complexities of E_MCP and HH_MCOP, we claim that E_MCP outperforms MH_MCOP in performance for the same amount of computational complexity.

Figure2. SR's of E_MCP and MH_MCOP with k-shortest paths of 50 nodes.

6-Related works

In wired networks MCP has been extensively studied. Jaffe [START_REF] Jaffe | Algorithms for finding paths with multiple constraints[END_REF] presented two heuristics. One takes pseudo-polynomial-time. The other takes polynomial-time. The latter is to minimize (w 1 (p) + d*w 2 (p)), where w 1 (p) is path length, w 2 (p) is path cost and d is a weight factor. Chen and Nahrstedt's approach [START_REF] Chen | On finding multi-constrained paths[END_REF] is to map the cost (or delay) of every link from an unbounded number of integers, which reduces the original NP-hard problem to a simpler problem that can be solved in polynomial time.

When Dijkstra's algorithm [START_REF] Bertsekas | Data Networks[END_REF] is used, the running time of their algorithm is O(x 2 |V| 2 ); when Bellman-Ford algorithm is used, the running time is O(x|V||E|), where x is an adjustable positive integer and its value can be as high as 10|V|, resulting in an overall running time of O(|V| 4 ). In [START_REF] Korkmaz | A randomized algorithm for finding a path subject to multiple QoS constraints[END_REF], Korkmaz et al.'s idea is firstly to prune all the links that cannot be on any feasible path connecting the source-destination pair. It then uses a randomized search to find a feasible path with a running time of O(|V| 3 ).

In [START_REF] Korkmaz | An efficient algorithm for finding a path subject to two additive constraints[END_REF], considering a combined weight αw 1 (p) +  βw 2 (p), Korkmaz et al. provided an approximation algorithm that searches for appropriate values of  α and β to determine a path subject to two additive constraints with a running time of a logarithmic number of Dijkstra's algorithm.

In [START_REF] Blokh | An approximation algorithm for combinatorial optimization problems with two parameters[END_REF] the authors proposed a similar dynamic algorithm for the MCP problem.

However, the computational complexity of this algorithm grows exponentially with the size of the network. In [START_REF] Chen | On finding multi-constrained paths[END_REF] the authors proposed a heuristic algorithm that modifies the MCP problem by scaling down the values of one of the two link weights to bounded integers. It was shown that the modified problem can be solved by using Dijkstra's shortest path algorithm and that the solution to the modified algorithm is also a solution to the original one.

7-Conclusions

Solving the MCP problem in ad hoc networks is a key issue for QoS application. In Extensive simulations demonstrate that E_MCP can achieve a higher probability of finding feasible paths than the original H_MCOP (or MH_MCOP) with better performance.

As future work, we will investigate how E_MCP performs in the presence of inaccurate state of information and what modifications need to be done.

Definition 1 :

 1 for the MCOP (Multi-Constrained Optimal Path) problem which only needs to run Dijkstra's algorithm (with slight modifications) twice to solve the MCOP problem with multiple constraints based on linear Lagrange relaxation. Multi-Constrained Optimal Path (MCOP) problem:

=

  represents the individually accumulated link weights along the path. The predecessor of this path is stored in label ] function of a foreseen sub-paths from s to u to d w.r.t. the nonlinear cost function λ g . accumulated cost of link weights and the primary cost along the already travelled segment of this path from s to u. The predecessor of u on this path is stored in the label ] The heuristic algorithm E_MCP for the MCP problem: As in H_MCOP, in E_MCP there are two directions in the algorithm: backward (from d to u) to estimate the cost of the remaining segment using the nonlinear cost function λ g and forward (from s to u) to find the most promising feasible path. In the backward direction, the E_MCP algorithm finds the path from every node u to d w.r.t. the cost function (.).

  , 2, …,K then return the path /*a feasible path is found*/ endif return failure Algorithm A3: The relaxation procedure of subroutine Reverse_Dijkstra in E_MCP Initially, Reverse_Dijkstra sets ∞ node u. It then starts at node d by setting ] [d r and ]

  K = to zeros. It explores the graph and eventually returns a path p from s to d. The main advantage of this modification is the enhancement of the probability of finding feasible paths by increasing λ .

  the ratios are very close to each other. The E_MCP executes Look_Ahead_Dijkstra to find the sub-path from s to an intermediate node v by evaluating (.) λ f w.r.t. the path itself. Moreover, the sub-path s-v must meet C1 and C2. On the other hand, E_MCP executes the Reverse_Dijkstra to choose the sub-path from an intermediate node v to d by evaluating (.) λ f which also must meet C1 and C2. In terms of weight, the sub-path from s to v has nearly the smallest ratio of weights accumulated at any node v to the upper bound, which is very close to that for any other weight.

Figure 1 .Figure 1

 11 Figure 1. SR of E_MCP and MH_MCOP used on random graphs with 10, 25 and 50

Figure 2

 2 bellow shows the SR's of both algorithms versus the number of considered shortest paths. E_MCP gives better SR using smaller values of k than MH_MCOP regardless of the number of nodes or correlation between the link weights. In other words, the SR of MH_MCOP increases with k and finally converges the SR that initially provided by E_MCP with smaller k. In essence, the SR provided by E_MCP with k = 1 can only achieved by MH_MCOP when k = 2.

1 f

 1 this article, we introduced a new heuristic algorithm E_MCP based on Lagrange relaxation technique for solving MCP problem in mobile ad hoc networks by introducing some efficient modifications to the H_MCOP algorithm, which proves a high success ratio in finding feasible paths. Hence, our modifications aimed to explore the full capability of the Lagrange relaxation and enhance the probability of finding feasible paths. We first, modified the relaxation procedure in Look_Ahead_Dijkstra's algorithm. Second, instead of using a linear cost function (.) , we employed the same nonlinear cost function (.) λ f in Reverse_Dijkstra.

Algorithm A5 :

 A5 The preference rule for Look_Ahead_Dijkstra in E_MCP The computational and space complexities of the resulting E_MCP algorithms are equal to that of Dijkstra's, since at most two modified versions of Dijkstra's algorithm

	are executed with the complexity of		O	(	n	log(	n	)	+	m	).	To improve the performance, the
	path algorithm is	O	(	km	log(	kn	)	+	k	2 m	).		Hence, the complexity of the E_MCP with k-
	shortest path is	O	(	n	log(	n	)	+	km	log(	kn	)	+	(	k	2	+	) 1	m	).

forward direction of the E_MCP can be used with the k-shortest path implementation of Dijkstra's algorithm presented in

[START_REF] Chong | On finding single-source single-destination k-shortest paths[END_REF]

. Note that k-shortest paths are considered w.r.t. the minimization of the nonlinear cost function. The complexity of the k-shortest

Table 1 .

 1 Ranges of link weights and the correlation between themThe constraints[START_REF] Korkmaz | Multi-constrained optimal path selection[END_REF] are also randomly generated, but their ranges determined based on the best paths w.r.t. 1

	1 w	( i	,	j	)	~] 10 , 1 [ uniform	1 w	( i	,	j	)	~	, 1 [ uniform	] 25	1 w	( i	,	j	)	~	] 10 , 1 [ uniform
	2 w	( i	,	j	)	~] 25 , 1 [ uniform	2 w	( i	,	j	)	~	] 50 , 1 [ uniform	2 w	( i	,	j	)	~	[ uniform	] 50 , 25
								OR														OR
	1 w	( i	,	j	)	~] 25 , 10 [ uniform											1 w	( i	,	j	)	~	, 10 [ uniform	] 25
	2 w	( i	,	j	)	~] 50 , 25 [ uniform											2 w	( i	,	j	)	~] 25 , 1 [ uniform