Marc Hoffmann
email: hoffmann@ceremade.dauphine.fr

Yating Liu
email: liu@ceremade.dauphine.fr

A statistical approach for simulating the density solution of a McKean-Vlasov equation

Keywords: Mathematics Subject Classification (2010): 60J60, 65C30, 65C35 Interacting particle systems, McKean-Vlasov models, Euler scheme, Oracle inequalities, Lepski's method

come     A statistical approach for simulating the density solution of a McKean-Vlasov equation

Introduction 1.Motivation

Let P 1 denote the set of probability distributions on R d , d ≥ 1, with at least one moment. Given a time horizon T > 0, an initial condition µ 0 ∈ P 1 and two functions

b : [0, T ] × R d × P 1 → R d and σ : [0, T ] × R d → R d ⊗ R d ,
a probability solution µ = (µ t (dx)) t∈[0,T ] of the nonlinear Fokker-Planck equation

∂ t µ + div b(•, •, µ)µ = 1 2 d k,k =1 ∂ 2 kk σσ µ µ t=0 = µ 0 (1)
exists, under suitable regularity assumptions on the data (µ 0 , b, σ), see in particular Assumptions 2.1, 2.2 and (2.3 or 2.14) in Section 2.1 and 2.4 below. Moreover, for each t > 0, µ t (dx) has a smooth density, also denoted by µ t (x), which is a classical solution to (1). The objective of the paper is to construct a probabilistic numerical method to compute µ T (x) for every T > 0 and x ∈ R d , with accurate statistical guarantees, that improves on previous works on the topic.

Model (1) includes in particular nonlinear transport terms of the form b(t, x, µ) = F f (x) + R d g(x -y)µ(dy) , for smooth functions F, f, g : R d → R d that account for evolution systems with common force f and mean-field interaction g under some nonlinear transformation F together with diffusion σ. Such models and their generalizations have inspired a myriad of application domains over the last decades, ranging from physics [START_REF] Martzel | Mean-field treatment of the many-body Fokker-Planck equation[END_REF][START_REF] Bossy | Instantaneous turbulent kinetic energy modelling based on Lagrangian stochastic approach in CFD and application to wind energy[END_REF] to neurosciences [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF][START_REF] Bossy | Clarification and complement to "Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF], structured models in population dynamics (like in e.g. swarming and flocking models [BFFT12, BFT15, MEK99, BCM07]), together with finance and mean-field games [LL18, CL18, CD18], social sciences and opinion dynamics [START_REF] Chazelle | Well-posedness of the limiting equation of a noisy consensus model in opinion dynamics[END_REF], to cite but a few.

The classical approach of Bossy and Talay [START_REF] Bossy | A stochastic particle method for the McKean-Vlasov and the Burgers equation[END_REF], subsequently improved by Antonelli and Kohatsu-Higa [START_REF] Antonelli | Rate of convergence of a particle method to the solution of the McKean-Vlasov equation[END_REF] paved the way: the probability µ t has a natural interpretation as the distribution L(X t ) of an R d -valued random variable X t , where the stochastic process (X t ) t∈[0,T ] solves a McKean-Vlasov equation

dX t = b(t, X t , L(X t ))dt + σ(t, X t )dB t L(X 0 ) = µ 0 , (2)
for some standard d-dimensional Brownian motion B = (B t ) t∈[0,T ] . The mean-field interpretation of (2) is the limit as N → ∞ of an interacting particle system

dX n t = b(t, X n t , N -1 N m=1 δ X m t )dt + σ(t, X n t )dB n t , 1 ≤ n ≤ N, L(X 1 0 , . . . , X N 0 ) = µ ⊗N 0 ,
where the (B n t ) t∈[0,T ] are independent Brownian motions. Indeed, under fairly general assumptions, see e.g. [START_REF] Henry P Mckean | Propagation of chaos for a class of non-linear parabolic equations[END_REF][START_REF] Sznitman | Topics in propagation of chaos[END_REF][START_REF] Méléard | Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models[END_REF] we have N -1 N n=1 δ X n t → µ t weakly as N → ∞. One then simulates an approximation of (X 1 t , . . . , X N t ) t∈[0,T ]

(3) by a stochastic system of particles ( X1,h t , . . . , XN,h t ) t∈[0,T ] (4) using an Euler scheme with mesh h, see (7) below for a rigorous definition. The synthesised data (4) are then considered as a proxy of the system (3) and a simulation of µ T (x) is obtained via the following nonparametric kernel estimator

µ N,h,η T (x) := N -1 N n=1 η -d K η -1 (x -Xn,h T ) . (5)
Here K : R d → R is kernel function that satisfies R d K(x)dx = 1 together with additional moment and localisation conditions, see Definition 2.4 below.

When K is a Gaussian kernel, [START_REF] Bossy | A stochastic particle method for the McKean-Vlasov and the Burgers equation[END_REF] and [START_REF] Antonelli | Rate of convergence of a particle method to the solution of the McKean-Vlasov equation[END_REF] obtain some convergence rates that depend on N, η, h and the data (µ 0 , b, σ). See also the comprehensive review of Bossy [START_REF] Bossy | Some stochastic particle methods for nonlinear parabolic PDEs[END_REF]. Our objective is somehow to simplify and improve on these results, up to some limitations of course, by relying on adaptive statistical methods. In order to control the strong error E [| µ N,h,η T (x)-µ T (x)| 2 ], a crucial step is to select the optimal bandwidth η that exactly balances the bias and the variance of the statistical error and that depends on the smoothness of the map x → µ T (x). To do so, we rely on recent data-driven selection procedures for η based on the Goldenshluger-Lepski's method in statistics [START_REF] Goldenshluger | Universal pointwise selection rule in multivariate function estimation[END_REF][START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF][START_REF] Goldenshluger | On adaptive minimax density estimation on R d[END_REF]. They do not need any prior knowledge of x → µ T (x), and this is a major advantage since the exact smoothness of the solution map (µ 0 , b, σ) → (x → µ T (x)) may be difficult to obtain, even when µ 0 , b, σ are given analytically. At an abstract level, if we assume that we can observe (X 1 t , . . . , X N t ) t∈[0,T ] exactly, a statistical theory has been recently

proposed to recover µ T (x) in [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF], a key reference for the paper. However, the fact that we need to simulate first a proxy of the idealised data via a system of interacting Euler schemes needs a significant adjustment in order to obtain precise error guarantees. For other recent and alternate approaches, we refer to [START_REF] Belomestny | Projected particle methods for solving McKean-Vlasov stochastic differential equations[END_REF][START_REF] Belomestny | Optimal stopping of McKean-Vlasov diffusions via regression on particle systems[END_REF] and the references therein.

Results and organisation of the paper

In Section 2.1 we give the precise conditions on the data (µ 0 , σ, b). We need a smooth and sub-Gaussian integrable initial condition µ 0 in order to guarantee the existence of a classical solution µ t (x) to the Fokker-Planck equation (1) and sub-Gaussianity of the process (X t ) t∈[0,T ] that solves the McKean-Vlasov equation (2). The coefficients b(t, x, µ t ), σ(t, x) are smooth functions of (t, x) and σ is uniformly elliptic, following the conditions of Gobet and Labart [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF].

In Section 2.2 we construct the estimator µ N,h,η T (x) of (5) via a regular kernel K of order , see Definition 2.4. It depends on the size N of the particle system, the step h > 0 of the Euler scheme and the statistical smoothing parameter η and the order of the kernel . Abbreviating K η = η -d K(η -1 •) for η > 0, we have the following decomposition for the simulation error:

µ N,h,η T (x) -µ T (x) = R d K η (x -•)d(μ N,h T -μh T ) + K η (μ h T -µ T )(x) + K η µ T (x) -µ T (x),
where denotes convolution, μN,h

T = N -1 N n=1 δ Xn,h
T and μh T is the probability distribution of the continuous Euler scheme at time t = T , constructed in (9) below. In turn, the study of the error splits into three terms, with only the first one being stochastic. We prove in Proposition 4.1 a Bernstein concentration inequality for the fluctuation of μN,h T (dy) -μh T (dy). In contrast to [START_REF] Malrieu | Concentration inequalities for Euler schemes[END_REF], it has the advantage to encompass test functions that may behave badly in Lipschitz norm, while they are stable in L 1 like a kernel K η for small values of η. It reads, for an arbitrary bounded test function ϕ:

P R d ϕ d(μ N,h T -μh T ) ≥ ε ≤ κ 1 exp - κ 2 N ε 2 |ϕ| 2 L 2 (µ T ) + |ϕ| ∞ ε
, for every ε ≥ 0, with explicitly computable constants κ 1 and κ 2 that only depend on T and the data (µ 0 , b, σ).

The proof is based on a change of probability argument via Girsanov's theorem. It builds on ideas developed in [START_REF] Lacker | On a strong form of propagation of chaos for McKean-Vlasov equations[END_REF] and [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF], but it requires nontrivial improvements in the case of approximating schemes, in particular fine estimates in Wasserstein distance between the distributions μh T and µ T , thanks to Liu [START_REF] Liu | Optimal Quantization: Limit Theorem, Clustering and Simulation of the McKean-Vlasov Equation[END_REF]. When the measure argument in the drift is nonlinear, we impose some smoothness in terms of linear differentiability, see Assumption 2.14. One limitation of this approach is that we can only encompass a measure term in the drift but not in the diffusion coefficient. The second approximation term K η (μ h T (x) -µ T (x) is managed via the uniform estimates of Gobet and Labart [GL08a] of μh T (x) -µ T (x) so that we can ignore the effect of K η as η → 0 by the L 1 -stability

|K η | L 1 = |K| L 1 . The final bias term K η µ T (x) -µ T (x) is controlled by standard kernel approximation.
In Theorem 2.9 we give our main result, namely a deviation probability for the error µ N,h,η T (x)µ T (x), provided h and η are small enough. In Theorem 2.10, we give a quantitative bound for the expected error of order p ≥ 1 that is further optimised in Corollary 2.11 when x → µ T (x) is k-times continuously differentiable. We obtain a (normalised) error of order N -k/(2k+d) + h, which combines the optimal estimation error in nonparametric statistics with the optimal strong error of the Euler scheme. Finally, we construct in Theorem 2.12 a variant of the Goldenshluger-Lepski's method that automatically selects an optimal data-driven smoothing parameter η without requiring any prior knowledge on the smoothness of the solution. As explained in Section 1.1, our approach suggests a natural way to take advantage of adaptive statistical methods in numerical simulation, when even qualitative information about the object to simulate are difficult to obtain analytically.

We numerically implement our method on several examples in Section 3. We show in particular that it seems beneficial to use high-order kernels (i.e. with many vanishing moments) rather than simply Gaussian ones or more generally even kernels that only have one vanishing moment. This is not a surprise from a statistical point of view, since x → µ T (x) is usually very smooth. Also, the data-driven choice of the bandwidth seems useful even in situations where our assumptions seem to fail (like for the Burgers equation, see in particular Section 3.3). The proofs are delayed until Section 4.

Main results

Notation and assumptions

For x ∈ R d , |x| 2 = x • x denotes the Euclidean norm. We endow P 1 with the 1-Wasserstein metric W 1 (µ, ν) = inf ρ∈Γ(µ,ν) R d ×R d |x -y|ρ(dx, dy) = sup |ϕ| Lip≤1 R d ϕ d(µ -ν),
where Γ(µ, ν) denotes the set of probability measures on R d × R d with marginals µ and ν.

All R d -valued functions f defined on [0, T ], R d , P 1 (or product of those) are implicitly measurable for the Borel σ-field induced by the (product) topology and written componentwise as f = (f 1 , . . . , f d ), where the f i are real-valued. Following [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF], we say that a function

f : [0, T ] × R d → R d is C k,l b for integers k, l ≥ 0 if the real-valued functions f i in the repre- sentation f = (f 1 , . . . , f d )
are all continuously differentiable with uniformly bounded derivatives w.r.t. t and x up to order k and l respectively. We write

C k b (or C l b ) if f is C k,l
b and only depends on the time (or space) variable.

For ϕ : R d → R, a σ-finite measure ν on R d and 1 ≤ p < ∞, we set

|ϕ| L p (ν) = R d |ϕ(x)| p ν(dx) 1/p , |ϕ| L p = R d |ϕ(x)| p dx 1/p , |ϕ| ∞ = sup x∈R d |ϕ(x)|. Assumption 2.1. The law L(X 0 ) = µ 0 (dx) in (2) or equivalently in (1) is sub-Gaussian, that is, R d exp(γ|x| 2 ) µ 0 (dx) < ∞ for some γ > 0.
The strong integrability assumption of µ 0 (together with the subsequent smoothness properties of b and σ) will guarantee that µ t is sub-Gaussian for every t > 0. This is the gateway to our change of probability argument in the Bernstein inequality of the Euler scheme in Proposition 4.1 below.

Assumption 2.2. The diffusion matrix

σ : [0, T ] × R d → R d ⊗ R d is uniformly elliptic, C 1,3 b and ∂ t σ is C 0,1 b .
By uniform ellipticity, we mean inf t,x ξ c(t, x)ξ ≥ δ|ξ| 2 for every ξ ∈ R d and some δ > 0, where c = σσ .

As for the drift b :

[0, T ] × R d × P 1 → R d ,
| b(t, 0, 0)| < ∞.
Assumption 2.3 implies that the function (t, x, µ) → b(t, x, µ) is Lipschitz continuous, i.e. for every (t, x, µ), (s, y, ν) ∈ [0, T ] × R d × P 1 , there exists some constant |b| Lip > 0 such that

|b(t, x, µ) -b(s, y, ν)| ≤ |b| Lip |t -s| + |x -y| + W 1 (µ, ν) (6)
where W 1 denotes the L 1 -Wasserstein distance on P 1 (see e.g. [CD18, (5.4) and Corollary 5.4] for the definition of the Wasserstein distance). As detailed in the proof of Proposition 4.4 below, Assumption 2.3 also implies that (t, x) → b(t, x, µ t ) is C 1,3 b . This together with Assumption 2.2 on the diffusion coefficient enables us to obtain a sharp approximation of the density of the Euler scheme of (2), a key intermediate result. This is possible thanks to the result of Gobet and Labart [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF], which, to our knowledge, is the best result in that direction. See also [START_REF] Bally | The law of the Euler scheme for stochastic differential equations: Convergence rate of the density[END_REF][START_REF] Konakov | Edgeworth type expansions for Euler schemes for differential equations[END_REF][START_REF] Guyon | Euler scheme and tempered distributions[END_REF] for similar results under stronger assumptions.

Construction of an estimator of µ T (x)

We pick an integer M ∈ N \ {0} for time discretisation and h := T M as time step. We set

t m := m • h, m ∈ {0, ..., M }. For 1 ≤ n ≤ N , let B n = (B n t ) t∈[0,T ] be N independent d-dimensional
       Xn,h tm+1 = Xn,h tm + h • b(t m , Xn,h tm , μN,h tm ) + √ h • σ(t m , Xn,h tm )Z n m+1 , 1 ≤ n ≤ N, Z n m+1 := 1 √ h B n tm+1 -B n tm , μN,h tm := 1 N N n=1 δ Xn,h tm , L X1,h 0 , . . . , Xh,N 0 = µ ⊗N 0 , (7)
with X1,h 0 , . . . , Xh,N 0 independent of (B n ) 1≤n≤N for safety. They are appended with their continuous version: for every 0 ≤ m ≤ M -1, t ∈ [t m , t m+1 ) and 1 ≤ n ≤ N :

Xn,h t = Xn,h tm + (t -t m )b(t m , Xn,h tm , μN,h tm ) + σ(t m , Xn,h tm )(B n t -B n tm ).
Definition 2.4. Let ≥ 0 be an integer. A -regular kernel is a bounded function K : R d → R such that for k = 0, . . . , ,

R d |x| k |K(x)|dx < +∞ and R d x k [1] K(x)dx = . . . = R d x k [d] K(x)dx = 1 {k=0} with x = (x [1] , . . . , x [d] ) ∈ R d . (8)
Remark 2.5. A standard construction of -order kernels is discussed in [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF], see also the classical textbook by Tsybakov [START_REF] Alexandre | Introduction to Nonparametric Estimation[END_REF]. In the numerical examples of the paper, we implement in dimension d = 1 the Gaussian high-order kernels described in Wand and Schucany [START_REF] Matthew | Gaussian-based kernels[END_REF]. Given a -regular kernel h in dimension 1, a multivariate extension is readily obtained by tensorisation:

for x = (x [1] , . . . , x [d] ) ∈ R d , set K(x) = h ⊗d (x) = d l=1 h(x [l]
), so that (8) is satisfied.

Finally, we pick a -regular kernel K, a bandwidth η > 0 and set

µ N,h,η T (x) := N -1 N n=1 η -d K η -1 (x -Xn,h T )
for our estimator of µ T (x). It is thus specified by the kernel K, the Euler scheme time step h, the number of particles N and the statistical smoothing parameter η.

Convergence results

Under Assumptions 2.1, 2.2 and (2.3 (or later 2.14)), (1) admits a unique classical solution (t, x) → µ t (x) see e.g. the classical textbook [START_REF] Vladimir I Bogachev | Fokker-Planck-Kolmogorov Equations[END_REF]. Our first result is a deviation inequality for the error µ N,h,η T (x) -µ T (x). We need some further notation.

Definition 2.6. The abstract Euler scheme relative to the McKean-Vlasov equation (2) is defined as

     Xh tm+1 = Xh tm + h • b(t m , Xh tm , µ tm ) + √ h • σ(t m , Xh tm )Z m+1 , Z m+1 = 1 √ h (B tm+1 -B tm ), L( Xh t0 ) = µ 0 , (9)
appended with its continuous version: for every 0 ≤ m ≤ M -1, t ∈ [t m , t m+1 ) and 1 ≤ n ≤ N :

Xh t = Xh tm + (t -t m )b(t m , Xh tm , µ tm ) + σ(t m , Xh tm )(B t -B tm ).
We set μh t = L( Xh t ). The term abstract for this version of the Euler scheme comes from its explicit dependence upon µ tm , that prohibits its use in practical simulation. It however appears as a natural approximating quantity.

Definition 2.7. For ε > 0, the ε-accuracy of the abstract Euler scheme is

h(ε) = sup h > 0, sup x∈R d |μ h T (x) -µ T (x)| ≤ ε , (10)
with the convention sup ∅ = ∞.

Proposition 4.4 below, relying on the estimates of [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF], implies that sup x∈R d |μ h T (x)-µ T (x)| is of order h. Therefore h(ε) is well defined and positive for every ε > 0.

We further abbreviate K η (x) = η -d K(η -1 x) for η > 0.

Definition 2.8. The bias (relative to a -regular kernel K) of a function ϕ : R d → R at x and scale η > 0 is

B η (ϕ, x) = sup 0<η ≤η R d K η (x -y)ϕ(y)dy -ϕ(x)
and for ε > 0, the ε-accuracy of the bias of ϕ is

η(ε) = sup η > 0, B η (ϕ, x) ≤ ε . (11)
Note that whenever ϕ is continuous in a vicinity of x, we have B η (ϕ, x) → 0 as η → 0, and η(ε) is well-defined and positive for every ε > 0.

Theorem 2.9. Work under Assumptions 2.1, 2.2 and 2.

3. Let x ∈ R d , N ≥ 2, ε > 0 and C > 0. For any 0 < h < min h( 1 3 |K| -1 L 1 ε ), CN -1 and 0 < η < η ( 1 3 ε), we have: P µ N,h,η T (x) -µ T (x) ≥ ε ≤ 2κ 1 exp -( 1 3 ) 2 κ 2 N ε 2 |K η (x -•)| 2 L 2 (µ T ) + |K η | ∞ ε , ( 12
)
where κ 1 , κ 2 are constants that only depend on C, x, T , the kernel K, the data (µ 0 , b, σ) and that are defined in Proposition 4.1 below.

Several remarks are in order:

1) We have |K η (x -•)| 2 L 2 (µ T ) ≤ sup y,x-y∈Supp(K) µ T (y)|K η | 2 L 2 = C(µ T , x, K)η -d and |K η | ∞ = η -d |K| ∞ , therefore (12) rather reads P µ N,h,η T (x) -µ T (x) ≥ ε ≤ c 1 exp -c 2 N η d ε 2 1 + ε ,
up to constants c 1 and c 2 that only depend on C, x, T , the kernel K and the data (µ 0 , b, σ).

2)

The factor 1 3 in the upper bound of h, η and the right-hand size of (12) can be replaced by an arbitrary constant in (0, 1) by modifying the union bound argument (31) in the proof. 3) The estimate (12) gives an exponential bound of the form exp(-cN η d ε 2 ) for the behaviour of the error µ N,h,η T (x) -µ T (x) for small ε, provided h and η are sufficiently small. This is quite satisfactory in terms of statistical accuracy, for instance if one wants to implement confidence bands: for any risk level α, with probability bigger than 1 -α, the error

| µ N,h,η T (x) -µ T (x)| is controlled by a constant times » N -1 η d log 1 α .
Theorem 2.9 does not give any quantitative information about the accuracy in terms of h and η. Our next result gives an explicit upper bound for the expected error of order p ≥ 1.

Theorem 2.10. Work under Assumptions 2.1, 2.2 and 2.3. For every p ≥ 1 and x ∈ R d , we have

E µ N,h,η T (x) -µ T (x) p ≤ κ 3 B η (µ T , x) p + (N η d ) -p/2 + h p . ( 13
)
for some κ 3 > 0 depending on p, x, T , |K| L 2 , |K| ∞ and the data (µ 0 , b, σ).

If z → µ T (z) is k-times continuously differentiable in a vicinity of x, then, for a regular kernel of order (k -1) + , we have B η (µ T , x) η k , see for instance the proof of Corollary 2.11. This enables one to optimise the choice of the bandwidth η:

Corollary 2.11. Assume that x → µ T (x) is C k
b for some k ≥ 1 and that K is -regular, with ≥ 0. In the setting of Theorem 2.10, the optimisation of the bandwidth choice η = N -1/(2 min(k, +1)+d) yields the explicit error rate

E µ N,h,η T (x) -µ T (x) p ≤ κ 4 N -min(k, +1)p/(2 min(k, +1)+d) + h p , (14)
for some κ 4 > 0 depending on p, , x, T, |K| L 2 , |K| ∞ and the data (µ 0 , b, σ).

Some remarks: 1) The dependence of κ 3 and κ 4 on (µ 0 , b, σ) is explicit via the bounds that appear in Assumptions 2.1, 2.2 and 2.3. 2) Corollary 2.11 improves the result of [AKH02, Theorem 3.1] for the strong error of the classical solution µ T (x) in the following way: in [START_REF] Antonelli | Rate of convergence of a particle method to the solution of the McKean-Vlasov equation[END_REF], the authors restrict themselves to the one-dimensional case d = 1 with p = 1 in the loss and take a Gaussian kernel for the density estimation, with bandwidth h = T M as in the time discretisation. They obtain the rate (in dimension d = 1)

N -1/2 h -1/2-+ h,
for arbitrary > 0 (up to an inflation in the constant when vanishes), and this has to be compared in our case to the rate (14) which gives

N -k/(2k+1) + h = N -1/2 N 1/2(2k+1) + h,
which yields possible improvement depending on the regime we pick for h, namely whenever N h -(2k+1) which is less demanding as k increases.

One defect of the upper bound ( 14) is that the optimal choice of η depends on the analysis of the bias B η (µ T , x) p , and more specifically on the smoothness k, and that quantity is usually unknown or difficult to compute: Indeed, the information that x → µ T (x) is k-times continuously differentiable with the best possible k is hardly tractable from the data (µ 0 , b, σ), essentially due to the nonlinearity of the Fokker-Planck equation (1). Finding an optimal η without prior knowledge on the bias is a long standing issue in nonparametric statistics. One way to circumvent this issue is to select η depending on the stochastic particle system ( X1,h T , . . . , XN,h T ) itself, in order to optimise the trade-off between the bias and the variance term. While this is a customary approach in statistics, to the best of our knowledge, this is not the case in numerical simulations. It becomes particularly suitable in the case of nonlinear McKean-Vlasov type models.

We adapt a variant of the classical Goldenshluger-Lepski's method [GL08b, GL11, GL14], and refer the reader to classical textbooks such as [START_REF] Giné | Mathematical foundations of infinite-dimensional statistical models[END_REF]. See also the illumating paper by [START_REF] Lacour | Estimator selection: a new method with applications to kernel density estimation[END_REF] that fives a good insight about ideas around bandwidth comparison and data driven selection. For simplicity, we focus on controlling the error for p = 2. Fix x ∈ R d . Pick a finite set H ⊂ [N -1/d (log N ) 2/d , 1], and define

A N η (T, x) = max η ≤η,η ∈H µ N,h,η T (x) -µ N,h,η T (x) 2 -(V N η + V N η ) + ,
where {x} + = max(x, 0) and

V N η = |K| 2 L 2 (log N )N -1 η -d , > 0. ( 15
) Let η N (T, x) ∈ argmin η∈H (A N η (T, x) + V N η ). ( 16
)
We then have the following so-called oracle inequality:

Theorem 2.12. Work under Assumptions 2.1, 2.2 and 2.3. For every x ∈ R d , we have

E µ N,h, η N (t,x) T (x) -µ T (x) 2 ≤ κ 5 min η∈H B η (µ T , x) 2 + N -1 η -d log N + h 2 , ( 17
)
for some κ 5 > 0 depending on x, T , |K| ∞ , Supp(K) and (µ 0 , b, σ), as soon as

Card(H) ≤ N , ≥ 12κ -1 2 max(sup y∈Supp(K) µ T (x -y), 4|K| -2 L 2 |K| 2 ∞ κ -1 2 ) and η d ≥ N -1 log N 3 log 2|K| 2 ∞ 2|K| 2 L 2
. Some remarks: 1) Up to an unavoidable log N -factor, known as the Lepski-Low phenomenon (see [START_REF] Lepskiȋ | A problem of adaptive estimation in Gaussian white noise[END_REF][START_REF] Mark | Nonexistence of an adaptive estimator for the value of an unknown probability density[END_REF]) in statistics, we thus have a data-driven smoothing parameter η N (T, x) that automatically achieves the optimal bias-variance balance, without affecting the effect of the Euler discretisation step h. 2) One limitation of the method is the choice of the pre-factor in the bandwidth selection procedure, which depends on upper bound of µ T locally around x, and more worryingly, on the constant κ -1 2 of Proposition 4.1 below which is quite large. In practice, and this is universal to all smoothing statistical methods, we have to adjust a specific numerical protocol, see Section 3 below.

When we translate this result in terms of number of derivatives for x → µ T (x), we obtain the following adaptive estimation result.

Corollary 2.13. Assume that x → µ T (x) is C k b for some k ≥ 1 and that K is -regular. Specify the oracle estimator µ N,h, η N (T,x) T (x) with H = {(N/ log N ) -1/(2m+d) , m = 1, . . . , + 1}.
In the setting of Theorem 2.12 for every x ∈ R d , we have

E µ N,h, η N (T,x) T (x) -µ T (x) 2 ≤ κ 6 log N N 2 min(k, +1) 2 min(k, +1)+d + h 2 ,
for some κ 6 > 0 depending on T , |K| ∞ , Supp(K) and the data (µ 0 , σ, b).

In practice, see in particular Section 3 below, if x → µ T (x) is very smooth (and this is the case in particular if x → b(t, x, µ t ) is smooth and σ constant) we are limited in the rate of convergence by the order of the kernel. This shows in particular that it is probably not advisable in such cases to pick a Gaussian kernel for which we have the restriction = 1.

Nonlinearity of the drift in the measure argument

In the cas of a drift (t, x, µ) → b(t, x, µ) with a nonlinear dependence in the measure argument µ, the assumptions are a bit more involved. For a smooth real-valued test function ϕ defined on R d , we set

A t ϕ(x) = b(t, x, µ t ) ∇ϕ(x) + 1 2 d k,k =1 c kk (t, x)∂ 2 kk ϕ(x), (18)
that can also be interpreted as the generator of the associated nonlinear Markov process of (X t ) t∈[0,T ] defined in (2). Following [START_REF] Carmona | Probabilistic theory of mean field games with applications[END_REF], we say that a mapping f : P 1 → R d has a linear functional derivative if there exists δ µ f :

R d × P 1 → R d such that f (µ ) -f (µ) = 1 0 dϑ R d δ µ f (y, (1 -ϑ)µ + ϑµ )(µ -µ)(dy), with |δ µ f (y , µ )-δ µ f (y, µ)| ≤ C(|y -y|+W 1 (µ , µ)) and |∂ y (δ µ f (y, µ )-δ µ f (y, µ))| ≤ CW 1 (µ , µ),
for some C ≥ 0. We can iterate the process via mappings

δ µ f : (R d ) × P 1 → R d for = 1, . . . , k defined recursively by δ µ f = δ µ • δ -1 µ f . Assumption 2.14. (Nonlinear representation of the drift.) For every µ ∈ P 1 , (t, x) → b(t, x, µ) is C 1,3
b . Moreover, for every x, (t, y) → A t δ µ b(t, x, y, µ t ) exists and is continuous and bounded.

Finally, µ → b(t, x, µ) admits a k-linear functional derivative (with k ≥ d for d = 1, 2 and k ≥ d/2 for d ≥ 3) that admits the following representation

δ k µ b(t, x, (y 1 , . . . , y k ), µ) = I,m j∈I (δ k µ b) I,j,m (t, x, y j , µ), (19)
where the sum in I ranges over subsets of {1, . . . , k}, the sum in m is finite and the mappings (x, y) → (δ k µ b) I,j,m (t, x, y j , µ) are Lipschitz continuous, uniformly in t and µ. Assumption 2.14 has two objectives: first comply with the property that (t, x) → b(t, x, µ t ) is C 1,3 b in order to apply the result of Gobet and Labart [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF] in Proposition 4.4 below and second, provide a sufficiently smooth structure in (19) in order to implement the change of probability argument of Proposition 4.1. While (19) appears a bit ad-hoc and technical, it is sufficiently general to encompass drift of the form

b(t, x, µ) = j F j t, x, R q 2 G j t, x, R d H j (t, x, z)µ(dz), z λ j (dz ) for smooth mappings F j (t, x, •) : R q3 → R d , G j (t, x, •) : R q1 × R q2 → R q3 , H j (t, x, •) : R d → R q1
and positive measures λ j on R q2 in some cases and combinations of these, see Jourdain and Tse [START_REF] Jourdain | Central limit theorem over non-linear functionals of empirical measures with applications to the mean-field fluctuation of interacting particle systems[END_REF]. Explicit examples where the structure of the drift is of the form 2.14 rather than 2.3 are given for instance in [CDFM18, Oel85, Mél96, JM98].

Theorem 2.15. Work under Assumptions 2.1, 2.2 and 2.14. The results of Theorems 2.9, 2.10, 2.12 and Corollary 2.11, 2.13 remain valid, up to a modification of the constants κ i , i = 1, . . . , 6.

Numerical implementation

We investigate our simulation method on three different examples. The simulation code is available via Google Colab.

1. A linear interaction as in Assumption 2.3 of the form b(t, x, y) = c(x-y). This case is central to several applications (see the references in the introduction for instance) and has moreover the advantage to yield an explicit solution for µ t (x), enabling us to accurately estimate the error of the simulation. See: bit.ly/3LID1rV.

2. A double layer potential with a possible singular shock in the common drift that "stresses" Assumption 2.3. Although the solution is not explicit, the singularity enables us to investigate the automated bandwidth choice on repeated samples and sheds some light on the effect of our statistical adaptive method. See: bit.ly/3yY332B and bit.ly/3ZnvhPj.

3. The Burgers equation in dimension d = 1. Although not formally within the reach Assumption 2.3, we may still implement our method. While we cannot provide with theoretical guarantees, we again have an explicit solution for µ t (x) and we can accurately measure the performance of our simulation method. See: bit.ly/3FEXMAM.

For simplicity, all the numerical experiments are conducted in dimension d = 1. The common findings of our numerical investigations can be summarised as follows: when the density solution x → µ t (x) is smooth, high-order kernels (with vanishing moments beyond order 1) outperform classical kernels such as the Epanechnikov kernel or a standard Gaussian kernel as shown in Cases 1 and 3. Specifically, we implement high-order Gaussian based kernels as developed for instance in Wand and Schucany [START_REF] Matthew | Gaussian-based kernels[END_REF]. Table 1 below gives explicit formulae for K and = 1, 3, 5, 7, 9.

Otherwise, for instance when the transport term has a common exactly Lipschitz continuous component (but no more) as in Case 1, the automated bandwidth method tends to adapt the window to prevent oversmoothing, even with high-order kernels. Overall, we find in all three examples that implementing high-order kernels is beneficial for the quality of the simulation.

order kernel K(x) 1 φ(x) 3 1 2 (3 -x 2 )φ(x) 5 1 8 (15 -10x 2 + x 4 )φ(x) 7 1 48 (105 -105x 2 + 21x 4 -x 6 )φ(x) 9 1 384 (945 -1260x 2 + 378x 4 -36x 6 + x 8 )φ(x)
Table 1: Construction of kernels K of order = 1, 3, 5, 7, 9, with φ(x) = (2π) -1/2 exp(-1 2 x 2 ), dimension d = 1, following [START_REF] Matthew | Gaussian-based kernels[END_REF].

The case of a linear interaction

We consider the simplest McKean-Vlasov SDE with linear interaction

dX t = - R (X t -x)µ t (dx)dt + dB t , L(X 0 ) = N (3, 1 2 )
where (B t ) t∈[0,T ] is a standard Brownian motion and N (ρ, σ 2 ) denotes the Gaussian distribution with mean ρ and variance σ 2 . We do not fulfill Assumption 2.3 here since b(t, x, y) = x -y is not bounded. We nevertheless choose this model for simulation since it has an explicit solution µ t as a stationary Ornstein-Uhlenbeck process, namely, abusing notation slightly,

µ t = L(X t ) = N (3, 1 2 ).
We implement the Euler scheme ( X1,h t , . . . , XN,h t ) t∈[0,T ] defined in (7) with b(t, x, µ) = -R (xy)µ(dy) and σ(t, x) = 1. We pick T = 1, h = 10 -2 T = 10 -2 , for several values of the system size N = 2 7 = 128, 2 8 = 256, . . . , 2 15 = 32768. We then compute

µ N,h,η T (x) := N -1 N n=1 η -d K η -1 (x -Xn,h T )
as our proxy of µ T (x) for = 1, 3, 5, 7, 9 and the kernels given in Table 1. We pick η = N -1/(2( +1)+1) according to Corollary 2.11 since µ T ∈ C k for every k ≥ 1. We repeat the experiment 30 times to obtain independent Monte-Carlo proxies ( µ N,h,η T ) j (x) for j = 1, . . . , 30 and we finally compute the Monte-Carlo strong error

E N = 1 30 30 j=1 max x∈D ( µ N,h,η T ) j (x) -µ T (x) 2 ,
where D is a uniform grid of 1000 points in [0, 6]. (The domain is dictated by our choice of initial condition µ 0 ∼ N (3, 1 2 ).)

Figure 1: Monte-Carlo (for 30 repeated samples) strong error for different kernel orders: log 2 E N as a function of log 2 N for = 1 (purple), = 3 (blue), = 5 (red), = 7 (orange), = 9 (green).

We see that a polynomial error in N is compatible with the data.

Figure 2: Least-square estimates of the slope α of log 2 E N = α log 2 N + noise in a linear model representation. We plot α as a function of the order of the kernel. We see that a higher order for the choice of the kernel systematically improves on the error rate, as predicted by the statistical bias-variance analysis.

In Figure 1, we display on a log-2 scale E N for = 1, 3, 5, 7, 9. In Figure 2 we display the least-square estimates of the slope of each curve of Figure 1 according to a linear model. We thus have a proxy of the rate α of the error E N ≈ N -α for different values of . We clearly see that higher-order kernels are better suited for estimating µ t . This is of course no surprise from a statistical point of view, but this may have been overlooked in numerical probability simulations.

A double layer potential

We consider an interaction consisting of a smooth long range attractive and small range repulsive force, obtained as the derivative of double-layer Morse potential. Such models are commonly used (in their kinetic version) in swarming modelling, see for instance [START_REF] Bolley | Stochastic mean-field limit: nonlipschitz forces and swarming[END_REF]. The corresponding McKean-Vlasov equation is

dX t = R U (X t -x)µ t (dx)dt + dB t , L(X 0 ) = N (0, 1),
where we pick U (x) = -exp(-x 2 ) + 2 exp(-2x 2 ). The potential U and its derivative U are displayed in Figure 3. η N (T, x) -d K η N (T, x) -1 (x -Xn,h T )

as our proxy of µ T (x) for = 1, 3, 5, 7, 9 according to K as in Table 1. The data-driven bandwidth η N (T, x) is computed via the minimisation ( 16), for which one still needs to set the penalty parameter arising in the Goldenschluger-Lepski method, see (15) in particular. The grid is set as H = {(N/ log N ) -1/(2m+1) , m = 1, . . . , + 1}, see in particular Corollary 2.13 in order to mimick the oracle.

In this setting, we do not have access to the exact solution µ T (x); we nevertheless explore several numerical aspects of our method via the following experiments:

• Investigate the effect of the order of the kernel (for an ad-hoc choice of the penalty parameter ). We know beforehand that the mapping x → µ T (x) is smooth, and we obtain numerical evidence that a higher-order kernel gives better results by comparing the obtained µ N,h, η N (T,x) T (x) for different values of as N increases in the following sense: for N = 2 5 particles, the estimates for high order kernels are closer to the estimates obtained for N = 2 10 than lower order kernels.

• Investigate the distribution of the data-driven bandwidth η N (T, x), for repeated samples as x varies and for different values of the penalty over the grid H = {(N/ log N ) -1/(2m+1) , m = 1, . . . , +1}. The estimator tends to pick the larger bandwidth with overwhelming probability, which is consistent with our prior knowledge that x → µ T (x) is smooth.

• In order to exclude an artifact from the preceding experiment, we conduct a cross-experiment by perturbing the drift with an additional Lipschitz common force (but not smoother) that saturates our Assumption 2.3. This extra transport term lowers down the smoothness of

x → µ T (x) and by repeating the preceding experiment, we obtain a different distribution of the data-driven bandwidth η N (T, x), advocating in favour of a coherent oracle procedure.

The effect of the order of the kernel

We display in Figure 4 the graph of x → µ N,h, η N (T,x) T (x) for T = 1, = 1, 3, 5, 7, 9 for N = 2 5 and N = 2 10 . The tuning parameter in the choice of η N (T, x) is set to [1] = 23. The same experiment is displayed in Figure 5 for N = 2 15 . The value N = 2 15 mimicks the asymptotic performance of the procedure as compared to N = 2 5 or N = 2 10 . We observe that the effect of the order of the kernel is less pronounced. A visual comparison of Figure 4 and 5 suggests that a computation with higher order kernels always perform better, since the shape obtained for = 9 for small values of N is closer to the asymptotic proxy N = 2 15 than the results obtained with smaller values of .

The distribution of the data-driven bandwidth

We pick several values for (namely = 20 -1 , 10 -1 , 1, 10, 20) and compute η N (T, x) accordingly for 1200 samples. Figure 6 displays the histogram of the η N (T, x) for x ∈ [-4, 4] (discrete grid with mesh 8 • 10 -2 ) for = 3 and = 5. We observe that the distribution is peaked around large bandwidths at the far right of the spectrum of the histogram, with comparable results for = 7 and = 9. In Figure 7, we repeat the experiment for = 5 over the restricted domain [-1 2 , 1 2 ] (discrete grid, mesh 10 -2 ) where we expect the solution µ T (x) to be more concentrated and see no significant difference. This result is in line of course with the statistical nonparametric estimation of a smooth signal (see e.g. the classical textbook [START_REF] Bernard W Silverman | Density estimation for statistics and data analysis[END_REF]).

[1] This choice is quite arbitrary, the values of the data-driven bandwidths showing stability as soon as ≥ 20.

Cross experiment by reducing the smoothness of x → µ T (x)

We repeat the previous experiment by adding a perturbation in the drift, via a common force V (x) = 2(1-|x|)1 {|x|≤1} . The drift now becomes b(x, µ) = V (x)+ R U (x-y)µ(dy). The common transport term V (x) is no smoother than Lipschitz continuous thus reducing the smoothness of the solution x → µ T (x). In the same experimental conditions as before, Figure 8 displays the distribution of η N (T, x) and is to be compared with Figure 7. We observe that the distribution is modified, in accordance with the behaviour of the oracle bandwidth of a signal with lower order of smoothness. This advocates as further empirical evidence of the coherence of the method. [-1 2 , 1 2 ] over a discrete grid of 100 points, when adding a common perturbation force V (x) = 2(1 -|x|)1 {|x|≤1} in the drift. The effect on the empirical distribution of η N (T, x) is in line with the behaviour of an oracle bandwidth that adjusts smaller bandwidths when the signal is less smooth.

The Burgers equation in dimension d = 1

We consider here the following McKean-Vlasov equation

dX t = R H(x -y)µ t (dy)dt + σdB t , L(X 0 ) = δ 0 , with σ = √ 0.
2, associated with the Burgers equation in dimension d = 1 for H(x -y) = 1 {y≤x} . Although the discontinuity at y = x rules out our Assumption 2.3, we may nevertheless implement our method, since a closed form formula is available for µ t (x). More specifically, for t > 0, the cumulative density function of µ t (x) is explicitly given by

M t (x) = ∞ 0 exp -1 σ 2 (x-y) 2 2t + y dy 0 -∞ exp -1 σ 2 (x-y) 2 2t dy + ∞ 0 exp -1 σ 2 (x-y) 2 2t
+ y dy . (x) for N = 2 10 for several values of . The values = 3 or = 5 provide with the best fit, showing again the benefit of higher-order kernels compared to a standard Gaussian kernel (with = 1). Similar results were obtained for other values of N . Yet, the distribution of η N (T, x) displayed in Figure 11 shows that our method tends to undersmooth the true density. The effect us less pronounced for N = 2 15 . In any case, our numerical results are comparable with [START_REF] Bossy | A stochastic particle method for the McKean-Vlasov and the Burgers equation[END_REF]. We first establish a Bernstein inequality for the fluctuations of μN,h T (dx) -μh T (dx). Proposition 4.1. Work under Assumptions 2.1, 2.2 and 2.3. Let N ≥ 2 and N h ≤ C. For any real-valued bounded function ϕ defined on R d and any ε ≥ 0, we have:

P R d ϕ d(μ N,h T -μh T ) ≥ ε ≤ κ 1 exp - κ 2 N ε 2 |ϕ| 2 L 2 (µ T ) + |ϕ| ∞ ε , ( 20
)
where κ 1 , κ 2 are constants that only depend on C, T and the data (µ 0 , b, σ) that are explicitly given in the proof.

The proof follows the strategy of Theorem 18 in [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF]. We repeat the main steps in order to remain self-contained and highlight the important modifications we need in the context of Euler scheme approximations.

Proof of Proposition 4.1. We work on a rich enough filtered probability space (Ω, F, (F t ) t≥0 , P) in order to accomodate all the random quantities needed. First, note that the abstract Euler scheme ( Xh t ) t∈[0,T ] defined by (9) solves

d Xh t = b(t, Xh t , µ t )dt + σ(t, Xh t )dB t , L( X0 ) = µ 0 , ( 21
)
where t is defined by

∀ 1 ≤ m ≤ M and t ∈ [t m , t m+1 ), t := t m .
Similarly, ( X1,h t , ..., XN,h t ) t∈[0,T ] defined by (7) solves

       d Xn,h t = b(t, Xn,h t , μN,h t )dt + σ(t, Xn,h t )dB n t , 1 ≤ n ≤ N μN,h t = 1 N N n=1 δ Xn,h t , L( X1,h 0 , . . . , XN,h 0 ) = µ ⊗N 0 . ( 22
)
Step 1. We construct a system (Y 1,h t , ..., Y N,h t

) t∈[0,T ] of N independent copies of the abstract Euler schemes (21) (thus without interaction) via

dY n,h t = b(t, Y n,h t , µ t )dt + σ(t, Y n,h t )dB n t , 1 ≤ n ≤ N (Y 1,h 0 , . . . , Y N,h 0 ) = ( X1,h 0 , . . . , XN,h 0 ), ( 23
) For t ∈ [0, T ], let L N,h t := N n=1 t 0 (c -1/2 b)(s, Y n,h s , 1 N N n=1 δ Y n,h s ) -(c -1/2 b)(s, Y n,h s , µ s ) dB n s and E t (L N,h • ) := exp L N,h t - 1 2 L N,h • t ,
where c -1/2 is any square root of c -1 = (σσ ) -1 and L N,h

• t denotes the predictable compensator of L N,h t . The following estimate is a key estimate to proceed to a change of probability. Its proof is delayed until the end of this section. Lemma 4.2. Work under Assumptions 2.1, 2.2 and 2.3. Assume that N h ≤ C for some constant C > 0. For every κ > 0, there exists δ(κ) > 0 such that

∀ δ ∈ [0, δ(κ)], sup t∈[0,T -δ(κ)] E P exp κ L N,h • t+δ -L N,h • t ≤ C ( 24
)
where C > 0 only depends on C, T and the data (µ 0 , b, σ).

By taking κ = 1 2 and by applying Novikov's criterion (see e.g. [KS91, Proposition 3.5.12, Corollary 3.5.13 and 3.5.14]), Lemma 4.2 implies that E t (L N,h

•

) t∈[0,T ] is a (F t , P) martingale as soon as N h ≤ C. We may then define another probability distribution Q on Ω, F T by setting

Q := E T (L N,h • ) • P. ( 25
)
By Girsanov's Theorem (see e.g. [KS91, Theorem 3.5.1]), we have

Q • (Y 1,h t , ..., Y N,h t ) -1 t∈[0,T ] = P • (X 1,h t , ..., X N,h t ) -1 t∈[0,T ] . ( 26
)
Step 2. We claim that for any subdivision 0 = t 0 < t 1 < ... < t K = T and for any F T -measurable event A N , we have

E P Q A N F tj-1 ≤ E P Q A N F tj 1 4 E P exp 2 L N,h • tj -L N,h • tj-1 1 4 , 1 ≤ j ≤ K. ( 27
)
The proof is the same as in Step 3 of the proof of Theorem 18 in [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF] and is inspired from the estimate (4.2) in Theorem 2.6 in [START_REF] Lacker | On a strong form of propagation of chaos for McKean-Vlasov equations[END_REF]. We repeat the argument: we have

E P Q A N F tj-1 = E P E Q Q A N F tj F tj-1 = E P E P E tj (L N,h • ) E tj-1 (L N,h • ) Q A N F tj F tj-1 = E P E tj (L N,h • ) E tj-1 (L N,h • ) Q A N F tj ,
where the second inequality follows from Bayes's rule in [KS91, Lemma 3.5.3]. Next, we have

E tj (L N,h • ) E tj-1 (L N,h • ) = E tj 2(L N,h • -L N,h tj-1 ) 1 2 • exp L N,h • tj -L N,h • tj-1 1 2 .
The process

E tj 2(L N,h • -L N,h tj-1 ) t∈[tj-1,T ] is a (F t , P) martingale if N h ≤ C. Hence E P E tj 2(L N,h • -L N,h tj-1 ) = 1 for every t ∈ [t j-1 , T ]. It follows that E P Q A N F tj-1 = E P E tj 2(L N,h • -L N,h tj-1 ) 1 2 • exp L N,h • tj -L N,h • tj-1 1 2 Q A N F tj ≤ E P î exp L N,h • tj -L N,h • tj-1 Q A N F tj 2 ó 1 2 ≤ E P exp 2 L N,h • tj -L N,h • tj-1 1 4 E P î Q A N F tj 4 ó 1 4 ≤ E P exp 2 L N,h • tj -L N,h • tj-1 1 4 E P Q A N F tj 1 4 ,
where the first two inequalities follows from Cauchy-Schwarz's inequality and the last inequality follows from Jensen's inequality. Thus ( 27) is established.

Step

3. Let A N ∈ F T . Since E t (L N,h • ) t∈[0,T ] is a (F t , P)
Q(A N ) = E Q Q(A N F 0 ) = E P Q(A N F 0 ) .
Now, let m ≥ 1 and take a subdivision 0 = t 0 < t 1 < . . . < t m = T such that

∀ 0 ≤ k ≤ m -1, t k+1 -t k ≤ δ(2)
where δ(2) is the constant in Lemma 4.2 for κ = 2. It follows by ( 27) that

E P Q(A N F 0 ) ≤ E P Q A N F T 4 -m m j=1 E P exp 2 L N,h • tj -L N,h • tj-1 j 4 ≤ P(A N ) 4 -m m j=1 E P exp 2 L N,h • tj -L N,h • tj-1 j 4 (since A N ∈ F T ) ≤ P(A N ) 4 -m sup t∈[0,T -δ(2)] E P exp 2 L N,h • t+ε -L N,h • t m(m+1) 8 ≤ P(A N ) 4 -m C m(m+1) 8 (28)
by applying Lemma 4.2 with κ = 2.

Step 4. We first recall the Bernstein's inequality: if Z 1 , . . . , Z N are real-valued centred independent random variables bounded by some constant Q, we have

∀ ε ≥ 0, P N n=1 Z n ≥ ε ≤ exp - ε 2 2 N n=1 E [Z 2 n ] + 2 3 Qε . (29)
Now, for ε ≥ 0, we pick

A N T := N n=1 ϕ Y n,h T - R d ϕ(y)μ h T (dy) ≥ N ε so that A N T ∈ F T . Since (Y 1,h T , . . . , Y N,h T
) are independent with common distribution μh T , we have

P(A N T ) ≤ exp - N ε 2 2|ϕ| 2 L 2 (μ T ) + 2 3 |ϕ| ∞ ε . (30)
by ( 29) with

Z n = ϕ(Y n,h T ) - R d ϕ(x)μ h T (dx), having E [Z 2 n ] ≤ |ϕ| 2 L 2 (μ h T ) and Q = |ϕ| ∞ . It follows by (26) that P R d ϕ d(μ N,h T -μh T ) ≥ ε = P N n=1 ϕ Xn,h T - R ϕ(y)μ h T (dy) ≥ N ε = Q(A N T ).
By

Step 3, we infer

Q(A N T ) ≤ P(A N T ) 4 -m C m(m+1) 8 ≤ C m(m+1) 8 exp -4 -m N ε 2 2|ϕ| 2 L 2 (μ h T ) + 2 3 |ϕ| ∞ ε ,
and we conclude by taking

κ 1 = C m(m+1) 8 and κ 2 = 2 -1 4 -m .
Completion of proof of Theorem 2.9

Recall the notation K η (x) = η -d K(η -1 x) and set ϕ ψ(x) = R d ϕ(x -y)ψ(y)dy for the convolution product between two integrable functions ϕ, ψ : R d → R. We have

P µ N,h,η T (x) -µ T (x) ≥ ε ≤ P µ N,h,η T (x) -K η μh T (x) ≥ ε/3 + 1 {|Kη μh T (x)-Kη µ T (x)|≥ε/3} + 1 {|Kη µ T (x)-µ T (x)|≥ε/3} . (31)
We have

K η μh T (x) -K η µ T (x) ≤ |K| L 1 sup x∈R d μh T (x) -µ T (x) since |K η | L 1 = |K| L 1 and this last term is strictly smaller than ε/3 for h < h |K| -1 L 1 1 3 ε , therefore the second indicator is 0. Likewise K η µ T (x) -µ T (x) ≤ B η (µ T , x) < ε/3 for η < η( 1 3 ε)
and the third indicator is 0 as well. Finally

P µ N,h,η T (x) -K η μh T (x) ≥ 1 3 ε = P R d K η (x -•)d(μ N,h T -μh T ) ≥ 1 3 ε ≤ 2κ 1 exp - κ 2 N ( 1 3 ε) 2 |K η (x -•)| 2 L 2 (µ T ) + |K η | ∞ ε
by Proposition 4.1 and Theorem 2.9 follows.

Proof of Lemma 4.2

Writing µ N,h s = 1 N N n=1 δ Y n,h s , for κ > 0, we have κ L N,h • t+δ -L N,h • t = κ N n=1 t+δ t b(s, Y n,h s , µ N,h s ) -b(s, Y n,h s , µ s ) c(s, Y n,h s ) -1 b(s, Y n,h s , µ N,h s ) -b(s, Y n,h s , µ s ) ds ≤ κ sup t∈[0,T ] |Tr(c(t, •) -1 )| ∞ N n=1 t+δ t b(s, Y n,h s , µ N,h s ) -b(s, Y n,h s , µ s ) 2 ds ≤ 2κ sup t∈[0,T ] |Tr(c(t, •) -1 )| ∞ N n=1 t+δ t b(s, Y n,h s , µ N,h s ) -b(s, Y n,h s , μh s ) 2 ds + N | b| 2 Lip t+δ t W 1 (μ h s , µ s ) 2 ds ,
∀ s, t ∈ [0, T ], s < t, E |X t -X s | 2 1/2 ≤ κ√ t -s
for some κ that depends on (µ 0 , b, σ) and T , see [

W 1 (μ h s , µ s ) ≤ κ 7 h 1/2 ,
where κ 7 depends on the data (µ 0 , b, σ) and T , by rewriting (2) as a Brownian diffusion The random variables (ξ n,n s ) 1≤n ≤N,n =n are centered, identically distributed and conditionally independent given Y n,h s . Moreover, for every integer m ≥ 1, we have the following estimate Lemma 4.3. We have

dX t = b µ (t, X t )dt + σ(t, X t
sup s∈[0,T ] E |Y n,h s | 2m ≤ κ m 8 m!
where κ 8 > 0 depends on T and the data (µ 0 , b, σ).

The proof is standard (see for instance [START_REF] Méléard | Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models[END_REF][START_REF] Sznitman | Topics in propagation of chaos[END_REF] or [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF] for a control of growth in the constant m and we omit it). From Lemma 4.3, we infer, for every m ≥ 1:

E ξ n,n s 2m | Y n,h s = E b(s, Y n,h s , Y n ,h s ) - R d b(s, Y n,h s , y)μ h s (dy) 2m Y n,h s ≤ E 2 2m-1 b(s, Y n,h s , Y n ,h s ) -b(s, Y n,h s , 0) 2m + b(s, Y n,h s , 0) - R d b(s, Y n,h s , y)μ h s (dy) 2m Y n,h s ≤ E 2 2m-1 | b| 2m Lip |Y n ,h s | 2m + | b| 2m Lip E |Y n ,h s | 2m Y n,h s = 2 2m | b| 2m Lip E |Y n ,h s | 2m ≤ 4| b| 2 Lip κ 8 m m!
E exp (N -1)(S n,N s ) 2 [k] 8κ 9 = E E exp (N -1)(S n,N s ) 2 [k] 8κ 9 Y n,h s ≤ 2,
and in turn

E exp (N -1)|S n,N s | 2 8dκ 9 ≤ d -1 d k=1 E exp (N -1)(S n,N s ) 2 [k] 8κ 9 ≤ 2. ( 32
) Likewise E exp |ξ n,n s | 2 8dκ 9 ≤ 2. ( 33
)
Abbreviating τ = 2κ sup t∈[0,T ] |Tr(c(t, •) -1 )| ∞ , it follows that 32) and (33). We obtain Lemma 4.2 with C = 2 exp(κ 10 C) with κ 10 = κ 2 7 384dκ8 .

E exp κ( L N,h • t+δ -L N,h • t ) ≤ e τ | b| 2 Lip κ 2 7 δN h E exp τ t+δ t N n=1 N -1 N n =1 ξ n,n s 2 ds ≤ e τ | b| 2 Lip κ 2 7 δN h 1 2δ t+δ t N -1 N n=1 E exp 4δτ N N -1 ξ n,n s 2 + E exp 4τ δN N -1 N S n,N s 2 ds ≤ 2e τ | b| 2 Lip κ 2 7 δN h as soon as δ ≤ (32τ dκ 9 ) -1 N N -1 = (64κ sup t∈[0,T ] |Tr(c(t, •) -1 )| ∞ dκ 9 ) -1 N N -1 by (

Proof of Theorem 2.10

We first need a crucial approximation result of the density of a diffusion process by the density of its Euler scheme counterpart. We heavily rely on the sharp results of Gobet and Labart [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF].

Proposition 4.4. Work under Assumptions 2.1, 2.2 and 2.3. For any 0 < t min ≤ T , we have

sup (t,x)∈[tmin,T ]×R d μh t (x) -µ t (x) ≤ κ 10 • h, (34)
for some κ 10 > 0 depending on t min , T, d and the data (µ 0 , b, σ).

Proof. For x 0 ∈ R d , let (ξ x0 t ) t∈[0,T ] be a diffusion process of the form

ξ x0 t = x 0 + t 0 b(s, ξ x0 s )ds + t 0 σ(s, ξ x0 s )dB s , (35)
where σ :

[0, T ] × R d → R d ⊗ R d and b : [0, T ] × R d → R d are both C 1,3
b and σ is uniformly elliptic and ∂ t σ is C 0,1 b . We associate its companion Euler scheme ( ξx0,h t ) t∈[0,T ] :

           ξx0,h tm+1 = ξx0,h tm + h • b(t m , ξx0,h tm ) + √ h σ(t m , ξx0,h tm )Z m+1 , Z n m+1 := 1 √ h B tm+1 -B tm , ξx0,h t0 = x 0 , ∀ t ∈ [t m , t m+1 ), ξx0,h t = ξx0,h tm + (t -t m ) b(t m , ξx0,h tm ) + σ(t m , ξx0,h tm )(B t -B tm ).
For t > 0, both L(ξ x0 t ) and L( ξx0,h t ) are absolutely continuous, with density p t (x 0 , x) and ph t (x 0 , x) w.r.t. the Lebesgue measure dx on R d . By Theorem 2.3 in [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF], there exist two constants c 1 and c 2 depending on T , d and the data ( σ, b) such that

ph t (x 0 , x) -p t (x 0 , x) ≤ c 1 h t -d+1 2 exp - c 2 |x -x 0 | 2 t . (36)
Taking σ(t, x) = σ(t, x) and b(t, x) = b(t, x, µ t ), we can identify μt (x) with R d ph t (x 0 , x)µ 0 (dx 0 ) and µ t (x) with R d p t (x 0 , x)µ 0 (dx 0 ). For every t ∈ [t min , T ], we have

μh t (x) -µ t (x) ≤ R d ph t (x 0 , x) -p t (x 0 , x) µ 0 (dx 0 ) ≤ c 1 h R d t -d+1 2 exp - c 2 |x -x 0 | 2 t µ 0 (dx 0 ) ≤ c 1 t -d+1 2 min • h =: κ 10 • h,
by Theorem 2.3 in [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF] since Assumption 2.3 implies (t, x) → b(t, x, µ t ) is C 1,3 b . More specifically, for (X t ) t∈[0,T ] solving (2), we have b(t, x, µ t ) = E b(t, x, X t ) .

By dominated convergence, x → b(t, x, µ t ) is C 3 b thanks to the regularity assumptions on x → b(t, x, y). For the regularity in time, by Itô's formula

E b(t, •, X t ) = E b(0, •, X 0 ) + t 0 E ∂ T b(s, •, X s ) ds + t 0 E A s b(s, •, X s ) ds, where (A s ) s∈[0,T ] is the family of generators defined in (18) applied to y → b(s, •, y) = ( b(s, •, y) [1] , . . . , b(s, •, y) [d] ) componentwise. Moreover, for each component 1 ≤ k ≤ d, A s b(s, •, X s ) [k] = E b(s, •, X s )] ∇ y b(s, •, X s ) [k] + 1 2 d l,l =1 c ll (s, X s )∂ y l y l b(s, •, X s ) [k] , hence E A s b(s, •, X s ) [k] = E b(s, •, X s )] E F k (s, X s ) + E G k (s, X s ) , with F k (s, y) = ∇ y b(s, •, y) [k] , and G k (s, y) = 1 2 d l,l =1 c ll (s, y)∂ y l y l b(s, •, y) [k] .
All three functions (s, y) → b(s, •, y), (s, y) → F k (s, y) and (s, y) → G k (s, y) are continuous and bounded by Assumption 2.3, and so are s → E b(s, •, X s )], s → E F k (s, X s ) and s → E G k (s, X s ) by dominated convergence using in particular that s → X s is stochastically continuous. Hence

s → E A s b(s, •, X s ) is continuous and s → E ∂ s b(s, •, X s ) too, using again Assumption 2.3. It follows that t → b(t, x, µ t ) is C 1 b on [0, T ].
Completion of proof of Theorem 2.10

We have

E µ N,h,η T (x) -µ T (x) p ≤ 3 p-1 E µ N,h,η T (x) -K η μh T (x) p + K η μh T (x) -K η µ T (x) p + K η µ T (x) -µ T (x) p ,
and we plan to bound each term separately. First, by Proposition 4.1

E µ N,h,η T (x) -K η μh T (x) p = ∞ 0 P µ N,h,η T (x) -K η μh T (x) ≥ z 1/p dz ≤ 2κ 1 ∞ 0 exp - κ 2 N z 2/p |K η (x -•)| 2 L 2 (µ T ) + |K η | ∞ z 1/p dz ≤ 2κ 1 ∞ 0 exp - κ 2 N η d z 2/p sup y∈Supp(K) µ T (x -y)|K| 2 L 2 + |K| ∞ z 1/p dz ≤ 2κ 1 c p κ -p/2 2 sup y∈Supp(K) µ T (x -y) p/2 |K| p L 2 (N η d ) -p/2 , stemming from the estimate ∞ 0 exp - az 2/p b + cz 1/p dz ≤ c p max a b -p/2 , a c -p , valid for a, b, c, p > 0, with c p = 2 ∞ 0 exp -1 2 (min(z, √ z)) 2/p dz. Next, K η μh T (x) -K η µ T (x) p ≤ sup y∈R d μh T (y) -µ T (y) p R d K η (x -y) dy p ≤ κ p 10 |K| p L 1 h p
by Proposition 4.4. Finally, by definition

K η µ T (x) -µ T (x) p ≤ B η (µ T , x) p ,
and we obtain Theorem 2.10 with

κ 3 = 3 p-1 max 2κ 1 c p κ -p/2 2 sup y∈Supp(K) µ T (x -y) p/2 |K| p L 2 , κ p 10 |K| p L 1 , 1 .

Proof of Corollary 2.11

By Taylor's formula, we have, for x, y ∈ R d , η > 0 and any 1 ≤ k ≤ k:

µ T (x + ηy) -µ T (x) = 1≤|α|≤k -1 ∂ α µ T (x) α! (ηy) α + k |α|=k (ηy) α α! 1 0 (1 -t) k -1 ∂ α µ T (x + tηy)dt, (37)
with multi-index notation α = (α 1 , . . . , α d ),

α i ∈ {0, 1, . . .}, α! = α 1 !α 2 ! • • • α d !, |α| = α 1 + α 2 + • • • + α d ,
and for x = (x 1 , . . . , x d ):

x α = x α1 1 x α2 2 • • • x α d d , and ∂ α f = ∂ |α| f ∂ α1 x 1 ∂ α2 x 2 • • • ∂ α d x d .
It follows that

K η µ T (x) -µ T (x) = η -d R d K η -1 (x -y) µ T (y) -µ T (x) dy = R d K(-y) µ T (x + ηy) -µ T (x) dy = R d K(-y) min(k, + 1) |α|=min(k, +1) (ηy) α α! 1 0 (1 -t) min(k, +1)-1 ∂ α µ T (x + tηy)dt dy,
thanks to (37) with k = min(k, + 1) and the cancellation property (8) of the kernel K that eliminates the polynomial term in the Taylor's expansion, hence

K η µ T (x) -µ T (x) ≤ |µ T | min(k, +1),K η min(k, +1) , with |f | min(k, +1),K = |α|=min(k, +1) |∂ α f | ∞ α! R d |y| α |K(y)|dy. (38)
It follows that

B η (µ T , x) p ≤ |µ T | p min(k, +1),K η min(k, +1)p . (39)
Applying now Theorem 2.10, we obtain

E µ N,h,η T (x) -µ T (x) p ≤ κ 3 |µ T | p min(k, +1),K η min(k, +1)p + (N η d ) -p/2 + h p ≤ κ 3 |µ T | p
min(k, +1),K + 1)N -min(k, +1)p/(2 min(k, +1)+1) + h p with the choice η = N -1/(2 min(k, +1)+d) hence the corollary with κ 4 = κ 3 (1 + |µ T | p min(k, +1),K ).

Proof of Theorem 2.12

This is an adaptation of the classical Goldenshluger-Lepski method [GL08b, GL11, GL14]. We repeat the main arguments and highlight the differences due to the stochastic approximation and the presence of the Euler scheme.

Abbreviating A N η (T, x) by A N η , we have, for every η ∈ H, E µ N,h, η N (T,x) T (x) -µ T (x) 2 ≤ 2E µ N,h, η N (T,x) T (x) -µ N,h,η T (x) 2 + 2E µ N,h,η T (x) -µ T (x) 2 ≤ 2E µ N,h, η N (T,x) T (x) -µ N,h,η T (x) 2 -V N η -V N η N (T,x) + + V N η + V N η N (T,x) + 2E µ N,h,η T (x) -µ T (x) 2 ≤ 2E A N max(η, η N (T,x)) + V N η + V N η N (T,x) + 2E µ N,h,η T (x) -µ T (x) 2 ≤ 2 E A N η + V N η + 2E A N η N (T,x) + V N η N (T,x) + 2E µ N,h,η T (x) -µ T (x) 2 ≤ 4E A N η + V N η + κ 3 (B η (x, µ T ) 2 + (N η d ) -1 + h 2 ) ≤ κ 11 E A N η + V N η + B η (x, µ T ) 2 + h 2 with κ 11 = 2 max(2 + κ3 |K| 2 L 2
, κ 3 ) by definition of η N (T, x) in ( 16), of V N η in (15) and using Theorem 2.10 with p = 2 to bound the last term.

We next bound the term E A N η . For η, η ∈ H, with η ≤ η, we start with the decomposition

µ N,h,η T (x) -µ N,h,η T (x) = µ N,h,η T (x) -K η µ h T (x) + K η µ h T (x) -K η µ T (x) + K η µ T (x) -µ T (x) -µ N,h,η T (x) -K η µ h T (x) -K η µ h T (x) -K η µ T (x) -K η µ T (x) -µ T (x) ,
and thus, by Proposition 4.4, we infer µ N,h,η T (x) -µ N,h,η T (x) ≤ µ N,h,η T (x) -K η µ h T (x) + µ N,h,η T (x) -K η µ h T (x) + 2κ 10 |K| L 1 h + 2 B η (µ T , x) using η ≤ η to bound the second bias term.

µ N,h,η T (x) -µ N,h,η T (x) 2 -V N η -V N η ≤ 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η + 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η + 3 2κ 10 |K| L 1 h + 2B η (x, µ T ) 2 .
Taking maximum for η ≤ η, we thus obtain max

η ≤η µ N,h,η T (x) -µ N,h,η T (x) 2 -V N η -V N η ≤ 3( µ N,h,η T (x) -K η µ T (x)) 2 -V N η + + max η ≤η 3( µ N,h,η T (x) -K η µ T (x)) 2 -V N η + + 24 κ 2 10 |K| 2 L 1 h 2 + B η (x, µ T ) 2 .
We finally bound the expectation of each term. First, we have, by Proposition 4.1

E 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η + = ∞ 0 P 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η ≥ z dz = ∞ 0 P µ N,h,η T (x) -K η µ h T (x) ≥ 3 -1/2 (z + V N η ) 1/2 dz ≤ 2κ 1 ∞ V N η exp - κ 2 N η d 1 3 z sup y∈Supp(K) µ T (x -y)|K| 2 L 2 + |K| ∞ 3 -1/2 z 1/2 dz ≤ 2κ 1 ∞ V N η exp - κ 2 N η d z 6 sup y∈Supp(K) µ T (x -y)|K| 2 L 2 dz + 2κ 1 ∞ V N η exp - κ 2 N η d z 1/2 2 √ 3|K| ∞ dz ≤ 12κ1 κ2 sup y∈Supp(K)
µ T (x -y)|K| 2 L 2 (N η d ) -1 N -κ2 /6 sup y∈Supp(K) µ T (x-y)

+ 48 √ 3κ1 κ2 |K| ∞ |K| L 2 1/2 (N η d ) -3/2 (log N ) 1/2 exp - 1/2 κ2|K| L 2 2 √ 3|K|∞ (N η d ) 1/2 (log N ) 1/2 ,
where we used ∞ ν exp(-z 1/2 )dz ≤ 4ν 1/2 exp(-ν 1/2 ) for ν ≥ 16. The specifications of and η d of Theorem 2.12 entail

E 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η + ≤ κ 12 (log N ) 1/2 N -2 ,
with κ 12 = 2 max( 12κ1 κ2 sup y∈Supp(K) µ T (x -y)|K| 2 L 2 , 48

√ 3κ1 κ2 |K| ∞ |K| L 2 1/2 ).
In the same way, we have the rough estimate

E max η ≤η 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η + ≤ η ∈H E 3 µ N,h,η T (x) -K η µ h T (x) 2 -V N η +
≤ Card(H)κ 12 (log N ) 1/2 N -2 ≤ κ 12 (log N ) 1/2 N -1 , using the previous bound. We thus have proved

E A N η ≤ 2κ 12 (log N ) 1/2 N -1 + 24 κ 2 10 |K| 2 L 1 h 2 + B η (x, µ T ) 2 ≤ κ 13 V N η + B η (x, µ T ) 2 + h 2 ,
with κ 13 = max 2κ 12 /|K| 2 L 2 , 24κ 2 10 |K| 2 L 1 , 24 . Back to the first step of the proof, we infer

E µ N,h, η N (T,x) T (x) -µ T (x) 2 ≤ κ 11 E A N η + V N η + B η (x, µ T ) 2 + h 2 ≤ κ 5 V N η + B η (x, µ T ) 2 + h 2 ,
where now κ 5 = κ 11 (1 + κ 13 ). Since η ∈ H is arbitrary, the proof of Theorem 2.12 is complete.

Proof of Corollary 2.13

For every η ∈ H, we have

B η (µ T , x) 2 ≤ |µ T | 2 k,K η 2 min(k, +1) ,
where |µ T | k,K is defined via (38), see (39) in the proof of Corollary 2.11. By Theorem 2.12, we thus obtain with κ 6 = 2κ 5 max(|µ T | 2 k,K , |K| 2 L 2 , 1), using the fact that the minimum in η ∈ H is attained for η = (N/ log N ) 1/(2 min(k, +1)+d) .

Proof of Theorem 2.15

We briefly outline the changes that are necessary in the previous proofs to extend the results to a nonlinear drift in the measure argument that satisfies Assumption 2.14. Theorem 2.9 under Assumption 2.14 Only Lemma 4.2 needs to be extended to the case of a nonlinear drift in order to obtain Proposition 4.1, the rest of the proof remains unchanged. This amounts to prove that the random variables b(s, Y n,h 2 + min W 1 ( µ N,h s , μs ) 2 , W 1 ( µ N,h s , μs ) 2k , following the proof of Lemma 21 in [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF]. Next, min W 1 ( µ N,h s , μs ) 2 , W 1 ( µ N,h s , μs ) 2k is sub-Gaussian with the right order in N , thanks to the sharp deviation estimates of Theorem 2 in Fournier and Guillin [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF]. As for the main part of the previous expansion, we use a bound of the form for some constant c depending on b, see Lemma 22 in [START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF], showing eventually that this term is sub-Gaussian with the right order in N . Lemma 4.2 follows.

E (R d ) l
Theorem

Remaining proofs

The proofs of Corollary 2.11, Theorem 2.12 and Corollary 2.13 remain unchanged under Assumption 2.14.

Figure 3 :

 3 Figure 3: Plot of the potential U (red) and its derivative U (blue). We implement the Euler scheme ( X1,h t , . . . , XN,h t ) t∈[0,T ] defined in (7) for coefficients b(x, µ) = R U (x -y)µ(dy) and σ(t, x) = 1. We pick T = 1, h = 10 -2 T = 10 -2 , for several values of the system size N = 2 5 = 32, 2 6 = 64, . . . , 2 16 = 65736. We then compute

Figure 4 :

 4 Figure 4: The graph of x → µ N,h, η N (t,x) t (x). The domain x ∈ [-4, 4] is computed over a discrete grid of 2000 points, i.e. mesh 4 • 10 -3 ) for N = 2 5 (Left) and N = 2 10 (Right).

Figure 5 :

 5 Figure 5: Same experiment as in Figure 4 mimicking an asymptotic behaviour of the procedure for N = 2 15 .

Figure 6 :

 6 Figure 6: Distribution of η N (T, x). The domain x ∈ [-4, 4] is computed over a discrete grid of 100 points, i.e. mesh 8 • 10 -2 ), N = 2 5 , 2 6 , • • • , 2 16 for = 3 (Top) and = 5 (Bottom).

Figure 7 :

 7 Figure 7: Same experiment as in Figure 6 (Bottom) on the restricted domain [-1 2 , 1 2 ] over a discrete grid of 100 points. The results are comparable with the experiment displayed in Figure 6 (Bottom).

Figure 8 :

 8 Figure8: Same experiment as in Figure7on the restricted domain [-1 2 , 1 2 ] over a discrete grid of 100 points, when adding a common perturbation force V (x) = 2(1 -|x|)1 {|x|≤1} in the drift. The effect on the empirical distribution of η N (T, x) is in line with the behaviour of an oracle bandwidth that adjusts smaller bandwidths when the signal is less smooth.

Figure 9 :

 9 Figure 9: Plots of x → M T (x) (Left) and x → µ T (x) (Right).

Figure 9

 9 Figure 9 displays the graph of x → M T (x) and x → µ T (x) for T = 1. We display in Figure 10 the reconstruction of x → µ T (x) via µ N,h, η N (T,x) T

Figure 10 :

 10 Figure 10: Reconstruction of µ T by µ N,h, η N (T,•) T with = 23, for different kernel orders for N = 2 10 (Top) and N = 2 15 (Bottom). Higher order kernels outperform the reconstruction provided with standard Gaussian kernels ( = 1) for N = 2 10 .

Figure 11 :

 11 Figure 11: Distribution of η N (T, x). The domain x ∈ [-3, 4] is computed for 1200 samples over a discrete grid of 100 points, i.e. mesh 8 • 10 -2 , N = 2 5 , 2 6 , . . . , 2 16 for = 7. The method tends pick the largest bandwidth.

  martingale, P and Q coincide on F 0 see e.g. [KS91, Section 3 -(5.5)] . It follows that

  )dB t and by applying the classical convergence result of the Euler scheme for a Brownian diffusion (see e.g. [Pag18, Theorem 7.2]) since for every s ∈ [0, T ], W 1 (μ h s , µ s ) ≤ E |X s -Xh s | . Writing b(s, Y n,h s , µ N,h s ) -b(s, Y n,

≤ exp κ9v 2 2 ,

 2 and in turn, by letting κ 9 = 24| b| 2 Lip κ 8 , we derive, for every v ∈ R,E exp(v(ξ n,n s ) [k] ) | Y n,h s k = 1, . . . , d.In other words, conditional onY n,h s , each component (ξ n,n s ) [k] of ξ n,n s = (ξ n,n s ) [1] , . . . , (ξ n,n s ) [d]is sub-Gaussian with variance proxy κ 9 (see e.g.[START_REF] Buldygin | Sub-Gaussian random variables[END_REF] and [Pau20, Theorem 2.1.1]). Consequently, conditional onY n,h s , each component (S n,N s ) [k] of S n,Ns is sub-Gaussian with variance proxy (N -1) -1 κ 9 , by the conditional independence of the random variables (ξ n,n s ) 1≤n ≤N,n =n . This implies in particular

  , η N (T,x) T (x) -µ T (x) 2 ≤ κ 5 max(|µ T | 2 k,K , |K| 2 L 2 , 1) × min η∈H (N η d ) -1 log N + η 2 min(k, +1) + h 2 ≤ κ 6 log N N 2 min(k, +1)2 min(k, +1)+d + h 2

s

  , µ N,h s ) -b(s, Y n,h s , μh s ) are sub-Gaussian, with the correct order in N . We may then repeat part of the proof of Proposition 19 of[START_REF] Della | Nonparametric estimation for interacting particle systems: McKean-Vlasov models[END_REF]. The smoothness assumption (19) for the nonlinear drift enables one to obtainb(s, Y n,h s , µ N,h s ) -b(s, Y n,h s , μh s ) 2 k-1 l=1 (R d ) l δ l µ b(s, Y n,h s , y l , μh s )( µ N,h s -μh s ) ⊗l (dy l )

  δ l µ b(s, Y n,h s , y l , μh s )( µ N,h s -μh s ) ⊗l (dy l ) 2p ≤ c p N -p p! for every p ≥ 1,

  by the uniform ellipticity of c thanks to Assumption 2.2 and by triangle inequality for the last inequality, with | b(t, x, y) -b(t, x, y ) ≤ | b| Lip |y -y | in particular. Recall that under Assumption 2.1, 2.2 and 2.3, there exists a unique strong solution X = (X t ) t∈[0,T ] of the McKean-Vlasov equation (2) satisfying

  LP23, Proposition 2.1]. Hence, the function (t, x) → b µ (t, x) := b(t, x, µ t ) with µ t = L(X t ) is 1 2 -Hölder continuous in t and Lipschitz continuous in x. Thus, one can obtain sup

	s∈[0,T ]

  2.10 under Assumption 2.14Again, only an extension of Proposition 4.4 is needed, the rest of the proof remains unchanged. This only amounts to show that (t,x) → b(t, x, µ t ) is C 1,3b . The smoothness in x is straightforward, and only the proof that t → b(t, x, µ t ) is C 1 b is required. By Itô's formula, we formally haved dt b(t, •, µ t ) = ∂ 1 b(t, •, µ t ) + E A t δ µ b(t, •, X t , µ t ) ,where (X t ) t∈[0,T ] is a solution to (2). Assumption 2.14 enables one to conclude that d dt b(t, •, µ t ) is continuous and bounded.


domain: shs sdv

text:
Marco Capocasa
email: marco.capocasa@uniroma1.it

Davide Venier

evolution, misinformation. 2023 Davide Venier Nutrition

Hal-03856652v2 Nutrition

*

Keywords: Nutrition, evolution, misinformation, evolution, misinformation eating habits, health, globalization, fake news, informed choices

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Human evolution is strongly connected with environmental conditions with which the subsistence strategies and the availability and variety of food are related and still largely dependent [START_REF] Boone | Subsistence Strategies and Early Human Population History: An Evolutionary Ecological Perspective[END_REF]Erlich and Erlich 2008). Over their history, human populations have occupied almost every place on the planet, experiencing different possibilities of survival which have inevitably defined various eating habits. The diversity of environmental conditions and food resources led to different selective pressures, favouring adaptations at the local level and determining genetic variation among populations [START_REF] James | Nutrition and its Role in Human Evolution[END_REF]. The relationships among environment, dietary habits and human genetic variation have been extensively investigated especially regarding the neolithic transition from a lifestyle of hunting and gathering to agro-pastoralism [START_REF] Ammerman | The Neolithic Transition and the Genetics of Populations in Europe[END_REF][START_REF] Bickle | Stable Isotopes and Dynamic Diets: The Mesolithic-Neolithic Dietary Transition in Terrestrial Central Europe[END_REF][START_REF] Schulting | Dietary Shifts at the Mesolithic-Neolithic Transition in Europe: An Overview of the Stable Isotope Data[END_REF][START_REF] Bergström | From Domestication Genomics Towards Molecular Ecology of Human Environments[END_REF]. This cultural change led to a variety of dietary changes in Neolithic populations. Particularly, the introduction of a great quantity of foods rich in starch in the daily diet of the early agropastoral communities. The digestion of starch depends on specific enzymes called amylases. One of them, alpha-amylase, is present in the saliva, while others are released at the pancreatic level. The gene responsible for the synthesis of this enzyme (AMY1) has undergone duplications during human evolution showing extensive copy number variation among populations [START_REF] Perry | Diet and the Evolution of Human Amylase Gene Copy Number Variation[END_REF]). This means that there are people having multiple copies of AMY1 in their genome, a condition contributing to a higher enzymatic availability for their digestion of starch. [START_REF] Perry | Diet and the Evolution of Human Amylase Gene Copy Number Variation[END_REF] highlighted a higher percentage of individuals carrying this genetic pattern among agro-pastoral societies than in communities of hunter-gatherers. More copies of this gene are significantly associated with higher levels of alpha-amylase [START_REF] Carpenter | Copy Number Variation of Human AMY1 is a Minor Contributor to Variation in Salivary Amylase Expression and Activity[END_REF], representing an improvement in the digestive process in populations with a starchy diet. The history of human eating habits does not depend only on the development of agricultural practices. Other most recent historical dynamics have contributed to its changing process. Firstly, the economic and social shifts linked to the transition of human societies from pre-eminently agricultural to industrial, which has taken place from the 18th century. Secondly, the increase in movements of people, information and goods linked to globalization, especially starting from the second half of the 20th century and, even more, in the new millennium. All these changes influenced lifestyles and, consequently, food choices, as they relate to the exponential growth in the productivity of the food industry [START_REF] Mendez | Globalization, Urbanization and Nutritional Change in the Developing World[END_REF].

Eating habits in the westernized world One of the most discussed effects that globalization has on eating habits regards its large-scale "westernization". This transition is not taking place at the same level everywhere, with local contexts where cultural and gastronomic identities are so marked as to curb the success of other external models [START_REF] Haile | Evolution of Human Diet and Effect of Globalization on Regional Diet with Emphasis to the Mediterranean Diet[END_REF][START_REF] García-Dorado | Economic Globalization, Nutrition and Health: A Review of Quantitative Evidence[END_REF][START_REF] Azzam | Is the World Converging to a 'Western Diet'?[END_REF]). However, a great variety of lowcost packaged food products, pivotal representatives of the "Western diet", rich in sugar, salt, and saturated fats, are available virtually everywhere [START_REF] Kopp | How Western Diet and Lifestyle Drive the Pandemic of Obesity and Civilization Diseases[END_REF]. This rapid change could play a crucial role in influencing people's health and quality of life [START_REF] Cordain | Origins and Evolution of the Western Diet: Health Implications for the 21st Century[END_REF][START_REF] García-Montero | Nutritional Components in Western Diet Versus Mediterranean Diet at the Gut Microbiota-Immune System Interplay. Implications For Health and Disease[END_REF]. Understanding these dynamics is useful to identify their effects on the human body at a biological level. As pointed out by a large body of research evidence, the Western diet is involved in determining imbalance in the gut microbiota. These investigations showed that the gut microbial diversity of huntergatherers is richer and varied than those of members of agropastoral groups and, even more, urban western communities (see [START_REF] Gupta | Geography, Ethnicity or Subsistence-Specific Variations in Human Microbiome Composition and Diversity[END_REF] and related citations therein). The main reason of this differences seems to be that a greater diversity of microbial species would confer a better resistance to parasites and efficiency in coping with the not optimal intake of nutrients due to the scarcity of food. The gut microbiota of hunter-gatherers greatly differs from that of the urban societies in terms of the richness of several genera of bacteria, such as Prevotella and Treponema, very useful for diets rich in fibres [START_REF] Obregon-Tito | Subsistence Strategies in Traditional Societies Distinguish Gut Microbiomes[END_REF]. Comparative studies on gut microbiota explain how different lifestyles and eating habits influence the bacterial composition of the human intestinal tract and the impact of its impoverishment on human health. Western diet is also strongly associated with obesity and other diseases, such as cardiovascular diseases and diabetes, which are among the ten with the highest mortality rate worldwide (GBD 2017 Causes of Death Collaborators 2018; [START_REF] Capocasa | A Light in the Dark: Open Access to Medical Literature and the COVID-19 Pandemic[END_REF][START_REF] Bisol | How Can We Get More Open Access to Medical Studies? Simple, Let's Take the Green Road[END_REF]. It is worth noting that more than half of the daily intake of calories in several countries of the western world such as Canada, the United States and the United Kingdom derives from these foods [START_REF] Monteiro | Ultra-Processed Foods: What They Are and How to Identify Them[END_REF]). In the last two decades, the sale of these products is growing dizzyingly also in low-and middle-income countries [START_REF] Monteiro | Ultra-Processed Products Are Becoming Dominant in the Global Food System[END_REF][START_REF] Moodie | Ultra-Processed Profits: The Political Economy of Countering the Global Spread of Ultra-Processed Foods -A Synthesis Review on the Market and Political Practices of Transnational Food Corporations and Strategic Public Health Responses[END_REF]. In 2004, the Food and Agriculture Organization (FAO) published a report entitled "Globalization of Food Systems in Developing Countries: Impact of Food Security and Nutrition", in which authors highlighted how unhealthy food market in developing countries was supported by the television commercials and by the spread of mobile phones and text messages. Emblematic is the case of Brazil, where 58% of TV commercials related to food products concerned those with high fat and sugar contents, only 9% regarding foods based on meat, beans, or eggs, and no commercials dedicated to fresh fruit and vegetables.

Nutrition science, between myths and lies

Television and the Internet have allowed amplification of the promotion of food products. At the same time, they represent the main tools for the dissemination of nutrition and healthy eating contents. However, scientists and nutrition professionals are still struggling to find the right way to communicate efficiently and effectively (throughout these channels) the positive implications for society that their studies could produce. This is because the "translation" of scientific works from the technical language to one suitable for a wider audience is an undertaking that hides several traps, such as oversimplification, sensationalism, and inaccuracy [START_REF] Secko | Four Models of Science Journalism: A Synthesis and Practical Assessment[END_REF]. This is particularly true in the case of nutrition science as it involves aspects of the daily life of everyone [START_REF] Rowe | Nutrition Science Communication: "If It's Transparent, Is That Why We Cannot See It?[END_REF]. Mass media not only have been used to advertise products, but also to disseminate fake news concerning the nutritional values of foods and their effects on human health (e.g. the case of the "milk consumption-hip fractures" controversial relationship; see [START_REF] Michaëlsson | Milk Intake and Risk of Mortality and Fractures in Women and Men: Cohort Studies[END_REF], but see also [START_REF] Hidayat | Systematic Review and Meta-Analysis of the Association Between Dairy Consumption and the Risk of Hip Fracture: Critical Interpretation of the Currently Available Evidence[END_REF]. However, even long before the invention of television and the worldwide web, word of mouth worked well for "food mythology" to reach a large audience. An example above all: "spinach is the richest source of iron". This is one of the longest-lived scientific legends, as it is now two centuries old. One version of the story tells about an error brought since the dawn of biochemistry, caused by confusing the iron content in the dry weight of spinach with that of fresh weight (see Rekdal 2014 and related studies cited therein). Other reconstructions refer to a misplaced decimal point in calculation by two different "guilty": Emil von Wolff and Gustav von Bunge [START_REF] Mielewczik | Spinach in Blunderland: How the Myth that Spinach is Rich in Iron Became an Urban Academic Legend[END_REF]. None of these stories would be true, however it is true that spinach is one of the vegetables that contain only a fair amount of iron (2.9 grams per 100 grams of fresh product) and that we can assimilate only a small quantity of this micronutrient due to its interactions with other dietary factors (i.e. polyphenols, phytates and calcium; [START_REF] Rodriguez-Ramiro | Estimation of the Iron Bioavailability in Green Vegetables Using an in Vitro Digestion/Caco-2 Cell Model[END_REF]). This is one of the many examples of food myths that the nutrition sciences have helped to unravel. A particularly arduous task in the globalized world, literally haunted by the proliferation of fake news mainly on the world wide web. The internet is widely used to find out information on nutrition and diets, although it cannot be view as the most reliable source [START_REF] Wangberg | Use of the Internet for Health Purposes: Trends in Norway 2000-2010[END_REF][START_REF] Goodman | Use of Nutritional Information in Canada: National Trends Between 2004 and 2008[END_REF]. [START_REF] Vosoughi | The Spread of True and False News Online[END_REF] pointed out the greater speed with which fake information spread on the web compared to truths. They also discussed the role of social networks noting how, through Twitter, fake news are 70% more likely to be shared compared to reliable news. This study clarifies how the web is a context in which false, or at least inaccurate news, can reach people needing information, leading them to make wrong choices. It is the case of the significant number of people who follow restrictive dietary patterns, such as lactose and gluten-free diets, without having performed the proper medical exams to verify if they have any intolerance (e.g. see [START_REF] Araya | Living with Gluten and Other Food Intolerances: Self-Reported Diagnoses and Management[END_REF].

Towards truly informed choices

The spread of fake news is due to the lack of scientific knowledge that is needed to distinguish which is a reliable information and which should be taken as a joke. As a matter of fact, people process information extrapolated from texts and videos differently, depending on their level of knowledge of the topics. Based on this elaboration, they may be vulnerable and become easily persuaded to make wrong and even harmful choices. The lower the knowledge, the greater the vulnerability and therefore the tendency to rely on the online pull of mixed contents, regardless of its quality. In an increasingly interconnected world, the popularity of social media contributed to make celebrities and influencers as role models. Followers may wonder how this well-known personality "has that physical aspect", desiring to resemble him/her, and thus wanting to know "what diet is he/she following", "what sporting activity is he/she practicing", or "which supplements are involved". In this way, popular personalities are seen as more reliable as experts, even if they do not have a specific scientific background. All these misconceptions bring people to follow them instead of searching for information from reliable sites, open access scientific papers or, even better, asking directly a nutritional or medical professional. How can we incentivize informed choices? First, we believe that nutrition professionals should play a primary and decisive role in providing reliable information to patients and their relatives. Particularly, we stress the importance of an accessible language to avoid technicalities that are perceived as hurdle by the patients. As pointed out by [START_REF] Mira | Barriers for an Effective Communication Around Clinical Decision Making: An Analysis of the Gaps Between Doctors' and Patients' Point of View[END_REF], communication barriers lead to a less involvement of patients and their reduced willingness to follow directions and suggestions coming from healthcare professionals. Furthermore, we claim that making a greater effort in the spread of scientific knowledge to the public, organizing public events and publishing easily readable articles on websites, blogs and magazines, could contribute to raise the awareness of the huge influence that eating habits have on human health.


Grigore Burloiu
email: grigore.burloiu@unatc.ro

Cinetic Unatc

Adaptive Drum Machine Microtiming with Transfer Learning and RNNs

We introduce rolypoly˜, the first drum machine for live performance that adapts its microtiming in relation to a human musician. We leverage state-of-the-art work in expressive performance modelling with recurrent nets, towards real-time application on the micro scale. Our models are pretrained on the Groove MIDI Dataset from Magenta, and then fine-tuned iteratively over several duet performances of a new piece. We propose a method for defining training targets based on previous performances, rather than a prior ground truth. The agent is able to adapt to human timing nuances, and can achieve effects such as morphing a rhythm from straight to swing.

INTRODUCTION

We are interested in the coupling between drum machine and human player, with a view to extending their mutual dynamics computationally [START_REF] Inverno | Heroic versus collaborative AI for the arts[END_REF]. We centre on the microtiming scale as a locus of expressive music interaction [START_REF] Leman | The expressive moment: How interaction (with music) shapes human empowerment[END_REF].

In studies on tempo and time-shift representation, [START_REF] Honing | Timing is tempo-specific[END_REF] posits that global tempo curves alone cannot account for the alterations observed in performances of the same material at different speeds. Nevertheless, score-driven automatic accompaniment has traditionally worked by computing such a curve to drive the warping of a backing track [START_REF] Raphael | Music plus one and machine learning[END_REF][START_REF] Cont | On the creative use of score following and its impact on research[END_REF]. Increasingly however, attention is also being paid to the interpretation of real-time accompaniment on the micro scale [START_REF] Xia | Spectral learning for expressive interactive ensemble music performance[END_REF][START_REF] Maezawa | Deep linear autoregressive model for interpretable prediction of expressive tempo[END_REF]. By leveraging advances in natural language processing, recurrent neural network (RNN)based machine learning architectures have been shown to produce state-of-the-art expressive outputs [START_REF] Jeong | Virtuosonet: A hierarchical attention rnn for generating expressive piano performance from music score[END_REF][START_REF] Gillick | Learning to groove with inverse sequence transformations[END_REF][START_REF] Oore | This time with feeling: Learning expressive musical performance[END_REF].

We delineate a set of design principles for a new musical agent. rolypoly˜is a score-driven drum machine with the following characteristics: lightweight (inference must be fast enough to run in real-time on an average, accessible computer), audio-interactive (must not only adapt to incoming sound from human musicians, but also "listen" to its own output, coupled with the former), progressive/scalable (must not rely on existing duet/ensemble ground truth; rather, it builds up a performance corpus, accruing learning along the way), and score-agnostic (does not require symbolic specification of the parts played by other musicians).

AGENT ARCHITECTURE

The piece to be interpreted is represented as a sequence of feature vector rows, encoding each drum hit and its context (tempo, time signature, beat phase) at timestep t, and the drum-target onset distance dif f t-1 (since this distance can only be measured in hindsight). The output is a drum microtiming offset, y t , with its corresponding estimation ŷt determining when in relation to the absolute notated time the t-th drum hit will actually be triggered.

Since we do not rely on an existing drum-human duet ground truth (no such dataset suitable for deep learning exists), we define y after and as a function of a performance, as follows.

We define a variable d t at timestep t as the cumulated realised offsets of score-to-drum and (variance-adjusted) drum-to-target:

d t = ŷt + A • dif f t σ dif f /σ ŷ , (1)
We then use d t to determine the ground truth drum offsets for training the next iteration of the model, by subtracting its mean (to keep outputs centred around zero) and again applying deviation normalisation:

y t = B • d t -µ d σ d /σ ŷ . (2)
A and B above are hyperparameters, controlling the weighting of the target offsets and the cumulated offsets, respectively. Their default setting is 1, allowing the timing output to adapt gradually, without major fluctuations. For rhythm morphing (see Section 3) they can be set to cancel out the effect of variance scaling.

We propose two RNN-derived architectures that predict each following drum timing based on the input features received in real time. Both are defined and trained using the PyTorch library, and communicate via OSC to Max, which handles the audio playback and analysis.

The first model is built on a 2-layer unidirectional LSTM network [START_REF] Hochreiter | Long short-term memory[END_REF] with 256 hidden units fed into tanh nonlinear activations. The result is passed through a linear layer with a single output, ŷ. The second model is a simplification of the Seq2Seq [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] architecture described in [START_REF] Gillick | Learning to groove with inverse sequence transformations[END_REF]. In our case the source sequence is the complete pre-performance input dataset (sans the unrealised audio-related features), and the target is the performance dataset up to the current timestep. The encoder is a bidirectional LSTM, and the decoder is a 2-layer unidirectional LSTM with 256 hidden units. As with the first model, a tanh nonlinearity and final linear layer project the decoder output to a one-dimensional activation, ŷ.

To pretrain our models we process a drums-only performance dataset, the Groove MIDI Dataset (GMD) from Magenta [START_REF] Gillick | Learning to groove with inverse sequence transformations[END_REF], the largest existing dataset of expressive drumming. We predict residual drum offsets, 1 resulting in audio-agnostic expressive interpreter models, ready to be fine-tuned with subsequent performances.

The source code and trained models 2 , notebook 3 and a demo video 4 are available

EXPERIMENTS

Figure 1 pictures the Seq2Seq model transitioning from a straight 4/4 beat to a "swing" shuffle, where offbeats are pushed slightly later-simply by "swinging" the respective notes on guitar for three iterations. Similarly we were able to morph the pattern x..x..x. to three equallydistanced triplets, as seen in the demo video.

The typical rolypoly˜use case consists in a song being performed multiple times and the model learning incrementally after each take. Heuristically we found the agent is able to limit the drum-target variance, and adapt to structural patterns in the piece.

FUTURE WORK

We are exploring improvements to e.g. the target audio representation, via the learning of a latent feature space of

Figure 1 .

 1 Figure 1. Rhythm morphing over three training iterations. Score: dotted lines. ŷ: blue stems. dif f : grey bars. Transition from straight (bottom) to swinging (top).

Thus, y measures the distance to the drum hit from its quantised position. While[START_REF] Gillick | Learning to groove with inverse sequence transformations[END_REF][START_REF] Makris | Conditional neural sequence learners for generating drums' rhythms[END_REF] used 16 steps per bar, we chose a quantisation step

of 24, to better account for triplets and swing.[START_REF] Leman | The expressive moment: How interaction (with music) shapes human empowerment[END_REF] See https://github.com/RVirmoors/rolypoly/ tree/master/py.