Nico Dellaert
email: n.p.dellaert@tue.nl

Jully Jeunet
email: jully.jeunet@dauphine.fr

Hospital admission planning to optimize major resources utilization under uncertainty

Keywords: operation theatre planning, intensive and medium care resources allocation, patient mix, tactical plan, operational schedule, integer linear programming

Admission policies for elective inpatient services mainly result in the management of a single resource: the operating theatre as it is commonly considered as the most critical and expensive resource in a hospital. However, other bottleneck resources may lead to surgery cancellations, such as bed capacity and nursing staff in Intensive Care (IC) units and bed occupancy in wards or medium care (MC) services. Our incentive is therefore to determine a master schedule of a given number of patients that are divided in several homogeneous categories in terms of the utilization of each resource: operating theatre, IC beds, IC nursing hours and MC beds. The objective is to minimize the weighted deviations of the resource use from their targets and probabilistic lengths of stay in each unit (IC and MC) are considered. We use a Mixed Integer Program model to determine the best admission policy. The resulting admission policy is a tactical plan, as it is based upon the expected number of patients with their expected characteristics. On the operational level, this tactical plan must be adapted to account for the actual arriving number of patients in each category. We develop several strategies to build an operational schedule that leans upon the tactical plan more or less closely. The strategies result from the combination of several options to create a feasible operational schedule from the tactical plan: overplanning, flexibility in selecting the patient groups to be operated and updating the tactical plan. The strategies were tested on real data from a Thoracic Surgery Centre over a 10-year simulation horizon. The performance was assessed by the average waiting time for patients, the weighted target deviations and some indicators of the plan changes between the tactical plan and the operational schedule. Simulation results show that the best strategies include overplanning, a limited flexibility and infrequent updates of the tactical plan.

Introduction

Booked hospital admission strategies for elective inpatient services most often rely on the optimization of the operating theatre, as it is commonly considered as the most critical and expensive resource. For instance, Hans et al. [START_REF] Hans | Robust surgery loading[END_REF] address the robust surgery loading problem which consists in assigning daily surgeries to operating rooms while minimizing the risk of overtime as durations of operation are uncertain. In addition to the operating room department, other references also consider the number of beds as a critical resource. Beliën and Demeulemeester [START_REF] Belien | Building cyclic master surgery schedules with leveled resulting bed occupancy[END_REF] design a stochastic model to minimize the beds shortages under capacity constraints related to the number or operating rooms and to the allocation of specific numbers of blocks to each surgeon. References addressing the bed capacity problem in itself generally consider emergency admissions together with elective patient admissions. Ridge et al. [START_REF] Ridge | Capacity planning for intensive care unit[END_REF] build a simulation model for bed capacity planning in Intensive Care, based on the queuing model logic. Utley et al. [START_REF] Utley | Analytical methods for calculating the capacity required to operate an effective booked admissions policy for elective inpatient services[END_REF] develop methods to estimate the bed capacity required to minimize the number of cancellations of booked elective patients.

Other references consider more resources like Guinet and Chaabane [START_REF] Guinet | Operating theatre planning[END_REF] who design a model for minimizing the hospitalizations costs including the resource overloads. The objective is to assign patients to operating rooms so as to minimize costs while satisfying the equipment and staffing constraints. Adan and Vissers [START_REF] Adan | Patient mix optimisation in hospital admission planning: a case study[END_REF] develop a model to generate a planning and mix of patients that minimize the deviations between resources consumption and their targets. The number of beds in the medium-care unit, the operating theatre capacity as well as the nursing hours and the number of beds available in the intensive care unit are all critical resources in this model. In the same vein, Vissers, Adan and Dellaert [START_REF] Vissers | Developing a platform for comparison of hospital admission systems: An Illustration[END_REF] consider a similar problem with additional restrictions in planning some combinations of patients and availability of resources. In this contribution, lengths of stay both in the medium and intensive care units are now assumed to be stochastic. The three above-mentioned references focus on planning issue at a tactical level, i.e. the mix of patients to be admitted within a medium term horizon (2 to 4 weeks, for instance).

The present paper pursues the work of Vissers et al. [START_REF] Vissers | Developing a platform for comparison of hospital admission systems: An Illustration[END_REF] by developing several strategies to determine the best scheduling of individual patients on an operational level. These strategies lean upon the admission policy on a tactical level that results from the application of the model of Vissers et al. [START_REF] Vissers | Developing a platform for comparison of hospital admission systems: An Illustration[END_REF]. The focus here is on the operational level, with the aim of assessing the usefulness of the planning at the tactical level to build the scheduling strategies.

The tactical plan is based on average arrivals of patients whereas operational strategies have to treat actual arrivals that obviously deviate from average. To deal with these random deviations some flexibility is necessary that can be created in different ways, for instance by regularly recalculating the patient mix, or by replacing patients from one group by patients from another group, or by creating some safety slack by initially overestimating the number of patients. We will again use the setting of the Thorax Centre Rotterdam to compare and to test the different operational strategies developed here.

The remainder of the paper is organized as follows. Section 2 describes the case study setting and ends with the research questions to be answered. Section 3 provides the mathematical model used to obtain the planning of patients on a tactical level. Section 4 explains the strategies implemented to obtain operational schedules of patients. The performance criteria used to assess each strategy also are described here. Numerical experiments are presented and discussed. Section 5 draws some conclusions and formulates recommendations for further work.

Case study setting

The patient flow of the Thorax Centre Rotterdam consists of scheduled patients (elective patients from the waiting list) and emergency patients requiring immediate surgery. We only take into account elective patients for a classical reservation policy is adopted for emergency patients.

Patients are usually admitted to the Medium Care unit (MC) one day before operation unless they come from another department of the University Hospital. After the operation they stay for some days in an Intensive Care unit (IC) and after recovery they may stay in the Medium Care unit (MC) for a few days.

The current tactical planning in the Thorax Centre has a strong focus on the operating theatre (OT) capacity. However, this planning may be improved by also taking into account all other critical resources involved. To this purpose, a cyclical schedule has been developed on a tactical level based upon average numbers of patients for the various groups. The corresponding mathematical model to obtain such a tactical schedule will be presented in the next section.

Patient groups, volumes and demand requirements

Patients have been grouped in several categories, each of these being relatively homogeneous in terms of consumption of OT and IC resources. Table 1 provides information on the patient groups considered, the expected duration of the operation for each group, the average length of stay at the IC (outliers excluded), and the average number of patients per group to be operated within a 4-week horizon.

Available resources

For each of the resources, there exists a maximum available capacity per day and a target utilization level. This data is displayed in Table 4 for each day and each resource. Defining a level of utilization lower than 100% allows for dealing with emergencies and fluctuations in number of patients. The data apply to every week in the 4-week planning cycle. Four operating theatres are available nine hours per day. From the total of 36 hours of capacity available per day, 29 hours are allocated to electives while the rest is reserved for emergencies. On fridays, the target utilization is lower. The IC unit has 10 beds available throughout the working week and 4 beds during the weekend. The target utilization level for the IC is set to 7 beds for electives patients throughout the working week and to 2 beds during the weekend. The MC has 36 beds available every day and the target utilization by electives is 27 beds throughout the whole week. The available IC nursing staff and target utilization of the IC nursing workload (in number of hours per day) is matched with the number of IC beds. The targets for IC-beds, MC-beds and IC-nursing are defined at a lower level compared to the target for OT hours, to deal with fluctuations in the number of patients.

OT hours IC beds MC beds IC nursing hours

The IC nursing hours (NH) required per day for patients of categories 1 to 7 are estimated to 12 hours throughout their stay in the IC unit unless for the second day in IC and for categories 5, 6 and 7 for which the needs are 24 hours. For category 7 there is also a need of 24 hours on day 3. Patients from category 8 only require 3 hours of IC nursing per day whatever their length of stay in IC.

The importance of each resource has been assessed by the stakeholders in the hospital who decided to allocate to each of them a certain absolute weight denoted by

€ a r , with € r = ot,ic,nh,mc { }.
Relative weights will be employed to evaluate the cost of deviating from the target utilization level. The relative weight € α r for resource r is defined as

€ α r = a r R r, j j=1 T ∑ a r R r, j j=1 T ∑ r= ot,ic,nh,mc { } ∑ , (1)
Where € R r, j is the target resource utilization for resource r in period j. These values can be read in Table 4. For instance € R ot,5 = 25.

Research questions

In the tactical planning we make a reservation of the resources for a fixed number of patients from the various groups to be treated during a cycle. However, if the average number of arriving patients is equal to the maximum number that can be treated, the system may become instable in the end as the variance is obviously not zero. For the operational planning we will consider several options to create feasibility: overplanning, flexibility in patient groups and periodic updates of the tactical plan. Our incentive is to determine which combinations of these strategies (overplanning, patient flexibility and updating) lead to the best schedule of electives in terms of waiting times for patients, in terms of deviations between the resources consumption and the target utilization and in terms of schedule changes for the operating specialists.

Problem formulation and mathematical model for the tactical planning

In this section we formulate the tactical planning problem as a mixed integer linear program. Let N be the number of patient categories and T the length of the cyclic operation schedule (N=8 and T=28 in the Thorax Centre problem). On each day of the operation schedule we have to decide on the number and mix of patients to be operated. Hence, the decision variables are X c,t denoting the number of patients from category c operated on day t of the operation schedule, where c = 1, 2, …, N and t = 1, 2, …, T. The objective is to determine the variables X c,t satisfying certain constraints and for which the expected utilization of all resources matches the target as close as possible. Over-utilizations and underutilizations relative to the targets are identically penalized, as decided by the Thorax Centre. Hereafter we first formulate the constraints on the decision variables X c,t and then the objective function.

The total number of patients of group c to be operated over the T-day cycle should be equal to the target patient throughput P c which values are displayed in the penultimate column of Table 1 (planned number of patients). Hence

€ X c,t t =1 T ∑ = P c , c = 1,...,N. (2)
To describe the constraints for the utilization of the resources we introduce the parameters C r,t and R r,t indicating the available capacity and target utilization, respectively, for resource r on day t, where € r = ot,ic,mc,nh

{

}. These values are given in Table 4. Let the auxiliary variables U r,t and O r,t denote the under-and over-utilization (with respect to the target). Then we get for the utilization of operating theatre

€ R ot,t -U ot,t ≤ s c c=1 N ∑ X c,t ≤ R ot,t + O ot,t , t = 1,...,T , (3)
where s c denotes the operation time dedicated to a patient of category c (see

€ R ic,t -U ic,t ≤ p ic,c, j j =0 L ic max ∑ c=1 N ∑ X c,t-j ≤ R ic,t + O ic,t , t = 1,...,T . (4)
In the above constraints we used the convention that the subscript t-j in X c,t-j should be treated modulo T: day 0 is the same as day T, day -1 is the same as day T-1 and so on. If w c,t denotes the IC nursing load (in hours) required for a patient of category c t days after operation, then we get for the expected nursing workload

€ R nh,t -U nh,t ≤ w c,t p ic,c, j j =0 L ic max ∑ c=1 N ∑ X c,t -j ≤ R nh,t + O nh,t , t = 1,...,T . (5)
Similarly, for the expected utilization of the MC unit we get

€ R mc,t -U mc,t ≤ X c,t+ j j=1 l c ∑ c=1 N ∑ + p mc,c, j j=0 L mc max ∑ c=1 N ∑ X c,t-j ≤ R mc,t + O mc,t , t = 1,...,T , (6)
where l c is the number of pre-operative days at the MC for category c and p mc,c,t is the probability that a patient from category c is at the MC unit t days after operation, t = 0, 1, 2, …, € L mc max (see Table 3). Further, for each of the resources, the available capacity should not be exceeded, so

€ R r,t + O r,t ≤ C r,t , r = ot,ic,nh,mc { }, t = 1,...,T . (7)
In addition to the constraints for the utilization of the resources we have to take into account restrictions valid for specific days of the operation schedule. If the number of operations for certain categories of patients is prescribed and fixed on certain days, then the corresponding variables X c,t are simply upper bounded accordingly. If the number of operations for certain combinations of patient categories is limited, then we have to require that

€ X c,t c∈S ∑ ≤ B t , t = 1,…,T , (8)
where S is a subset of the patient categories and B t denotes the maximum number of patients from the categories S c ∈ that can be operated on day t of the operation schedule. The objective is to minimize the weighted sum of under-and over-utilization,

€ α r r∈Ω ∑ U r,t + O r,t ( ) t=1 T ∑ , (9)
where the relative weight € α r for resource r is defined as in equation (1). Our planning problem therefore consists in minimizing the objective function in (9) subject to constraints (2) to (8) and the integrality constraint

€ X c,t ∈ 0,1,2,… { } , c = 1,…,N, t = 1,…,T , (10)
The resultant mixed integer program is implemented in C language and linked with ILOG CPLEX 9.0 as a callable optimization library. Solving the deterministic version of the model to optimality did not take much computation time, but the stochastic model could not be solved to optimality. After 24 hours of computation time, there is still a 12 percent optimality gap.

The operational schedule

If the average number of arriving patients is close to the maximum number that can be treated, the system will become instable in the end. One solution is to create a security slack by reserving capacity for more patients than the average. This overplanning strategy consists in increasing the target throughput of patients on the basis of which the tactical plan is derived. As an additional measure, we can allow the replacement of patients from one group by patients from another group to avoid unused capacity in some days. An alternative is to make a new tactical plan regularly, either every 3 months or every year, and build the operational schedule on the basis of a number of required patients depending partly on the expected number of patients and partly on the waiting list of patients. These strategies are further described and illustrated in the next section.

The strategies

The strategies we will consider to build the operational schedule consist of different options for overplanning, patient flexibility and rescheduling. For the sake of illustration, let us consider 3 categories of patient (for instance, categories 3, 4 and 5 of the Thorax Centre problem) over a cycle of 7 days. Table 6 provides the optimal tactical plan corresponding to that problem. Overplanning. We consider two possibilities for the amount of planned patients in the tactical plan: P1. Plan the minimum amount of patients per group (no overplanning). P2. Plan the amount of patients per group in such a way that less than 5 percent of the patients have to wait more than one cycle. For determining this amount, we calculate the steady state probabilities in a simple queuing model.

Planned admission over a 7-day cycle

In the Thorax Centre problem, option (P1) simply consists in calculating the tactical plan on the basis of the initial P-value, with P={8,10,67,13,3,2,1,7} (see Table 1, penultimate column). Option (P2) results in the computation of the tactical plan using P={9,11,70,15,4,3,2,9}.

Flexibility. For a better operational use of the reserved capacities emanating from the tactical plan, we consider three options: F1. No flexibility. We follow the tactical plan unless the number of actual arriving patients is inferior to the planned number. In the latter case, the operation is cancelled. F2. Full flexibility. If some patients are planned on a certain day (i.e. the number of planned patients in the tactical plan is not zero), we use the reserved capacities to schedule the actual arriving patients with the longest waiting times. This means that some planned categories can be cancelled and replaced by others with a longest queue. F3. Partial flexibility. Amongst the categories that are planned, we replace each category with no actual arriving patient with another planned category for which there are actual arriving patients. The selected category for replacement is the one for which the product of the number of planned patients by the number of arrivals is maximum. This avoids the cancellation of some operations that would result from a lack of patients. Contrary to the previous flexibility option, a category that is not planned in the tactical plan can not be scheduled.

To illustrate the no-flexibility option (F1), let us consider the planned numbers of patients on day 3 in the tactical plan of Table 6. These planned numbers are {0,2,3} for the 3 categories. An actual arriving number of patients equal to {1,0,4} would lead to a stream of scheduled patients: {0,0,3}. The arriving patient of the first category is not considered as no patient of that category was planned on that day. No patient of the second category is scheduled although we planned to operate 2 of these patients, since there is no actual arriving patient of this category. Only 3 of the 4 arriving patients are scheduled as 3 patients are planned so we follow the tactical plan to avoid exceeding capacities.

To illustrate the full flexibility option (F2), let us consider the planned number in the tactical plan on day 1: {5,1,0}. While treating each category, the waiting time of each group is measured and the category with the highest waiting time is selected, leading to a stream of scheduled patients: {2,2,2}. Such a schedule is definitely different from the tactical plan that was notably suggesting not to operate patients of the third category.

Consider again the planned number of patients on the first day: {5,1,0} and assume the actual arriving patients are equal to {10,0,1}. Under the partial flexibility option (F3), the second category for which there is no actual patient is replaced with the first category but could not have been replaced with the third group as it is not in the tactical plan. Thus, (F3) leads to a scheduled stream of patients: {6,0,0} whereas the no flexibility option would have resulted in the schedule: {5,0,0}. The full flexibility option could produce the stream: {4,1,1}.

Updating. To update the tactical plan, we replace part of the expected number of patients with the actual number of patients in queue. Updated values of target throughputs for each category,

€ P c ' , are computed according to (11)
where € λ c is the average number of patients (see last column in Table 1) and € {P c 0 } c=1..N is the initial stream of target throughput values (overplanned values under option P2 or regular values under option P1). We thus compute a tactical plan on the basis of these new target throughput values for patients. We consider 3 frequencies of updating: U1. Compute a new tactical plan with updated values € P c ' every 3 cycles U2. Update the numbers and the tactical plan every year U3. Do not update

€ P c ' = round 11 12 λ c + 1 6 round(λ c ) + 1 3 (P c 0 -λ c )       , c = 1,...,N,
Combining the two options for the overplanning, the three options for flexibility and the 3 updating possibilities leads to 18 different strategies for which we will compute several performance criteria.

Performance criteria

The performance of an operational schedule will be assessed by the weighted sum of deviations as defined in Eq. (9). To compute the consumption of resources associated to such a schedule, we will use the actual length of stay of each patient, rather than using the expected length of stay we considered in the tactical plan.

The average waiting time for the patients will also be used as a performance indicator of the operational schedule. Under a FIFO assumption the average waiting time over a cycle is defined as the difference between the total number of arriving patients and the total number of scheduled patients divided by the number of periods in a cycle (cycle length). Note that in a hospital usually an additional component of waiting time has to be introduced to account for the period that patients have to be notified. This period is usually one week. In our performance measure we will not include this period.

We consider 4 indicators of the scheduling differences between the operational plan and the tactical plan. Recall that € X c,t designates the number of patients of category c in the tactical plan on day t. We let € Y c,t be the equivalent in the operational schedule. Indicator C measures the number of cancelled operations of any category over a cycle. If A straightforward indicator to measure the differences between two tactical plans if we use the updating options is the number of new days in which operations are programmed. If, in a previous plan, we had € X c,t = 0 and if we have € X c,t > 0 in the new plan, then the number of plan changes is increased by 1.

Numerical experiment

We ran the 18 combinations of strategies over a simulation horizon of 10 years. Arrival of patients was simulated through a Poisson process. For each scheduled patient, we drew an actual IC stay and MC stay according to the empirical distributions that were obtained from the hospital and described in Tables 2 and3. Each tactical plan was obtained by running Cplex and we stopped the optimization process after the 7 th best integer solution found so as to limit the computation time. Simulation results are displayed in Table 7. For each combination of strategies, Table 7 provides the associated values of the performance indicators as described in the previous subsection. In addition, we compute a global indicator of the plan changes which is the weighted sum of the 5 indicators of plan changes and the target deviations with the following weights {2,10,1,5,4,10} for IS, I, CS, the average number of plan changes and the target deviations respectively. The resultant index is called "global volatility" and appears in the penultimate column of Table 7 From Table 7, it appears that the average waiting time clearly improves with more flexibility. Under the full flexibility option, the waiting times are minimum, since this option precisely consists in scheduling in priority the categories with the highest waiting times. Introducing more flexibility however has a negative impact on the weighted target deviations. Without any flexibility the lowest target deviations values are obtained, as we stick as closely as possible to the tactical plan that seeks to minimize the deviations. An increasing flexibility is obviously associated to higher values for the plan changes indicators, as exhibited by the numbers in the columns 4-8 of Table 7.

Overplanning the number of patients leads to a clear reduction in the average waiting times, as the volume of planned patients is increased whereas the same number of patients still arrives. On a global basis, overplanning has little influence on the target deviations, unless under the partial flexibility option for which overplanning is detrimental to the deviations. The indicators of plan changes are roughly higher under overplanning, possibly because with more patients there also are more possibilities of canceling and adding unplanned patients.

The simulation results also shows that, at least in the stationary demand situation that we considered, updating the tactical plan does not contribute to a better performance and only makes the schedule of the operating specialists more uncertain as indicated by the values of the average number of plan changes (see the penultimate column in Table 7).

From these results, some dominant strategies appear in bold numbers in Table 7, with respect to the global volatility and the average waiting time. These strategies are plotted according to these 2 dimensions in Figure 1. The dominant strategies correspond to a limited flexibility and infrequent updating as most of them include no or partial flexibility and no update or a yearly one. The overplanning option is included in almost all the dominant strategies. By increasing the number of planned patients we can decrease the average waiting time in a more efficient way than by calculating a different tactical plan.

Conclusions and recommendations

In many organizations, the capacity planning is based upon standard durations for the different process phases: surgery, IC and MC stay. Vissers et al [START_REF] Vissers | Patient mix optimization in cardiothoracic surgery planning: a case study[END_REF] have considered such a deterministic approach to determine the optimal patient mix for a cardiothoracic surgery department. In this paper we have extended their model, by considering the operational planning for the operations theatre and the subsequent stay in the IC unit and in the MC unit. Based upon a big sample of patients of a Dutch cardiothoracic surgery department, we created an empirical distribution for the durations of the IC phase and the MC phase and used this in our mixed integer linear programming model, leading to a cyclic master operation schedule minimizing weighted deviations between realized and targeted resource use. This master schedule was used at the operational level to assign patients to operating slots. In order to make the operational plan feasible we considered options like overplanning, patient flexibility and rescheduling the tactical plan based upon the waiting list. The simulation results showed that updating the tactical plan was not very efficient. It also showed that deviating too much from the tactical plan leads to inefficient resource use. This enhances the need of a good tactical plan as a basis to build an effective operational schedule. It should also be noted that the improvement of the average waiting time for patients by implementing flexibility strategies also degrades the weighted target deviation.

€

  X c,t > 0 and if € Y c,t < X c,t , the number of cancellations of patients of category c on day t equals € X c,t -Y c,t . Indicator C sums up the daily cancellations if any, over a cycle. Indicator CS gives the number of cancelled categories over a cycle. If € X c,t > 0 and if € Y c,t = 0 then CS is increased by 1. Conversely, indicator I measures the increase number of patients for a category over the planned amount. If € X c,t > 0 and € Y c,t > X c,t then I is increased by € Y c,t -X c,t . At last, indicator IS provides the number of unplanned categories. If € X c,t = 0 and € Y c,t > 0 then IS is increased by € Y c,t .

Figure 1 .

 1 Figure 1. Weighted target deviation versus average waiting time

Patient group Example procedures OT duration (hours) IC-stay (days) Planned # patients Average # patients

	(1) Child simple	Closure ventricular septal defect	4	1.1	8	7.36
	(2) Child complex	Arterial switch	8	1.1	10	9.36
	(3) Adult, short OT, short IC	Coronary bypass (CABG)	4	1.3	67	66.00
	(4) Adult, long OT, short IC	Mitral valve plasty	8	1.5	13	12.73
	(5) Adult, short OT, middle IC	CBAG with expected medium IC stay	4	1.6	3	2.64
	(6) Adult, long OT, middle IC	Heart transplant	8	4.0	2	1.55
	(7) Adult, long OT, long IC	Thoraco-abdominal aneurysm, ELVAD	8	7.0	1	0.36
	(8) Adult, very short OT, no IC Cervical mediastinoscopy	2	0.2	7	6.91
	Table 1. Patient groups, operation duration, length of stay and 4-week numbers
	2.2 Length of stay
	We use a stochastic length of stay for IC and MC, based on empirical data of 2006
	including a sample of 576 patients. Table 2 provides information on the length of
	stay distribution at the IC for each patient group. For instance, any patient of
	category 8 has 79% chance to stay zero day in IC.

Probability of length of stay IC (days)

	Patient group	0	1	2	3	4	5	6	7	8	9	10
	(1) Child simple	0.07	0.87 0.02 0.02 0.02	0	0	0	0	0	0
	(2) Child complex	0	0.90 0.08 0.02	0	0	0	0	0	0	0
	(3) Adult, short OT, short IC	0.01	0.83 0.11 0.03 0.01 0.01	0	0	0	0	0
	(4) Adult, long OT, short IC	0	0.83 0.10 0.04	0	0.01 0.01	0	0	0	0.01
	(5) Adult, short OT, middle IC	0	0.79 0.07 0.07	0	0	0	0.07	0	0	0
	(6) Adult, long OT, middle IC	0	0	0.14 0.44 0.14 0.14	0	0	0.14	0	0
	(7) Adult, long OT, long IC	0	0	0	0	0	0	0	1	0	0	0
	(8) Adult, very short OT, no IC	0.79	0.21	0	0	0	0	0	0	0	0	0

Table 2 .

 2 Length of stay distribution at IC per patient group

Table 3

 3 displays the same information type for the MC.

	Probability

of length of stay post-op MC (days)

	Patient group	0	1	2	3	4	5	6	7	8	9	10	>10
	(1) Child simple	0.74	0	0	0	0.02 0.1	0.07 0.05 0.02	0	0	0
	(2) Child complex	0.83	0	0	0	0	0	0.04 0.04 0.02 0.02	0	0.05
	(3) Adult, short OT, short IC	0	0.01 0.01 0.04 0.32 0.24 0.12 0.09 0.05 0.03 0.04 0.05
	(4) Adult, long OT, short IC	0.03	0	0	0.01 0.12 0.16 0.18 0.15 0.10 0.04 0.04 0.17
	(5) Adult, short OT, middle IC	0	0	0	0	0.07 0.07 0.07 0.20	0	0.20 0.20 0.19
	(6) Adult, long OT, middle IC	0	0	0	0	0	0	0	0.14	0	0	0.14 0.72
	(7) Adult, long OT, long IC	0	0	0	0	0	0	0	0	0	0	1	0
	(8) Adult, very short OT, no IC	0.21	0.3 0.08 0.15 0.13 0.05	0	0.05	0	0.03	0	0

Table 3 .

 3 Length of stay distribution at MC per patient group

Table 4 .

 4 Available resources

		Capacity	Target	Capacity Target Capacity	Target	Capacity	Target
	Monday	36	29	10	7	36	27	133
	Tuesday	36	29	10	7	36	27	133
	Wednesday	36	29	10	7	36	27	133
	Thursday	36	29	10	7	36	27	133
	Friday	36	25	10	7	36	27	133
	Saturday	0	0	4	2	36	27	52
	Sunday	0	0	4	2	36	27	52

Table 5 .

 5 Table 5 provides the absolute and relative weights for each resource. Absolute and relative weights per resource

	Resource r	Absolute weight	a r	Relative weight	α r
	OT hours	8		0.167
	IC beds	10		0.756
	MC beds IC nursing	€ 3 5		0.047 € 0.029

Table 1 ,

 1 third column). To formulate the constraints for the expected utilization of the IC unit we introduce the notation p ic,c,t denoting the probability that a patient from category c is (still) at the IC unit t days after operation, t = 0, 1, 2, …,

	L ic max where
	€

Table 6 .

 6 Example of an optimal tactical plan

									Planned #	Average #
	Patient group	1	2	3	4	5	6	7	patients	patients
	(3) Adult, short OT, short IC	5	4	0	4	0	0	7	20	19.85
	(4) Adult, long OT, short IC	1	1	2	1	0	0	0	5	4.72
	(5) Adult, short OT, middle IC	0	1	3	0	0	0	0	4	3.64

Table 7 .

 7 . Performance of the 18 strategies

								Average	Weighted		Average
			Update	Indicators of plan change		# plan	target	Global	waiting
	Over-planning	Flex. Freq.	IS	I	CS	C	changes deviation	volatility time
	1. No Overplan. No	No	0	0	3.37	2.05	0	53.51	548.72	25.11
	2. No Overplan.	No	Year	0	0	4.44	2.34	22.78	54.25	649.76	18.73
	3. No Overplan.	No	Quarter	0	0	4.4	2.42	20.58	56.85	667.32	14.28
	4. No Overplan.	Partial No	0	6.1	4.01	5.37	0	54.77	639.56	9.66
	5. No Overplan.	Partial Year	0	4.19	5.06	2.98	21.22	54.43	691.04	12.43
	6. No Overplan.	Partial Quarter	0	4.97	4.31	3.14	19.61	54.63	694.45	11.66
	7. No Overplan.	Full	No	36.17	10.35	15.63 20.32	0	71.51	1008.17	6
	8. No Overplan.	Full	Year	36.84	10.25	17.2 19.15	19.78	72.03	1088.55	8.66
	9. No Overplan.	Full	Quarter	36.42	10.68	17.89 19.24	18.48	69.65	1064.15	8.81
	10. Overplan.	No	No	0	0	8.89	5.5	0	59.97	636.09	8.89
	11. Overplan.	No	Year	0	0	9.2	3.56	15.89	58.54	675.96	10.65
	12. Overplan.	No	Quarter	0	0	8.36	3.45	14.12	57.8	660.09	9.92
	13. Overplan.	Partial No	0	16.69	14.12 15.14	0	69.3	949.72	3.71
	14. Overplan.	Partial Year	0	13.04	10.96	9.01	16	60.65	856.91	4.96
	15. Overplan.	Partial Quarter	0	14.08	12.58 10.66	15.48	66.93	937.9	4.21
	16. Overplan.	Full	No	25.15	18.56	28.66 18.77	0	74.32	1101.61	2.53
	17. Overplan.	Full	Year	30.42	16.82	26.1 17.76	13.78	69.94	1098.46	3.02
	18. Overplan.	Full	Quarter	29.58	16.53	25.45 18.49	11.06	68.49	1071.5	3.02

The trade-off between these 2 dimensions must therefore be considered. Further future work is to look at the reservations policies for emergency admissions.


S Val Ère

Bitseki Penda

Angelina Roche
email: roche@ceremade.dauphine.fr

LOCAL BANDWIDTH SELECTION FOR KERNEL DENSITY ESTIMATION IN BIFURCATING MARKOV CHAIN MODEL

Keywords: Bandwidth selection, Bifurcating autoregressive process, Bifurcating Markov chains, Binary trees, Nonparametric kernel estimation. Mathematics Subject Classification (2010): 62G05, 62G10, 62G20, 60J80, 60F05, 60F10, 60J20, 92D25

. Drawing inspiration from dimension jump methods for model selection, we also provide an algorithm to select the best constant in the penalty.

Introduction

First introduced by [START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF] [START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF], Bifurcating Markov chains models (BMCM for short) has recently received particular attention for its application to cell lineage study. [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF] [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], have proposed such a model to detect cellular aging in Escherichia Coli and proved laws of large numbers and central limit theorem for this class of stochastic process. Bitseki-Penda et al. (2014) [START_REF] Bitseki Penda | Deviation inequalities, moderate deviations and some limit theorems for bifurcating Markov chains with application[END_REF] have then completed these asymptotic results and proved concentration inequalities.

To the best of our knowledge, kernel density estimation for the BMCM were considered first by Doumic & al. [17], where they estimate the division rate of population of cells reproducing by symmetric division, i.e. cell reproduction where each fission produces two equal daughter cells. After this work, Bitseki & al. [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF], have used the wavelets methodology to study the nonparametric estimation of the density of BMCM. They propose an adaptive estimator in dimension 1. Recently, Bitseki and Olivier [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] have studied the Nadaraya-Watson type estimators of a BMCM that they called nonlinear bifurcating autoregressive process. The latter model can be seen as an adaptation of nonlinear autoregressive process on binary regular tree. We mention that, except in [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF], all the estimations done in the previous works are non adaptive. In particular, the question of data-driven bandwidth selection was not addressed in [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] and [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF]. The main objective of this work is then to propose a data-driven method for choosing the bandwidth for the kernel estimator of the invariant measure in the multi-dimensional BMCM, following ideas from the works of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF].

The idea of the method is to select the bandwidth minimizing an empirical criterion imitating the bias-variance decomposition of the risk of the kernel estimator. More precisely, let H be a collection of bandwidths and let ( ν h ) h∈H be a family of kernel estimators of an unknown density ν. Then, we select the bandwidth h as h = arg min h∈H { A(h) + bV (h)} 1 with A(h) an empirical version of the bias of the estimator ν h and V (h) a penalty term with the same order than the variance. The now so-called Goldenshuger-Lepski methodology ( [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF], but also [START_REF] Goldenshluger | On adaptive minimax density estimation on R d[END_REF][START_REF] Rebelles | Lp adaptive estimation of an anisotropic density under independence hypothesis[END_REF][START_REF]Rebelles Pointwise adaptive estimation of a multivariate density under independence hypothesis[END_REF][START_REF] Patschkowski | Adaptation to lowest density regions with application to support recovery[END_REF][START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF]), initially developed for density estimation, has been applied in many contexts such as deconvolution problems [START_REF] Comte | Anisotropic adaptive kernel deconvolution[END_REF], conditional cumulative distribution function estimation [START_REF] Chagny | Adaptive and minimax estimation of the cumulative distribution function given a functional covariate[END_REF][START_REF] Chagny | Adaptive warped kernel estimators[END_REF], regression problems [START_REF] Chagny | Adaptive warped kernel estimators[END_REF][START_REF] Chagny | Adaptive estimation in the functional nonparametric regression model[END_REF], conditional density estimation [START_REF] Chagny | Adaptive warped kernel estimators[END_REF][START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF], hazard rate estimation [START_REF] Chagny | Adaptive warped kernel estimators[END_REF], white noise model [START_REF] Lepski | Adaptive estimation over anisotropic functional classes via oracle approach[END_REF], kernel empirical risk minimization (including robust regression) [START_REF] Chichignoud | Bandwidth selection in kernel empirical risk minimization via the gradient[END_REF], Lévy processes [START_REF] Bec | Adaptive pointwise estimation for pure jump Lévy processes[END_REF], Cox model [START_REF] Guilloux | Adaptive kernel estimation of the baseline function in the Cox model with high-dimensional covariates[END_REF], stochastic differential model [START_REF] Dion | Bidimensional random effect estimation in mixed stochastic differential model[END_REF]. Under suitable assumptions on the kernel, it is shown in [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] that this selection rule leads to a minimax adaptive estimator on a general class of regular functions, for a general class of L s -risks. Pointwise versions of the Goldenshluger-Lepski selection rule have been less considered [START_REF] Comte | Anisotropic adaptive kernel deconvolution[END_REF] [START_REF] Comte | Anisotropic adaptive kernel deconvolution[END_REF], [START_REF] Rebelles | Lp adaptive estimation of an anisotropic density under independence hypothesis[END_REF] [START_REF]Rebelles Pointwise adaptive estimation of a multivariate density under independence hypothesis[END_REF] and Chagny and Roche (2016) [START_REF] Chagny | Adaptive estimation in the functional nonparametric regression model[END_REF]). The interest of such an approach is that the bandwidth is selected with a local criterion which realizes the best bias-variance compromise at the point where the estimator is calculated. On the contrary, integrated versions of the Goldenshluger-Lepski selection rule select the same bandwidth at all points. Let us mention that, in all the articles cited above, the theoretical results rely on concentration inequalities for sums of i.i.d. random variables such as Bernstein Inequality or Talagrand Inequality. In our context, such results are not applicable. Hence, we prove a Bernstein-type Inequality for functionals of BMC, where the functions are kernels and convolution of kernels. Compare to those obtained in [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF], our inequalities are more complete in the sense that the deviation parameter can take all the positive values. More precisely, its values do not depend on the size of the samples, with is essential for our theoretical results.

Ideally, the penalty term V , called "the majorant" by Goldenshluger and Lepski, depends entirely on the kernel and the observations (it does not depend on the density ν). The selected estimator is then ν h . However, for the BMCM, the variance term V in the previous selection rule contains a term which may depend on the unknown density ν. Moreover, this term is generally not estimable from observations. To resolve this problem, we propose a modification of Goldenshluger-Lepski rule's selection, inspired by the works of [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF] [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF]. As suggested in [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF], the constant term in bV (x, h) is then selected automatically from the data with an algorithm inspired by the works of [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF].

The paper is organized as follows. The model is defined in Section 2. Section 3 is devoted to the definition of the estimator. In Section 4, we provide a numerical study of our estimator. The proofs are given in Section 5.

Definitions

We are now going to give a precise definition of a BMC. First we note that this class of stochastic processes has been introduced by Guyon [21] in order to understand the mechanisms of cell division. Indeed, these stochastic processes are well adapted to study a population (or more generally, any dynamic system) where each individual (or more generally, each particle) in one generation gives birth to two individuals in the next one. In the sequel, we will then use the language of the population dynamic to define the sets of interest.

Let (Ω, F, (F m , m ∈ N), P) be a filtered probability space. Let (X u , u ∈ T) be a sequence of random variables defined on (Ω, P), taking values in R d , where d ≥ 1, and indexed by the infinite binary tree T = ∞ m=0 {0, 1} m , with the convention that {0, 1} 0 = ∅. We equip R d with its usual Borel σ-field. Now, we will see T as a given population. Then each individual u of this population is represented by a sequence of 0's and 1's, and has two descendants, u0 and u1. The initial individual of the population is ∅. For all m ∈ N, let G m be the set of individuals belonging to the m-th generation, and T m the set of individuals belonging to the m first generations. We have:

G m = {0, 1} m , T m = m q=0 G q and T = m≥0 G m .
For an individual u ∈ G m , we set |u| := m its length (i.e. the generation to which it belongs).

2.1. Bifurcating Markov chain [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF].

Definition 1 (T-transition probability). Let

P : R d × B((R d ) 2 ) → [0, 1] (with B(R d ) 2 = B(R d ) ⊗ B(R d ) the usual product σ-field on (R d ) 2 ). Then P is a T-transition probability if • x → P(x, A) is measurable for all A ∈ B(R d ) 2 . • A → P(x, A) is a probability measure on ((R d ) 2 , B(R d ) 2 ) for all x ∈ R d . For a B(R d ) 3 -measurable function f : (R d ) 3 → R, we denote (when it is defined) by Pf the B(R d )-measurable function x ∈ R d → (R d ) 2
f (x, y, z)P(x, dy, dz).

Definition 2 (Bifurcating Markov chain). Let µ be a probability measure on (R d , B(R d )) and P a T-transition probability. We say that (X u , u ∈ T) is a (F m )-bifurcating Markov chain with initial distribution µ and T-transition probability P (denoted T-BMC in the sequel) if

• X u is F m measurable for all u ∈ G m . • X ∅ has distribution µ.
• For all m ∈ N, and for all family

(f u , u ∈ G m ) of B(R d )-measurable functions from (R d ) 3 to R, E u∈Gm f u (X u , X u0 , X u1 )|F m = u∈Gm Pf u (X u ),
where u0 := (u, 0) ∈ G m+1 and u1 := (u, 1) ∈ G m+1 .

2.2. Tagged-branched chain [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF][START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF]. Let (X u , u ∈ T) be a T-BMC with initial distribution µ and T-transition probability P. We denote by P 0 and P 1 respectively the first and second marginals of P. More precisely

P 0 (x, B) = P(x, B × R d ) and P 1 (x, B) = P(x, R d × B), for all x ∈ R d and all B ∈ B(R d ).
Let Q be the mixture of P 0 and P 1 with equal weights

Q = 1 2 P 0 + 1 2 P 1 .
The Markov chain Y := (Y m ) m∈N on R d with initial value Y 0 := X ∅ and transition probability Q is called the tagged-branch chain.

In all the paper, we will denote by Q m the mth iterated of Q recursively defined by the formulas

Q 0 (x, •) = δ x and Q m+1 (x, B) = S Q(x, dy)Q m (y, B) ∀B ∈ B(R d ).
It is well known that Q m is a transition probability in (R d , B(R d )). In particular, we have

E[f (Y m )] = µQ m f for all measurable function f : R d → R.
In the sequel, we will assume that the Markov chain Y is ergodic-that is say, that there exists a unique distribution ν on (R d , B(R d ) such that, for all measurable function f : R d → R,

lim m→∞ E[f (Y m )] = R d f dν.
We will also assume that the distribution ν has a density, that we also denote by ν, with respect to the Lebesgue measure.

As previous works have shown (see for example [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF][START_REF] Bitseki Penda | Deviation inequalities, moderate deviations and some limit theorems for bifurcating Markov chains with application[END_REF]), the analysis of a BMC (X u , u ∈ T) is strongly related to the asymptotic behavior of the tagged-branched chain (Y m , m ∈ N), and therefore to the knowledge of the invariant distribution ν. We stress that this distribution is unknown and it is not directly observable, in such a way that its estimation from the data is of great interest. The aim is to estimate ν from the observation of a subpopulation (X u , u ∈ T n ).

3. Estimation of the stationary distribution ν 3.1. Definition of the estimator. We suppose that we observe the process (X u , u ∈ T) up to the n-th generation. We denote by |T n | = 2 n+1 -1 the cardinality of T n . Based on the observation of (X u , u ∈ T n ), we propose the following estimator of ν

(1) ν h (x) = 1 |T n | u∈Tn K h (x -X u ) ∀x ∈ R d , where for all x = (x 1 , ..., x d ) t ∈ R d and for all h = (h 1 , ..., h d ) t ∈]0, +∞[ d K h (x) = K(x 1 /h 1 , ..., x d /h d )/ d j=1 h j
and K is a kernel, that is to say a function K : R d → R which verifies R d K(u)du = 1. The vector h is usually called a bandwidth. In the sequel, |h| will denote the product

h 1 ו • •×h d .
It is well known in the kernel estimation theory that the choice of h is of great interest. Indeed the amount of smoothing is controlled by a judicious choice of the bandwidth. Now, in order to tackle the issue of the choice of this bandwidth, we will need the following assumptions.

Assumption 3 (Assumptions on the kernel).

K 1 = R d |K(t)|dt < +∞, K 2 2 = R d |K(t)| 2 dt < +∞ and K ∞ = sup t∈R d |K(t)| < +∞.

Assumption 4 (Uniform geometric ergodicity condition).

There exists two constants ρ ∈ (0, 1/2) and M > 0 such that for all bounded ν-integrable function g, for all x ∈ R d and m ≥ 0

|Q m g(x) - R d gdν| ≤ M g ∞ ρ m .
This assumption is verified for severals models of BMC. We refer for example to [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] where Bitseki and Olivier have shown it for NBAR processes (see [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] Lemma 20). For a precise definition of NBAR processes, we refer to section 4.

Assumption 5 (Assumption on P, P 0 , P 1 , Q and ν). We assume that the transitions P, P 0 , P 1 and Q admit densities with respect to the Lebesgue measure that we denote with the same notations. Moreover, we assume that

P ∞ < +∞, P 0 ∞ < +∞, P 1 ∞ < +∞, ν ∞ < +∞, and Q ∞ < +∞.
3.2. Bias-variance decomposition. We consider a pointwise quadratic risk

E[( ν h (x) -ν(x)) 2 ],
where x ∈ R d is a given point.

The results of [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] (see the proof of Proposition 21) allows to obtain the following upper-bound on the risk, where we have make explicit the constants that appear. Proposition 6. Under Assumption 3 to Assumption 5, we have

(2) E[( ν h (x) -ν(x)) 2 ] ≤ 2(K h * ν(x) -ν(x)) 2 + 2 C(P, ν) |T n ||h| , where * denotes the convolution product f * g(x) = R d f (x-t)g(t)
dt for all functions f, g integrable over R d and

C(P, ν) = C I ( √ 2 -1) 2
with

C I = (1 + 1 1 -2ρ 2 )( Q ∞ + ν ∞ ) 2 + M 2 + C P and C P = 2 K 2 2 ( Q ∞ + ν ∞ ) + P ∞ + ν ∞ ( K 2 2 + P 0 ∞ + P 1 ∞ )
Inequality ( 2) can be seen as a bias-variance decomposition in the sense that

E ν [ ν h (x)] = 1 |T n | u∈Tn E ν [K h (x -X u )] = R d K h (x -t)dν(t) = K h * ν(x),
where E ν is the expectation with respect to the measure ν and C(P, µ)/|T n ||h| is an upper-bound on the variance term E ( ν h (x) -K h * ν(x)) 2 .

Selection rule. Given a family of bandwidths

H n ⊂ [0, +∞[ d ,
we have a family of estimators ( ν h ) h∈Hn and the aim is to select an estimator in this family with risk close to the unknown oracle risk

E ( ν h * (x) -ν(x)) 2 with h * = arg min h∈Hn E ( ν h (x) -ν(x)) 2 .
Imitating the decomposition given in Equation (2) we could consider the following selection rule inspired by the work of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] (3)

h = arg min h∈Hn {A(x, h) + bV (x, h)}
where

• H n ⊂ [0, +∞[ d is a finite collection of bandwidths; • A(x, h) = max h ∈Hn ( ν h (x) -K h * ν h (x)) 2 -aV (x, h ) + with b ≥ a ≥ 1; • V (x, h) = C(P, µ) log(|T n |)/|T n ||h|.
The constants a and b can be different, as suggested in Section 5 of [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF]. In Section 3.4, we provide a method to choose both a and b, as well as the quantity C(P, µ) appearing in the variance term V (x, h). The lower bound a = 1 is a critical minimal value, in the sense that the procedure fails for lower values of a (see [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF]).

We prove the following oracle-inequality on the selected estimator ν h .

E ( ν h (x) -ν(x)) 2 ≤ C 1 min h∈Hn B h (x) + log(|T n |) |T n ||h| + C 2 |T n | , (4)
where C 1 , C 2 > 0 do not depend on n nor x and

B h (x) = max h ∈Hn (K h * ν(x) -K h * K h * ν(x))
= (x 1 , ..., x d ) ∈ D, x j ∈ R such that (x 1 , ..., x j-1 , x j , x j+1 , ..., x d ) ∈ D, ∂ βj f ∂x βj j (x 1 , , ..., x j-1 , x j , x j+1 , ..., x d ) - ∂ βj f ∂x βj j (x) ≤ L|x j -x j | βj -βj .
Moreover, suppose that K is a kernel of order ∈ N (with ≥ max j=1,...,d {β j }) that is to say, R K(t)dt = 1, R x j K(x) = 0 for all j = 1, ..., and that

R |x| |K(u)| < +∞.
Suppose also that, for all n, there exists h * = (h where β = d/(1/β 1 + ... + 1/β d ) is the harmonic mean of β. Then, there exists a constant C > 0 such that

sup x∈D E ( ν h (x) -ν(x)) 2 ≤ C |T n | log(|T n |) -2 β/(2 β+1) .
This corollary is a direct consequence of Theorem 7 and [28] Proposition 2.1. Remark that the assumption on the bandwidth collection is verified e.g. by

H n = h max k -α , k = 1, ...., (|T n |h max / log(|T n |)) 1/α d
for all constant h max > 0 and all α > 1.

3.4. Estimation of the constant appearing in the variance. Due to the term C(P, µ) which is hardly estimable, the variance term V (x, h) is not calculable in practice. We propose an algorithm to estimate it, based on slope estimation, as developed in [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF][START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF].

As suggested in [START_REF] Lacour | Minimal penalty for the Goldenshluger-Lepski method[END_REF], we take b = 2a, hence in order to calculate the estimator, it is sufficient to have a good value of κ = aC(P, µ). The following algorithm is inspired by the procedure described in [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF], for model selection purposes. However, note that, in our case, both terms A(x, h) and V (x, h) depends on the constant κ, which is not the case in model selection contexts where only the penalty term depends on the constant. The selection of the grid of κ is then different here. 2 and κ 1 = 0.

κ m = |T n | log(|T n |) max h,h ∈Hn |h |( ν h (x) -K h * ν h (x))
2. While s ≤ s max (i) Generate a sequence (κ j ) 1≤j≤m such that and

κ j = κ 1 + j-1 m-1 (κ m -κ 1 )
, for all j = 1, ..., m. (ii) Calculate h j := h(κ j ) as the minimizer of the criterion (3) with a = κ/C(P, µ) and b = 2a. (iii) Set

j jump = arg max j=1,...,m-1 1 | h j | - 1 | h j+1 | , κ 1 = κ jjump , κ 1 = κ jjump+1 and s = s + 1. 3. Return h jjump+1 .
The algorithm search the value of κ for which the variance of the estimator increases significantly and select a slightly larger value. This allows to select an estimator with a reasonable variance. The same reasoning has given rise to the so-called dimension jump method for model selection purposes (see [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF]). The chosen value for the initialization of κ m comes from the following reasoning. Setting κ ≥ κ m as suggested: by definition, for all h, h ∈ H n ,

|T n |h log(|T n |) ( ν h (x) -K h * ν h (x)) 2 ≤ κ
which implies that A(x, h) = 0 for all h ∈ H n and that the criterion (3) will select the smaller bandwidth in H n . On the contrary, if κ > κ m , let h, h ∈ H n for which the minimum is attained in (5), we have

|T n | log(|T n |h ) ( ν h (x) -K h * ν h (x)) 2 > κ
which implies that A(x, h) > 0 and the criterion may select a bandwidth which is not the smaller one. Hence, the values of κ for which the criterion (3) may select suitable values of the bandwidth can not be greater than κ m . That is the reason why we consider an initial grid in the interval [0, κ m ]. However, the initial value κ m may be very large compared to the optimal value of κ. The loop 2. allows to search among small values of κ while avoiding the choice of a too large m which could be very expensive in terms of computation time. In practice s max = 2 and m = 20 seems to be a reasonable choice.

Simulations

We shall now illustrate the previous results on simulated data coming from various models of bifurcating Markov chains.

Bifurcating autoregressive processes.

Bifurcating autoregressive processes (BAR, for short) were first introduced by Cowan and Staudte [START_REF] Cowan | The bifurcating autoregressive model in cell lineage studies[END_REF] in order to study the data from cell division, where each individual in one generation gives birth to two children in the next generation. This model has been widely studied over the last thirty years (see for e.g. [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] and references therein). Recently, Bitseki and Olivier in [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] have proposed an extension of BAR process initially introduced by Cowan and Staudte. Their model is defined as follows.

Let X u ∈ R be a quantitative data associated to the cell u ∈ T, for example the growth rate of E. Coli. Then the quantities X u0 and X u1 associated to u0 and u1 the two children of u are linked to X u through the following autoregressive equations ( 6)

L(X ∅ ) = µ, and for u ∈ T,

   X u0 = f 0 (X u ) + ε u0 , X u1 = f 1 (X u ) + ε u1 ,
where µ is a distribution probability on R and f 0 , f 1 : R ; R. The noise (ε u0 , ε u1 ), u ∈ T forms a sequence of independent and identically distributed bivariate centered random variables with common density g ε on R × R. The process (X u , u ∈ T) defined by ( 6) is a bifurcating Markov chain with T-transition probability

P(x, dy, dz) = g ε (y -f 0 (x), z -f 1 (x))dydz.
Under some assumptions on µ, f 0 , f 1 and g ε , it has been shown in [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] that the process X satisfies all the good properties needed for our theoretical results (we refer to [START_REF] Bitseki Penda | Autoregressive Functions Estimation in Nonlinear Bifurcating Autoregressive Models[END_REF] for more details). We note that the previous model can be seen as an adaptation of nonlinear autoregressive model when the data have a binary tree structure. Furthermore, the original BAR process in [START_REF] Cowan | The bifurcating autoregressive model in cell lineage studies[END_REF] is defined for linear link functions f 0 and f 1 with f 0 = f 1 . Now, for our numerical illustrations, we build a BAR process living in S := [0, 1] as follows. First, we choose X ∅ such that L(X ∅ ) = Beta(2, 2), where Beta(2, 2) is the standard Beta distribution with shape parameters [START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF][START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF]. Then for u ∈ T and conditionally on X u = x, we construct X u0 and X u1 independently in such a way that P(X u0 ∈ dy, X u1 ∈ dz) = P(x, y, z)dydz, where P(x, •, •) := P(x, •) ⊗ P(x, •) and

P(x, y) := (1 -x) y(1 -y) 2 B(2, 3) + x y 2 (1 -y) B(3, 2) , x, y ∈ [0, 1]
with B(α, β) the normalizing constant of a standard Beta distribution with shape parameters α and β. Now, one can prove that this process is stationary, it has an explicit invariant density, which is crucial to evaluate the quality of estimation of our method: this is a standard Beta distribution with shape parameters [START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF][START_REF] Basawa | Non-Gaussian bifurcating models and quasi-likelihood estimation[END_REF]. One can also prove that

E [X u0 |X u ] = E [X u1 |X u ] = 1/5X u + 2/5,
in such a way that the equations ( 6) are satisfied with f 0 (x) = f 1 (x) = 1/5x + 2/5 (for more details, we refer for e.g. to [START_REF] Pitt | Constructing First Order Stationary Autoregressive Models via Latent Processes[END_REF]). Now, it is no hard to verify that this process satisfies our required assumptions.

We simulate the n first generations of the process X, with n = 10 (hence the size of T n is |T n | = 2 10 = 1024). We consider the Gaussian kernel K(t) = 1 √ 2π e -t 2 /2 . The results are given in Figure 1. We are now interested in growth-fragmentation models. Theses models describes (for e.g.) the evolution of cells which grow and divide randomly over time (see for e.g. [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] and references therein). The model we are going to study is a simplification of the one studied in [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF]; it is defined as follows. Let S be a subset of [0, ∞) and let B : S → [0, ∞) be a continuous function (the splitting rate). Each cell u ∈ T grows exponentially with a common rate τ > 0 and when it reaches a certain size x, it splits at rate B(x), and gives birth to two offspring (u0 and u1) of size x/2. Next, this two offspring, u0 and u1, start a new life independently of each other. Clearly, the process (X u , u ∈ T), where X u is the size of the cell u at birth, is a BMC. It is proved in [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] that the T-transition probability of this BMC is P(x, •, •) := P(x, •) ⊗ P(x, •), where the density of P(x, •) is given by

P(x, y) := B(2y) τ y exp - y x/2 B(2z) τ z dz 1 {y≥x/2} ,
for x ∈ S and y ∈ S/2. It can also be seen that the probability transition of the tagged-branch chain is Q = P. Doumic et al. [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] have proved that P admits an invariant probability measure ν having a density, that we still denote by ν(•), with respect to the Lebesgue measure. It is also known from [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] that the rate function B(•) and the invariant density ν(•) verify

B(x) = τ x 2 ν B (x/2)
x x/2 ν B (z)dz in such a way that a natural estimator for B(x), based on the observation of (X u , u ∈ T n ) is

B n (x) = τ x 2 ν h (x/2) 1 |Tn| u∈Tn 1 {x/2≤Xu<x} ∨ n
, where ν h is the kernel estimator defined in (1) and n is a threshold which ideally go to 0. Moreover, Bitseki et al. [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF] have proved that under suitable assumptions on the splitting rate B(•), the process X satisfies all the good properties needed for our theoretical results. For all the previous assertions, we refer to [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF][START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF] for more details. The strategy is then to use our results for the estimation of ν B (x/2) for all x ∈ S. Now, for our numerical illustrations, we will work with the splitting rate using in [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF]. We choose τ = 2, S = (0, 5) and for all x ∈ S, B has the form

B(x) = x x -5 + 3T 2(x - 7 2 )
where

T (x) = (1 + x)1 {-1≤x<0} + (1 -x)1 {0≤x≤1} is a tent shaped function.
With this choice of B, the required assumptions for our theoretical results are satisfied (see [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF]). We simulate the n first generations of the process X, with n = 15. The results are given in Figure 2.

Proofs

The proof rely on the lemma below, which is a Bernstein-type inequality.

Lemma 12. Let (X u , u ∈ T n ) be a bifurcating Markov chain on R d with initial distribution µ and T-transition probability P . Under the assumption of uniform geometric ergodicity, we have for all δ > 0

(7) P 1 |T n | u∈Tn (K h * K h (x -X u ) -E ν [K h * K h (x -X u )]) > δ ≤ 2 exp δc K,Q,ν,M c ρ 4c ρ,M K 1 K ∞ δ 3 + c K,Q,ν,M c ρ exp   - δ 2 |T n ||h | 2 c K,Q,ν,M c ρ + 4c ρ,M K 1 K ∞δ 3   ,
where

c ρ,M = M (1 + ρ) 1 -2ρ , c ρ = 3 + 2 1 -2ρ , c K,Q,ν,M = 8 max{2 K 2 1 K 2 2 ( Q ∞,∞ + ν ∞ ); max{ Q ∞,∞ + ν ∞ ; M K 1 K ∞ } 2 } We also have for all δ > 0, (8) P 1 |T n | u∈Tn (K h (x -X u ) -E ν [K h (x -X u )]) > δ ≤ 2 exp δc K,Q,ν,M c ρ 4c ρ,M K ∞δ 3 + c K,Q,ν,M c ρ exp   - δ 2 |T n ||h| 2 c K,Q,ν,M c ρ + 4c ρ,M K ∞δ 3   where c K,Q,ν,M = 8 max{M K ∞ ; ( Q ∞,∞ + ν ∞ ) K 1 ; ( Q ∞,∞ + ν ∞ ) K 2 2 }.
Remark 13. As mentioned above, these inequalities are more complete than those obtained in [START_REF] Bitseki Penda | Adaptive estimation for bifurcating Markov chains[END_REF], since the deviation parameter δ does not depend on the size of the data. We stress that this fact is essential for our theoretical results.

Proof. We will do the proof of [START_REF] Bitseki Penda | Deviation inequalities, moderate deviations and some limit theorems for bifurcating Markov chains with application[END_REF]. The proof of (8) follows the same lines. Let λ > 0 and δ > 0. By Chernoff inequality, we have

P 1 |T n | u∈Tn (K h * K h (x -X u ) -E ν [K h * K h (x -X u )]) > δ ≤ exp (-δλ|T n |) E exp λ u∈Tn g(X u ) ,
where the function g is defined by

g(y) = K h * K h (x -y) -E ν [K h * K h (x -X ∅ )] .
For all u ∈ G n-1 , we have, on the one hand

|g(X u0 ) + g(X u1 ) -2Qg(X u )| ≤ 2M (1 + ρ) K h * K h ∞ .
Using the Young's inequality, we have

K h * K h ∞ ≤ K 1 K ∞ /|h | and therefore |g(X u0 ) + g(X u1 ) -2Qg(X u ))| ≤ 2c ρ,M K 1 K ∞ /|h |.
On the other hand, we have

E (g(X u0 ) + g(X u1 ) -2Qg(X u )) 2 |F n-1 ≤ 4 Q ∞,∞ K h * K h 2 2 ≤ c K,Q,ν,M /|h |). Now for all λ ∈ (0, 3/(2c ρ,M K 1 K ∞ (|h |) -1
)), we have from Bennett inequality

E [exp (λ (g (X u0 ) + g (X u1 ) -2Qg (X u ))) |F n-1 ] ≤ exp   c K,Q,ν,M (|h |) -1 λ 2 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3  
and using the Markov property, this leads us to

E exp λ u∈Tn g(X u ) ≤ exp   c K,Q,ν,M (|h |) -1 λ 2 |G n-1 | 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3   × E   exp   λ u∈Tn-3 g(X u )   × exp   λ u∈Gn-2 I (2) (X u )   × u∈Gn-2 E exp J (2) (X u , X u0 , X u1 ) F n-2  
where

I (2) (X u ) = (g + 2Qg + 2 2 Q 2 g)(X u ) J (2) (X u , X u0 , X u1 ) = (g + 2Qg)(X u0 ) + (g + 2Qg)(X u1 ) -2(Qg + 2 2 Q 2 g)(X u ).
Now for the second step, we will do the same thing with J (2) (X u , X u0 , X u1 ) instead of g(X u0 ) + g(X u1 ) -2Qg(X u ). For all u ∈ G n-2 we have

|J (2) (X u , X u0 , X u1 )| ≤ 2c ρ,M K 1 K ∞ (|h |) -1 .
We also have E[(J (2) (X u , X u0 , X u1 )) 2 |F n-2 ] ≤ 4Q((g + 2Qg) 2 )(X u ). Now we will control the right hand side of the previous inequality. On the one hand, we have for all y ∈ R d ,

Qg(y) = K h (t) K h (x -z -t) Q(y, z)dtdz -E ν [K h * K h (x -X ∅ ].
Using the change of variables t = hu and z = x -hu -h v, where for x = (x 1 , ..., x d ) t , y = (y 1 , ..., y d ) t ∈ R d , xy denotes the vector (x 1 y 1 , ..., x d y d ) t , we obtain

K h (t) K h (x -z -t) Q(y, z)dtdz| ≤ Q ∞,∞ .
In the same way we prove that

|E ν [K h * K h (x -X ∅ ]| ≤ ν ∞ . This leads us to |Qg(y)| ≤ Q ∞,∞ + ν ∞ and therefore, for all m ≥ 1 |Q m g(y)| ≤ Q ∞,∞ + ν ∞ .
On the other hand, uniform geometric ergodicity assumption implies that for all m ≥ 1 and for all y ∈ R d ,

Q m g(y) ≤ M ρ m K h * K h ≤ M K 1 K ∞ ρ m (|h |) -1 ,
and therefore, for all m ≥ 1 and for all y ∈ R d , we have

Q m g(y) ≤ max{ Q ∞,∞ + ν ∞ , M K 1 K ∞ } × inf{1, ρ m (|h |) -1 }.
For all y ∈ R d , we also have

Qg 2 (y) ≤ 2Q(K h * K h ) + 2(E[K h * K h (x -X ∅ )]) 2 ≤ 2 K 2 1 K 2 2 ( Q ∞,∞ + ν ∞ ), in such a way that for all u ∈ G n-2 , 4Q((g + 2Qg) 2 )(X u ) ≤ 8Qg 2 (X u ) + 8Q((2Qg) 2 )(X u ) ≤ c K,Q,ν,M ((|h |) -1 + 2 inf{1, ρ(|h |) -1 } 2 )
and thus, we obtain

E[J (2) (X u , X u0 , X u1 )|F n-2 ] ≤ c K,Q,ν,M ((|h |) -1 + 2 inf{1, ρ(|h |) -1 } 2 ). Once again, for all λ ∈ (0, 3/(2c ρ,M K 1 K ∞ (|h |) -1
)) and for all u ∈ G n-2 , we have from Bennett's inequality

E[J (2) (X u , X u0 , X u1 )|F n-2 ] ≤ exp   c K,Q,ν,M ((|h |) -1 + 2 inf{1, ρ(|h |) -1 } 2 )λ 2 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3   .
It then follows that

E exp λ u∈Tn g(X u ) ≤ exp   c K,Q,ν,M (|h |) -1 λ 2 |G n-1 | 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3   × exp   c K,Q,ν,M ((|h |) -1 + 2 inf{1, ρ(|h |) -1 } 2 )λ 2 |G n-2 | 2 1 - 2c ρ,M K 1 K ∞ (|h |) -1 λ 3   × E   exp   λ u∈Tn-4 g(X u )   × exp   λ u∈Gn-3 I (3) (X u )   × u∈Gn-3 E exp J (3) (X u , X u0 , X u1 ) F n-2   where
I (3) (X u ) = (g + 2Qg + 2 2 Q 2 g + 2 3 Q 3 g)(X u ); J (3) (X u , X u0 , X u1 ) = (g + 2Qg + 2 2 Q 2 g)(X u0 ) + (g + 2Qg + 2 2 Q 2 g)(X u1 ) -2(Qg + 2Q 2 g + 2 2 Q 3 g)(X u ).
Now, iterating this method, we are led to

E exp λ u∈Tn g(X u ) ≤ exp       n m=1 c K,Q,ν,M (|h |) -1 + m-1 l=1 2 l ((|h |) -1 ρ l ∧ 1) 2 |G n-m| λ 2 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3       × exp c K,Q,ν,M (|h |) -1 + n m=1 2 m 1 ∧ ρ m (|h |) -1 . Set m * = log |h |/ log ρ . Then we have n m=1   (|h |) -1 + m-1 l=1 2 l ((|h |) -1 ρ l ∧ 1) 2   |G n-m| = n m=m * +1   (|h |) -1 + m * l=1 2 l + m-1 l=m * 2 l h -1 ρ l 2   |G n-m| + m * m=1   (|h |) -1 + m-1 l=1 2 l 2   |G n-m| ≤ 6 + (1 + 1 1 -2ρ ) 2 (|h |) -1 |T n |.
We also have

(|h |) -1 + n m=1 2 m 1 ∧ ρ m (|h |) -1 ≤ c ρ (|h |) -1 .
In view of the above, for all λ ∈ (0, 3/(2c ρ,M K 1 K ∞ (|h |) -1 )) we have

P 1 |T n | u∈Tn (K h * K h (x -X u ) -E ν [K h * K h (x -X u )]) > δ ≤ exp   -λδ|T n | + c K,Q,ν,M c ρ (|h |) -1 |T n |λ 2 2 1 - 2c ρ,M K 1 K ∞(|h |) -1 λ 3   × exp λc K,Q,ν,M c ρ (|h |) -1 . Taking λ = (δ(|h |))/(c K,Q,ν,M c ρ + (4c ρ,M K 1 K ∞ δ)/3), we obtain P 1 |T n | u∈Tn (K h * K h (x -X u ) -E ν [K h * K h (x -X u )]) > δ ≤ exp δc K,Q,ν,M c ρ 4c ρ,M K 1 K ∞δ 3 + c K,Q,ν,M c ρ exp   - δ 2 |T n |(|h |) 2 c K,Q,ν,M c ρ + 4c ρ,M K 1 K ∞δ 3   .
The result follows since we can do the same thing for -g instead of g. Now, the proof of (8) follows the same lines and this ends the proof.

Proof of Theorem 7. We start from the following decomposition, true for all h ∈ H n

( ν h (x) -ν(x)) 2 ≤ 3( ν h (x) -K h * ν h (x)) 2 + 3(K h * ν h (x) -ν h (x)) 2 + 3( ν h (x) -ν(x)) 2 . Hence, since K h * ν h (x) = K h * ν h (x)
, by definition of A(x, h) and then by definition of h and the fact that a ≤ b

( ν h (x) -ν(x)) 2 ≤ 3(( ν h (x) -K h * ν h (x)) 2 -aV (x, h) + aV (x, h)) h * ν h (x) -ν h (x)) 2 -aV (x, h) + aV (x, h)) + 3( ν h (x) -ν(x)) 2 ≤ 3(A(x, h) + bV (x, h)) +3(A(x, h) + bV (x, h)) + 3( ν h (x) -ν(x)) 2 ≤ 6(A(x, h) + bV (x, h)) + 3( ν h (x) -ν(x)) 2 . (9) Now it remains to upper-bound E[A(x, h)]. We have ( ν h (x) -K h * ν h (x)) 2 ≤ 3( ν h (x) -K h * ν(x)) 2 + 3(K h * ν h (x) -K h * K h * ν(x)) 2 +3(K h * ν(x) -K h * K h * ν(x)) 2 .
With a rough upper-bound of the max h∈Hn by the h∈Hn we get

E[A(x, h)] ≤ 3 E max h ∈Hn ( ν h (x) -K h * ν(x)) 2 -a V (x, h ) 6 + + 3 E max h ∈Hn K h * ν h (x) -K h * K h * ν(x)) 2 -a V (x, h ) 6 + + 3 max h ∈Hn (K h * ν(x) -K h * K h * ν(x)) 2 ≤ 3 h ∈Hn E     1 |T n | u∈Tn K h (x -X u ) -E ν [K h (x -X u )] 2 -a V (x, h ) 6   +   +3 h ∈Hn E     1 |T n | u∈Tn K h * K h (x -X u ) -E ν [K h * K h (x -X u )] 2 -a V (x, h ) 6   +   +3 max h ∈Hn (K h * ν(x) -K h * K h * ν(x)) 2 ≤ T 1 + T 2 + B h (x).
We first give an upper-bound for T 1 . Let h ∈ H n fixed, now remark that, by Lemma 12,

T 1 = +∞ 0 P     1 |T n | u∈Tn K h (x -X u ) -E ν [K h (x -X u )] 2 -a V (x, h ) 6   + ≥ t   dt ≤ +∞ 0 P 1 |T n | u∈Tn K h (x -X u ) -E ν [K h (x -X u )] ≥ t + a V (x, h ) 6 dt ≤ +∞ aV (x,h )/6 P 1 |T n | u∈Tn K h (x -X u ) -E ν [K h (x -X u )] ≥ √ u du ≤ +∞ aV (x,h )/6 exp √ uc K c ρ 4cρ K ∞ √ u 3 + c K c ρ exp - u|T n ||h| 2(c K c ρ + (4/3)c ρ K ∞ √ u) du ≤ I 1 + I 2 ,
where Hence, gathering [START_REF] Chagny | Adaptive warped kernel estimators[END_REF] and [START_REF] Chagny | Adaptive and minimax estimation of the cumulative distribution function given a functional covariate[END_REF], there exists C > 0 depending only on C(P, µ), K, c K and ρ, such that

I
T 1 ≤ C card(H n ) log 2 (|T n |)|T n | - √ ac * ≤ C |T n | -1
with C = 3 max{C 1 , C 2 , C 2 }, c * = min{c * 1 , c * 2 } as soon as a > 4/(c * ) 2 .

Acknowledgement

Remark 11 .

 11 The choice b = 2a is arbitrary since, in theory, any value b such that b ≥ a might work as soon as a sufficiently large. Other choices may be made such as b = a which is the usual choice of Lepski's method, b = a + ε, with ε > 0 or b = (1 + ε)a.

Figure 1 .

 1 Figure 1. Plot of 10 estimators (red dashed lines) obtained from independent copies of X. The black solid line represents the function ν to estimate.

Figure 2 .

 2 Figure 2. Plot of 1 estimator (red solid line) obtained from a copy of X. The yellow solid line represents the function B to estimate.

3 + 3 +I 1 ≤ C 1 aC 2 √ 2 / 3 √ 3 ) 3 +c 2 (√ 6 C 6 √

 331122333266 c K c ρ exp -u|T n ||h| 2(c K c ρ + (4/3)c ρ K ∞ c K c ρ exp -u|T n ||h| 2(c K c ρ + (4/3)c ρ K ∞ √ u) du,where we recall that, since for all h ∈ H n , |h | ≥ log(|Tn| |Tn| , we have V (x, h ) ≤ C(P, µ). We first upper-bound I 1 , ρ K ∞ aC(P, µ)) du≤ C 1 |T n ||h | exp -√ ac * 1 |T n ||h | ≤ C 1 log(|T n |) exp -√ ac * 1 log(|T n |) = C 1 |T n | - √ ac * 1 ,(10)withC 1 = exp( aC(P, µ)/6), C 1 = C 1 2(c K c ρ +(2 √ 2/3 √ 3)c ρ K ∞ aC(P, µ)), c * 1 = C(P, µ)/(12(c K c ρ + (c ρ K ∞ C(P, µ))), using the fact that a ≥ 1 (remark that c * 1 does not depend on a). We turn now to I 2 , remark that the function u→ √ uc K c ρ /( 4cρ K ∞ √ u K c ρ ) is non decreasing and converges to 3c K c ρ /(4c ρ K ∞ ) when u → ∞,hence it is bounded by this quantity. We have then, using again a ≥ 1, |T n ||h |) 2 exp(-√ ac * 2 |T n ||h |) ≤ C 2 log(|T n |)|T n | - √ ac * C 2 = exp(3c K c ρ /(4c ρ K ∞ ), C 2 = C 2 √ a(4c K c ρ +163 (P, µ)c ρ K ∞ ), C 2 = 8C 2 √ 6c K c ρ / C(P, µ) 6c K c ρ / C(P, µ) + 4 3 c ρ K ∞ .

  Theorem 7. Under Assumption 4, if min h∈Hn |h| ≥ log(|T n |)/|T n |, then there exists a minimal value a min > 0 independent of n, such that, for all a > a min ,

  Remark 8. The form of the bias term B h (x) in Inequality (4) is very similar to the one obtained for pointwise adaptive kernel density estimation in [27,Theorem 1]. It can be replaced by an upper-bound, e.g. K 1 ν -K h * ν ∞ , coming from the Young's inequality, as in [14,Theorem 2].

2

.

Once the previous theorem is proved, an immediate corollary follows Corollary 9. Suppose that the assumptions of Theorem 7 are verified and that ν ∈ Σ(β, L, D) where β = (β 1 , ..., β d ) ∈]0, +∞[ d and Σ(β, L, D) is the set of all β-Hölder densities on the open set D ⊆ R d i.e. the set of all functions f : D → R which admits, for all j = 1, ..., d, partial derivatives with respect to x j up to the order β j and verifies, for all x

+ C 2 log 2 (|T n |)|T n | - √ ac * 2 ,(11) with

The authors want to thank Claire Lacour for her helpful advices on constant calibration.