N Zoroa

E Lesigne

M J Fernández-Sáez

P Zoroa

& Casas

The coupon collector urn model with unequal probabilities in ecology and evolution

Keywords: Coupon collector's problem, parasitoid, stochastic dominance, strong dominance, ecology. 2010 Mathematics Subjet Classication: 60C05, 60E15, 92B05 1

The sequential sampling of populations with unequal probabilities and with replacement in a closed population is a recurrent problem in ecology and evolution. Many of these questions can be reformulated as urn problems, often as special cases of the coupon collector problem, most simply expressed as the number of coupons that must be collected to have a complete set. We aimed to apply the coupon collector model in a comprehensive manner to one example -hosts (balls) being searched (draws) and parasitized (ball color change) by parasitic wasps -to evaluate the inuence of dierences in sampling probabilities between items on collection speed.

Based on the model of a complete multinomial process over time, we dene the distribution, distribution function, expectation and variance of the number of hosts parasitized after a given time, as well as the inverse problem, estimating the sampling eort. We develop the relationship between the risk distribution on the set of hosts and the speed of parasitization and propose a more elegant proof of the weak stochastic dominance among speed of parasitization, using the concept of Schur convexity and the Robin Hood transfer numerical operation.

Numerical examples are provided and a conjecture about strong dominance is proposed. The speed at which new items are discovered is a function of the entire shape of the sampling probability distribution. The sole comparison of values of variances is not sucient to compare speeds associated to two dierent distributions, as generally assumed in ecological studies.

Introduction

The description of sequential sampling of a population of individuals for which the probability of being selected does not vary until a specic event, such as the collection of all or some types of individuals or a specic subgroup of the population, occurs is a common problem in ecology and evolution studies. In probability theory, such problems are often treated as urn problems, generally as the coupon collector problem (CCP). The CCP is a mathematical model that belongs to the family of urn problems that can be formulated as follows:

A company issues coupons of dierent types, each with a particular probability of being issued. The object of interest is the number of coupons that must be collected to obtain a full collection. This problem has been widely studied. The rst ndings concerned the classical problem in which all coupons are equally likely to be obtained [START_REF] Feller | An Introduction to Probability Theory and Its Applications: Volume One[END_REF]). Rapid advances have been made in this eld [START_REF] Boneh | The coupon-colector problem revisited[END_REF][START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF][START_REF] Anceaume | New results on a generalized coupon collector problem using Markov chains[END_REF], but they have gone largely unnoticed by most scientists working in ecology and evolution. This is partly due to diculties in making the correct analogies, partly due to a lack of worked examples and partly because each eld devises its own vocabulary, procedures and formalism. In ecological sciences, for example, a vibrant eld of theoretical and applied ecological statistics developed in the 1950s from the repeated sampling of populations to estimate biodiversity richness (McArthur, 1957;[START_REF] Simpson | Measurement of diversity[END_REF]. This eld could greatly benet from the latest advances in the CCP [START_REF] Bunge | Estimating the number of species: a review[END_REF][START_REF] Huillet | Sampling from Dirichlet partitions: estimating the number of species[END_REF].

Related problems deal with the relative abundance of species from a community containing many species [START_REF] Dennehy | Bacteriophages as model organisms for virus emergence research[END_REF], or the sampling eort required to achieve a particular level of coverage [START_REF] Neal | Sampling Eciency and Biodiversity[END_REF]. Increases in the number of new hosts being infected or superinfected are a topic of great importance in population dynamics and epidemiology [START_REF] Daley | Epidemic modelling: an introduction[END_REF][START_REF] Lloyd-Smith | Superspreading and the eect of individual variation on disease emergence[END_REF][START_REF] Keeling | Modelling infectious diseases in humans and animals[END_REF]. Several of the questions posed in capture-recapture studies relate to the coupon collector problem. Occupancy problems and related capture-recapture techniques are, indeed, dened as problems in which the probability of a given species occupying a given state at a given time must be determined (see the review [START_REF] Arnold | Special feature. Modelling demographic processes in marked populations: Proceedings of the Euring 2013 Analytical Meeting[END_REF]the paper Hernandez-Suarez andHiebeler, 2011). In ethological sciences, the estimation of a repertoire of signals in animal communication is considered as a form of the CCP, because vocal repertoire size is a key behavioral indicator of the complexity of the vocal communication system in birds and mammals [START_REF] Kershenbaum | Estimating vocal repertoire size is like collecting coupons: A theoretical framework with heterogeneity in signal abundance[END_REF]. In genetics and evolution, the coupon collector problem has been recognized as such only occasionally, despite these elds having generated some of the most elegant theorems and uses of other urn processes [START_REF] Ewens | The sampling theory of selectively neutral alleles[END_REF]Donelly, 1986;[START_REF] Doumas | How many trials does it take to collect all dierent types of a population with probability p?[END_REF]. Indeed, the coupon collector problem has been used in the context of exhaustive haplotype sampling in phylogeography, [START_REF] Dixon | A means of estimating the completeness of haplotype sampling using the Stirling probability distribution[END_REF], determining the number of benecial mutations as a function of sequence lines [START_REF] Tenaillon | The molecular diversity of adaptive convergence[END_REF], and estimations of the size of the library required to target a particular percentage of the non-essential genome displaying a given property [START_REF] Vandewalle | Characterization of genome-wide ordered sequence-tagged Mycobacterium mutant libraries by Cartesian Pooling-Coordinate Sequencing[END_REF], for example.

Urn models have been much more widely used for modelling host-parasitoid systems than in other topics of ecology. We therefore used the biological context and formalism of parasitism by parasitic wasps, as the results obtained with this system can easily be extended to other ecological and evolutionary contexts. Parasitic wasps search for insects hosts, such as caterpillars, in which they lay a single, or multiple eggs. In solitary wasp species, only a single wasp develops fully in a given host. Parasitism can thus be formalized as a probabilistic dynamic process with hosts as `balls' and parasitoids changing their `color' by parasitizing them. In work beginning more than a century ago [START_REF] Fiske | Superparasitism: an important factor in the natural control of insects[END_REF][START_REF] Thompson | La théorie mathématique de l'action des parasites entomophages et le facteur du hasard[END_REF], the pioneering population dynamicists assumed that hosts were found and attacked on successive occasions governed by exponential laws in continuous time. The number of draws was thus considered to be random and the number of eggs for a given host was assumed to follow a Poisson law [START_REF] Montovan | The puzzle of partial resource use by a parasitoid wasp[END_REF]. If we assume that the number of draws is xed, then the distribution of the number of eggs for a given host is binomial, but closely approximates a Poisson distribution in large host populations. The proportion of the population without eggs (the zero class) is of particular interest, because these hosts survive parasitism and produce ospring for the next generation.

In eld studies however, observed distributions are generally more aggregated than would be expected under the assumed Poisson distribution (Hemerick et al., 2002). Aggregation is interpreted as the result of heterogeneity in the risk of being found, due to dierences in location, accessibility, appearance, color, developmental stage or any other trait [START_REF] Hassell | The spatial and temporal dynamics of host-parasitoid interactions[END_REF][START_REF] Murdoch | Consumer-Resource Dynamics (MPB-36[END_REF]. The risk distribution greatly inuences the stability of the host-parasitoid system and has been widely studied [START_REF] May | Host-parasitoid systems in patchy environments: a phenomenological model[END_REF][START_REF] Ives | Variability and parasitoid foraging eciency: a case study of pea aphids and Aphidius ervi[END_REF][START_REF] Singh | Skewed attacks, stability, and host suppression[END_REF].

All these models focus on the distribution of eggs over the entire population of hosts, after a given time or a given number of draws (Figure1). However, the use of this distribution greatly decreases the amount of information available, as it collapses individual host histories. Parasitism is a multinomial process (Fig- ure1), in which time corresponds to host draws. Its dynamics determines, for example, the percentage of hosts parasitized at the end of the season, the opportunity and time at which alternative pest control methods need to be deployed in supplement in biological control with parasitoids, and the time required to achieve a given degree of control by parasitic wasps. In the present paper, we aimed to model parasitism as a multinomial urn process over time and we study the speed of parasitization (Figure1). We consider host encounters followed by oviposition without discrimination. The parasitism process described above can be considered as a coupon collector problem. In this case, there is a nite population of hosts diering in appearance, location, developmental stage or other factors. This heterogeneity results in dierent probabilities of hosts being found by parasitoids. These probabilities, p h for host h, do not change over time.

Our work therefore entails describing in depth the coupon collector problem, highlighting unnoticed analogies among previous works within the probability literature, and comparing the inuence of the degree of heterogeneity among hosts on the speed of infection. We give a compact and hopefully more elegant proof than previously known of the following fact : the more the distribution p on the set of hosts is heterogeneous, the more the (random) number Y of parasitized hosts after a given number of draws is small; in other terms, there is a monotonic relationship between the majorization relation on the set of probability distributions p with the stochastic dominance on the set of random numbers Y .

This paper is structured as follows. In Section 2, we dene a succession S n , n = 1, 2, ..., of N -dimensional random variables describing the state of the host population over time, in which time, n, is given by the number of attacks on the set of hosts. Each marginal distribution of S n provides us information about a subset of hosts, including, in particular, the h-th component representing the number of times that host h has been attacked by a parasitoid between times 1 and n. In Section 3 we dene the random variables Y n , n = 1, 2, ..., representing the number of parasitized hosts after n draws. We also compute the distribution, the distribution function, the expectation and the variance of Y n . We found no examples of calculations of this value in previous studies and therefore believe this aspect to be novel. We obtain the expected number of draws required for all hosts in a given subset to be parasitized and provide upper and lower bounds for this value in Section 4. In Section 5, we apply the results developed in previous sections to two particular risk distributions on the set of hosts. We rst use the uniform distribution, and then a distribution corresponding to a host population with two dierent kinds of hosts. We calculate the most relevant values for each of these cases. In Section 6, we develop the relationship between the speed of parasitization and the risk distribution in the set of hosts. A narrower risk distribution is associated with faster parasitization. Thus, parasitization is fastest when the risk distribution is uniform. We highlight this nding with numerical examples in Section 7 and propose a conjecture on strong stochastic dominance in Section 8.

Modelling parasitism as an urn process

We assume a nite population of hosts, constant for the entire duration of the experiment. The parasitoid population is irrelevant, but we assume that the number of eggs that can be laid in the host population is not limiting.

The situation is developed in successive stages or draws. At each stage, a parasitoid attacks a host, in which it lays an egg. The model is based on the fundamental assumption that successive draws are mutually independent. The hosts dier in appearance due to intrinsic qualities, and these dierences modify their probability of being attacked by a parasitoid. If the hosts are named 1, 2, 3, . . . , N , then host h has a probability p h ≥ 0 ( p h = 1) of being attacked by a parasitoid in a draw. This probability does not change during the process.

We will say that p 1 ,p 2 ,...,p N or (p 1 , p 2 , ..., p N ) is the risk distribution for the set of hosts H = {1, 2, ..., N }.

The underlying probability space of our model is (Ω, S, P ), where the elements of Ω are all the possible histories of parasitism, that is Ω = H N , equipped with its product σ-algebra S and the probability P given by Kolmogorov´s Theorem: if we therefore x i 1 , i 2 , i 3 , . . . , i n in H, the probability of the event {ω = i 1 i 2 i 3 . . . i n j n+1 j n+2 • • • : for some j k in H, k > n} is p i1 p i2 . . . p in . When necessary, the vector (p 1 , p 2 , ..., p N ) will be denoted by a single letter p, and the probability P will be denoted by P p .

We can describe this situation by dening a succession of random variables, S n = (S n1 , S n2 , ..., S nN ), n = 1, 2, ...

(2.1)
where S nj denotes the number of eggs in host j after n draws.

Variable S n represents the state of the host population after n draws, that is, the distribution of eggs over the total population of hosts. If host i was visited r i times between stages 1 and n, for i = 1, 2, ..., N , then S n takes the value (r 1 , r 2 , ..., r N ). This variable has a multinomial distribution with parameters n, p 1 ,p 2 ,...,p N , that is, for every integers r 1 , r 2 ,..., r N , 0 ≤ r i ≤ n, r 1 +r 2 +...+r N = n,

P (S n1 = r 1 , S n2 = r 2 , ..., S nN = r N ) = n! r 1 !r 2 !...r N ! p r1 1 p r2 2 ...p r N N .
(2.2)

The marginal distribution of S n , for i 1 , i 2 , ..., i h distinct elements of {1, 2, ..., N } is given by

P (S ni1 = r i1 , S ni2 = r i2 , ..., S ni h = r i h ) = n! r i1 !r i2 !...r i h !(n -r ij )! p ri 1 i1 p ri 2 i2 ...p ri h i h q n-ri j i1i2...i h , (2.3) with 0 ≤ r ij ≤ n, j = 1, 2, ..., h, r i1 + r i2 + ... + r i h ≤ n, q i1i2...i h = 1 - h j=1 p ij .
This is the probability that, after n draws host i 1 has been visited r i1 times by the parasitoids, host i 2 r i2 times and host i h r i h times, without considering the rest of the hosts.

In particular, the component S nh of S n has a binomial distribution with parameters n, p h ,

P (S nh = r) = n! r!(n -r)! p r h q n-r h , r = 0, 1, 2, ..., n (2.4)
where

q h = 1 -p h = i =h p i .
This variable represents the state of host h after n draws. Thus, P (S nh = r) is the probability that host h has been attaked r times during the n draws.

The expected value and variance of this random variable are

E(S nh ) = np h , V ar(S nh ) = np h (1 -p h ).
Let (e 1 , e 2 , ..., e N ) denote the canonical base of the space R N . We emphasize that the process (S n ) n≥1 is the random walk on Z N + with independent increments obeying the following law: S n+1 -S n = e k with probability p k . The statistical behavior of this process is also very well known.

Note that, in this model, the sequence of random subsets of H, describing the set of parasitized hosts over time, is a Markov chain, and it is not dicult to give a precise description of its probability transitions. However, it is not straightforward to study this Markov chain directly.

Number of parasitized hosts after n draws

Let Y n be the random variable representing the number of parasitized hosts after n draws, that is Y n = k if there are exactly k parasitized hosts after n draws. In this section we study this random variable obtaining expressions for: From now on, for any integer h > 0 and real x, we write

x (h) = x(x -1)(x -2)...(x -h + 1), x (0) = 1, x h = x (h) h! = x(x -1)...(x -h + 1) h! , and x 0 = 1. (3.1)
The distribution and the distribution function of Y n have been obtained in previous studies, see [START_REF] Anceaume | New results on a generalized coupon collector problem using Markov chains[END_REF] .

By applying the general inclusion and exclusion principle, we nd that, for any distinct elements j 1 , j 2 , ..., j k of H, and denoting p j1 j2...j k = p j1 + p j2 + ... + p j k , P (the set of parasitized hosts after n draws is {j 1 , j 2 , ...,

j k }) = p n j1 j2...j k - {i1,i2,...,i k-1 }⊂{j1,j2,...,j k } p n i1 i2...i k-1 + {i1,i2,...,i k-2 }⊂{j1,j2,...,j k } p n i1 i2...i k-2 -... + (-1) k-1 i∈{j1,j2,...,j k } p n i ,
from which we deduce that

P (Y n = k) = {j1,j2,...,j k }⊂{1,2,...,N } p n j1 j2...j k - N -k + 1 N -k {j1,j2,...,j k-1 }⊂{1,2,...,N } p n j1 j2...j k-1 + N -k + 2 N -k {j1,j2,...,j k-2 }⊂{1,2,...,N } p n j1 j2...j k-2 -... +(-1) k-1 N -1 N -k {j}⊂{1,2,...,N } p n j .
(3.2)

Using the notation p A = i∈A p i for any A ⊂ H, this can be written in a more compact form

P (Y n = k) = A⊂H,|A|≤k (-1) k-|A| N -|A| k -|A| p n A , for 0 ≤ k ≤ N, and k ≤ n. (3.3)
where |A| denotes the number of elements of the set A.

Let us now consider the distribution function of Y n ,

P (Y n ≤ k) = k j=1 P (Y n = j) = k j=1 A⊂H,|A|≤j (-1) j-|A| N -|A| j -|A| p n A = A⊂H,|A|≤k   k-|A| i=0 (-1) i N -|A| i   p n A .
(3.4)

As, for any integers K and k ≥ 0 the equality

k i=0 (-1) i K i = (-1) k K -1 k holds, we obtain P (Y n ≤ k) = A⊂H,|A|≤k (-1) k-|A| N -|A| -1 k -|A| p n A , k = 1, 2, ..., N. (3.5)
A similar expression can be seen in [START_REF] Anceaume | New results on a generalized coupon collector problem using Markov chains[END_REF] . From (3.5) we can calculate the moments of Y n . Let

m [n] k = A⊂H,|A|=k p n A , for every k ≤ N k j=1 P (Y n ≤ j) = k l=1 (-1) k-l N -l -2 k -l m [n] l , this gives, with k = N -1 and k = N , N -1 j=1 P (Y n ≤ j) = m [n] N -1 , N j=1 P (Y n ≤ j) = m [n] N + m [n] N -1 = 1 + m [n] N -1 .
And, bearing in mind that

E(Y n ) = N j=1 P (Y n ≥ j) = 1 + N j=1 P (Y n > j) = 1 + N - N j=1 P (Y n ≤ j),
we obtain the well known formula:

E(Y n ) = N -m [n] N -1 = N - N i=1
(1 -p i ) n .

(3.6)

We were unable to nd any expression for E(Y 2 n ) and the variance of Y n , in previous studies. These two quantities can be obtained as follows. We compute,

for k ≤ N k t=1 t j=1 P (Y n ≤ j) = k l=1 (-1) k-l N -l -3 k -l m [n] l , then, for k = N -2, N -1 and N we obtain N -2 t=1 t j=1 P (Y n ≤ j) = m [n] N -2 , N -1 t=1 t j=1 P (Y n ≤ j) = m [n] N -2 + m [n] N -1 ,
and

N t=1 t j=1 P (Y n ≤ j) = m [n] N -2 + 2m [n] N -1 + m [n]
N .

The last identity can be written:

N j=1 (N -j + 1)(N -j + 2) 2 P (Y n = j) = m [n] N -2 + 2m [n] N -1 + m [n] N , this gives (N + 1)(N + 2) 2 - 2N + 3 2 E(Y n ) + 1 2 E(Y 2 n ) = m [n] N -2 + 2m [n] N -1 + m [n] N , therefore E(Y 2 n ) = 2m [n] N -2 -(2N -1)m [n] N -1 + N 2 = 2 1≤i<j≤N (1 -p i -p j ) n -(2N -1) N i=1 (1 -p i ) n + N 2 and V ar(Y n ) = 2 1≤i<j≤N (1 -p i -p j ) n + N i=1 (1 -p i ) n 1 - N i=1 (1 -p i ) n . (3.7)
4

The number of draws required to reach a given level of parasitism

The expected number of draws required for the parasitization of k unparasitized hosts may be of considerable interest. For example, we might want to know the expected number of draws required for k of the hosts occupying a determined region, or with probabilities of parasitization greater (or less) than a given value, etc., are parasitized. We dene below a random variable representing the number of draws required for the event of interest to happen and we obtain its expectation. We also describe the relationship between the random variables dened here and the variables Y n dened in Section 3. Let us consider that, at a given stage of the process, there is a set K ⊂ H of unparasitized hosts, this is our set of interest, and the remaining hosts H -K are or are not parasitized. Let us use X to denote the number of hosts in the set H -K attacked by the parasitoids before one of the hosts in K is attacked.

As this process involves the repeating of independent trials, the random variable X follows a geometric distribution with parameter p = i∈K p i , (or a degenerate distribution if K = H). Thus,

E(X) = i∈H-K p i i∈K p i . (4.1) Now, let k and N 1 be integers 1 ≤ k ≤ N 1 ≤ N . Let H 1 be a subset of the set of hosts, H, and H 2 = H -H 1 , |H 1 | = N 1 . We can assume that H 1 = {1, 2, ...N 1 } without lost of generality.
If we consider the hosts of set H 1 to be unparasitized, then we can dene T k,N1 as the random number of draws required to ensure that k hosts of set H 1 have been parasitized. Its expectation is the expected number of draws required for k hosts of set H 1 be parasitized. The case H 1 = H has been studied before and dierent expressions for E(T k,N ) have been obtained. We include these at the end of this section. In [START_REF] Boneh | The coupon-colector problem revisited[END_REF], an expression is proposed for the particular case in which k = N 1 = N . The expression obtained here is more general.

Let i 1 , i 2 , ..., i k be distinct elements of H 1 . Let D i1i2...i k be the event dened by the fact that the rst k hosts of set H 1 parasitized (i.e. attacked by a parasitoid for rst time) are hosts i 1 , i 2 , ..., i k and are parasitized in the precise order i 1 , i 2 , ..., i k . In other words, some of the hosts of set H 2 may be attacked rst, followed by host i 1 . Next, some hosts of H 2 ∪ {i 1 } may be attacked, followed by host i 2 , etc. Let p = i∈H1 p i , q = 1 -p = i∈H2 p i . Then

P (D i1i2...i k ) = P (first host of H 1 parasitized is i 1 ) P (second host of H 1 parasitized is i 2 | first host of H 1 parasitized was i 1 ) ... P (k-th host of H 1 parasitized isi k | first (k -1) hosts of H 1 parasitized were i 1 , i 2 , ..., i k-1 ).
Both in the case q = 0 (H 1 = {1, 2, ..., N }) and the case q > 0

P (first parasitized host of H 1 is i 1 ) = p i1 1 -q , where q = i∈H2 p i > 0.
For the rest of the factors

P (h -th parasitized host of H 1 is i h | first parasitized hosts of H 1 were i 1 , i 2 , ..., i h-1 ) = ∞ r=0 p i h (q + h-1 j=1 p ij ) r = p i h 1 -q - h-1 j=1 p ij , h = 1, 2, ..., k, where q = i∈H2 p i , therefore, P (D i1i2...i k ) = k j=1 p ij p(p -p i1 )(p -p i1 -p i2 )...(p - k-1 j=1 p ij ) . (4.2) Let Π k be the set of all k-permutations of 1, 2, ..., N 1 . Then the events D i1i2...i k with (i 1 i 2 ...i k ) ∈ Π k constitute a partition of Ω, i.e. D i1i2...i k ∩ D j1j2...j k = if i 1 i 2 ...i N1 = j 1 j 2 ...j N1 and (i1i2...i k )∈Π k P (D i1i2...i k ) = 1.
We can then write E(T k,N1 ) as follows,

E(T k,N1 ) = (i1i2...i k )∈Π k E(T k,N1 |D i1i2...i k )P (D i1i2...i k ). (4.3)
To compute the conditional expectations E(T k,N1 |D i1i2...i k ), let us denote by X h the random variable representing the number of draws elapsed after h -1 hosts of the set H 1 being parasitized and before a new host of the set

H 1 is parasitized, 1 ≤ h ≤ k. We can then write T k,N1 = X 1 + 1 + X 2 + 1 + ... + X k + 1 = X 1 + X 2 + ... + X k + k. (4.4)
and therefore

E(T k,N1 |D i1i2...i k ) = k h=1 E(X h |D i1i2...i k ) + k (4.5) but E(X h |D i1i2...i k ) = E (X h |already parasitized hosts are those of H 2 and i 1 i 2 ...i h-1 ) (4.6)
One direct application of (4.1) would then be:

E(X h |D i1i2...i k ) = q + h-1 j=1 p ij p - h-1 j=1 p ij for h = 1, 2, ..., k. (4.7)
From (4.5) and (4.7)

E(T k,N1 |D i1i2...i k ) = k h=1 q + h-1 j=1 p ij p - h-1 j=1 p ij + k = k h=1 ( q + h-1 j=1 p ij p - h-1 j=1 p ij + 1) = 1 p + 1 p -p i1 + 1 p -p i1 -p i2 + ... + 1 p - k-1 j=1 p ij = 1 1 -q + 1 1 -q -p i1 + ... + 1 1 -q - k-1 j=1 p ij , (4.8)
where

q = i∈H2 p i = N i=N1+1 p i and p = i∈H1 p i = N1 i=1 p i .
Bearing in mind (4.3), (4.2) and (4.8) we can state the following:

Proposition 4.1. The expected value of T k,N1 is

E(T k,N1 ) = (i1i2...i k )∈Π k 1 p + 1 p -p i1 + 1 p -p i1 -p i2 + ... + 1 p - k-1 j=1 p ij k j=1 p ij p(p -p i1 )(p - 2 j=1 p ij )...(p - k-1 j=1 p ij ) , (4.9)
where Π k is the set of all k-permutations of set {1, 2, ..., N 1 }, i.e. the arrangements of length k of dierent elements of {1, 2, ..., N 1 }.

Thus, E(T k,N1 ) given by (4.9) is the expected number of draws required for k hosts of a set of unparasitized hosts H 1 ⊂ H with cardinality N 1 , to be parasitized. This value is generally dicult to obtain because the number of terms required for its computation is the number of k-permutations of 1, 2, ..., N 1 , that is N

(k) 1 = N 1 (N 1 -1)...(N 1 -k + 1
). This value is huge when N 1 and k are large. It is therefore important to obtain upper and lower bounds for this value.

Proposition 4.2. Let k be given and p 1 , p 2 , ...,p N1 be real numbers satisfying p 1 ≥ p 2 ≥ ... ≥ p N1 . Then, the maximum of E(T k,N1 |D i1i2...i k ) dened by (4.8) over all possible choices of the k-subsets {i

1 , i 2 , ..., i k } of H 1 is E(T k,N1 |D 1,2,...,k ) = 1 N1 i=1 p i + 1 N1 i=2 p i + ... + 1 N1 i=k p i , (4.10)
and the minimum is

E(T k,N1 |D N1,N1-1,...,N1-k+1 ) = 1 N1 i=1 p i + 1 N1-1 i=1 p i +...+ 1 N1-k+1 i=1 p i . (4.11) Proof . From hypothesis p 1 ≥ p 2 ≥ ... ≥ p N1 , it follows directly that N1 i=h p i ≤ N1 j=h p ij ≤ N1-h+1 i=1 p i , h = 1, 2..., N 1 , (4.12) then E(T k,N1 |D 1,2,...,k ) = 1 p + 1 p -p 1 + 1 p - 2 i=1 p i + ... + 1 p - k-1 i=1 p i ≥ 1 p + 1 p -p i1 + 1 p -p i1 -p i2 + ... + 1 p - k-1 j=1 p ij ≥ 1 p + 1 p -p N1 + 1 p - N1 i=N1-1 p i + ... + 1 p - N1 i=N1-k+2 p i = E(T k,N1 |D N1,N1-1,...,N1-k+1 )
and the proof is complete.

Proposition 4.3. Let p 1 , p 2 , ..., p N1 be real numbers satisfying 0 ≤ p i ≤ 1, for i = 1, 2, ..., N 1 and p 1 ≥ p 2 ≥ ... ≥ p N1 . It is then true that

E(T k,N1 |D 1,2,...,k ) ≥ E(T k,N1 ) ≥ E(T k,N1 |D N1,N1-1,...N1-k+1 )
In other words, E(T k,N1 |D 1,2,...,k ) and E(T k,N1 |D N1,N1-1,...N1-k+1 ) are upper and lower bounds, respectively, for the expected number of draws required for k hosts of the set H 1 to be parasitized.

Furthermore, the mode of the distribution on the events D i1i2...i k , (i 1 , i 2 , ..., i k ) ∈ Π k , is D 1,2,...,k , i.e. the order of parasitism of k hosts mostly likely to occur is 1, 2, ..., k.

Proof . The rst part of this proposition is a straightforward consequence of the previous proposition.

The second part comes directly from the fact that

P (D 1,2,...,k ) ≥ P (D i1i2...i k ) for (i 1 i 2 ...i k ) ∈ Π k .
which follows from (4.2) and (4.12). Propositions 4.2 and 4.3 prove that, if p 1 ≥ p 2 ≥ ... ≥ p N1 , then the most likely order of parasitization of k hosts in H 1 is the preferential order 1, 2, ..., k.

Moreover the shortest scenario (in terms of expectation) for the parasitization of k hosts of H 1 is the sequence extending from the least likely host, N 1 , to the most likely host, N 1 -k + 1, in the correct order. The longest scenario (in terms of expectation) for the parasitization of k hosts of H 1 extends from the most likely, 1, to the least likely host, k, in the correct order. These results can be intuitively explained as follows; let us suppose that host 1 is parasitized in the rst place. The probability of a new host of the set H 1 -{1} being parasitized is then q -p 1 . This value is less than any other value q -p j with j = 1. It is therefore more dicult for a host of the set H 1 -{1} to be parasitized than for a host of the set H 1 -{j}, j = 1, to be parasitized.

The repeated application of this reasoning explains the rst inequality of the proposition. The second inequality can be explained in a similar manner.

For simplicity, we denote T k,N by T k in the particular case in which N 1 = N . Recalling the denitions of these random variables and the random variables Y n , we obtain the following relations

P (Y n ≤ k -1) = P (T k > n), then P (Y n ≤ k -1) = 1 -P (T k ≤ n) and P (T k = n) = P (Y n-1 ≤ k -1) -P (Y n ≤ k -1).
From (3.3), (3.5) and above equalities we see that

P (T k > n) = A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n A , (4.13) P (T k ≤ n) = 1 - A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n A and P (T k = n) = A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n-1 A - A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n A = A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n-1 A (1 -p A ). E(T k ) = E(T k,N
) is the expected number of draws required for k hosts are parasitized. Dierent expressions have been described for this expectation [START_REF] Boneh | The coupon-colector problem revisited[END_REF][START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF]. From (4.13) it follows immediately that

E(T k ) = ∞ n=0 P (T k > n) = ∞ n=0   A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| p n A ,   = A⊂H,|A|≤k-1 (-1) k-1-|A| N -|A| -1 k -1 -|A| 1 1 -p A .
This expression was obtained in [START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF]. In [START_REF] Boneh | The coupon-colector problem revisited[END_REF] the following expression was obtained,

E(T k ) = k-1 r=0 u r ˆt≥0 N i=1
(1 + u(e pit -1))e -t dt, where x r f (x) is the coecient of x r in the power series development of f (x).

If k = N 1 = N , then E(T N ) = E(T N,N
) is the expected number of draws required to obtain complete parasitism. From (4.9)

E(T N ) = (i1i2...i N )∈Π N N -1 r=0 1 1 - r j=1 p ij N i=1 p i N k=1 N j=k p ij (4.14)
where Π N is the group of permutations of {1, 2, ...N } . This expression for E(T N ) is proposed in [START_REF] Boneh | The coupon-colector problem revisited[END_REF]. The authors provide no proof for this formula, and we have found no proof elsewhere.

Applications to various risk distributions

In this section, we consider two dierent risk distributions on the set of hosts and compute the most relevant values for every each.

The uniform distribution. The situation in which risk is distributed uniformly, i.e. all the hosts have the same probability of being parasitized, with:

p 1 = p 2 = ... = p N = 1 N (5.1)
has been widely studied. In this case, the expectation and variance of the random variable Y n representing the number of parasitized hosts after n draws are

E(Y n ) = N - (N -1) n N n-1 , V ar(Y n ) = (N -1)(N -2) n N n-1 + (N -1) n (N n-1 -(N -1) n ) N 2n-3 .
and the expected number of draws for k new hosts to be parasitized (4.9) is

E(T k,N1 ) = N 1 N 1 + 1 N 1 -1 + ... + 1 N 1 -k + 1 ,
which, in the case in which k = N , can be written as the following well-known formula

E(T N ) = N 1 + 1 2 + 1 3 ... + 1 N .
It is clear that in this case the upper and lower bounds for E(T k,N1 ) obtained in Proposition 4.3, are both equal to E(T k,N1 ), and all the probabilities

P (D i1i2...i k ) are equal to 1 N (k) 1 .
Two kinds of hosts. The two types of host situation is an idealization of the following cases. Hosts which are dead, either because they were previously parasitized or because they produced artifacts such as mines and galls, remain in the ecosystem for much longer than the existence of the host. They can make up to 90% of the host population. They can be still attractive to parasitoids long after the host death. Parasitoids will not lay eggs in them, but they will be checked carefully, implying a waste of time of up to 20% [START_REF] Casas | Foraging behaviour of a leafminer parasitoid in the eld[END_REF][START_REF] Casas | Parasitoid behaviour: predicting eld from laboratory[END_REF]). In such cases, it is possible to envision two categories, living and dead hosts, while being interested in the rate of parasitism of the living ones only.

Let us now consider the situation in which there are two kinds of hosts and, therefore, two dierent probabilities of being detected by a parasitoid.

In a population of N hosts, each of the hosts 1, 2, ..., m has a probability α of being parasitized, and each hosts m + 1, m + 2, ..., N has a probability β of being parasitized, such that

p 1 = p 2 = ... = p m = α, p m+1 = p m+2 = ... = p N = β.
(5.2)

The probability of host 1 being visited r 1 times, host 2 r 2 times, etc, for r 1 + r 2 + ... + r N = n, given by (2.2) is in this case

P (S n1 = r 1 , S n2 = r 2 , ..., S nN = r N ) = n! r 1 !r 2 !...r N ! α i≤m ri β i>m ri 0 ≤ r 1 ≤ n, 0 ≤ r 2 ≤ n, ..., 0 ≤ r N ≤ n, r 1 + r 2 + ... + r N = n.
The probability that, after n draws host i 1 had been chosen r i1 times by the parasitoids, host i 2 r i2 times and host i h r i h times, without taking the other hosts into account, is given by (2.3). It is equal to

P (S ni1 = r i1 , S ni2 = r i2 , ..., S ni h = r i h ) = n! r i1 !r i2 !...r i h !(n -r ij )! α i j ≤m ri j β i j >m ri j   1 - ij ≤m α - ij >m r ij β   n-ri j
.

We will now calculate the expected number of parasitized hosts after n draws with this risk distribution, using the results obtained in Section 2.

Let Y n be the random variable representing the number of parasitized hosts after n draws. From (3.3) it follows that

P (Y n = k) = k j=1 (-1) k-j N -j k -j j i=0 m i N -m j -i (iα -(j -i)β) n
and the expected value of Y n , (3.6), is equal to

E(Y n ) = N -m(1 -α) n -(N -m)(1 -β) n .
To compute the expected number of draws for k hosts of a set H 1 ⊂ H of unparasitized hosts to be parasitized, we will name the hosts of the set H 1 , hosts 1, 2, ..., N 1 . Without any loss of generality, we can assume p 1 = p 2 = ... = p m1 = α and p m1+1 = p m1+2 = ... = p N1 = β. Let Π k be the set of all k-permutations of the integers 1, 2, ..., N 1 . For every I = (i 1 , i 2 , ..., i k ) ∈ Π k , let A I ⊂ {1, 2, ..., k} be the set dened by j ∈ A I if i j ≤ m 1 . It is clear that the probability P (D i1,i2,...,i k ) = P (D I ) given by (4.2) is, in this case,

P (D I ) = α |A I | β k-|A I | p(p -γ 1 )(p - 2 j=1 γ j )...(p - k-1 j=1 γ j )
, where

γ j = α if j ∈ A I β if j / ∈ A I (5.3) Then, if A I = A I for I ∈ Π k and I ∈ Π k , it follows directly that P (D I ) = P (D I ).
We can therefore dene an equivalence relation on Π k as follows: I is related to I if A I = A I . We denote by I the equivalence class of I, and by Π k the set whose elements are the equivalence classes of the elements of Π k , that is

Π k = I : I ∈ Π k .
There are as many equivalence classes as subsets of {1, 2, ..., k} with cardinalities greater than or equal to max {0, k -n 1 }, where n 1 = N 1 -m 1 , and less than or equal to min {k, m 1 }, and the cardinalities of these equivalence classes are

Ī = m (h) 1 n (k-h) 1 if |A I | = h.
Given the above considerations, it is clear that E(T k,N1 ) can be written in this case as:

E(T k,N1 ) = I∈Π k 1 p + 1 (p -γ 1 ) + ... + 1 (p - k-1 j=1 γ j ) = α |A I | β k-|A I | p(p -γ 1 )(p - 2 j=1 γ j )...(p - k-1 j=1 γ j ) = Ī∈ Πk I∈ Ī 1 p + 1 (p -γ 1 ) + ... + 1 (p - k-1 j=1 γ j ) α |A I | β k-|A I | p(p -γ 1 )(p - 2 j=1 γ j )...(p - k-1 j=1 γ j ) = Ī∈ Πk m (|A I |) 1 n (k-|A I |) 1 1 p + 1 (p -γ 1 ) + ... + 1 (p - k-1 j=1 γ j ) α |A I | β k-|A I | p(p -γ 1 )(p - 2 j=1 γ j )...(p - k-1 j=1 γ j ) = min{k,m1} h=max{0,k-n1} |A I |=h m (h) 1 n (k-h) 1 1 p + 1 (p -γ 1 ) + ... + 1 (p - k-1 j=1 γ j ) α h β k-h p(p -γ 1 )(p - 2 j=1 γ j )...(p - k-1 j=1 γ j )
.

where γ j is dened by (5.3).

Let us suppose that α > β.

To obtain an upper bound for E(T k,N1 ),

we distinguish two cases, k ≤ m 1 and k > m 1 . If k ≤ m 1 then E(T k,N1 |D 1,2,...,k ) = 1 m 1 α + n 1 β + 1 (m 1 -1)α + n 1 β + ... + 1 (m 1 -k + 1)α + n 1 β , if k > m 1 , this upper bound is E(T k,N1 |D 1,2,...,k ) = 1 m 1 α + n 1 β + 1 (m 1 -1)α + n 1 β + ... + 1 n 1 β + 1 (n 1 -1)β + ... + 1 (n 1 + m 1 -k + 1)β .
Similarly, to obtain a lower bound for E(T k,N1 ) we distinguish the cases k ≤ n 1 and k > n 1 . If k ≤ n 1 this lower bound is

E(T k,N1 |D N1,N1-1,...,N1-k+1 ) = 1 m 1 α + n 1 β + 1 m 1 α + (n 1 -1)β + ... + 1 m 1 α + (n 1 -k + 1)β , and if k > n 1 , a lower bound for E(T k,N1 ) is E(T k,N1 |D N1,N1-1,...,N1-k+1 ) = 1 m 1 α + n 1 β + 1 m 1 α + (n 1 -1)β + ... + 1 m 1 α + 1 (m 1 -1)α + 1 (n 1 + m 1 -k + 1)α .
The maximum of the values P (D i1,i2,...,i k ) is

P (D 1,2,...,k ) =          α k k-1 h=0 ((m 1 -h)α + n 1 β) , if k ≤ m 1 α m1 β k-m1 m1 h=0 ((m 1 -h)α + n 1 β) k-m1-1 l=1 (n 1 -l)β , if k > m 1
In the extreme case that there is only one host with a probability α of being parasitized and the others have probability β of being parasitized, we obtain the following expressions for E(T k,N1 ).

If the host with probability α of being parasitized does not belong to set H 1 , then

E(T k,N1 ) = 1 β k-1 j=0 1 N 1 -h .
If the host with probability α of being parasitized belongs to set H 1 , then

E(T k,N1 ) = (N 1 -1) (k) 1 α + (N 1 -1)β + 1 α + (N 1 -2)β + ... + 1 α + (N 1 -k)β β k (α + (N 1 -1)β)(α + (N 1 -2)β)...(α + (N 1 -k)β) + (N 1 -1) (k-1) k-1 j=1 1 α + (N 1 -1)β + 1 α + (N 1 -2)β + ... + 1 α + (N 1 -j)β + 1 (N 1 -j)β + 1 (N 1 -j -1)β + ... + 1 (N 1 -k + 1)β αβ k-1 (α + (N 1 -1)β)(α + (N 1 -2)β)...(α + (N 1 -j)β)(N 1 -j)β(N 1 -j -1)β...(N 1 -k + 1)β + (N 1 -1) (k-1) 1 α + (N 1 -1)β + 1 α + (N 1 -2)β + ... + 1 α + (N 1 -k)β αβ k-1 (α + (N 1 -1)β)(α + (N 1 -2)β)...(α + (N 1 -k)β) .
6 Relationship between the risk distribution and the speed of parasitization

In the preceding sections, we studied the process of parasitization for a given risk distribution in the set of hosts. In this section we compare this process for dierent risk distributions. We show how parasitization speed depends on the risk distribution, and its scatter in particular. We use the concept of majorization to formalize the idea that risk distributions have dierent degrees of spread. This notion dates from the start of the 20th century. A comprehensive review of the theory can be found in [START_REF] Marshall | Inequalities: Theory of majorization and its applications[END_REF].

Less spread distributions are associated with faster parasitization. In other words, the more spread out the risk distribution, the larger the number of draws required for a given number of hosts to be parasitized. Thus the distribution function for the rst time parasitization of a given number of hosts, viewed as a function of the vector p, is Schur convex (see the denition at the end of this section). The mathematical community studying the coupon collector problem seems to be largely unaware of it, but this result is not new and can be found in [START_REF] Wong | A majorization theorem for the number of distinct outcomes in N independent trials[END_REF]. This result constitutes the rst part of Theorem 6.1.

We give a proof more concise and clearer than previous proposal. Moreover, our method provides a precise result for strict Schur convexity. This renement constitutes the second part of Theorem 6.1. We make use in our proof of the relationship between the concept of majorization and the numerical operation known as "Robin Hood transfer", described below.

is satised, where Then, we obtain, with a = N -2,

n = k -1, b = N -k 0≤r≤k-1 (-1) k-1-r N -2 r N -2 -r k -1 -r = k -2 k -1 = 0,
and the lemma follows.

Lemma 6.2. Let q 1 , q 2 , ...,q M be non-negative real numbers and I = {1, 2, ..., M }. For every A ⊂ I let q A = i∈A q i . Then, for any integer m ≥ 0 , A⊂I,|A|≤r

(-1) r-|A| M -|A| r -|A| q m A ≥ 0.
Moreover, if m ≥ r and at least r of the values q 1 , q 2 , ...,q M are greater than zero, then A⊂I,|A|≤r

(-1) r-|A| M -|A| r -|A| q m A > 0.
Proof . If all the q i are zero, there is nothing to prove. Let us suppose that s = M i=1 q i > 0. Let p i = q i /s, i = 1, 2, ..., M . These values dene the probability distribution p = (p 1 , p 2 , ..., p M ) on I. From (3.3) it follows that A⊂I,|A|≤r

(-1) r-|A| M -|A| r -|A| q m A = s m A⊂I,|A|≤r (-1) r-|A| M -|A| r -|A| p m A = s m P p (Y m = r),
which proves the lemma.

Let p = (p 1 , p 2 , ..., p N ) denote a probability distribution p over the set H. Suppose that p is not uniform. We can assume p 1 < p 2 without loss of generality. Let 0 < h ≤ p2-p1

2 , α = 1 -h p2-p1 . We then dene a new risk distribution p by applying a Robin Hood transfer as follows p = (p 1 + h, p 2 -h, p 3 , p 4 , ..., p N ) = (αp 1 + (1 -α)p 2 , αp 2 + (1 -α)p 1 , p 3 , ..., p N ).

(6.1)

We indeed have p ≺ p.

Theorem 6.1. Let p be a non uniform probability distribution over H. Without loss of generality, we can assume that p 1 < p 2 . Let p be dened by (6.1). Then, for all k between 1 and N -1,

P p (Y n ≤ k) ≥ P p (Y n ≤ k), (6.2) which is equivalent to P p (T k+1 ≤ n) ≤ P p (T k+1 ≤ n) (6.3)
Moreover, if at least k -1 of the values p 3 , p 4 , ..., p N are non-zero, then

P p (Y n ≤ k) > P p (Y n ≤ k), n = k + 1, k + 2, k + 3... (6.4) which is equivalent to P p (T k+1 ≤ n) < P p (T k+1 ≤ n), n = k + 1, k + 2, k + 3... (6.5)
where p is dened by (6.1).

Proof . Let H = {3, 4, ..., N }. According to (3.5) we have:

P p (Y n ≤ k) = A⊂H,|A|≤k (-1) k-|A| N -|A| -1 k -|A| p n A = A⊂H ,|A|≤k (-1) k-|A| N -|A| -1 k -|A| p n A + A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 1 ) n + A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 2 ) n + A⊂H ,|A|≤k-2 (-1) k-2-|A|) N -3 -|A| k -2 -|A| (p A + p 1 + p 2 ) n . Then P p (Y n ≤ k) -P p (Y n ≤ k) = A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 1 ) n + A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 2 ) n - A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 1 + h) n - A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + p 2 -h) n .
Let f be the real function dened by

f (x) = A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| (p A + x) n , x ∈ R
This function is a polynomial of degree less than or equal to n. The coecient of x n is equal to

A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| = 0≤r≤k-1 (-1) k-1-r N -2 r N -2 -r k -1 -r
and this is equal to 0 by Lemma 6.1. The coecient of x n-j for j = 1, 2, ..., n is n j

A⊂H ,|A|≤k-1 (-1) k-1-|A| N -2 -|A| k -1 -|A| p j A ,
by the rst part of Lemma 6.2 with I = H = {3, 4, ..., N }, M = |I| = N -2, m = j and r = k -1, it follows that these coecients are greater than or equal to zero. This polynomial function is then convex on [0, +∞), so that

f (p 1 ) + f (p 2 ) ≥ f (αp 1 + (1 -α)p 2 ) + f (αp 2 + (1 -α)p 1 ) = f (p 1 + h) + f (p 2 -h), f (p 1 ) + f (p 2 ) -f (p 1 + h) -f (p 2 -h) ≥ 0.
However, this inequality is the same as

P p (Y n ≤ k) -P p (Y n ≤ k) ≥ 0,
which gives (6.2). Recalling the relationship between the random variables Y i and the random variables T j , we also obtain

P p (T k+1 ≤ n) ≤ P p (T k+1 ≤ n),
which is (6.3).

Moreover, from the second part of Lemma 6.2. it follows that if at least k -1 of the values p 3 , p 4 , ..., p N are greater than zero and n ≥ k + 1, then the coecient of x n-k+1 is greater than zero, where n -k + 1 ≥ 2. So, at least one monomial of degree greater than or equal to 2 appears in the polynomial. The We can state the following corollaries.

Corollary 6.1. Let p = (p 1 , p 2 , ..., p N ) and q = (q 1 , q 2 , ..., q N ) be risk distributions on H = {1, 2, ..., N }. If p ≺ q then, for every n ≥ 1 and every k ≥ 1 P p (Y n ≤ k) ≤ P q (Y n ≤ k), (6.6) is satised and P p (T k+1 ≤ n) ≥ P q (T k+1 ≤ n).

Furthermore, if the distributions p and q are actually dierent, meaning that they do not dier only by a permutation, then the preceding inequalities are strict, except in trivial cases. More precisely, denoting by j the number of non zero p i values (and remarking that the number of non-zero q i values is at most j), we have: If k ≥ n or k ≥ j then P p (Y n ≤ k) = P q (Y n ≤ k) = 1 and P p (T k+1 ≤ n) = P q (T k+1 ≤ n) = 0;

If n ≥ 2, k < n and k < j, then P p (Y n ≤ k) < P q (Y n ≤ k) and P p (T k+1 ≤ n) > P q (T k+1 ≤ n).

Proof . As it is possible to go from vector q to vector p by a nite sequence of Robin Hood transfers, the corollary follows directly from Theorem 6.1, which proves that each transfer decreases the quantity P p (Y n ≤ k). We just have to consider the cases in which this quantity is strictly decreased. Remark 6.1. We can interpret the results obtained above in terms of the theory of Schur-convex functions. A real-valued function φ dened on a set A ⊂ R N is said to be Schur-convex on A if, for every x and y pair of elements in A such that x ≺ y the inequality φ(x) ≤ φ(y) holds. The rst part of Corollary 6.1 states that the map p → P p (Y n ≤ k) is Schur-convex. This was dominates the random variable X´if there is a critical value c ∈ R such that, for any d ∈ D

•if d ≤ c, then P (X = d) ≤ P ´(X´= d),

•if d > c, then P (X = d) ≥ P ´(X´= d).

It is easy to show that strong dominance implies weak dominance, but that the converse is not true. Coming back to our CCP model, we propose the following:

Conjecture. If p ≺ q, then the random variable Y n dened on the probability space (Ω, P p ) strongly stochastically dominates the random variable Y n dened on the probability space (Ω, P q ) . This conjecture has been tested on various examples, but we have been able to prove it formally for only a few values of the pair (n, N ), namely for n = 2 or 3 and any N , and for n = 4 and N ≤ 5.

In applications, strong dominance reinforces weak dominance. It gives more precise statements concerning the relative probabilities that a given number of hosts are parasitized after a given number of eggs laid, for two risk distributions.

Figure 1 .

 1 Figure 1. Attacks of caterpillar larvae hosts by parasitic wasps as an urn process in discrete time n. A wasp oviposit 10 eggs among ve hosts (n = 10) (A). The outcome of the fundamental multinomial process (B) is summarized in the marginal distribution of the number of eggs per individual host at a given time S 10 = (3, 1, 3, 0, 3) (C), in the frequency distribution of eggs among hosts (D) and in the number of hosts attacked over time Y n (E)

  its probability mass function (3.3), distribution function (3.5), expectation (3.6) and variance (3.7).

  convexity is then strict, and we can writef (p 1 ) + f (p 2 ) -f (p 1 + h) -f (p 2 -h) > 0, and P p (Y n ≤ k) -P p (Y n ≤ k) > 0, n = k + 1, k + 2, k + 3... which is equivalent to P p (Y n ≤ k) > P p (Y n ≤ k), n = k + 1, k + 2, k + 3... and therefore to P p (T k+1 ≤ n) < P p (T k+1 ≤ n), n = k + 1, k + 2, k + 3...This completes the proof.

Figure 2 .

 2 Figure 2. Graphic A: Distribution functions of ve variables Y 12 corresponding to ve dierent risk distributions, p 1 , p 2 ,..., p 5 , satisfying p 1 ≺ p 2 ≺ p 3 ≺ p 4 ≺p 5 . It can be observed that P pi (Y 12 ≤ k) < P pi+1 (Y 12 ≤ k), for k = 1, 2, ..., 11, i = 1, 2, 3, 4. Graphic B: Distribution functions of ten variables T 6 corresponding to ten risk distributions, p 1 , p 2 , ..., p 10 satisfying p 1 ≺ p 2 ≺ ... p 9 ≺p 10 . Graphic C: Distribution functions of ten variables T 9 corresponding to the same previous risk distributions. In graphics B and C it can be observed that P pi (T k ≤ n) > P pi+1 (T k ≤ n), for n = k, k + 1, ... i = 1, 2, ...9

Figure 3 .

 3 Figure 3. Comparison of distribution functions of random variables T k corresponding to two unrelated risk distributions, i.e. neither p ≺ q nor q ≺ p, to show how these distribution functions act in dierent ways depending on the value of k.

Acknowledgments N.Z. and M.J.F.S acknowledge the nancial support of the Fundación Seneca of the Comunidad Autónoma de la Región de Murcia, project 19320/IP/14. N.Z. is also grateful to the University François-Rabelais of Tours, for its support and hospitality.

In this section, we work with dierent risk distributions, requiring further notation and denitions. Given a risk distribution p = (p 1 , p 2 , ..., p N ), we denote by P p the probability distribution induced by p on the σ-eld over the space of the all the possible incidences of parasitization.

Given (p 1 , p 2 , ..., p N ) in R N , we denote by (p1, p2, ..., p N ) the N -uple obtained by permutation of p i such that p1 ≥ p2 ≥ . . . ≥ p N .

The following denitions are given in [START_REF] Marshall | Inequalities: Theory of majorization and its applications[END_REF]. Definition 6.1. Let p 1 , p 2 , ..., p N , q 1 , q 2 , ..., q N , be real numbers. We say that p = (p 1 , p 2 , ..., p N ) is majorized by q = (q 1 , q 2 , ..., q N ), and we write p ≺ q,

It is clear that when we apply this denition to the comparison of two risk distributions, the last equality is trivially satised.

Let q = (q 1 , q 2 , ..., q N ) ∈ R N . If q h < q k we can transfer an amount ∆, 0 < ∆ < q k -q h from q k to q h to obtain the following new risk distribution q = (q 1 , q 2 , ...,q N ), where q h = q h + ∆, q k = q k -∆ and q i = q i for i = h, k. Then, q is less spread out than the initial distribution, that is, q ≺ q. Such operations involving the shifting of some income from one individual to a poorer individual, are described, somewhat poetically, as Robin Hood transfers [START_REF] Arnold | Majorization and the Lorenz Order: A Brief Introduction[END_REF]. If we dene α = 1 -∆ q k -q h then we can write q h = q h + ∆ = αq h + (1 -α)q k and q k = q k -∆ = αq k + (1 -α)q h . Proposition 6.1. The following conditions are equivalent: a) p ≺ q, b) p can be derived from q by successive applications of a nite number of Robin Hood transfers.

It is not dicult to prove this equivalence. It was proved for the rst time in [START_REF] Muirhead | Some methods applicable to identities and inequalities of symmetric algebraic functions of n letters[END_REF] for vectors of non-negative integer components. Lemma 6.1. Let k and N be integers satisfying

already proved in [START_REF] Wong | A majorization theorem for the number of distinct outcomes in N independent trials[END_REF], and was stated as a conjecture in [START_REF] Anceaume | New results on a generalized coupon collector problem using Markov chains[END_REF].

Corollary 6.2. Let u = (1/N, 1/N, ..., 1/N ) be the uniform distribution on H = {1, 2, ..., N } and p = (p 1 , p 2 , ..., p N ) any other risk distribution on H.

Proof . It can be clearly seen that u = (1/N, 1/N, ..., 1/N ) is majorized by any other distribution on H and the corollary follows.

Remark 6.2. The results obtained in Corollary 6.1 and Corollary 6.2 can be expressed in terms of a comparison of probability distributions as follows.

If p ≺ q, then relation (6.6) proves that the random variable Y n dened on the probability space determined by p on the space of the random sets of H = {1, 2, ..., N } is weakly stochastically dominated by the random variable Y n dened on the probability space determined by q. Corollary 6.2 proves that the random variable Y n dened on the probability space determined by the uniform distribution u = (1/N, 1/N, ..., 1/N ) is always weakly stochastically dominated by the random variable Y n dened on the probability space determined by any other probability distribution on H. Remark 6.3. After the redaction of this section, we have seen a similar study in [START_REF] Anceaume | Optimization results for a generalized coupon collector problem[END_REF]. In particular, they prove inequalities (30) and (31) of Theorem 6.1. However, our contribution still presents a real interest, thanks to the quality of the argument based on use of fundamental formulas ( 6) and ( 7) in dierent contexts, and because we obtain cases of strict inequalities.

Illustrative examples

In this section we show graphically the relationships satised among the distribution functions of random variables Y n as well as the distribution functions of random variables T k , when their corresponding risk distributions are able to be compared by majorization.

The distribution functions of ve variables Y n are represented in graphic A of Figure 2. They correspond to ve dierent risk distributions, p 1 , p 2 , p 3 , p 4 , and p 5 , satisfying p 1 ≺ p 2 ≺ p 3 ≺ p 4 ≺p 5 . These are distributions on the set {1, 2, ..., 12} (so N = 12), p 1 is the uniform distribution, p i = ( 1 /10i+2, ..., 1 /10i+2, 10(i-1)+1 /10i+2) for i = 2 and 3, and p i = ( 1 /45(i-3)+12, ..., 1 /45(i-3)+12, 1 /45(i-3)+12, 45(i-3)+1 /45(i-3)+12) for i = 4 and 5. We have also used n = 12, and it can be observed that P pi (Y 12 ≤ k) < P pi+1 (Y 12 ≤ k), with k = 1, 2, ..., 11, i = 1, 2, 3, 4.

The distribution functions of ten variables T k are represented in every one of the graphics B and C in Figure 2. N = 10 and the risk distributions are the same in both cases; p 1 is the uniform distribution and p i = ( 1 /5(i+1), ..., 1 /5(i+1), i /5(i+1), 4(i-1)+1 /5(i+1)), for i = 2, 3, ..., 10. For these risk distributions p 1 ≺ p 2 ≺ ... p 9 ≺p 10 is satised. In graphic B of Figure 2, k = 6 and the values of n lie between 6 to 50. In graphic C of Figure 2, k = 9 and the values of n lie from 9 to 100. It can be seen that P pi (T k ≤ n) > P pi+1 (T k ≤ n), for n = k, k + 1, ... i = 1, 2, ...9, in both graphics.

Figure 3 compares distribution functions of random variables T k corresponding to two unrelated risk distributions p and q, i.e. neither p ≺ q nor q ≺ p.

Thus, these distribution functions act in dierent ways depending on the value of k. We include three dierent graphics, each bearing two curves. These curves are the distribution functions of two random variables T k . The risk distributions associated with these random variables are, in the three graphics, p = ( 3 /85, 3 /85, 3 /85, 3 /85, 3 /85, 12 /85, 13 /85, 13 /85, 13 /85, 19 /85) and q = ( 3 /81, 4 /81, 4 /81, 4 /81, 4 /81, 5 /81, 5 /81, 5 /81, 15 /81, 32 /81). In the rst graphic k = 5, in the second k = 8 and in the third k = 9. In the last two cases the distribution functions cross.

They do not cross in the rst.

A conjecture on strong dominance

In Section 6 we used an order relationship between random variables (or more precisely between their distributions) that can be dened formally as follows.

Definition 8.1. Let X and X´be two real random variables, dened on probability spaces (Ω, P ) and (Ω , P ), respectively. We say that the random variable X weakly stochastically dominates the random variable X´if the cumulative distribution function of X´dominates the cumulative distribution function of X, that is, for any t ∈ R, P (X ≤ t) ≤ P ´(X´≤ t) .

The main result of Section 6 is that if p ≺ q, then the random variable Y n dened on the probability space (Ω, P p ) weakly stochastically dominates the random variable Y n dened on the probability space (Ω, P q ) .

A particular case of weak dominance is that one in which inequalities apply not only to the cumulative distribution functions, but also to the distributions themselves. We will refer to this situation as strong dominance, and we provide a formal denition of strong dominance below, for the case of discrete random variables. (A similar denition can be given for continuous random variables with densities). In short, X strongly dominates X´if, for any small enough value d, P (X = d) ≤ P ´(X´= d), and if for any other possible value e, P (X = e) ≥ P ´(X´= e) . Definition 8.2. Let X and X´be two real random variables, dened on probabilities spaces (Ω, P ) and (Ω , P ), respectively, and taking values in a denumerable set D. We say that the random variable X strongly stochastically


Pedro Machado Manhães De Castro

Olivier Devillers

Expected Length of the Voronoi Path in a High Dimensional Poisson-Delaunay Triangulation

Keywords:

Let X n be a d dimensional Poisson point process of intensity n. We prove that the expected length of the Voronoi path between two points at distance 1 in the Delaunay triangulation associated with X n is

2 ) for all n ∈ N and d → ∞. In any dimension, we provide a precise interval containing the exact value, in 3D the expected length is between 1.4977 and 1.50007.

Introduction

Finding paths in a Delaunay triangulation is a classical problem in computational geometry [START_REF] Devillers | Walking in a Triangulation[END_REF]. In the context of random points, several kind of paths have been studied in two dimensions such as straight walk [START_REF] Bose | On the stabbing number of a random Delaunay triangulation[END_REF][START_REF] Devroye | Expected time analysis for Delaunay point location[END_REF], cone walk [START_REF] Broutin | Efficiently navigating a random Delaunay triangulation[END_REF], visibility walk [START_REF] Devillers | The worst visibility walk in a random Delaunay triangulation is O( √ n)[END_REF], the shortest path [START_REF] Chenavier | Stretch Factor of Long Paths in a planar Poisson-Delaunay Triangulation[END_REF], or the Voronoi path [START_REF] Baccelli | Markov paths on the Poisson-Delaunay graph with applications to routing in mobile networks[END_REF].

In this paper we take interest in the stretch ratio of a particular path in the Delaunay triangulation -the Voronoi path -and study its expected length in dimension d when the point set is a Poisson point process of density n with n going to infinity. The Voronoi path links the seeds of the Voronoi regions intersected by a line segment. An illustration for dimension 2 is given in Figure 1. The main result of this paper is the computation of upper and lower bound on the expected length of the Voronoi path. These bounds show that the asymptotic behavior of the length is 2 ). Table 1 provides the values of our bounds for small dimensions aside the exact value obtained from numerical integration.

» 2d π + O(d - 1
Obviously, the length of the Voronoi path give an upper bound on the length of the shortest path and provide an upper bound on the expected stretch ratio of long walk in the Delaunay triangulation of a random point set.

Previous results provide values only for dimension two, where the expected length of the Voronoi path is 4 π 1.27 and its variance is O(n -1 2 ) [START_REF] Baccelli | Markov paths on the Poisson-Delaunay graph with applications to routing in mobile networks[END_REF][START_REF] Devillers | Shortcuts in the Voronoi path in a planar Poisson-Delaunay triangulation[END_REF].

Voronoi Path

We define the Voronoi path V P χ as the list of closest neighbors in a set of points χ of a point moving linearly from s = (0, 0, 0) to t = (1, 0, 0). This path is x-monotone, it starts at the closest neighbor of s and reaches the closest neighbor of t. It uses a sequence of edges of DT χ the Delaunay triangulation of χ.

Notice that if somebody wants a path that goes actually from s to t, he can consider V P χ∪{s,t} . This path differs from V P χ only by few edges around s and around t.

An edge p 1 p 2 belongs to V P χ iff the unique ball B x (p 1 , p 2 ) centered on the x axis passing through p 1 and p 2 is centered on the segment [st] and does not contain any other points of χ (see Figure 1). Thus, denoting M (p 1 , p 2 ) the center of B x (p 1 , p 2 ), we can write the length of the t s

M (p1, p2) Bx(p1, p2) p 1 p 2

Voronoi path as

(V P χ ) = 1 2 (p1,p2)∈χ 2 1 [Bx(p1,p2)∩χ=∅] 1 [M (p1,p2)∈[st]] p 1 p 2
where the 1 2 arise because each edge is counted twice, once for each orientation. Now we turn our interest to the case where the point set is X n a Poisson point process of intensity n. Using Slivnyak-Mecke formula [9, Theorem 3.3.5], we have:

E [ (V P Xn )] = n 2 2 (R 3 ) 2 P [B x (p 1 , p 2 ) ∩ X n = ∅] 1 [M (p1,p2)∈[st]] p 1 p 2 dp 1 dp 2 . (1)
The path V P Xn∪{s,t} has about the same length, actually E (V P Xn ) -

(V P Xn∪{s,t} ) = O Ä n -1 d ä
can be proven as an easy generalization from dimension two [START_REF] Devillers | Shortcuts in the Voronoi path in a planar Poisson-Delaunay triangulation[END_REF].

Voronoi Path in Dimension 3

We start by illustrating our method in three dimension.

Lemma 1. In dimension 3

1.497706663 788984278470257640690697143 745000536337515228912680960 √ 2 ≤ E [ (V P Xn )] ≤ 4523370364712510658076963509 4264485828690604413776035840 √ 2 1.500066356.
Proof. The integral at Equation ( 1) is computed by substitution. The points p 1 and p 2 are defined by their sphere B x (p 1 , p 2 ) and their spherical coordinates on that sphere. Let Φ be the function

Φ : R × R + × ([0, π) × [0, 2π)) 2 -→ (R 3 ) 2 (x, r, α 1 , β 1 , α 2 , β 2 ) -→ (p 1 , p 2 ),
with

p i = Ñ x + r cos α i r sin α i cos β i r sin α i sin β i é = (x, 0, 0) + ru i with u i = Ñ cos α i sin α i cos β i sin α i sin β i é .
Φ is a C 1 -diffeomorphism up to a null set. Its Jacobian J Φ has as determinant:

det(J Φ ) = 1 cos α 1 -r sin α 1 0 0 0 0 sin α 1 cos β 1 r cos α 1 cos β 1 -r sin α 1 sin β 1 0 0 0 sin α 1 sin β 1 r cos α 1 sin β 1 r sin α 1 cos β 1 0 0 1 cos α 2 0 0 -r sin α 2 0 0 sin α 2 cos β 2 0 0 r cos α 2 cos β 2 -r sin α 2 sin β 2 0 sin α 2 sin β 2 0 0 r cos α 2 sin β 2 r sin α 2 cos β 2 = r 4 sin α 1 sin α 2 (cos α 1 -cos α 2 )
Inria Then we substitute the new variables:

1 E [ (V P Xn )] = n 2 2 ∞ -∞ ∞ 0 π 0 π 0 2π 0 2π 0 P [B((x,0,0), r) ∩ Xn = ∅] 1 [(x,0,0)∈[st]] •r u 1 u 2 | det(J Φ )|dβ 1 dβ 2 dα 1 dα 2 drdx = n 2 2 ∞ -∞ ∞ 0 π 0 π 0 2π 0 2π 0 e -n 4 3 πr 3 1 [x∈[0,1]] • r u 1 u 2 •r 4 sin α 1 sin α 2 | cos α 1 -cos α 2 |dβ 1 dβ 2 dα 1 dα 2 drdx = n 2 2 Å 1 0 dx ã Å ∞ 0 e -n 4 3 πr 3 r 5 dr ã E q. ( 7
) • π 0 π 0 2π 0 2π 0 sin α 1 sin α 2 | cos α 1 -cos α 2 | • u 1 u 2 dβ 1 dβ 2 dα 1 dα 2 = n 2 2 3 16n 2 π 2 π 0 π 0 2π 0 2π 0 sin α 1 sin α 2 | cos α 1 -cos α 2 |• u 1 u 2 dβ 1 dβ 2 dα 1 dα 2 (2)
Using the trivial bound u 1 u 2 ≤ 2 we get

π 0 π 0 2π 0 2π 0 sin α 1 sin α 2 | cos α 1 -cos α 2 | • 2dβ 1 dβ 2 dα 1 dα 2 = Å 2 π 0 α 2 0 sin α 1 sin α 2 (cos α 1 -cos α 2 ) dα 1 dα 2 ã • Å 2π 0 2π 0 dβ 1 dβ 2 ã = 8 3 E q. ( 9
) • 4π 2 = 32π 2 3 .
Plugging this bound in Equation ( 2) gives E [ (V P Xn )] ≤ 2.

We can improve on this result. Expressing u 1 u 2 in terms of the spherical coordinates we have

u 1 u 2 = (u 1 -u 2 ) 2 = u 2 1 + u 2 2 -2u 1 • u 2 = √ 2 -2u 1 • u 2 = √ 2 1 -cos α 1 cos α 2 -sin α 1 sin α 2 cos(β 1 -β 2 ) (3) = √ 2 1 -δ with δ = (cos α 1 cos α 2 + sin α 1 sin α 2 cos(β 1 -β 2 )) .
We can use Taylor expansion to bound u 1 u 2 . If T k (y) is the Taylor expansion up to degree k of √ 1 -y, for k odd and y ∈ [-1, 1] we have

T k (y) -T k (1)y k+1 ≤ 1 -y ≤ T k (y). ( 4
)
The side figure illustrates this for k = 3:

1 -y 2 -y 2 8 -y 3 16 -5y 4 16 ≤ 1 -y ≤ 1 -y 2 -y 2 8 -y 3 16 .
One can compute the following integrals (e.g., using Maple):

2 π 0 α 2 0 2π 0 2π 0 sin α 1 sin α 2 | cos α 1 -cos α 2 |T 41 (δ)dβ 1 dβ 2 dα 1 dα 2 = 262994759490085880230232381 23281266760547350903521280 π 2 2 π 0 α 2 0 2π 0 2π 0 sin α 1 sin α 2 | cos α 1 -cos α 2 | T 41 (δ)-T 41 (1)δ 42 dβ 1 dβ 2 dα 1 dα 2 = 1507790121570836886025654503 133265182146581387930501120 π 2 .
Plugging these values in Equation (2) using Equation ( 3) and ( 4) yields the bounds announced in the lemma statement.

A numerical evaluation of the integral at Equation ( 2) gives a value in that interval, pretty close to 3 2 and we conjecture that 3 2 is the correct value.

Voronoi Path in Dimension d

Theorem 2. In dimension d, the expected length of the Voronoi path is bounded by:

Γ d 2 4 2 4d-5 d π 2 (2d -2)! Å 1 - d -1 4d 2 -1 ã √ 2 ≤ E [ (V P Xn )] ≤ Γ d 2 4 2 4d-5 d π 2 (2d -2)! Å 1 + 1 4d -2 ã √ 2
and the behavior when d → ∞ is:

… 2d π - 1 4 √ 2dπ + O Ä d -3 2 ä ≤ E [ (V P Xn )] ≤ … 2d π + 3 4 √ 2dπ + O Ä d -3 2 ä
The rest of the section is devoted to the proof of Theorem 2.

As in dimension 3, we compute by substitution the integral in Equation (1) defining the points p 1 and p 2 by their sphere B x (p 1 , p 2 ) and their spherical coordinates on that sphere. Let Φ be the function

Φ : R × R + × [0, π] d-1 × [0, 2π) 2 -→ (R d ) 2 (x, r, α 1,1 , . . . , α 1,d-1 , α 2,1 , . . . , α 2,d-1 ) -→ (p 1 , p 2 ), with p i =            x 0 0 . . . 0 . . . 0 0            + ru i with u i =              cos α i,1 sin α i,1 cos α i,2 sin α i,1 sin α i,2 cos α i,3 . . . j-1 l=1 sin α i,l cos α i,j . . . d-2 l=1 sin α i,l cos α i,d-1 d-2 l=1 sin α i,l sin α i,d-1             
.

Φ is a C 1 -diffeomorphism up to a null set. Its Jacobian J Φ has as determinant: 7)

det(JΦ) = r 2(d-1) d-2 i=2 sin d-i-1 (α1,i) sin d-i-1 (α2,i) • sin d-2 (α1,1) sin d-2 (α2,1)| cos α1,1 -cos α2,1|. Inria E [ (V P Xn )] = n 2 2 ∞ -∞ ∞ 0 (S d-1 ) 2 P [B((x,0,...,0), r) ∩ Xn = ∅] 1 [(x,0,...,0)∈[st]] •r u 1 u 2 | det(J Φ )|dα 1,1 . . . dα 1,d-1 dα 2,1 . . . dα 2,d-1 drdx = n 2 2 1 0 dx × ∞ 0 e -n π d/2 Γ( d 2 +1) r d vo lu m e in S d-1 r 2d-1 dr × (S d-1 ) 2 u 1 u 2 d-2 i=2 sin d-i-1 (α 1,i ) sin d-i-1 (α 2,i ) • sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1 = Γ d 2 + 1 2 2dπ d E q. (
× (S d-1 ) 2 u 1 u 2 d-2 i=2 sin d-i-1 (α 1,i ) sin d-i-1 (α 2,i ) (5)
• sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1
We use the Taylor expansion of Equation 4with k = 1 to bound u 1 u 2 :

√ 2 Å 1 - δ 2 - δ 2 2 ã ≤ u 1 u 2 ≤ √ 2 Å 1 - δ 2 ã (6)
with

δ = u 1 • u 2 = d-2 i=1 i-1 j=1 sin α 1,j sin α 2,j cos α 1,i cos α 2,i + d-2 j=1 sin α 1,j sin α 2,j • cos(α 1,d-1 -α 2,d-1 ).
We write δ as a sum of three terms

δ = δ A + δ B + δ C with δ A = d-2 j=1 sin α 1,j sin α 2,j • cos(α 1,d-1 -α 2,d-1 ), δ B = d-2 i=2 i-1 j=1 sin α 1,j sin α 2,j cos α 1,i cos α 2,i , δ C = cos α 1,1 cos α 2,1 .
Replacing u 1 u 2 by its bounds in Equation ( 5) yields the computation of the following integral for k ∈ {0, 1, 2}:

I k = (S d-1 ) 2 δ k d-2 i=2 sin d-i-1 (α 1,i ) sin d-i-1 (α 2,i ) • sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1 4.1 Computation of I 0 I 0 = 2π 0 2π 0 dα 1,d-1 dα 2,d-1 × d-2 i=2 π 0 sin d-i-1 (α 1,i )dα 1,i × d-2 i=2 π 0 sin d-i-1 (α 2,i )dα 2,i × π 0 π 0 sin d-2 (α 1,1 ) sin d-2 (α 2,1 )| cos α 1,1 -cos α 2,1 |dα 1,1 dα 2,1 = 4π 2 × d-2 i=2 Γ d-i 2 Γ d-i+1 2 2 π d-3 E q. (8 ) × π 0 π 0 sin d-2 (α 1,1 ) sin d-2 (α 2,1 )| cos α 1,1 -cos α 2,1 |dα 1,1 dα 2,1 = 4π 2 × 1 Γ d-1 2 2
te le sc op in g

× π d-3 × 2 2d ((d -2)!) 2 (2d -2))! E q. (9 ) = 2 2d+2 π d-1 ((d -2)!) 2 Γ d-1 2 2 (2d -2))!

Computation of I 1

Decomposing on the different terms of δ, we can write

I 1 = I A 1 + I B 1 + I C 1 . I A 1 = (S d-1 ) 2 d-2 j=1 sin α1,j sin α2,j • cos(α 1,d-1 -α 2,d-1 ) d-2 j=2
sin d-j-1 (α1,j) sin d-j-1 (α2,j)

• sin d-2 (α1,1) sin d-2 (α2,1) • | cos α1,1 -cos α2,1|dα 1:2,1:d-1 I B 1 = d-2 i=2 (S d-1 ) 2 i-1 j=1
sin α1,j sin α2,j cos α1,i cos α2,i

d-2 j=2
sin d-j-1 (α1,j) sin d-j-1 (α2,j)

• sin d-2 (α1,1) sin d-2 (α2,1) • | cos α1,1 -cos α2,1|dα 1:2,1:d-1 I C 1 = (S d-1 ) 2 cos α1,1 cos α2,1 d-2 j=2
sin d-j-1 (α1,j) sin d-j-1 (α2,j)

• sin d-2 (α1,1) sin d-2 (α2,1) • | cos α1,1 -cos α2,1|dα 1:2,1:d-1
We have I A 1 = 0 since integrating over α 1,d-1 and α 2,d-1 create a null factor according to Equation (10). We also have I B 1 = 0 since integrating over α 1,i create a null factor in each term of the sum using Equation (11).

To compute I C 1 , we can integrate all variables different from cos α 1,1 cos α 2,1 in the same way as we have done for computing I 0 and get

I1 = I C 1 = 4π d-1 Γ d-1 2 2 × π 0 π 0 sin d-2 (α1,1
) sin d-2 (α2,1) cos(α1,1) cos(α2,1)| cos α1,1 -cos α2,1|dα1,1dα2,1

Inria

Then using Equation (12) one can finally compute I 1 . Observing that the result of Equation (12) is the same as Equation ( 9) up to a factor -(2d -1) we get

I 1 = - 1 2d -1 I 0

The Upper Bound

Using the values of I 0 and I 1 in Equation ( 5), we get

E [ (V P Xn )] ≤ Γ d 2 + 1 2 2dπ d √ 2 I 0 -1 2 I 1 = Γ d 2 + 1 2 2dπ d I 0 Å 1 + 1 4d -2 ã √ 2 = Γ d 2 + 1 2 2dπ d 4π d-1 Γ d-1 2 2 2 2d ((d -2)!) 2 (2d -2))! Å 1 + 1 4d -2 ã √ 2 = Γ d 2 4 2 4d-5 d π 2 (2d -2)! Å 1 + 1 4d -2 ã √ 2 = … 2d π + 3 4 √ 2dπ + O Ä d -3 2 ä when d → ∞ 4.4 Computation of I 2
We can easily get an upper bound of integral I 2 by noticing that δ 2 ≤ 1 and thus I 2 ≤ I 0 , but this yields an unsatisfactory lower bound for the length of the shortest path. So we compute I 2 exactly.

The integral I 2 can be split in 6 terms according to the development of δ 2 :

δ 2 = δ 2 A + δ 2 B + δ 2 C + 2δ A δ B + 2δ A δ C + 2δ B δ C ,
as for the computation of I 1 , because of Equations 10 and 11, the three terms corresponding to δ A δ B , δ A δ C , and δ B δ C yields null integrals. Thus

I 2 = I A 2 + I B 2 + I C 2 ,
where

I A 2 = (S d-1 ) 2 d-2 j=1 sin α 1,j sin α 2,j • cos(α 1,d-1 -α 2,d-1 ) 2 d-2 j=2
sin d-j-1 (α 1,j ) sin d-j-1 (α 2,j )

• sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1 = (S d-1 ) 2 d-2 i=2 sin d-i+1 (α 1,i ) sin d-i+1 (α 2,i ) cos 2 (α 1,d-1 -α 2,d-1 ) • sin d (α 1,1 ) sin d (α 2,1 )| cos(α 1,1 ) -cos(α 2,1 )| dα 1:2,1:d-1 = d-2 i=2 π 0 sin d-i+1 (α 1,i )dα 1,i 2 • 2π 0 2π 0 cos 2 (α 1,d-1 -α 2,d-1 )dα 2,d-1 dα 1,d-1 • π 0 π 0 sin d (α 1,1 ) sin d (α 2,1 )| cos(α 1,1 ) -cos(α 2,1 )|dα 2,1 dα 1,1 = d-2 i=2 Γ d-i 2 + 1 Γ d-i+1 2 + 1 2 π d-3 E q. (8 ) • 2π 2 E q. (1 3) • 2 2d+4 d! 2 (2d + 2)! E q. (9 ) = 2π d-1 Γ d+1 2 2
te le sc op in g

• 2 2d+4 (d!) 2 (2d + 2)! = 2 2d+5 π d-1 (d!) 2 Γ d+1 I B 2 . = (S d-1 ) 2 d-2 i=2 i-1 j=1 sin α 1,j sin α 2,j cos α 1,i cos α 2,i 2 d-2 j=2
sin d-j-1 (α 1,j ) sin d-j-1 (α 2,j )

• sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1 = (S d-1 ) 2 d-2 i=2 sin d-i-1 (α 1,i ) sin d-i-1 (α 2,i ) × d-2 i=2 i-1 j=2 sin 2 (α 1,j ) sin 2 (α 2,j ) cos 2 (α 1,i ) cos 2 (α 2,i ) × sin d (α 1,1 ) sin d (α 2,1 )| cos(α 1,1 ) -cos(α 2,1 )| dα 1:2,1:d-1
since all terms where cos α 1,i is not squared have null integral by Equation (11)

= d-2 i=2 i-1 j=2 π 0 sin d-j+1 (α 1,j )dα 1,j 2 Å π 0 sin d-i-1 (α 1,i ) cos 2 (α 1,i )dα 1,i ã 2 d-2 j=i+1 π 0 sin d-j-1 (α 1,j )dα 1,j 2 × 2π 0 2π 0 dα 2,d-1 dα 1,d-1 × π 0 π 0 sin d (α 1,1 ) sin d (α 2,1 )| cos(α 1,1 ) -cos(α 2,1 )|dα 2,1 dα 1,1 = á d-2 i=2 i-1 j=2 Γ d-j+2 2 Γ d-j+3 2 2 E q. (8 ) Ç 1 d -i + 1 • Γ d-i 2 Γ d-i+1 2 å 2 E q. ( 1 4
)

d-2 j=i+1 Γ d-j 2 Γ d-j+1 2 2 E q. (8 ) ë π d-3 × 4π 2 × 2 2d+4 d! 2 (2d + 2)! E q. (1 2) = 1 Γ d+1 2 2 d-1 i=3 Ç Γ 1 + i 2 iΓ i 2 å 2
te le sc op in g

× 4π d-1 2 2d+4 d! 2 (2d + 2)! = 1 Γ d+1 2 2 d-1 i=3 1 2 2 Γ pr op er ti es × 4π d-1 2 2d+4 d! 2 (2d + 2)! = Ç (d -3)/4 Γ d+1 2 2 å × 4π d-1 2 2d+4 d! 2 (2d + 2)! = π d-1 (d -3)2 2d+4 (d!) 2 Γ d+1 I C 2 = (S d-1 ) 2 (cos α 1,1 cos α 2,1 ) 2 d-2 i=2 sin d-i-1 (α 1,j ) sin d-i-1 (α 2,i ) • sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) • | cos α 1,1 -cos α 2,1 |dα 1:2,1:d-1 = d-2 i=2 π 0 sin d-i-1 α 1,i dα 1,i 2 × 2π 0 2π 0 dα 2,d-1 dα 1,d-1 × π 0 π 0 sin d-2 (α 1,1 ) sin d-2 (α 2,1 ) cos 2 (α 1,1 ) cos 2 (α 2,1 )| cos(α 1,1 ) -cos(α 2,1 )|dα 2,1 dα 1,1 = d-2 i=2 Γ d-i 2 Γ d-i+1 2 2 π d-3 E q. ( 8
) × 4π 2 × 2 2d+2 (7d -1)d! 2 (d -1) 2 d(2d + 2)! E q. (1 5) = 4π d-1 Γ d+1 2 2
te le sc op in g

• 2 2d+4 (d!) 2 (2d + 2)! • 7d -1 16d = 2 2d+2 π d-1 (d!) (7d -1) d Γ d+1 2 2 (2d + 2)! .
Summing the three terms together and simplifying gives:

I 2 = 2 2d+5 π d-1 (d!) 2 Γ d+1 2 2 (2d + 2)! + π d-1 (d -3)2 2d+4 (d!) 2 Γ d+1 2 2 (2d + 2)! + 2 2d+2 π d-1 (d!) 2 (7d -1) d Γ d+1 2 2 (2d + 2)! = 2 2d+2 π d-1 ((d -2)!) 2 Γ d-1 2 2 (2d -2))! • d 2 (d -1) 2 d-1 2 2 (2d + 2)(2d + 1)2d(2d -1) 8 + 4(d -3) + 7d -1 d = 2 2d+2 π d-1 ((d -2)!) 2 Γ d-1 2 2 (2d -2))! • 4d 2 (2d + 2)(2d + 1)2d(2d -1) • (d + 1)(4d -1) d = I 0 • 4d -1 4d 2 -1

The Lower Bound

Using the values of I 0 , I 1 , and I 2 in Equation ( 5), we get

E [ (V P Xn )] ≥ Γ d 2 + 1 2 2dπ d √ 2 I 0 -1 2 I 1 -1 2 I 2 = Γ d 2 + 1 2 2dπ d I 0 Å 1 + 1 4d -2 - 4d -1 4d 2 -1 ã √ 2 = Γ d 2 4 2 4d-5 d π 2 (2d -2)! Å 1 - d -1 4d 2 -1 ã √ 2 = … 2d π - 1 4 √ 2dπ + O Ä d -3 2 ä when d → ∞ Inria d k

2.3635

Table 2: Lower and upper bounds for the expected length of the Voronoi path using Taylor expansion of order k and exact values from numerical integration.

Small Dimensions

For d small and similarly to what we have done in dimension 3, we can use symbolic computation to get better bounds using higher order Taylor expansions, and to get idea of the exact value using numerical integration. Results are in Table 2, Maple code is available with this report.

A Useful Integrals

Proofs are below.

∞ 0 e -cr d r 2d-1 dr = 1 d • c 2 . ( 7
) π 0 sin d α dα = Γ d+1 2 Γ 1 + d 2 • √ π. ( 8
) π 0 π 0 sin d α 1 sin d α 2 | cos α 1 -cos α 2 |dα 1 dα 2 = 2 2d+4 d! 2 (2d + 2)! . ( 9
) 2π 0 2π 0 cos (α 1 -α 2 )dα 2 dα 1 = 0. ( 10
) π 0 sin d α cos αdα = 0. ( 11
) π 0 π 0 sin d α 1 sin d α 2 | cos α 1 -cos α 2 | cos α 1 cos α 2 dα 1 dα 2 = - 2 2d+4 d! 2 (2d + 3)! . ( 12
) 2π 0 2π 0 cos 2 (α 1 -α 2 )dα 2 dα 1 = 2π 2 . ( 13
) π 0 sin d (α) cos 2 (α)dα = Å 1 d + 2 ã Γ d+1 2 Γ 1 + d 2 • √ π. ( 14
) π 0 π 0 sin d α 1 sin d α 2 cos 2 α 1 cos 2 α 2 | cos α 1 -cos α 2 |dα 1 dα 2 (15) = 2 2d+6 (7d + 13)(d + 2)! 2 (d + 1) 2 (d + 2)(2d + 6)! . b a sin d α dα =              - 1≤i≤d i odd 1 i Å Γ( d+1 2 )Γ(1 + i 2 ) Γ(1 + d 2 )Γ( 1+i 2 ) ã sin i-1 (α) cos(α) b a , if d is odd, - 2≤i≤d i even 1 i Å Γ( d+1 2 )Γ(1 + i 2 ) Γ(1 + d 2 )Γ( 1+i 2 ) ã sin i-1 (α) cos(α) b a + Γ( d+1 2 ) Γ(1 + d 2 ) α √ π b a , if d is even. ( 16
) π 0 α sin d α cos α dα = - Å 1 d + 1 ã Γ d 2 + 1 Γ 1 + d+1 2 • √ π. (17)

Proofs

Proof of Equation [START_REF] Devillers | Walking in a Triangulation[END_REF]. Classic.

Proof of Equation (16). Integration by parts gives

b a sin d α dα = -cos(α) sin d-1 (α) b a + b a cos 2 α(d -1) sin d-2 (α) dα.
Using the identity cos 2 (α) = 1 -sin 2 (α) and simplifying gives b a

sin d α dα = - Å 1 d ã cos(α) sin d-1 (α) b a + Å d -1 d ã b a sin d-2 (α) dα.
1 i Ç Γ( d+1 2 )/ Γ(1 + d 2 ) √ π Γ( i+1 2 )/ Γ(1 + i 2 ) √ π å sin i-1 (α) cos(α) b a + Γ( d+1 2 ) Γ(1 + d 2 ) α √ π b a .
Analogously, for d odd, we use

k k=1 2i 2i + 1 = Γ(d + 1) √ π Γ(d + 3 2 ) , to get by induction b a sin d α dα = - 1≤i≤d i odd 1 i Ç Γ( d+1 2 ) √ π /Γ(1 + d 2 ) Γ( i+1 2 ) √ π /Γ(1 + i 2 ) å sin i-1 (α) cos(α) b a .
Grouping the two expressions above completes the proof.

Inria

Proof of Equation [START_REF] Devroye | Expected time analysis for Delaunay point location[END_REF]. It is a direct corollary of Equation 16.

Proof of Equation (17). Integration by parts gives

π 0 α sin d α cos α dα = 1 d + 1 α sin d+1 (α) π 0 - 1 d + 1 π 0 sin d+1 α dα = - 1 d + 1 π 0 sin d+1 α dα = - Å 1 d + 1 ã Γ d 2 + 1 Γ 1 + d+1 2 • √ π E q. ( 8
)
.

Proof of Equation [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]. First, we use the symmetry with respect to the line α 1 = α 2 to split the integral in two equal parts without absolute values.

integral = 2 α 1 ,α 2 ∈[0,π] 2 α 1 <α 2 sin d α 1 sin d α 2 (cos α 1 -cos α 2 )dα 1 dα 2 = 2 Ç π 0 sin d α 1 cos α 1 π α 1 sin d α 2 dα 2 dα 1 - π 0 sin d α 2 cos α 2 α 2 0 sin d α 1 dα 1 dα 2 å = 2 π 0 sin d x cos x Å π x sin d ydy - x 0 sin d ydy ã . Let f d (x) = Ä π x sin d ydy - x 0 sin d ydy ä
. Then, by using Equation 16, we get

f d (x) =                2 1≤i≤d i odd 1 i Ç Γ( d+1 2 )Γ(1 + i 2 ) Γ(1 + d 2 )Γ( 1+i 2 ) å sin i-1 x cos x, if d is odd, 2 2≤i≤d i even 1 i Ç Γ( d+1 2 )Γ(1 + i 2 ) Γ(1 + d 2 )Γ( 1+i 2 ) å sin i-1 x cos x + Γ( d+1 2 ) Γ(1 + d 2 ) √ π Å 1 - 2x π ã , if d is even.
First, suppose d is odd. Then, replacing the equation above in the simplified integral gives

2 π 0 sin d x cos x Ö 2 1≤i≤d i odd Γ d+1 2 Γ(1 + i/2) Γ(1 + d/2)Γ i+1 2 1 i sin i-1 x cos x è dx,
which simplifies to

4 1≤i≤d i odd Γ d+1 2 Γ(1 + i/2) Γ(1 + d/2)Γ i+1 2 1 i π 0 sin d+i-1 x cos 2 xdx.
Using Equation 14 gives

4 1≤i≤d i odd 1 i Γ d+1 2 Γ(1 + i/2) Γ(1 + d/2)Γ i+1 2 Å 1 d + i + 1 ã Γ d+i 2 Γ d+i+1 2 √ π.
Taking d = 2k -1, i = 2j -1, replacing in the equation above and summing, gives

4 (4k + 1) √ πΓ(k)Γ(2k) Γ(2k + 3 2 )Γ(1 + k) .
Finally, replacing k by d+1 2 and simplifying gives

2 √ π Γ(d + 1) d+1 2 Γ(d + 1 + 1 2 )
. Now, suppose d is even. Then, again, replacing f d in the simplified integral gives

2 π 0 sin d x cos x         2 2≤i≤d i even Γ d+1 2 Γ(1 + i/2) Γ(1 + d/2)Γ i+1 2 1 i sin i-1 x cos x A + Γ( d+1 2 ) Γ(1 + d 2 ) √ π 1 - 2x π B         dx,
We handle the A term first, which is pretty similar to the d odd case. Again, for the A term, expanding the sum and using Equation 14gives

4 2≤i≤d i even 1 i Γ d+1 2 Γ(1 + i/2) Γ(1 + d/2)Γ i+1 2 Å 1 d + i + 1 ã Γ d+i 2 Γ d+i+1 2 √ π.
Taking d = 2k, i = 2j, replacing in the equation above and summing, gives

4 Ç √ πΓ(2k + 1) (2k + 1)Γ(2k + 3 2 ) - 2 (2k + 1) 2 å ,
Finally, replacing k by d/2 and simplifying gives

2 √ π Γ(d + 1) d+1 2 Γ(d + 1 + 1 2 ) -2 1 d+1 2
2 .

Now, we handle the B term, which is

2 π 0 sin d x cos x Ç Γ d+1 2 Γ(1 + d/2) å √ π Å 1 - 2x π ã dx.
Expanding it, gives

2 π 0 sin d x cos x Ç Γ d+1 2 Γ(1 + d/2) å √ πdx -2 • 2 π • π 0 x sin d x cos x Ç Γ d+1 2 Γ(1 + d/2) å √ πdx.
The left term is null because of Equation 11. And plugging Equation 17 in the right term above gives

2 • Γ( d+1 2 ) Γ(1 + d/2) √ π • Å 1 d + 1 ã • Γ(1 + d/2) Γ(1 + d+1 2 ) • √ π • 2 π = 2 1 d+1 2 2 .
Summing A term with B term gives

2 √ π Γ(d + 1) d+1 2 Γ(d + 1 + 1 2 )
.

Inria

Therefore, the expression above holds for both odd and even d. Moreover, it can be simplified to the following nicer expression 2 2d+4 (d!) 2 (2d + 2)! .

Proof of Equation (10). Symmetry on the circle gives 0.

Proof of Equation (11). Symmetry with respect to π/2 gives 0.

Proof of Equation (12). As in the proof of Equation 9, we simplify this integral to get rid of the absolute value and isolate each of its variable.

integral = 2 α 1 ,α 2 ∈[0,π] 2 α 1 <α 2 sin d α 1 sin d α 2 (cos α 1 -cos α 2 ) cos α 1 cos α 2 dα 1 dα 2 = 2 Ç - π 0 sin d α 1 cos α 1 π α 1 sin d α 2 cos 2 α 2 dα 2 dα 1 + π 0 sin d α 2 cos α 2 α 2 0 sin d α 1 cos 2 α 1 dα 1 dα 2 å = 2 π 0 sin d x cos x Å - π x sin d y cos 2 ydy + x 0 sin d y cos 2 ydy ã . = -2 π 0 sin d x cos xf d (x)dx + 2 π 0 sin d x cos xf d+2 (x)dx
where f d is defined in the proof of Equation 9. The first integral is exactly the opposite as the one in Equation 9, the second integral is quite similar, for d odd it gives: Substracting the two integrals yields the claimed result. The case with d even can be solved similarly.

Proof of Equation (13). Simple integration.

Proof of Equation (14). Use cos 2 x = 1 -sin 2 x and Equation 8, then simplify.

Proof of Equation (15). As in the proof of Equations 9 and 12, we simplify this integral to get rid of the absolute value and isolate each of its variable. where f d is defined in the proof of Equation 9. The first and fourth integrals are the same as the one in Equation 9(with d replaced by d + 2 in the fourth integral), the second integral is the same as the one in Equation 12, the third integral is quite similar, for d odd it gives: Adding the four integrals yields the claimed result. The case with d even can be solved similarly.

  d, we get that, for d even, b a sin d α dα = -2≤i≤d i even

b a sin 1 α

 1 dα = -cos(α)

integral = 2 α 1 ,α 2 ∈[0,π] 2 α 1 <α 2 sin d α 1 sin d α 2 (cos α 1 -cos α 2 ) cos 2 α 1 cos 2 α 2 dα 1 dα 2 = 2 Ç π 0 sin d α 1 cos 3 α 1 π α 1 sin d α 2 cos 2 α 2 dα 2 dα 1 - π 0 sin d α 2 cos 3 α 2 α 2 0sin d α 1 cos 2 α 1 dα 1 dα 2 å = 2 π 0 sin d x cos x( 1 -sin 2 x) Å π x sin d y( 1 - 0 sin d y( 1 - 2 π 0 sin

 21222212220102220110120 sin 2 y)dy -x sin 2 y)dy ã . cos xf d (x)dx -2 π 0 sin d x cos xf d+2 (x)dx -d+2 x cos xf d (x)dx + 2 π 0 sin d+2 x cos xf d+2 (x)dx

Table 1 :

 1 Figure 1: The Voronoi path. Lower and upper bounds for the expected length of the Voronoi path.

	RR n°8947

Useful integrals on exponential or trigonometric functions are given in appendix.RR n°8947

RR n°8947

(2d + 2)! , Inria

(2d + 2)! , RR n°8947

rue du Jardin Botanique CS20101 54603 Villers-lès-Nancy Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Sylvain Delchini

Abdeltif Lahfid

Brice Lacroix

Alexis Plunder

Tectonic-metamorphic evolution of the Jebilet massif (Morocco) in the context of the Variscan orogeny

come

The Jebilet massif belongs to the western Moroccan Meseta characterized by nearly complete Paleozoic sedimentary sequence folded and metamorphosed under low-grade greenschist-to amphibolite-facies during the Variscan orogeny, and intruded by widespread syn-to late-orogenic Carboniferous granitoids (Michard et al., 2008, 2010 andreferences therein). It has been previously proposed that Jebilet massif first underwent a regional, greenschist facies metamorphic event (D1 phase), followed by a high-T regional and contact metamorphism that reached the hornfels/amphibolite facies conditions (D2 and D2/D3 phases). However our recent observations of staurolite likely belonging to D1 assemblage obscure the previously proposed metamorphic conditions. In order to refine the metamorphic model we propose to revise the entire metamorphic conditions of the Jebilet massif.

To address this issue, detailed structural, mineralogical, thermobarometric and Raman Spectroscopy on Carbonaceous Material (RSCM) methods have been employed to study the tectono-metamorphic evolution of the Jebilet. The results obtained for the metapelitic rocks that underwent D2/D3 higher metamorphism grades (hornfels/amphibolite facies), show four dominant mineral assemblages: (1) Chlorite-Biotite, (2) Cordierite-Biotite, (3) Andalusite-Garnet-Biotite, and (4) Andalusite-Cordierite-Biotite. The corresponding RSCM temperatures vary between 474 ± 50 °C and 628 ± 50 °C. The computed pseudo-sections for samples from the hornfels/amphibolite facies confirm the peak temperatures measured by the RSCM method. Our structural and mineralogical results support the occurrence of the Garnet-Staurolite assemblage during the D1 regional metamorphism event. These results bring new constrains on the evolution of the Jebilet massif during the Variscan orogeny of the western Moroccan Meseta and strongly suggest that the D1 phase reached higher conditions than those initially proposed, refining the metamorphic history of the Jebilet massif. As the Jebilet massif contains many ore deposit, this work will contribute to determine the main thermal event responsible of concentration the mineralization.