Zakia Mediouni

Guillaume Potherat

Xavier Barrere

Alain Debure

Alexis Descatha

Renal failure

Keywords: Occupational disease, Kidney disease (Failure), (Organic)Solvent, (Etiologic diagnosis)

Introduction

The etiological work-up of chronic kidney disease with renal failure often includes a search for toxic causes, including environmental or occupational causes . However, although many investigations are now frequently performed in patients as a result of recent progress, [1 ] the individual and collective benefit-risk balance of the aetiological work-up must be carefully considered.

The advantages and disadvantages of complementary investigations and subsequent management will be discussed in the light of a case report concerning patient consulting an occupational health unit.

Results/Case report

Medical data

A 47-year-old man was referred by his occupational physician to our occupational medicine unit for an advice concerning occupational management following the discovery of advanced chronic renal failure (declared as occupational disease, fitness for work).

Laboratory tests performed in a context of deteriorating asthenia and recently discovered severe hypertension showed elevated serum creatinine of 285 mol.L and 248 mol.L when repeated 6 months later. Creatinine clearance was 26 ml.min and 28 ml.min 6 months

1 - 1 - 1 - 1 -
later. Investigation of this renal failure did not reveal any urinary tract obstruction, nephrotic syndrome, or electrolyte abnormalities.

Haemoglobin and serum calcium were at the lower limit of normal. The aetiological work-up did not reveal any significant cardiovascular disease. Ultrasound showed no obstruction and a small, regular right kidney with no signs of polycystic disease, no vascular abnormalities and no thrombosis on Doppler studies. He had no particular personal or family medical history, with no history of infection or any particular treatment. The haematological, immunological and metabolic assessment was normal. In view of the patient s occupation, an ' occupational cause was therefore considered.

Occupational data

The patient had worked as a painting assistant in several companies since the age of 18 in Algeria and then in France for the last 7 years. In his current job, for the last 6 years, the patient prepares cars painting 5 days a week from 8:30 a.m. to 12:00 a.m. and from 1:30 p.m. to 5:00 p.m. Paint preparation is performed in a workshop without windows or mechanical ventilation. The two doors of the workshop open onto a courtyard, but are kept closed during working hours. The patient s personal protection equipment comprises gloves, ' overalls and FFP1 masks that he does not use. The kitchen where the patient eats his meals is situated in the same workshop. Safety data forms for the paints prepared by the patient indicate that he was concomitantly exposed to various organic solvents while occupying his current job (especially unsubstituted aliphatic hydrocarbons such as naphtha and higher benzene analogues) in proportions ranging from 0.1 to 80 . No heavy metals were found in the painting.

%

In practice, evaluation of the patient s exposure to organic solvents can be reasonably considered to be greater than 20 hours a week ' for 6 years, which corresponds to a high level of exposure to organic solvents. The patient s exposure to organic solvents during his ' previous jobs could not be determined, but was very likely.

Review of the literature

A review of the literature confirms that some of these solvents have been incriminated in the development of renal failure with varying levels of evidence , as some studies appear to link organic solvents with certain forms of glomerulonephritis (immunological

[2 -6 ]
mechanism), such as membranoproliferative glomerulonephritis and IgA glomerulonephritis . The odds ratios linking exposure to [6 -9 ] organic solvents with glomerulonephritis reported in the literature range between 2.8 to 8.9 with a probable dose-effect relationship between the degree of exposure and severity of the disease , , .

[3 4 6 ] A number of arguments also indicate toxicity of organic solvents on proximal tubular cells causing necrosis via oxidative stress (experimental studies and studies on glue sniffers , . Several studies provide concordant results concerning the aggravating role of [1 10 ] solvent-induced renal failure , , . For example, a large cohort study demonstrated significant deterioration of membrane [1 10 11 ] nephropathy towards end-stage renal disease due to exposure to organic solvents, especially toluene/xylene and gasoline, fuel, gas-oil .

[12 ]

Discussion

Various investigations could be proposed in this case to elucidate the associations between exposure to organic solvents and renal failure.

General investigations have been ever made to exclude non-occupational causes of renal failure. A possible lead exposure before the 1980 s could have been hypothesised, especially by releasing of lead ions from bones. Nevertheless, even though we could not exclude this ' aetiology, the haematological investigations were normal.

Kidney needle biopsy is theoretically essential to identify the site of the kidney lesion and was discussed in this case to eliminate a lesion with no known occupational cause. However, in view of the advanced degree of renal failure of our patient, biopsy would very probably have shown only nonspecific and inconclusive abnormalities attributed to the invasive nature of the technique.

Salivary gland histology was also envisaged to exclude amyloidosis, but the absence of renal hypertrophy made such a diagnosis unlikely and, in any case, no specific treatment is available for this disease.

Confirmation and quantification of occupational exposure were also considered in order to corroborate solvent exposure (metrological samples), but such investigations only provide confirmation of exposure when they are positive and do not provide any information about the link between exposure and the patient s disease (especially as the exposure responsible for the disease may date back many years). ' Furthermore, negative results of this type of metrological study do not exclude possible exposure (solvents not assayed, solvents not handled on the day of sampling, etc.). Consequently, in view of the possible direct link between exposure to organic solvents and renal failure, and the possibly essential nature of the disease (no aetiology demonstrated), the patient was given a medical certificate in support of an occupational disease. The advantage for the patient of exclusion of the role of exposure would be to allow him to continue his current occupation. Nevertheless, in view of the products handled and the available data concerning this type of occupation, the probability of exposure is relatively high with a low negative predictive value for metrology. It therefore appears reasonable to suspend the patient s job ' aptitude and request a change of job to slow the progression to dialysis or renal transplantation.

The possible collective benefit of further investigations needed however to be discussed after consulting several specialists, as precise quantification of the risk in the workshop by environmental (metrological) and laboratory measures of exposure or its effects, for example by screening other employees, could provide a benefit in terms of collective prevention. This would allow primary prevention of unexposed employees, but also earlier detection of any renal abnormalities (secondary prevention). The disadvantage of this type of approach is the economic risk for small financially fragile companies. Apart from providing better individual medical, social and financial cover, declaration of an occupational disease can also contribute to improvement of overall working conditions for other employees.

Conclusion

This case report illustrates the importance of limiting certain invasive investigations that only provide a limited benefit to the patient. However, other investigations with a collective benefit can be considered. A multidisciplinary approach (general practitioner, nephrologist, occupational health physician and specialist in toxicology) is therefore necessary to discuss the most appropriate work-up, based on assessment of the individual and collective benefit-risk balance.

References:


MD Alexis Descatha
email: alexis.descatha@rpc.aphp.fr

Unusual delayed reaction after H1N1 vaccine

Thimerosal is a traditional high-prevalence allergen in the general population, which has been used as an antimicrobial agent in contact lens and vaccines [START_REF] Nielsen | Allergic contact sensitization in an unselected Danish population. The Glostrup Allergy Study, Denmark[END_REF]. However, prevalence in the general population is decreasing in many countries, such as reactions to thimerosal, because of the exclusion of thimerosal in vaccines [START_REF] Thyssen | The epidemiology of contact allergy in the general population--prevalence and main findings[END_REF]. We reported a probable case of allergic contact dermatitis three weeks after the H1N1 vaccine, involving thimerosal additive.

A 55 year-old woman, who is a medical doctor, noticed an important skin reaction on her left arm (figure 1). The lesion was characterized by erythema and edema, followed by the appearance of papules and vesicles. The skin reaction was exactly at the same site of injection of the Influenza A H1N1 vaccine (Pandemrix© vaccine, Influenza A H1N1 2009 Monovalent AS03-Adjuvanted Vaccine, GSK), but 22 days after the shot. No history of atopy or documented allergic contact dermatitis was found. Only moderate pruritus occurred, with only one application of topical steroid. She had no other significant medical condition of any kind and no treatment at this time. The patch test of thimerosal (thiomersal, additive of the vaccine) was weakly positive a week after the initial visit ('+', erythema, infiltration and papules). The evolution was favorable with disappearance in three weeks. A complete testing, including most of the component of the vaccine was performed two month later, in order to avoid recurrence and to confirm the association with thiomerosal. The patch test only found a positivity for thimerosal (but important this time '+++'), and was involved in the allergic contact dermatitis presented by our colleague.

There is a debate about thimerosal reactions: some studies raised the possibility that most of the positive reactions to it are irrelevant. However, even though reactions are to it are decreasing, there a consistency about thimerosal sensitivity does exist in contact dermatitis patients [START_REF] Slodownik | Thimerosal--is it really irrelevant?[END_REF]. Immediate generalized reaction to thimerosal from an influenza vaccine have been previously reported [START_REF] Mcmahon | Inactivated influenza vaccine (IIV) in children <2 years of age: examination of selected adverse events reported to the Vaccine Adverse Event Reporting System (VAERS) after thimerosalfree or thimerosal-containing vaccine[END_REF][START_REF] Lee-Wong | A generalized reaction to thimerosal from an influenza vaccine[END_REF][START_REF] Zheng | Thimerosal in influenza vaccine: an immediate hypersensitivity reaction[END_REF], but this was the first reported case of allergic contact dermatitis three weeks after the H1N1 vaccine, probably caused by thimerosal additive.

Allergic contact dermatitis presented by our colleague come from delayed-type hypersensitivity reaction [START_REF] Usatine | Diagnosis and management of contact dermatitis[END_REF], mediated by hapten-specific T cells [START_REF] Saint-Mezard | Allergic contact dermatitis[END_REF]. Even though T-cell mediated reactions usually starting from 2-8 hours up to 2 days after vaccination [START_REF] Fritsche | Vaccine hypersensitivity-update and overview[END_REF], very late reactions have been described for other medications [up one month [START_REF] Rajabian | Severe facial dermatitis as a late complication of aesthetic rhinoplasty; a case report[END_REF]]. Other differential cause of acute dermatitis has been excluded by anamnesis and patch tests. No cause of immunodeficiency has been found. It is however possible, the lesion appeared before our colleague noticed it, but only a few days before.

In conclusion, we report an allergic contact dermatitis involving thimerosal additive. Patient should be aware of this possible problem with this unusual delayed reaction.

Figure 1

 1 Figure 1 Eczema lesion of left shoulder, three days after the beginning of the rash (25 days after the H1N1 vaccine shot).

Acknowledgments

We would like to thank: Dr AM Jonathan (allergologist), Dr MN Crepy (dermatologist), Dr Lillo-Le-Louet (phrarmacovigilance center), Mr H Leroux (MD student) and our colleague who accepted the publication of her case.

Abstract (if necessary)

We reported a probable case of allergic contact dermatitis three weeks after the H1N1 vaccine, involving thimerosal additive. Patients should be aware of this possible problem with this unusual delayed reaction.


J Lévy Véhel

Wavelet packet based digital watermarking

We present a method for digital image watermarking based on the modification of certain subsets of the wavelet packet decomposition. These subsets are determined both from a secret key and an image dependent procedure that chooses a best basis from an energy criterion. The mark is set by imposing a parity constraint at each level of the decomposition. We elaborate on the choice of some of the parameters of the model, showing how they can be tuned so as to obtain good resistance to attacks. Examples are displayed to assess the validity of our approach.

Introduction

Digital image watermarking has attracted a lot of interest in recent years, due in particular to the development of Internet and the World Wide Web. The aim is to protect ownership by including in the image a copyright information. This information, or mark, has to be set in such a way that it is invisible: indeed, it must not alter the viewing content and, in addition, it should not be easy to remove. Furthermore, the mark must be resistant to attacks directed at erasing it. Such attacks are of two kinds, whether one uses cryptographic or image processing methods. We shall be concerned in this paper only with the second type of attacks.

A number of methods have been proposed to insert robust and invisible watermarks. Some operate directly in pixel space [START_REF] Nikolaidis | Copyright protection of images using robust digital signatures[END_REF], other in a transform domain, such as Fourier [START_REF] Bordland | Phase watermarking of digital images[END_REF] or DCT [START_REF] Koch | Towards robust and hidden image copyright labeling[END_REF]. We propose here to study a wavelet packet based watermarking procedure. Working in the wavelet domain yields a number of advantages: First, it allows to control in a precise way the location both in space and scale of the mark. Second, wavelet coefficients give a structured way of representing the information: as is well known for instance from studies in image compression, in most cases, only a few coefficients are large, indicating where the information lies in a given image in terms of scale and space. This is useful both for invisibility and robustness concerns.

Using wavelet packets adds another degree of freedom because it allows to select frequency independently of scale. Finally, wavelet based algorithms are fast and allow to reconstruct the image. Section 2 briefly recalls some basic facts about wavelet packets decomposition. Section 3 presents our watermarking algorithm. Some details on the choice of various parameters are discussed in section 4. Finally, section 5 presents numerical experiments.

Recalls on Wavelet Packets

Wavelet packet decomposition (WPD) is a generalization of the dyadic wavelet transform (DWT) where the low pass parts are further analyzed. Figure 1 shows the decomposition tree of an image: the coefficients of each packet are obtained by successive filtering and decimation along lines and columns. As is apparent from figure 1, this representation is redundant, and it is possible to extract a basis by selecting coefficients in an appropriate way. Figure 2 shows three levels of decomposition of an image, along with one of the bases corresponding to a particular tiling of the space/frequency domain. Usually, one defines a "best basis" as a basis that optimizes a certain criterion. A popular choice is to minimize the entropy of the representation [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF].

A word on notation: C p;i;j denotes the packet at resolution p in the frequency region indexed by i; j. Individual values of the coefficients in this packet are denoted C p;i;j (k), where k codes for the spatial translation.

Watermarking based on the parity of certain subsets associated to a best basis

Our method consists in inserting a mark by modifying the best basis associated to an image so that it respects some parity constraints. More precisely, the mark will be a sequence of 0-s and 1-s ; we first compute the best basis using an energy criterion detailed below (the classical entropy criterion is not fitted to our needs because it is not robust to compression). We then extract certain subsets from this basis, using a private information coded in a secret key. The image is modified so that, at each level p of the decomposition, the number of packets in the best basis which are selected by the secret key is odd or even according to whether the p th element in the mark is 0 or 1. This procedure inserts in some sense a "virtual" mark, since it is set on the structure of the best basis. We now make the different steps precise (see figure 3). Since the mark is set on the best basis, the basis should be both robust to attacks and contain enough vectors so that sufficiently rich marks may be implemented. For instance, we do not wish to select only packets at the highest or lowest resolution levels because they would not fulfill these constraints. The following criterion for best basis selection helps to attain our objectives. In short, we shall elect a packet to be in the basis B if it has sufficient energy and if its offsprings do not share this property. Formally: C p;i;j 2 B iff C p;i;j 2 F and C p+1;2i;2j 6 2 F , where the set F is defined in a recursive way: F = fC 0;0;0 g fC d;x;y = 8( 1 ; 2 ) 2 f0; 1g 2 ; X k jC d;x+ 1r1;y+ 2r2 (k)j 2 ; r 1 = 1 2 (x 2]); r 2 = 1 2 (y 2]); C d To each sub-base B i , we then associate a scalar e i : e i is 0 if the number of elements in B i is even and is 1 otherwise (see figure 4). We first modify SB according to the watermark that identifies the owner of the image. This is modeled by an operator : P(B) W 7 ! P(B) (SB; W) 7 ! SB = B 1 ; B 2 ; :::; B m where the set of watermarks W is simply f0; 1g m . SB is computed in the following way:

SB V 1 0 K V* B B* W= 0
1. if w i = e i , then the elements of SB at are equal to those of SB 2. otherwise, we suppress the element of SB that has highest energy and add at scale below its four offsprings to fill B. SB is filled by the offsprings selected by E.

3. the same procedure is performed iteratively for all the subset.

Once SB is completely defined, we modify the image to get the corresponding best basis B . To this end, the coefficients of the offsprings of the elements suppressed from SB are increased in absolute value so that the energy in their packet becomes larger than . This procedure is iterated for each subset starting from the lowest resolution (see figure 4). Finally, to make the process more resistant to attacks, we slightly modify the whole image in the following way : we multiply those coefficients which are "close" to the thresthold by a value that move them away from the thresthold.

Finally, the decoding process simply follows the same steps 1,2,3 as above, and then replaces step 4 by the computation of the parity of the number of vectors at each level of the m secret sub-bases. The presence of the mark is then classicaly assessed through computation of the Hamming distance between the observed parity vector and the theoretical one. Note that the decoding process does not require the knowledge of the original image.

Fine tuning of the parameters

Our method requires that a number of choices be made regarding the various parameters. In general, it is a delicate matter to derive values which are universally optimal. This would imply the knowledge of all possible images and all possible attacks, which is clearly impossible. A way out is to restrict the considered range of images and/or attacks. For some applications, for instance photographic reproduction of artworks, it is quite natural to focus on a given type of images. In this paper, we will rather specify a set of possible attacks and keep all freedom on the kind of marked images.

We have chosen to use as a model for the set of all possible attacks the Stirmark software, which has served as a benchmark for testing the robustness of watermarks by various authors [START_REF] Kutter | A fair benchmark for image watermarking systems[END_REF]. Our aim in this section is thus to try and set some of the parameters of our method so as to maximize its resistance to the attacks implemented in Stirmark. More precisely, we will consider the choice of the energy threshold in the best basis selection. To simplify notation, we denote from now on c i the energy of the packet at "location" i, where i an index that subsumes the information of scale and frequency. Obvious requirements concerning the energy threshold are the following : on the one hand, should be set small enough so that sufficiently many packets have energy larger than , and that a best basis may be found. On the other hand, a large is desirable, because selection of low energy packets should be forbidden (they are not robust to attacks such as compression or filtering). If we consider specifically the problem of robustness to attacks, something more precise may be said. Let I be the set of the indices of all wavelet packets coefficients, J be the subset of all indices in the marked image such that the corresponding coefficients are eligible to be in the best basis, i.e. i 2 J if c i > , and K be the complement of J in I. Since the secret key is chosen independently of the best basis, an ideal choice for (as far as robustness is concerned) would be one that makes both the sets J and K stable under the action of any attack. In other words, an ideal situation would be one where c i > ) d i > and c i < ) d i < , where c i denotes the energy of a packet in the marked image and d i the corresponding energy in the image after an attack. Indeed, such a choice would insure that the same coefficients are considered whether the image has been attacked or not, yielding in turn the same parity at each level, and making the mark insensitive to attacks.

While it is not possible to find such a in general, we may try to be as close as possible to this ideal case. In this view, a reasonable approach is to use an asymmetric procedure and to look for two reals ( ; 0 ) so that the number of coefficients for which c i and d i 0 have opposite signs is minimum in some sense. Denoting e and e the maximum and minimum of the energy of all packets, we consider the following stochastic modeling of our problem.

For a fixed image D, we consider the set C of all its wavelet packets energy coefficients c j . The set of elementary events is composed of "all identified attacks", for instance all attacks implemented in Stirmark. Our probability triplet is thus ( ; B; P), where B is the natural algebra associated with and P is the uniform probability. For each index i 2 I, the wavelet packet energy coefficient d i is a random variable whose probability distribution F i is determined by the result of the action of all possible attacks on the particular coefficient c i . We consider the following random function :

S( ; 0 ) = X i2I (c i )sgn(d i 0 )
where sgn(x) = x jxj ; sgn(0) = 0, and is a C 1 function that approximates the sgn function, such as the arctangent (this will simplify the optimization step). We seek to maximize B( ; 0 ) = E(S( ; 0 )) w.r.t. ( ; 0 ) 2 = e; e],

where E denotes expectation and the bounds set on ( , 0 ), are a weak version of the requirements described at the beginning of this subsection, namely that enough coefficients should be eligible, but packets with too low energy should not be included. Such restrictions also allow our problem to be well posed, otherwise ( , 0 ) = (0; 0) and ( , 0 ) = (1; 1) would be degenerate solutions independent of the input image. It is easy to check that:

B( ; 0 ) = X i2I (c i )(1 2F i ( 0 )) (1)
If the F i are C 2 , so is B( ; 0 )1 , and its maximum may be found by classical gradient descent once the F i are known. To learn the F i -s, we adopted the following strategy: for each image in a learning set, we computed for each c i the empirical distribution corresponding to all attacks. We checked that, for a given position i, the empirical distributions F I i were roughly similar for all images I. We then approximated the common F i using a shifted Gamma distribution with density of the form: are free parameters to be estimated from empirical distribution. The choice of the Gamma distribution and the form of the shift were made on the basis of the observed empirical distributions and on "semi-heuristic"considerations. In particular, if we assume that one can model the action of attacks as the addition of Gaussian white noise to the wavelet coefficients, then the energies should be distributed as a chisquare, which is a particular case of Gamma distribution. In addition, the fact Gamma distributions are stable by linear combinations is also useful for our purpose. Figure 5 shows examples of empirical distributions of the wavelet packets energy which have been submitted to the 89 attacks implemented in Stirmak. Figure 6 shows a fitted Gamma distribution using a maximum likelihood estimate of a; b; c. Once the parameters and thus the F i are known, it is an easy task to compute ( ; 0 ) that maximize B( ; 0 ). While the above modeling works fine for attacks such as compression or small angle rotations, other ones based on cropping or low-pass filtering are not well controlled by this approach. These attacks can however be taken care of efficiently through adequate restrictions on the secret key. This will be presented elsewhere.

Numerical Experiments

We show an experiment on a boat image of size 256 256 pixels coded on 8 bits. Figure 7 shows the original signal along with a marked image containing a 32 bits watermark with redundancy equal to 71 . The PSNR between the two images is 49,27 dB. Figure 8 shows the watermark detector response to 1000 randomly generated watermarks after an attack consisting of JPEG compression with 30 % quality. This and other experiments made on several images with various attacks show that our method is an efficient and robust watermarking technique.

Figure 1 .Figure 2 .

 12 Figure 1. Quaternary wavelet packet decomposition tree, obtained by successive filtering and decimation along the lines and the columns. h and g denotes respectively the high pass and low pass filters of the DWT.

Figure 3 .

 3 Figure 3. Different steps in the watermarking process.

step 1 :

 1 Wavelet packet decomposition. step 2: Best basis selection.

Figure 4 .

 4 Figure 4. Encoding of the mark on a sub-base. step 4: Encoding of the mark.

  f i (x) = (x a c i ) b 1 e c(x a ci) H(x a c i ) where = c b (b) and H is the Heaviside function. a; b; c

Figure 5 .

 5 Figure 5. Empirical distributions of the attacked wavelet packet energies, for packet C8;203;201 in the Lena (left) and Barbara (right) images.

Figure 6 .

 6 Figure 6. Empirical distributions of the attacked wavelet packet energies and the fitted version.

Figure 7 .

 7 Figure 7. Original (left) and watermarked image (right).

Figure 8 .

 8 Figure 8. Watermark detector response to 1000 randomly generated watermarks after JPEG compression.

  x 2] means x modulo 2. The best basis is thus of maximal depth and contains only packets which have "sufficient" energy: a packet is in the best basis if it belongs to F but its offsprings do not. F is the set of packets which have energy larger than and such that their brothers, father and uncles also have energy larger than . The important question of the optimal choice of the threshold is addressed in the next section. step 3: Extraction of the sub-bases.

	1;E(x=2);E(y=2) 2 F g

where Let B be the set of all admissible bases corresponding to the image. We define an operator E : B 7 ! P(B) B 7 ! SB = B 1 ; B 2 ; :::; B m where B 1 ; B 2 ; :::; B m are disjoint subsets of B on which the mark will be implemented. SB is the secret key associated to the process and we have : SB B. Note that, of course, the extraction operator E should not easily be invertible.

While there are a finite number of attacks in Stirmark, we may still assume that F i has a density by considering that the various parameters (e.g. the compression rate in JPEG) can be tuned in a continuous way.


Marie Doumic-Jauffret
email: marie.doumic-jauffret@inria.fr

Marc Hoffmann
email: marc.hoffmann@ensae.fr

Patricia Reynaud-Bouret
email: patricia.reynaud-bouret@unice.fr

Vincent Rivoirard
email: vin-cent.rivoirard@dauphine.fr

M Doumic

Nonparametric estimation of the division rate of a size-structured population

Keywords: Lepski method, Oracle inequalities, Adaptation, Aggregation-fragmentation equations, Statistical inverse problems, Nonparametric density estimation, Cell-division equation. Mathematical Subject Classification: 35A05, 35B40, 45C05, 45K05, 82D60, 92D25

1 Introduction

Motivation

Structured models have long served as a representative deterministic model used to describe the evolution of biological systems, see for instance [START_REF] Perthame | Transport equations arising in biology[END_REF] or [START_REF] Metz | Formulating models for structured populations[END_REF] and references therein. In their simplest form, structured models describe the temporal evolution of a population structured by a biological parameter such as size, age or any significant trait, by means of an evolution law, which is a mass balance at the macroscopic scale. A paradigmatic example is given by the transport-fragmentation equation in cell division, that reads

             ∂ ∂t n(t, x) + ∂ ∂x g 0 (x)n(t, x) + B(x)n(t, x) = 4B(2x)n(t, 2x
), t ≥ 0, x ≥ 0, gn(t, x = 0) = 0, t > 0, n(t = 0, x) = n 0 (x), x ≥ 0.

(1.1)

The mechanism captured by Equation (1.1) can be described as a mass balance equation (see [START_REF] Banks | Estimation of Cell Proliferation Dynamics Using CFSE Data[END_REF][START_REF] Metz | Formulating models for structured populations[END_REF]): the quantity of cells n(t, x) of size x at time t is fed by a transport term g 0 (x) that accounts for growth by nutrient uptake, and each cell can split into two offsprings of the same size according to a division rate B(x). Supposing g 0 (x) = κg(x), where we suppose a given model for the growth rate g(x) known up to a multiplicative constant κ > 0, and experimental data for n(t, x), the problem we consider here is to recover the division rate B(x) and the constant κ.

In [START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF], Perthame and Zubelli proposed a deterministic method based on the asymptotic behavior of the cell amount n(t, x) : indeed, it is known (see e.g. [START_REF] Michel | General relative entropy inequality: an illustration on growth models[END_REF][START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF]) that under suitable assumptions on g and B, by the use of the general relative entropy principle (see [START_REF] Perthame | Transport equations arising in biology[END_REF]), one has ∞ 0 n(t, x)e -λt -n 0 , φ N (x) φ(x) dx -→ t→∞ 0

(1.2)

where n 0 , φ = n 0 (y)φ(y)dy and φ is the adjoint eigenvector (see [START_REF] Michel | General relative entropy inequality: an illustration on growth models[END_REF]). The density N is the first eigenvector, and (λ, N ) the unique solution of the following eigenvalue problem

     κ ∂ ∂x g(x)N (x) + λN (x) = 4BN (2x) -BN (x), x > 0, B(0)N (0) = 0, N (x)dx = 1, N (x) ≥ 0, λ > 0.
(1.3) Moreover, under some supplementary conditions, this convergence occurs exponentially fast (see [START_REF] Perthame | Exponential decay for the fragmentation or cell-division equation[END_REF]). Hence, in the rest of this article, we work under the following analytical assumptions.

Assumption 1 (Analytical assumptions).

1. For the considered nonnegative functions g and B and for κ > 0, there exists a unique eigenpair (λ, N ) solution of Problem (1.3).

2. This solution satisfies, for all p ≥ 0, x p N (x)dx < ∞ and 0 < g(x)N (x)dx < ∞.

3. The functions N and gN belong to W s+1 with s ≥ 1, and in particular N ∞ < ∞ and (gN ) ′ 2 < ∞. (W s+1 denotes the Sobolev space of regularity s + 1 measured in L 2 -norm.) [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF]. We have g ∈ L ∞ (R + ) with R + = [0, ∞).

Hereafter • 2 and • ∞ denote the usual L 2 and L ∞ norms on R + . Assertions 1 and 2 are true under the assumptions on g and B stated in Theorem 1.1 of [START_REF] Doumic | Eigenelements of a General Aggregation-Fragmentation Model[END_REF], under which we also have N ∈ L ∞ . Assertion 3 is a (presumably reasonable) regularity assumption, necessary to obtain rates of convergence together with the convergence of the numerical scheme. Assertion 4 is restrictive, but mandatory in order to apply our statistical approach.

Thanks to this asymptotic behavior provided by the entropy principle (1.2), instead of requiring time-dependent data n(t, x), which is experimentally less precise and more difficult to obtain, the inverse problem becomes: How to recover (κ, B) from observations on (λ, N ) ? In [START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF][START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF], as generally done in deterministic inverse problems (see [START_REF] Engl | Regularization of Inverse Problems[END_REF]), it was supposed that experimental data were pre-processed into an approximation N ε of N with an a priori estimate of the form N -N ε ≤ ε for a suitable norm • . Then, recovering B from N ε becomes an inverse problem with a certain degree of ill-posedness. From a modelling point of view, this approach suffers from the limitation that knowledge on N is postulated in an abstract and somewhat arbitrary sense, that is not genuinely related to experimental measurements.

The statistical approach

In this paper, we propose to overcome the limitation of the deterministic inverse problems approach by assuming that we have n data, each data being obtained from the measurement of an individual cell picked at random, after the system has evolved for a long time so that the approximation n(t, x) ≈ N (x)e λt is valid. This is actually what happens if one observes cell cultures in laboratory after a few hours, a typical situation for E. Coli cultures for instance, provided, of course, that the underlying aggregation-fragmentation equation is valid.

Each data is viewed as the outcome of a random variable X i , each X i having probability distribution N (x)dx. We thus observe (X 1 , . . . , X n ), with

P(X 1 ∈ dx 1 , . . . , X n ∈ dx n ) = n i=1 N (x i )dx i ,
and where P(•) hereafter denotes probability1 . We assume for simplicity that the random variables X i are defined on a common probability space (Ω, F, P) and that they are stochastically independent. Our aim is to build an estimator of B(x), that is a function x Bn (x, X 1 , . . . , X n ) that approximates the true B(x) with optimal accuracy and nonasymptotic estimates. To that end, consider the operator

(λ, N ) T(λ, N )(x) := κ ∂ ∂x g(x)N (x) + λN (x), x ≥ 0. (1.4)
From representation (1.3), we wish to find B, solution to T(λ, N ) = L(BN ), where

L ϕ (x) := 4ϕ(2x) -ϕ(x), (1.5)
based on statistical knowledge of (λ, N ) only. Suppose that we have preliminary estimators L and N of respectively T(λ, N ) and N, and an approximation L -1 k of L -1 . Then we can reconstruct B in principle by setting formally

B := L -1 k ( L) N .
This leads us to distinguish three steps that we briefly describe here. The whole method is fully detailed in Section 2.

The first and principal step is to find an optimal estimator L for T(λ, N ). To do so, the main part consists in applying twice the Goldenschluger and Lepski's method [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] (GL for short). This method is a new version of the classical Lepski method [START_REF] Lepski | One problem of adaptive estimation in Gaussian white noise[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 1. Upper bounds[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 2. Statistical models without optimal adaptation[END_REF][START_REF] Lepski | On problems of adaptive estimation in white Gaussian noise[END_REF]. Both methods are adaptive to the regularity of the unknown signal and the GL method furthermore provides with an oracle inequality. For the unfamiliar reader, we discuss adaptive properties later on, and explain in details the GL method and the oracle point of view in Section 2.

1. First, we estimate the density N by a kernel method, based on a kernel function K.

We define N = Nĥ where Nh is defined by (2.1) and the bandwidth ĥ is selected automatically by (2.3) from a properly-chosen set H. (see Section 2.1 for more details).

A so-called oracle inequality is obtained in Proposition 2 measuring the quality of estimation of N by N . Notice that this result, which is just a simplified version of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF], is valid for estimating any density, since we have only assumed to observe an n-sample of N, so that this result can be considered per se.

2. Second, we estimate the density derivative (up to g) D = ∂ ∂x (gN ), again by a kernel method with the same kernel K as before, and select an optimal bandwidth h given by Formula (2.6) similarly. This defines an estimator D := Dh where Dh is specified by (2.4), and yields an oracle inequality for D stated in Proposition 3. In the saemway as for N, this result has an interest per se and is not a direct consequence of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF].

From there, it only remains to find estimators of λ and κ. To that end, we make the following a priori (but presumably reasonable) Assumption 2 on the existence of an estimator λn of λ.

Assumption 2 (Assumption on λn ).

There exists some q > 1 such that

ε λ,n = E[| λn -λ| q ] 1/q < ∞, R λ,n = E[ λ2q n ] < ∞.
Indeed, in practical cell culture experiments, one can track n individual cells that have been picked at random through time. By looking at their evolution, it is possible to infer λ in a classical parametric way, via an estimator λn that we shall assume to possess from now on2 . Based on the following simple equality

κ = λρ g (N ) where ρ g (N ) = R + xN (x)dx R + g(x)N (x)dx , (1.6)
obtained by multiplying (1.3) by x and integrating by part, we then define an estimator κn by (2.8). Finally, defining L = κn D + λn N ends this first step. The second step consists in the formal inversion of L and its numerical approximation: For this purpose, we follow the method proposed in [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] and recalled in Section 2.4. To estimate H := BN , we state

Ĥ := L -1 k ( L) (1.7)
where L -1 k is defined by (2.10) on a given interval [0, T ]. A new approximation result between L -1 and L -1 k is given by Proposition 4. The third and final step consists in setting B := Ĥ N , clipping this estimator in order to avoid explosion when N becomes too small, finally obtaining

B(x) := max(min( B(x), √ n), - √ n).
(1.8)

Rates of convergence

Because of the approximated inversion of L on [0, T ], we will have access to error bounds only on [0, T ]. We set f 2 2,T = T 0 f 2 (x)dx. for the L 2 -norm restricted to the interval [0, T ]. If the fundamental (yet technical) statistical result is the oracle inequality for Ĥ stated in Theorem 1 (see Section 2.5), the relevant part with respect to existing works in the non-stochastic setting [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF][START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF] is its consequence in terms of rates of convergence. For presenting them, we need to assume that the kernel K has regularity and vanishing moments properties.

Assumption 3 (Assumptions on K). The kernel K is differentiable with derivative K ′ . Furthermore, K(x)dx = 1 and K 2 and K ′ 2 are finite. Finally, there exists a positive integer m 0 such that K(x)x p dx = 0 for p = 1, . . . , m 0 -1 and I(m 0 ) := |x| m 0 K(x)dx is finite.

Then our proposed estimators satisfy the following properties.

Proposition 1. Under Assumptions 1, 2 and 3, let us assume that R λ,n and √ nǫ λ,n are bounded uniformly in n and specify L -1 k with k = n f. Assume further that the family of bandwidth

H = H = {D -1 : D = D min , ..., D max } depends on n is such that 1 ≤ D min ≤ n 1/(2m 0 +1) and n 1/5 ≤ D max ≤ n 1/2 for all n. Then Ĥ satisfies, for all s ∈ [1; m 0 -1] E Ĥ -H q 2,T = O n -qs 2s+3 , (1.9)
Furthermore, if the kernel K is Lipschitz-regular, if there exists an interval [a, b] in (0, T ) such that

[m, M ] := [ inf x∈[a,b] N (x), sup x∈[a,b] N (x)] ⊂ (0, ∞), Q := sup x∈[a,b] |H(x)| < ∞,
and if ln(n) ≤ D min ≤ n 1/(2m 0 +1) and n 1/5 ≤ D max ≤ (n/ ln(n)) 1/(4+η) for some η > 0, then B satisfies, for all s ∈ [1, m 0 -1],

E ( B -B)1 [a,b] q 2 = O n -qs 2s+3 .
(1.10)

Remarks and comparison to other works

1) Let us establish formal correspondences between the methodology and results when recovering B from (1.3) from the point of view of statistics or PDE analysis. ¡ After renormalization, we obtain the rate n -s/(2s+3) for estimating B, and this corresponds to ill-posed inverse problems of order 1 in nonparametric statistics. We can make a parallel with additive deterministic noise following Nussbaum and Pereverzev [START_REF] Nussbaum | The degrees of ill-posedness in stochastic and deterministic noise models[END_REF] (see also [START_REF] Mair | Statistical inverse estimation in Hilbert scales[END_REF] and the references therein). Suppose we have an approximate knowledge of N and λ up to deterministic errors ζ 1 ∈ L 2 and ζ 2 ∈ R with noise level ε > 0: we observe

N ε = N + εζ 1 , ζ 1 2 ≤ 1, (1.11)
and

λ ε = λ + εζ 2 , |ζ 2 | ≤ 1. From the representation B = L -1 T(N, λ) N ,
where T(N, λ) is defined in (1.4), we have that the recovery of T(N, λ) is ill-posed in the terminology of Wahba [START_REF] Wahba | Practical approximate solutions to linear operator equations when the data are noisy[END_REF] for it involves the computation of the derivative of N . Since L is bounded with an inverse bounded in L 2 and the dependence in λ is continuous, the overall inversion problem is ill-posed of degree a = 1. By classical inverse problem theory for linear cases 3 , this means that if N ∈ W s , the optimal recovery rate in L 2 -error norm should be ε s/(s+a) = ε s/(s+1) (see also the work of Doumic, Perthame and collaborators [START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF][START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF]).

Suppose now that we replace the deterministic noise ζ 1 by a random Gaussian white noise: we observe

N ε = N + εB (1.12)
where B is a Gaussian white noise, i.e. a random distribution in W -1/2 that operates on test functions ϕ ∈ L 2 and such that B(ϕ) is a centered Gaussian variable with variance ϕ 2 2 . Model (1.12) serves as a representative toy model for most stochastic error models such as density estimation or signal recovery in the presence of noise. Let us formally introduce the α-fold integration operator I α and the derivation operator ∂. We can rewrite (1.12) as N ε = I 1 (∂N ) + εB and applying I 1/2 to both side, we (still formally) equivalently observe

Z ε := I 1/2 N ε = I 3/2 (∂N ) + εI 1/2 B.
We are back to a deterministic setting, since in this representation, we have that the noise εI 1/2 B is in L 2 . In order to recover ∂N from Z ε , we have to invert the operator I 3/2 , which has degree of ill-posedness 3/2. We thus obtain the rate

ε s/(s+3/2) = ε 2s/(2s+3) = n -s/(2s+3)
for the calibration ε = n -1/2 dictated by (1.12) when we compare our statistical model with the deterministic perturbation (see for instance [START_REF] Nussbaum | Asymptotic equivalence of density estimation and white noise[END_REF] for establishing formally the correspondence ε = n -1/2 is a general setting). This is exactly the rate we find in Proposition 1: the deterministic error model and the statistical error model coincide to that extent4 .

2) The estimators Ĥ and B do not need the exact knowledge of s as an input to recover this optimal rate of convergence. We just need to know an upper bound m 0 -1 to choose the regularity of the kernel K. This capacity to obtain the optimal rate without knowing the precise regularity is known in statistics as adaptivity in the minimax sense (see [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF] for instance for more details). It is close in spirit to what the discrepancy principle can do in deterministic inverse problems [START_REF] Engl | Regularization of Inverse Problems[END_REF]. However, in the deterministic framework, one needs to know the level of noise ε, which is not realistic in practice. In our statistical framework, this level of noise is linked to the size sample n through the correspondence ε = n -1/2 .

3) Finally, note that the rate is polynomial and no extra-logarithmic terms appear, as it is often the case when adaptive estimation is considered (see [START_REF] Lepski | One problem of adaptive estimation in Gaussian white noise[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 1. Upper bounds[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 2. Statistical models without optimal adaptation[END_REF][START_REF] Lepski | On problems of adaptive estimation in white Gaussian noise[END_REF]).

The next section explains in more details the GL approach and presents our estimators to a full extent, including the fundamental oracle inequalities. It also elaborates on the methodology related to oracle inequality. The main advantage of oracle inequalities is that they hold nonasymptotically (in n) and that they guarantee an optimal choice of bandwidth with respect to the selected risk. Section 3 is devoted to numerical simulations that illustrate the performance of our method. Proofs are delayed until Section 4.

Construction and properties of the estimators 2.1 Estimation of N by the GL method

We first construct an estimator of N . A natural approach is a kernel method, which is all the more appropriate for comparisons with analytical methods (see [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] for the deterministic analogue). The kernel function K should satisfy the following assumption, in force in the sequel.

Assumption 4 (Assumption on the kernel density estimator). K : R → R is a continuous function such that K(x)dx = 1 and K 2 (x)dx < ∞.

For h > 0 and x ∈ R, define Nh (x) := 1 n n i=1 K h (x -X i ), (2.1)
where

K h (x) = h -1 K(h -1 x). Note in particular that E( Nh ) = K h ⋆ N
, where ⋆ denotes convolution. We measure the performance of Nh via its squared integrated error, i.e. the average L 2 distance between N and Nh . It is easy to see that

E[ N -Nh 2 ] ≤ N -K h ⋆ N 2 + E[ K h ⋆ N -Nh 2 ],
with

E[ K h ⋆ N -Nh 2 2 ] = 1 n 2 E n i=1 K h (x -X i ) -E K h (x -X i ) 2 dx = 1 n 2 n i=1 E K h (x -X i ) -E K h (x -X i ) 2 dx ≤ 1 n E K 2 h (x -X 1 )dx = K h 2 2 n = K 2 2
nh .

Applying the Cauchy-Schwarz inequality, we obtain

E[ N -Nh 2 ] ≤ N -K h ⋆ N 2 + 1 √ nh K 2 .
The first term corresponds to a bias term, it decreases when h → 0. The second term corresponds to a variance term, which increases when h → 0. If one has to choose h in a family H of possible bandwidths, the best choice is h where

h := argmin h∈H N -K h ⋆ N 2 + 1 √ nh K 2 . (2.2)
This ideal compromise h is called the "oracle": it depends on N and then cannot be used in practice. Hence one wants to find an automatic (data-driven) method for selecting this bandwidth. The Lepski method [START_REF] Lepski | One problem of adaptive estimation in Gaussian white noise[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 1. Upper bounds[END_REF][START_REF] Lepski | Asymptotic minimax adaptive estimation. 2. Statistical models without optimal adaptation[END_REF][START_REF] Lepski | On problems of adaptive estimation in white Gaussian noise[END_REF] is one of the various theoretical adaptive methods available for selecting a density estimator. In particular it is the only known method able to select a bandwidth for kernel estimators. However the method do not usually provide a non asymptotic oracle inequality. Recently, Goldenschluger and Lepski [START_REF] Goldenshluger | Uniform bounds for norms of sums of independent random functions[END_REF] developed powerful probabilistic tools that enable to overcome this weakness and that can provide with a fully data-driven bandwidth selection method. We give here a practical illustration of their work: how should one select the bandwidth for a given kernel in dimension 1?

The main idea is to estimate the bias term by looking at several estimators. The method consists in setting first, for any x and any h, h ′ > 0,

Nh,h ′ (x) := 1 n n i=1 (K h ⋆ K h ′ )(x -X i ) = (K h ⋆ Nh ′ )(x).
Next, for any h ∈ H, define

A(h) := sup h ′ ∈H Nh,h ′ -Nh ′ 2 - χ √ nh ′ K 2 + = sup h ′ ∈H max 0, Nh,h ′ -Nh ′ 2 - χ √ nh ′ K 2 ,
where, given ε > 0, we set χ := (1 + ε)(1 + K 1 ). The quantity A(h) is actually a good estimator of N -K h ⋆ N 2 up to the term K 1 (see (4.2) and (4.3) in Section 4). The next step consists then in setting ĥ := arg min

h∈H A(h) + χ √ nh K 2 , (2.3)
and our final estimator of N is obtained by putting N := Nĥ . Let us specify what we are able to prove at this stage.

Proposition 2. Assume N ∈ L ∞ and work under Assumption 4. If H ⊂ {D -1 , D = 1, . . . , D max } with D max = δn for δ > 0, then, for any q ≥ 1, E N -N 2q 2 ≤ C(q)χ 2q inf h∈H K h ⋆ N -N 2q 2 + K 2q 2 (hn) q + C 1 n -q ,
where C(q) is a constant depending on q and C 1 is a constant depending on q, ε, δ, K 2 , K 1 and N ∞ .

The previous inequality is called an oracle inequality, for we have 2q) and ĥ is performing as well as the oracle h up to some multiplicative constant. In that sense, we are able to select the best bandwidth within our family H. [START_REF] Goldenshluger | Uniform bounds for norms of sums of independent random functions[END_REF], we do not consider the case where H is an interval and we do not specify K except for Assumption 4. This simpler method is more reasonable from a numerical point of view, since estimating N is only a preliminary step. The probabilistic tool we use here is classical in model selection theory (see Section 4 and [16]) and actually, we do not use directly [START_REF] Goldenshluger | Uniform bounds for norms of sums of independent random functions[END_REF]. In particular the main difference is that, in our specific case, we are able to get max(H) fixed whereas Goldenschluger and Lepski [START_REF] Goldenshluger | Uniform bounds for norms of sums of independent random functions[END_REF] require max(H) to tend to 0 with n. The price to pay is that we obtain a uniform bound (see Lemma 1 in Section 4.3) which is less tight, but that will be sufficient for our purpose.

E[ N -N 2 ] ≤ (E[ N -N 2q 2 ]) 1/(

Remark 1. As compared to the results of Goldenschluger and Lepski in

Estimation of ∂ ∂x g(x)N(x) by the GL method

The previous method can of course be adapted to estimate

D(x) := ∂ ∂x g(x)N (x) .
We adjust here the work of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] to the setting of estimating a derivative. We again use kernel estimators with more stringent assumptions5 on K.

Assumption 5 (Assumption on the kernel of the derivative estimator). The function

K is differentiable, K(x)dx = 1 and (K ′ (x)) 2 dx < ∞.
For any bandwidth h > 0, we define the kernel estimator of D as

Dh (x) := 1 n n i=1 g(X i )K ′ h (x -X i ) = 1 nh 2 n i=1 g(X i )K ′ x-X i h .
(2.4)

Indeed E( Dh (x)) = K ′ h (x -u)g(u)N (u)du = K ′ h ⋆ (gN ) (x) = K h ⋆ (gN ) ′ (x).
Again we can look at the integrated squared error of D. We obtain the following upper bound:

E[ Dh -D 2 ] ≤ D -K h ⋆ D 2 + E[ K h ⋆ D -Dh 2 ],
with

E[ K h ⋆ D -Dh 2 2 = 1 n 2 E n i=1 g(X i )K ′ h (x -X i ) -E g(X i )K ′ h (x -X i ) 2 dx = 1 n 2 n i=1 E g(X i )K ′ h (x -X i ) -E g(X i )K ′ h (x -X i ) 2 dx ≤ 1 n E g 2 (X 1 )K ′ h 2 (x -X 1 )dx ≤ g 2 ∞ K ′ h 2 2 n = g 2 ∞ K ′ 2 2 nh 3 .
Hence, by Cauchy-Schwarz inequality

E[ D -Dh 2 ] ≤ D -K h ⋆ D 2 + 1 √ nh 3 g ∞ K ′ 2 .
Once again, there is a bias-variance decomposition, but now the variance term is of order

1 √ nh 3 K ′ 2 g ∞ .
We therefore define the oracle by

h := argmin h∈ H D -K h ⋆ D 2 + 1 √ nh 3 g ∞ K ′ 2 . (2.5)
Now let us apply the GL method in this case. Let H be a family of bandwidths. We set for any h, h ′ > 0,

Dh,h ′ (x) := 1 n n i=1 g(X i )(K h ⋆ K h ′ ) ′ (x -X i ) and Ã(h) := sup h ′ ∈ H Dh,h ′ -Dh ′ 2 -χ √ nh ′3 g ∞ K ′ 2 + ,
where, given ε > 0, we put χ := (1 + ε)(1 + K 1 ). Finally, we estimate D by using D := Dh with

h := argmin h∈ H Ã(h) + χ √ nh 3 g ∞ K ′ 2 . (2.6)
As before, we are able to prove an oracle inequality for D.

Proposition 3. Assume N ∈ L ∞ . Work under Assumption 5. If H = {D -1 , D = 1 .
. . , Dmax }, with Dmax = δn for δ > 0, then for any q ≥ 1,

E D -D 2q 2 ≤ C(q) χ2q inf h∈ H K h ⋆ D -D 2q 2 + g ∞ K ′ 2 √ nh 3 2q + C1 n -q ,
where C(q) is a constant depending on q and C1 is a constant depending on q, ε, δ,

K ′ 2 , K ′ 1 , g ∞ and N ∞ .

Estimation of κ (and λ)

As mentioned in the introduction, we will not consider the problem of estimating λ and we work under Assumption 2: an estimator λn of λ is furnished by the practitioner prior to the data processing for estimating B. It becomes subsequently straightforward to obtain an estimator of κ by estimating ρ g (N ), see the form of (1.6). We estimate ρ g (N ) by ρn :

= n i=1 X i n i=1 g(X i ) + c , (2.7)
where c > 0 is a (small) tuning constant 6 . Next we simply put κn = λn ρn .

(2.8)

Approximated inversion of L

From (2.7), the right-hand side of (1.3) is consequently estimated by κn D + λn N .

It remains to formally apply the inverse operator L -1 . However Given ϕ, the dilation equation

L ψ (x) = 4ψ(2x) -ψ(x) = ϕ(x), x ∈ R + (2.9)
admits in general infinitely many solutions, see Doumic et al. [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF], Appendix A. Nevertheless, if ϕ ∈ L 2 , there is a unique solution ψ ∈ L 2 to (2.9), see Proposition A.1. in [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF], and moreover it defines a continuous operator L -1 from L 2 to L 2 . Since K h and gK ′ h belong to L 2 , one can define a unique solution to (2.9) when ϕ = κn D + λn N . This inverse is not analytically known but we can only approximate it via the fast algorithm described below.

Given T > 0 and an integer k ≥ 1, we construct a linear operator L -1 k that maps a function ϕ ∈ L 2 into a function with compact support in [0, T ] as follows. Consider the regular grid on [0, T ] with mesh k -1 T defined by

0 = x 0,k < x 1,k < • • • < x i,k := i k T < . . . < x k,k = T.
We set

ϕ i,k := k T x i+1,k x i ,k ϕ(x)dx for i = 0, . . . , k -1,
and define by induction the sequence 7

H i,k (ϕ) := 1 4 (H i/2,k (ϕ) + ϕ i/2,k ), i = 0, . . . , k -1,
what gives, for i = 0 and i = 1

H 0 (ϕ) := 1 3 ϕ 0,k , H 1 (ϕ) := 4 21 ϕ 0,k + 1 7 ϕ 1,k .
Finally, we define

L -1 k (ϕ)(x) := k-1 i=0 H i,k (ϕ)1 [x i,k ,x i+1,k ) (x). (2.10)
As stated in the introduction, we eventually estimate H = BN by

Ĥ = L -1 k (κ n D + λn N ).
The stability of the inversion is given by the fact that L -1 k : L 2 → L 2 is continuous, see Lemma 4 in Section 4.3, and by the following approximation result between L -1 and L -1 k .

7 for any sequence ui, i = 1, 2, . . ., we define

u i/2 := u i/2 if i is even 1 2 (u (i-1)/2 + u (i+1)/2 ) otherwise.
Proposition 4. Let T > 0 and ϕ ∈ W 1 . Let L -1 (ϕ) denote the unique solution of (2.9) belonging to L 2 . We have for k ≥ 1:

L -1 k (ϕ) -L -1 (ϕ) 2,T ≤ C T √ k ϕ W 1 , with C < 1 √ 6 .
Hence, L -1 k behaves nicely over sufficiently smooth functions. Moreover the estimation of N and the estimators κn and λn are essentially regular. Finally we estimate B as stated in (1.8). The overall behaviour of the estimator is finally governed by the quality of estimation of the derivative D, which determines the accuracy of the whole inverse problem in all these successive steps.

Oracle inequalities for Ĥ

We are ready to state our main result, namely the oracle inequality fulfilled by Ĥ. where Nĥ and Dh are kernel estimators defined respectively by (2.1) and (2.4), and where we have selected ĥ and h by (2.3) and (2.6). Moreover take κn as defined by (2.7) and (2.8) for some c > 0.

The following upper bound holds for any n :

E Ĥ -H q 2,T ≤ C 1 R λ,n inf h∈ H K h ⋆ D -D q 2 + g ∞ K ′ 2 √ nh 3 q + inf h∈H K h ⋆ N -N q 2 + K 2 √ nh q + ε q λ,n + ( N W 1 + gN W 2 ) T √ k q + C 2 n -q 2 ,
where C 1 is a constant depending on q, g, N , ε, ε, K 1 and c; and

C 2 is a constant depending on q, g, N, ε, ε, δ, δ, K 2 , K 1 , K ′ 2 , K ′ 1 .
1) Note that the upper bound quantifies the additive cost of each step used in the estimation method. The first part is an oracle bound on the estimation of D times the size of the estimator of λ. The second part is the oracle bound for N . Of course, the results are sharp only if λ is good, which can be seen through ε λ,n . Finally, the bound is also governed by the approximated inversion through the only term where k appears. The last term is just a residual term, that will be in most of the cases negligible with respect to the other terms. In particular since all the previous errors are somehow unavoidable, this means that, as far as our method is concerned, our upper bound is the best possible that can be achieved in order to select the different bandwidths h, up to multiplicative constants. Moreover, one can see how the limitation in k influences the method and how large k needs to be chosen to guarantee that the main error comes from the fact that one estimates a derivative.

2) The result holds for any n. In particular, we expect the method to perform well for small sample size, see the numerical illustration below. This also shows that we are able to select a good bandwidth as far as the kernel K is fixed, even if there is no assumption on the moments of K and consequently on the approximation properties of K. In the next simulation section, we focus on a Gaussian kernel which has only one vanishing moment (hence one cannot really consider minimax adaptation for regular function with it) but for which the bandwidth choice is still important in practice. The previous result guarantees an optimal bandwidth choice even for this kernel, up to some multiplicative constant.

3) From this oracle result, we can easily deduce the rates of convergence of Proposition 1 at the price of further assumptions on the kernel K, i.e. Assumption 3 defined in the Introduction. Section 1.2.2 of [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF] recalls how to build compactly supported kernels satisfying Assumption 3. If m 0 satisfies Assumption 3, then for any s ≤ m 0 , for any f ∈ W s ,

K h ⋆ f -f 2 ≤ C f W s h s ,
where C is a constant that can be expressed by using K and m 0 (see Theorem 8.1 of [START_REF] Härdle | Wavelets, Approximation and Statistical Applications[END_REF]). Now, it is sufficient to choose h = D -1 of order n -1/(2s+3) to obtain (1.9). The complete proof of Proposition 1 is delayed until the last section.

Numerical illustration

Let us illustrate our method through some simulations.

The numerical protocol

First, we need to build up simulated data: to do so, we depart from given g, κ and B on a regular grid [0, dx, . . . , X M ] and solve the direct problem by the use of the so-called power algorithm to find the corresponding density N and the principal eigenvalue λ (see for instance [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] for more details). We check that N almost vanishes outside the interval [0, X M ]; else, X M has to be increased. The density N being given on this grid, we approximate it by a spline function, and build a n-sample by the rejection sampling algorithm. For the sake of simplicity, we do not simulate an approximation on λ and keep the exact value, thus leading, in the estimate of Theorem 1, to R λ,n = λ 2q and λn = λ.

We then follow step by step the method proposed here and detailed in Section 2.

1. The GL method for the choice of N = Nh. We take the classical Gaussian kernel K(x) = (2π) -1/2 exp -x 2 /2 , set D max = n and limit ourselves to a logarithmic sampling H = {1, 1/2, . . . , 1/9, 1/10, 1/20, . . . , 1/100, 1/200, . . . , 1/n} in order to reduce the cost of computations (The GL method is indeed the most time-consuming step in the numerical protocol).

2. The GL method for the choice of Dh . The procedure is similar except that we choose here D max = √ n. The selected bandwidths h and h can be different. We check that the GL method does not select an extremal point of H.

3. The choice of κn , as defined by (2.8).

4. The numerical scheme described in Section 2.4 and in [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] for the inversion of L.

The division by N and definition of B as described in (1.8).

At each step, we compare, in L 2 -norm, the reconstructed function and the original one: N vs N, ∂ ∂x (g N ) vs ∂ ∂x (gN ), Ĥ vs BN and finally B vs B.

Results on simulated data

We first test the three cases simulated in [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] in which the numerical analysis approach was dealt with. Namely on the interval [0, 4], we consider the cases where g ≡ 1 and first B = B 1 ≡ 1, second B(x) = B 2 (x) = 1 for x ≤ 1.5, then linear to B 2 (x) = 5 for x ≥ 1.7. This particular form is interesting because due to this fast increase on B, the solution N is not that regular and exhibits a 2-peaks distribution (see Figure 1). Finally, we test

B(x) = B 3 (x) = exp(-8(x -2) 2 ) + 1.
In Figures 1 and2, we show the simulation results with n = 5.10 4 (a realistic value for in vitro experiments on E. Coli for instance) for the reconstruction of N, ∂ ∂x (gN ), BN and B.

One notes that the solution can well capture the global behavior of the division rate B, but, as expected, has more difficulties in recovering fine details (for instance, the difference between B 1 and B 3 ) and also gives much more error when B is less regular (case of B 2 ). One also notes that even if the reconstruction of N is very satisfactory, the critical point is the reconstruction of its derivative. Moreover, for large values of x, even if N and its derivative are correctly reconstructed, the method fails in finding a proper division rate B. This is due to two facts: first, N vanishes, so the division by N leads to error amplification. Second, the values taken by B(x) for large x have little influence on the solutions N of the direct problem: whatever the values of B, the solutions N will not vary much, as shown by Figure 1 (left). A similar phenomenon occurred indeed when solving the deterministic problem in [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF] (for instance, we refer to Fig. 10 of this article for a comparison of the results).

We also test a case closer to biological true data, namely the case B(x) = x 2 and g(x) = x. The results are shown on Figures 3 and4 for n-samples of size 10 3 , 5.10 3 , 10 4 and 5.10 4 .

One notes that reconstruction is already very good for N when n = 10 3 , unlike the reconstruction of ∂ ∂x (gN ) that requires much more data. Finally, in Table 3.2 we give average error results on 50 simulations, for n = 1000, g ≡ B ≡ 1. We display the relative errors in L 2 norms, (defined by ||φ -φ|| L 2 /||φ|| L 2 ), and their empirical variances. In Table 3.2, for the case g(x) = x and B(x) = x 2 , we give some results on standard errors for various values of n, and compare them to n -1/5 , which is the order of magnitude of the expected final error on BN, since with a Gaussian kernel we have s = 1 in Proposition 1. We see that our numerical results are in line with the theoretical estimates: indeed, the error on H is roughly twice as large as n -1/5 .

Proofs

In Section 4.1, we first give the proofs of the main results of Section 2. This allows us, in Section 4.2, to prove the results of Section 2.5, which require the collection of all the results of Section 2, i.e. the oracle-type inequalities on the one hand and a numerical analysis result on the other hand. This illustrates the subject of our paper that lies at the frontier between these fields. Finally, we state and prove the technical lemmas used in Section 4.1. These technical tools are concerned with probabilistic results, namely concentration and Rosenthal-type inequalities that are often the main bricks to establish oracle inequalities,

Proofs of the main results of Section 2

Proof of Proposition 2

For any h * ∈ H, we have: We obtain

N -N 2 ≤ Nĥ -N ĥ,h * 2 + N ĥ,h * -Nh * 2 + Nh * -N 2 ≤ A 1 + A 2 + A 3 , with A 1 := Nĥ -N ĥ,h * 2 ≤ A(h * ) + χ n ĥ K 2 , A 2 := N ĥ,h * -Nh * 2 ≤ A( ĥ) + χ √ nh * K 2 and A 3 := Nh * -N 2 .
N -N 2 ≤ A(h * ) + χ n ĥ K 2 + A( ĥ) + χ √ nh * K 2 + Nh * -N 2 ≤ 2A(h * ) + 2χ √ nh * K 2 + Nh * -N 2 . (4.1)
Since we have

A(h * ) = sup h ′ ∈H Nh * ,h ′ -Nh ′ 2 - χ √ nh ′ K 2 + ≤ sup h ′ ∈H Nh * ,h ′ -E[ Nh * ,h ′ ] -Nh ′ -E[ Nh ′ ] 2 - χ √ nh ′ K 2 + + E[ Nh * ,h ′ ] -E[ Nh ′ ] 2 (4.2)
and for any x and any

h ′ E Nh * ,h ′ (x)) -E( Nh ′ (x) = (K h * ⋆ K h ′ )(x -u)N (u)du -K h ′ (x -v)N (v)dv = K h * (x -u -t)K h ′ (t)N (u)dtdu -K h ′ (x -v)N (v)dv = K h * (v -u)K h ′ (x -v)N (u)dudv -K h ′ (x -v)N (v)dv = K h ′ (x -v) K h * (v -u)N (u)du -N (v) dv,
we derive

E( Nh * ,h ′ ) -E( Nh ′ ) 2 ≤ K 1 E h * 2 , (4.3)
where

E h * (x) := (K h * ⋆ N )(x) -N (x), x ∈ R +
represents the approximation term. Combining (4.1), (4.2) and ( 4.3) entails

N -N 2 ≤ Nh * -N 2 + 2 K 1 E h * 2 + 2χ √ nh * K 2 + 2ζ n , with
ζ n := sup h ′ ∈H Nh * ,h ′ -E[ Nh * ,h ′ ] -Nh ′ -E[ Nh ′ ] 2 - χ √ nh ′ K 2 + = sup h ′ ∈H K h * ⋆ Nh ′ -E[ Nh ′ ] -( Nh ′ -E[ Nh ′ ]) 2 - (1 + ε)(1 + K 1 ) √ nh ′ K 2 + ≤ (1 + K 1 ) sup h ′ ∈H Nh ′ -E[ Nh ′ ] 2 - (1 + ε) √ nh ′ K 2 + . Hence E N -N 2q 2 ≤ q E Nh * -N 2q 2 + K 2q 1 E h * 2q 2 +χ 2q K 2q 2 (nh * ) q +(1+ K 1 ) 2q E[ξ 2q n ] ,
where

ξ n = sup h ′ ∈H Nh ′ -E( Nh ′ ) 2 - (1 + ε) √ nh ′ K 2 + .
Now, we have:

E Nh * -N 2q 2 ≤ 2 2q-1 E Nh * -E[ Nh * ] 2q 2 + E[ Nh * ) -N 2q 2 ≤ 2 2q-1 E Nh * -E[ Nh * ] 2q 2 + E h * 2q 2 .

Then, by setting

Kc h * (X i , x) := K h * (x -X i ) -E(K h * (x -X 1 )),
we obtain

E Nh * -E[ Nh * ] 2q 2 = E 1 n n i=1 Kc h * (X i , x) 2 dx q ≤ 2 q-1 n 2q E n i=1 Kc 2 h * (X i , x)dx q + E 1≤i,j≤n i =j Kc h * (X i , x)Kc h * (X j , x)dx q . Since Kc 2 h * (X i , x)dx = K h * (x -X i ) -E K h * (x -X 1 ) 2 dx ≤ 2 K 2 h * (x -X i )dx + E K h * (x -X 1 ) 2 dx ≤ 2 K h * 2 2 + E K 2 h * (x -X 1 ) dx ≤ 4 K h * 2 2 = 4 h * K 2 2 ,
the first term can be bounded as follows

E n i=1 Kc 2 h * (X i , x)dx q ≤ 4n h * K 2 2 q .
For the second term, we apply Theorem 8.1.6 of de la Peña and Giné (1999) (with 2q ≥ 2) combined with the Cauchy-Schwarz inequality:

E 1≤i,j≤n i =j Kc h * (X i , x)Kc h * (X j , x)dx q ≤ E 1≤i,j≤n i =j Kc h * (X i , x)Kc h * (X j , x)dx 2q 1 2 ≤ q n q E Kc h * (X 1 , x)Kc h * (X 2 , x)dx 2q 1 2 ≤ q n q E Kc 2 h * (X 1 , x)dx 2q 1 2 ≤ q 4n h * K 2 2 q .
It remains to deal with the term E(ξ 2q n ). By Lemma 1 below, we obtain

E[ξ 2q n ] ≤ q,η,δ K 2 , K 1 , N ∞ n -q
and the conclusion follows.

Proof of Proposition 3

The proof is similar to the previous one and we avoid most of the computations for simplicity. For any

h 0 ∈ H, Dh -D 2 ≤ Dh -Dh ,h 0 2 + Dh ,h 0 -Dh 0 2 + Dh 0 -D 2 ≤ Ã1 + Ã2 + Ã3 , with Ã1 := Dh -Dh ,h 0 2 ≤ Ã(h 0 ) + χ n h3 g ∞ K ′ 2 , Ã2 := Dh ,h 0 -Dh 0 2 ≤ Ã( h) + χ nh 3 0 g ∞ K ′ 2 and Ã3 := Dh 0 -D 2 .
Then,

Dh -D 2 ≤ 2 Ã(h 0 ) + 2 χ nh 3 0 g ∞ K ′ 2 + Dh 0 -D 2 .
To study Ã(h 0 ), we first evaluate

E[ Dh 1 ,h 2 (x)] -E[ Dh 2 (x)]. = (K h 1 ⋆ K h 2 ⋆ (gN ) ′ )(x) -(K h 2 ⋆ (gN ) ′ )(x) = D(t)(K h 1 ⋆ K h 2 )(x -t)dt -D(t)K h 2 (x -t)dt = D(t) K h 1 (x -t -u)K h 2 (u)dudt -D(t)K h 2 (x -t)dt = D(t) K h 1 (v -t)K h 2 (x -v)dvdt -D(v)K h 2 (x -v)dv = K h 2 (x -v) D(t)K h 1 (v -t)dt -D(v) dv = (K h 2 ⋆ Ẽh 1 )(x),
where we set, for any real number x

Ẽh 1 (x) := D(t)K h 1 (x -t)dt -D(x) = (K h 1 ⋆ D)(x) -D(x). (4.4)
It follows that

Ã(h 0 ) = sup h∈ H Dh 0 ,h -Dh 2 -χ √ nh 3 g ∞ K ′ 2 + ≤ sup h∈ H Dh 0 ,h -E[ Dh 0 ,h ] -Dh -E[ Dh ] 2 -χ √ nh 3 g ∞ K ′ 2 + + E[ Dh 0 ,h ] -E[ Dh ] 2 ≤ sup h∈ H Dh 0 ,h -E[ Dh 0 ,h ] -( Dh -E[ Dh ]) 2 -χ √ nh 3 g ∞ K ′ 2 + + K 1 Ẽh 0 2 ≤ (1 + K 1 ) sup h∈ H Dh -E[ Dh ] 2 - (1 + ε) √ nh 3 g ∞ K ′ 2 + + K 1 Ẽh 0 2 , (4.5)
In order to obtain the last line, we use the following chain of arguments:

Dh 0 ,h (x) = 1 n n i=1 g(X i ) K ′ h (x -X i -t)K h 0 (t)dt = K h 0 (t) 1 n n i=1 g(X i )K ′ h (x -X i -t) dt and E Dh 0 ,h (x) = K h 0 (t) g(u)K ′ h (x -u -t)N (u)du dt, therefore Dh 0 ,h (x) -E Dh 0 ,h (x) = K h 0 (t)G(x -t)dt = K h 0 ⋆ G(x), with G(x) = 1 n n i=1 g(X i )K ′ h (x -X i ) -g(u)K ′ h (x -u)N (u)du = Dh (x) -E Dh (x) . Therefore Dh 0 ,h -E[ Dh 0 ,h ] 2 ≤ K h 0 1 G 2 ≤ K 1 Dh -E[ Dh ] 2 ,
which justifies (4.5). In the same way as in the proof of Proposition 2, we can establish the following:

E Ẽh 0 2q 2 = E Dh 0 -D 2q 2 ≤ q Ẽh 0 2q 2 + g ∞ K ′ 2 nh 3 0 2q
.

Finally, we successively apply (4.4), (4.5) and Lemma 1 in order to conclude the proof.

Proof of Proposition 4

We use the notation and definitions of Section 2.4. We have

L -1 k (ϕ) -L -1 (ϕ) 2 2,T = k-1 i=0 x i+1,k x i,k H i,k -L -1 (ϕ)(x) 2 dx := k-1 i=0 L i,k .
We prove by induction that for all i, one has

L i,k ≤ C 2 T 2 k 2 ϕ 2 W 1 .
The result follows by summation over i. We first prove the two following estimates:

x i+1,k x i,k ϕ i,k -ϕ(x) 2 dx ≤ T 2 4π 2 k 2 ϕ 2 W 1 , (4.6
)

|ϕ i+1,k -ϕ i,k | 2 ≤ T k ϕ 2 W 1 . (4.7)
By definition, ϕ i,k is the average of the function ϕ on the interval [x i,k , x i+1,k ] of size T k . Thus (4.7) is simply Wirtinger inequality applied to ϕ ∈ W 1 on the interval [x i,k , x i+1,k ].F or(4.7), we use the Cauchy-Schwarz inequality:

|ϕ i+1,k -ϕ i,k | 2 = k 2 T 2 x i+1,k x i,k ϕ(x + T k ) -ϕ(x) dx 2 = k 2 T 2 x i+1,k x i,k x+ T k x ϕ ′ (z)dz dx 2 ≤ k 2 T 2 x i+1,k x i,k T k ϕ W 1 dx 2 = T k ϕ 2 W 1 .
We are ready to prove by induction the two following inequalities:

L i,k ≤ C 2 1 T 2 k 2 ϕ 2 W 1 , (4.8)
|H i+1,k -H i,k | 2 ≤ C 2 2 T k ϕ 2 W 1 . (4.9)
for two constants C 1 and C 2 specified later on. First, for i = 0, we have

L 0,k = T k 0 |H 0,k (ϕ) -L -1 (ϕ)(x)| 2 dx = T k 0 1 3 ϕ 0,k -L -1 (ϕ)(x) 2 dx.
We recall (see Proposition A.1. of [START_REF] Doumic | Numerical Solution of an Inverse Problem in Size-Structured Population Dynamics[END_REF]) that

L -1 (ϕ)(x) = ∞ n=1 2 -2n ϕ(2 -n
x), and we use the

fact that 1 3 = ∞ n=1 2 -2n and for a, b > 0, ab ≤ 1 2 (a 2 + b 2 ) in order to write L 0,k = T k 0 ∞ n=1 2 -2n ϕ 0,k -ϕ(2 -n x) 2 dx ≤ ∞ n,n ′ =1 2 -2n-2n ′ T k 0 |ϕ 0,k -ϕ(2 -n x)| 2 dx ≤ 1 3 ∞ n=1 2 -n 2 -n T k 0 |ϕ 0,k -ϕ(y)| 2 dy ≤ 1 3 T 2 4π 2 k 2 ϕ 2 W 1 .
This proves the first induction assumption for i = 0, and

|H 1,k -H 0,k | 2 = 1 7 (ϕ 1,k -ϕ 0,k ) 2 ≤ 1 7 2 T k ϕ 2 W 1 ,
proves the second one. Let us now suppose that the two induction assumptions are true for all j ≤ i -1, and take i ≥ 1. Let us first evaluate

L i,k = x i+1,k x i,k H i,k -L -1 (ϕ)(x) 2 dx = 1 16 x i+1,k x i,k H i 2 ,k + ϕ i 2 ,k -L -1 (ϕ)( x 2 ) -ϕ( x 2 ) 2 dx.
We distinguish the case when i is even and when i is odd. Let i be even: then, by definition

L i,k ≤ 1 8 x i+1,k x i,k H i 2 ,k -L -1 (ϕ)( x 2 ) 2 dx+ 1 8 x i+1,k x i,k ϕ i 2 ,k -ϕ( x 2 ) 2 dx ≤ 1 4 (C 2 1 + 1 4π 2 ) T 2 k 2 ϕ 2 W 1
by the induction assumption and Assertion 4.7 on ϕ for j = i 2 . If i is odd, we write by definition

L i,k = 1 16 x i+1,k x i,k H i-1 2 ,k +ϕ i-1 2 ,k -L -1 (ϕ)( x 2 )-ϕ( x 2 )+ 1 2 (H i+1 2 ,k -H i-1 2 ,k )+ 1 2 (ϕ i+1 2 ,k -ϕ i- 1 2 ,k ) 2 dx.
Hence, re-organizing terms, we can write

L i,k ≤ 1 2 x i-1 2 +1,k x i-1 2 ,k H i-1 2 ,k -L -1 (ϕ)(y) 2 dy + 1 2 x i-1 2 +1,k x i-1 2 ,k ϕ i-1 2 ,k -ϕ(y) 2 dy + 1 16 T k (H i+1 2 ,k -H i-1 2 ,k ) 2 + 1 16 T k (ϕ i+1 2 ,k -ϕ i-1 2 ,k ) 2
Putting together the four inequalities above (the estimates for ϕ and the induction assumptions), we obtain

L i,k ≤ T 2 k 2 ϕ 2 W 1 C 2 1 2 + 1 8π 2 + C 2 2 16 + 1 16
and (4.8) is proved. It remains to establishe (4.9). Let us write it for i even (the case of an odd i is similar):

|H i+1,k -H i,k | 2 = 1 16 |H i+1 2 -H i 2 + ϕ i+1 2 -ϕ i 2 | 2 = 1 32 |H i 2 +1 -H i 2 + ϕ i 2 +1 -ϕ i 2 | 2 .
Hence, as previously, we obtain

|H i+1,k -H i,k | 2 ≤ 1 16 T k ϕ 2 W 1 (C 2 2 + 1).
To complete the proof, we remark that C

Proof of Theorem 1 and Proposition 1

Proof of Theorem 1

It is easy to see that

H -H 2,T = L -1 k (κ D + λn N ) -L -1 (L(BN )) 2,T ≤ L -1 k (κ D + λn N ) -L -1 k (L(BN )) 2,T + L -1 k (L(BN )) -L -1 (L(BN )) 2,T ≤ L -1 k (κ D + λn N -(κD + λN )) 2,T + 1 3 T √ k L(BN ) W 1 ,
thanks to Proposition 4. Note that L(BN ) = κ(gN ) ′ + λN so that we can write

L(BN ) W 1 ≤ C( N W 1 + gN W 2 ).
We obtain, thanks to Lemma 4 that gives the boundedness of the operator L -1 k :

H -H 2,T ≤ κn D -κD 2,T + λn N -λN 2,T + ( N W 1 + gN W 2 ) T √ k ≤ |κ n | D -D 2 + | λn | N -N 2 + |κ n -κ| D 2 + | λn -λ| N 2 + ( N W 1 + gN W 2 ) T √ k ≤ | λn ||ρ n | D -D 2 + | λn | N -N 2 + |ρ n -ρ g (N )| D 2 + N 2 + |ρ g (N )| D 2 | λn -λ| + ( N W 1 + gN W 2 ) T √ k .
Taking expectation and using Cauchy-Schwarz inequality, we obtain for any q ≥ 1,

E[ H -H q 2,T ] ≤ q E[ λ2q n ] 1/2 E[ρ 4q n ] 1/4 E[ D -D 4q 2 ] 1/4 + E[ N -N 2q 2 ] 1/2 + D q 2 E[|ρ n -ρ g (N )| 2q ] 1/2 + N 2 + ρ g (N ) D 2 q E[| λn -λ| q ] + ( N W 1 + gN W 2 ) T √ k q .
Now, Lemma 3 gives the behaviour of E[|ρ n -ρ g (N )| 2q ]. In particular, we obtain E[ρ 4q n ] ≤ q,g,N,c . We finally apply successively Propositions 2 and 3 to obtain the proof of Theorem 1.

Proof of Proposition 1

We have already proved (1.9). It remains to prove (1.10). We introduce the event

Ω n = {2 N (x) ≥ m for any x ∈ [a, b]}. Then, for n larger that Q 2 , E b a B(x) -B(x) 2 dx q 2 = E b a B(x) -B(x) 2 dx × 1 Ωn q 2 + E b a B(x) -B(x) 2 dx × 1 Ω c n q 2 ≤ E b a B(x) -B(x) 2 dx × 1 Ωn q 2 + 2(b -a)(n + Q 2 ) q 2 P(Ω c n ) ≤ E b a Ĥ N - H N 2 × 1 Ωn q 2 + 4n(b -a) q 2 P(Ω c n ) ≤ E b a ĤN -N H N N 2 × 1 Ωn q 2 + 4n(b -a) q 2 P(Ω c n ) ≤ q,m,M,Q E Ĥ -H q 2 + E N -N q 2
+ 4n(b -a) q 2 P(Ω c n ). The first term of the right hand side is handled by (1.9) and Proposition 2. The second term is handled by Lemma 2 that establishes that P(Ω c n ) = O(n -q ).

Technical lemmas

Concentration inequalities

We first state the following concentration result. Note that a more general version of this result can be found in [START_REF] Goldenshluger | Uniform bounds for norms of sums of independent random functions[END_REF]. We nevertheless give a proof for the sake of completeness.

Lemma 1. We have the following estimates

• Assume that K 2 , K 1 , g ∞ and N ∞ are finite. For every q > 0, introduce the grid H ⊂ {1, 1/2, ..., 1/D max } and D max = δn for some δ > 0. Then, for every η > 0

E sup h∈H Nh -E[ Nh ] 2 - (1 + η) √ nh K 2 2q + ≤ q,η,δ K 2 , K 1 , N ∞ n -q .
• Assume that K ′ 2 , K ′ 1 , g ∞ and N ∞ are finite. For every q > 0, introduce the grid H ⊂ {1, 1/2, ..., 1/ Dmax } and Dmax = δn for some δ > 0. Then for every η > 0

E sup h∈ H Dh -E[ Dh ] 2 - (1 + η) √ nh 3 g ∞ K ′ 2 2q + ≤ q,η, δ, K ′ 2 , K ′ 1 , g ∞, N ∞ n -q .
Proof. Let X be a real random variable and let us consider the random process

∀ t ∈ R, w(t, X) = ϕ(X)Ψ(t -X),
where ϕ and Ψ are measurable real-valued functions. Let X 1 , ..., X n be n independent and identically distributed random variables with the same distribution as X and let us consider the process

∀ t ∈ R, ξ ϕ,Ψ (t) = n i=1 w(t, X i ) -E w(t, X) .
First, let us study the behavior of Hence one can apply Talagrand's inequality (see the version of Bousquet in the independent and identically distributed case [16, p 170]). For all ε, x > 0, one has:

ξ ϕ,Ψ 2 = ξ 2 ϕ,Ψ (t)
P ξ ϕ,Ψ 2 ≥ (1 + ε) E[ ξ ϕ,Ψ 2 ] + √ 2vx + c(ε)bx ≤ e -x ,
where c(ε)

= 1/3 + ε -1 , v = n sup a∈A E a(t) w(t, X) -E w(t, X) dt 2 ,
and b = sup y∈R,a∈A a(t) w(t, y) -E w(t, X) dt.

We study each term of the right hand term within the expectation.

• Obviously, one has:

E ξ ϕ,Ψ 2 ≤ E n i=1 w(t, X i ) -E w(t, X) 2 dt 1/2 = n i=1 E w(t, X i ) -E w(t, X) 2 dt 1/2 ≤ n E w(t, X) 2 dt 1/2
.

But we easily see that for all y,

w 2 (t, y)dt ≤ ϕ 2 ∞ Ψ 2 2 , (4.10)
hence

E[ ξ ϕ,Ψ 2 ] ≤ √ n ϕ ∞ Ψ 2 .
• Since v is a supremum of variance terms,

v ≤ n sup a∈A E a(t)w(t, X)dt 2 ≤ n sup a∈A E |w(t, X)|dt a 2 (t)|w(t, X)|dt ≤ n sup y∈R |w(t, y)|dt × sup t∈R E[|w(t, X)|] ≤ n ϕ 2 ∞ N ∞ Ψ 2 1 .
• The Cauchy-Schwarz inequality and (4.10) give

b = sup y∈R w(•, y) -E[w(., X)] 2 ≤ sup y∈R w(•, y) 2 + E w(t, X) 2 dt 1/2 ≤ 2 ϕ ∞ Ψ 2 .
The main point here is that √ v may be much smaller than E[ ξ ϕ,Ψ 2 ]. So, for all ε, x > 0,

P ξ ϕ,Ψ 2 ≥ (1 + ε) √ n ϕ ∞ Ψ 2 + ϕ ∞ N 1/2 ∞ Ψ 1 √ 2nx + 2c(ε) ϕ ∞ Ψ 2 x ≤ e -x .
Now we consider a family M of possible functions ϕ and Ψ. Let us introduce some strictly positive weights L ϕ,Ψ and let us apply the previous inequality to x = L ϕ,Ψ + u for u > 0.

Hence with probability larger than 1 -(ϕ,Ψ)∈M e -L ϕ,Ψ e -u , for all (ϕ, Ψ) ∈ M, one has

ξ ϕ,Ψ 2 ≤ (1 + ε) √ n ϕ ∞ Ψ 2 + ϕ ∞ N 1/2 ∞ Ψ 1 2nL ϕ,Ψ + 2c(ε) ϕ ∞ Ψ 2 L ϕ,Ψ + ϕ ∞ N 1/2 ∞ Ψ 1 √ 2nu + 2c(ε) ϕ ∞ Ψ 2 u. (4.11) Let M ϕ,Ψ = (1 + ε) √ n ϕ ∞ Ψ 2 + ϕ ∞ N 1/2 ∞ Ψ 1 2nL ϕ,Ψ + 2c(ε) ϕ ∞ Ψ 2 L ϕ,Ψ .
It is also easy to obtain an upper bound of R q for any q ≥ 1 with

R q = E sup (ϕ,Ψ)∈M ξ ϕ,Ψ 2 -M ϕ,Ψ 2q + = +∞ 0 P sup (ϕ,Ψ)∈M ξ ϕ,Ψ 2 -M ϕ,Ψ 2q + ≥ x dx. Indeed R q ≤ (ϕ,Ψ)∈M +∞ 0 P ξ ϕ,Ψ 2 -M ϕ,Ψ 2q + ≥ x dx.
Then, let us take u such that

x = f (u) 2q := ϕ ∞ N 1/2 ∞ Ψ 1 √ 2nu + 2c(ε) ϕ ∞ Ψ 2 u 2q , so dx = 2q f (u) 2q-1 √ 2n ϕ ∞ N 1/2 ∞ Ψ 1 1 2 √ u + 2c(ε) ϕ ∞ Ψ 2 du. Hence R q ≤ (ϕ,Ψ)∈M +∞ 0 e -(L ϕ,Ψ +u) 2q f (u) 2q-1 √ 2n ϕ ∞ N 1/2 ∞ Ψ 1 1 2 √ u + 2c(ε) ϕ ∞ Ψ 2 du ≤ 2q (ϕ,Ψ)∈M e -L ϕ,Ψ +∞ 0 f (u) 2q e -u u -1 du, ≤ q,ε (ϕ,Ψ)∈M e -L ϕ,Ψ n q N q ∞ ϕ 2q ∞ Ψ 2q 1 +∞ 0 e -u u q-1 du + ϕ 2q ∞ Ψ 2q 2 +∞ 0 e -u u 2q-1 du .
Finally, we have proved that • First, when ϕ = 1/n and Ψ = 1/hK(•/h), the family M corresponds to the family H. In that case M ϕ,Ψ and L ϕ,Ψ will respectively be denoted by M h and L h . The upper bound given in (4.12) becomes

R q ≤ q,ε (ϕ,Ψ)∈M e -L ϕ,Ψ n q N q ∞ ϕ 2q ∞ Ψ 2q 1 + ϕ 2q ∞ Ψ 2q 2 . ( 4
R q ≤ q,ε h∈H e -L h n -q N q ∞ K 2q 1 + n -2q h -q K 2q 2 . (4.13)
Now it remains to choose L h . But

M h = (1 + ε) 1 √ nh K 2 + N 1/2 ∞ K 1 2L h n + 2c(ε) 1 n √ h K 2 L h .
Let θ > 0 and let

L h = θ 2 K 2 2 2 N ∞ K 2 1 √ h .
Obviously the series in (4.13) is finite and for any h ∈ H, since h ≤ 1, we have:

M h ≤ (1 + ε + θ) K 2 √ nh + c(ε)θ 2 K 3 2 N ∞ K 2 1 1 nh .
Since D max = δn, one obtains that

M h ≤ 1 + ε + θ + c(ε)θ 2 K 3 2 √ δ N ∞ K 2 1 K 2 √ nh .
It remains to choose ε = η/2 and θ small enough such that θ +

c(ε)θ 2 K 3 2 √ δ N ∞ K 2 1 = η/2
to obtain the desired inequality.

• Secondly, if ϕ = g/n and Ψ = 1/h 2 K ′ (•/h)the family M corresponds to the family H. So, M ϕ,Ψ and L ϕ,Ψ will be denoted by M ′ h and L ′ h respectively. The upper bound given by (4.12) now becomes

R q ≤ q,ε h∈ H e -L ′ h n -q h -2q N q ∞ g 2q ∞ K ′ 2q 1 + n -2q h -3q g 2q ∞ K ′ 2q 2 . (4.14) But M ′ h = (1+ε) 1 √ nh 3/2 g ∞ K ′ 2 + g ∞ N 1/2 ∞ K ′ 1 2L ′ h nh 2 +2c(ε) g ∞ K ′ 2 L ′ h nh 3/2 .
Let θ > 0 and let

L ′ h = θ 2 K ′ 2 2 2 N ∞ K ′ 2 1 h
.

Obviously the series in (4.14) is finite and we have:

M ′ h ≤ (1 + ε + θ) g ∞ K ′ 2 √ nh 3 + c(ε)θ 2 K ′ 3 2 g ∞ N ∞ K ′ 2 1 1 nh 5/2 . But h 2 ≥ ( δn) -1 . Hence M ′ h ≤ 1 + ε + θ + c(ε)θ 2 K ′ 3 2 g ∞ δ N ∞ K ′ 2 1 K ′ 2 √ nh 3 .
As previously, it remains to choose ε and θ accordingly to conclude.

The second result is based on probabilistic arguments as well.

Lemma 2. Under Assumptions and notations of Proposition 1, if there exists an interval

[a, b] in (0, T ) such that [m, M ] := [ inf x∈[a,b] N (x), sup x∈[a,b] N (x)] ⊂ (0, ∞), Q := sup x∈[a,b] |H(x)| < ∞,
and if ln(n) ≤ D min ≤ n 1/(2m 0 +1) and n 1/5 ≤ D max ≤ (n/ ln(n)) 1/(4+η) for some η > 0 fixed, then there exists C η , a constant depending on η, such that for n large enough,

P(Ω c n ) ≤ C η n -q .
Proof. Let x, y be two fixed point of [a, b] and for h ∈ H ′ , first let us look at

Z h (x, y) = Nh (x) -Nh (y) -[K h ⋆ N (x) -K h ⋆ N (y)].
One can apply Bernstein inequality to Z h (x, y) (see e.g. (2.10) and (2.21) of [START_REF] Massart | Concentration inequalities and model selection[END_REF]) to obtain:

E e λZ h (x,y) ≤ exp λ 2 v 2 (x, y) 2(1 -λc(x, y)) , ∀λ ∈ (0, 1/c(x, y)), with c(x, y) = 1 3n sup u∈R |K h (x -u) -K h (y -u)| and v 2 (x, y) = 1 n R |K h (x -u) -K h (y -u)| 2 N (u)du. But, with K the Lipschitz constant of K, c(x, y) ≤ K 3nh 2 |x -y| and v 2 (x, y) ≤ K 2 nh 4 |x -y| 2 . Let x 0 be a fixed point of [a, b].
One can apply Theorem 2.1 of [START_REF] Baraud | A Bernstein-type inequality for suprema of random processes with applications to model selection in non-Gaussian regression[END_REF] which gives that, if

Z h = sup x∈[a,b] |Z h (x, x 0 )| then for all positive x, P(Z h ≥ 18( v 2 (x + 1) + c(x + 1)) ≤ 2e -x , with v 2 = K 2 nh 4 |b -a| 2 and c = K 3nh 2 |b -a|.
But one can similarly prove, using Bernstein's inequality, that

P | Nh (x 0 ) -K h ⋆ N (x 0 )| ≥ N,K x nh + K x nh ≤ 2e -x . Hence P ∃x ∈ [a, b], | Nh (x) -K h ⋆ N (x)| ≥ N,K x nh 4 + K
x nh 2 ≤ 4e -x . We apply this inequality with x = q log(n) + D η , D = D ′ min , . . . , D ′ max . We obtain

P ∃h ∈ H ′ , ∃x ∈ [a, b], | Nh (x) -K h ⋆ N (x)| ≥ N,K,q log(n) nh 4+η + K,q log(n) nh 2+η ≤ 4
Dmax D=1 e -D η n -q = η n -q .

Consequently, since log(n) nh 4+η → 0 uniformly in h ∈ H ′ , for n large enough,

P(∃h ∈ H ′ , ∃x ∈ [a, b], | Nh (x) -K h ⋆ N (x)| ≥ m/4) ≤ η n -q .
But since N ∈ W s+1 , s ≥ 1, then N ′ ∈ W s so N ′ is bounded on R. Then

Rosenthal-type inequality

The following result studies the behavior of the moments of ρn . We recall that Assumption 1 ensures that E[X p i ] = x p N (x)dx < ∞, for any p ≥ 1. Hence, for the first term A n , using (4.15), we have: Let us turn to the term B n : (recall that Assumption 1 states that g ∞ < ∞, which also implies E[g(X 1 ) 2 ] < ∞). Since g is positive, the Bernstein inequality (see Section 2.2.3 of [START_REF] Massart | Concentration inequalities and model selection[END_REF]) gives: The boundedness of L -1 k Lemma 4. For any function ϕ, we have:

A n = E
B n = E n i=1 X i n i=1 g(X i ) + c - n i=1 X i n g(x)N (x)dx + c p ≤ E 1 n n i=1 X i 2p 1/2 × E 1 n n i=1 g(X i ) -g(x)N (x)dx
P(Ω n ) ≥ 1 -2n -γ . ( 4
L -1 k (ϕ) 2,T ≤ 1 3 ϕ 2,T .
Proof. We have:

|L -1 k (ϕ)(x)| 2 dx = T k k-1 i=0 |H i,k (ϕ)| 2 = T k k-1 i=0 1 4 (H i 2 ,k (ϕ) + ϕ i 2 ,k ) 2 = T k [ k-1 2 ] j=0 1 16 H j,k + ϕ j,k 2 + [ k 2 ] j=1 1 64 (H j,k + ϕ j,k + H j-1,k + ϕ j-1,k ) 2 ≤ T k [ k-1 2 ] j=0 1 8 H 2 j,k + ϕ 2 j,k + [ k 2 ] j=1 1
16 (H 2 j,k + ϕ 2 j,k + H 2 j-1,k + ϕ 2 j-1,k

≤ 1 4 T k k-1 i=0 H 2 j,k + ϕ 2 j,k = 1 4 |L -1 k (ϕ)(x)| 2 dx + T k k-1 i=0 ϕ 2 j,k .
At the second line, we have distinguished the i's that are even and the i's that are odd. At the third line, we have used the inequalities (a + b) 2 ≤ 2a 2 + 2b 2 and (a + b + c + d) 2 ≤ 4(a 2 + b 2 + c 2 + d 2 ). By substraction, we obtain

|L -1 k (ϕ)(x)| 2 dx ≤ 1 3 T k k-1 i=0 ϕ 2 j,k .
The Cauchy-Schwarz inequality gives:

x i+1,k x i,k ϕ 2 (x)dx 2 ≤ T k x i+1,k
x i,k ϕ 2 (x)dx, so that

T k k-1 i=0 ϕ 2 j,k = T k k-1 i=0 k 2 T 2 x i+1,k x i,k ϕ(x)dx 2 ≤ T 0 ϕ 2 (x)dx
and finally we obtain the desired result:

|L -1 k (ϕ)(x)| 2 dx ≤ 1 3 ϕ 2 (x)dx.

Theorem 1 .

 1 Work under Assumptions 1 and 2 and let K a kernel satisfying Assumptions 4 and 5. Define H ⊂ {D -1 , D = 1 . . . , D max } with D max = δn for δ > 0 and H ⊂ {D -1 : D = 1, ..., Dmax } with Dmax = δn for δ > 0. For k ≥ 1 and T > 0, let us define L -1 k by (2.10) on the interval [0, T ]. Finally, define the estimator Ĥ = L -1 k ( λn Nĥ + κn Dh ),

Figure 1 :

 1 Figure 1: Reconstruction of N (left) and of ∂ ∂x (gN ) (right) obtained with a sample of n = 5.10 4 data, for three different cases of division rates B.

Figure 2 :

 2 Figure 2: Reconstruction of BN (left) and of B (right) obtained with a sample of n = 5.10 4 data, for three different cases of division rates B.

Figure 3 :

 3 Figure 3: Reconstruction of N (left) and of ∂ ∂x (gN ) (right) obtained for g(x) = x and B(x) = x 2 , for various sample sizes.

Figure 4 :

 4 Figure 4: Reconstruction of BN (left) and of B (right) obtained for g(x) = x and B(x) = x 2 , for various sample sizes.

. 12 )

 12 Now let us evaluate what this inequality means for each set-up.

KN= h tK t 1 0N

 1 h ⋆ N (x) -N (x) = K h (x -y)(N (y) -N (x)′ (x + u(y -x))dudy ′ (x -uht)dudt. So, |N (x) -K h ⋆ N (x)| ≤ N,K h.Because of the definition of H ′ this term tends to 0, uniformly in x and h. Hence for n large enough, for all x ∈ [a, b] and all h ∈ H ′K h ⋆ N (x) ≥ N (x) -m/4 ≥ 3m/4.Consequently, for n large enough,P(∃h ∈ H ′ , ∃x ∈ [a, b], 2 Nh (x) < m) ≤ η n -q .

Lemma 3 .. 1 p

 31 For any p ≥ 2,E |ρ n -ρ g (N )| p ≤ p,g,N,c n -p/2 .Proof. We have:E |ρ n -ρ g (N )| p ≤ p (A n + B n ),withA n := E n i=1 X i n g(x)N (x)dx + c -R + xN (x)dx R + g(x)N (x)dx p and B n := E n i=1 X i n i=1 g(X i ) + c -n i=1 X i n g(x)N (x)dx + c pWe use the Rosenthal inequality (see e.g. the textbook[START_REF] Hall | Martingale Limit Theory and Its Applications[END_REF]): if R 1 , . . . , R n are independent centered variables such that E[|R 1 | p ] < ∞, with p ≥ 2 then ≤ p n -p n E |R 1 | p + n E R 2

  N,c n -p ≤ p,g,N,c n -p/2 .

  X i ) -g(x)N (x)dx i ) -g(x)N (x)dx ≤ 2γ Var(g(X 1 )) log n n + γ g ∞ log n 3n

. 16 ) 2 ≤

 162 Therefore, we bound from above the term B n by a constant p,g,N timesE p,g,N,c n -p/2 + p,g,N √ 2(nc -1 g ∞ ) p n -γ/2 ≤ p,g,N,c n -p/2 ,where we have used (4.15) for the first term and (4.16) for the second one. This concludes the proof of the lemma.

  dt If B denotes the unit ball in L 2 and A is a countable subset of B, we have

			1/2	.
	a∈B ξ ϕ,Ψ 2 = sup	a(t)ξ ϕ,Ψ (t)dt
	= sup	a(t)ξ ϕ,Ψ (t)dt
	a∈A
		n
	= sup a∈A	i=1	a(t) w(t, X i ) -E w(t, X) dt.

In the sequel, we denote by E(•) the expectation operator with respect to P(•) likewise.

Mathematically sepaking, this only amounts to enlarge the probability space to a rich enough structure that captures this estimator. We do not pursue that here.

although here the problem is nonlinear, but that will not affect the argument.

The statistician reader willl note that the rate n -s/(2s+3) is also the minimax rate of convergence when estimating the derivative of a density, see[START_REF] Gasser | Optimal convergence properties of kernel estimates of derivatives of a density function[END_REF].

For sake of simplicity we use the same kernel to estimate N and D but this choice is not mandatory.

In practice, one can take c = 0.

Acknowledgment: We warmly thank Oleg Lepski for fruitful discussions that eventually led to the choice of his method for the purpose of this article. The research of M. Hoffmann is supported by the Agence Nationale de la Recherche, Grant No. ANR-08-BLAN-0220-01. The research of P. Reynaud-Bouret and V. Rivoirard is partly supported by the Agence Nationale de la Recherche, Grant No. ANR-09-BLAN-0128 PARCIMONIE. The research of M. Doumic is partly supported by the Agence Nationale de la Recherche, Grant No. ANR-09-BLAN-0218 TOPPAZ.