J-B Brunel

Quentin Clairon

Florence D'alché-Buc

Parametric Estimation of Ordinary Dierential Equations with Orthogonality Conditions

Keywords: Gradient Matching, Nonparametric statistics, Methods of Moments, Plug-in Property, 1 Variational formulation, Sobolev Space

Dierential equations are commonly used to model dynamical deterministic systems in applications. When statistical parameter estimation is required to calibrate theoretical models to data, classical statistical estimators are often confronted to complex and potentially ill-posed optimization problem. As a consequence, alternative estimators to classical parametric estimators are needed for obtaining reliable estimates. We propose a gradient matching approach for the estimation of parametric Ordinary Dierential Equations observed with noise. Starting from a nonparametric proxy of a true solution of the ODE, we build a parametric estimator based on a variational characterization of the solution. As a Generalized Moment Estimator, our estimator must satisfy a set of orthogonal conditions that are solved in the least squares sense. Despite the use of a nonparametric estimator, we prove the root-n consistency and asymptotic normality of the Orthogonal Conditions estimator. We can derive condence sets thanks to a closed-form expression for the asymptotic variance. Finally, the OC estimator is compared to classical estimators in several (simulated and real) experiments and ODE models in order to show its versatility and relevance with respect to classical Gradient Matching and Nonlinear Least Squares estimators. In particular, we show on a real dataset of inuenza infection that the approach gives reliable estimates. Moreover, we show that our approach can deal directly with more elaborated models such as Delay Dierential Equation (DDE).

1 Introduction

Problem position and motivations

Dierential Equations are a standard mathematical framework for modeling dynamics in physics, chemistry, biology, engineering sciences, etc and have proved their eciency in describing the real world. Such models are dened thanks to a time-dependent vector eld f , dened on the state-space X ⊂ R d and that depends on a parameter θ ∈ Θ ⊂ R p , d, p ≥ 1. The vector eld is then a function from [0, 1]×X ×Θ to R d . If φ(t) is the current state of the system, the time evolution is given by the following Ordinary Dierential Equation, dened for t ∈ [0, 1] by: φ(t) = f (t, φ(t), θ) (1.1) where dot indicates derivative with respect to time. An important task is then the estimation of the parameter θ from real data. [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF] proposed a signicant improvement to this statistical problem, and gave motivations for further statistical studies. We are interested in the denition and in the optimality of a statistical procedure for the estimation of the parameter θ from noisy observations y 1 , . . . , y n ∈ R d of a solution at times t 1 < • • • < t n .

Most works deal with Initial Value Problems (IVP), i.e. with ODE models having a given (possibly unknown) initial value φ(0) = φ 0 . There exists then a unique solution φ(•, φ 0 , θ) to the ODE (1.1) dened on the interval [0, 1], that depends smoothly on φ 0 and θ.

The estimation of θ is a classical problem of nonlinear regression, where we regress y on the time t.

If φ 0 is known, the Nonlinear Least Square Estimator θNLS (NLSE) is obtained by minimizing

Q LS n (θ) = n i=1 |y i -φ(t i , φ 0 , θ)| 2 (1.2)
where |•| is the classical Euclidean norm. The NLSE, Maximum Likelihood Estimator or more general M-estimators [START_REF] Van De Geer | Empirical processes in M-estimation[END_REF] are commonly used because of their good statistical properties (root-n consistency, asymptotic eciency), but they come with important computational diculties (repeated ODE integrations and presence of multiple local minima) that can decrease their interest. We refer to [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF] for a detailed overview of the previous works in this eld. An adapted NLS estimator (dedicated the specic diculties of ODEs) is also introduced and studied in [START_REF] Xue | Sieve estimation of constant and time-varying coecients in nonlinear ordinary dierential equation models by considering both numerical error and measurement error[END_REF].

Global optimization methods are then often used, such as simulated annealing, evolutionary algorithms ( [START_REF] Moles | Parameter estimation in biochemical pathways: a comparison of global optimization methods[END_REF] for a comparison of such methods). Other classical estimators are obtained by interpreting noisy ODEs as state-space models: ltering and smoothing techniques can be used for parameter inference [START_REF] Cappé | Inference in Hidden Markov Models[END_REF], which can provide estimates with reduced computational complexity [START_REF] Quach | Estimating parameters and hidden variables in non-linear state-space models based on odes for biological networks inference[END_REF][START_REF] Ionides | Inference for nonlinear dynamical systems[END_REF][START_REF] Ionides | Iterated ltering[END_REF].

Nevertheless, the diculty of the optimization problem is the outward sign of the ill-posedness of the inverse problem of ODE parameter estimation, [START_REF] Hein W Engl | Inverse problems in systems biology[END_REF]. Hence some improvements on classical estimation have been proposed by adding regularization constraints in an appropriate way.

Starting from dierent methods used for solving ODEs, dierent estimators can be developed based on a mixture of nonparametric estimation and collocation approximation. This gives rise to Gradient Matching (or Two-Step) estimators that consists in approximating the solution φ with a basis expansion {B 1 , . . . , B K }, such as cubic splines. The rationale is to estimate nonparametrically the solution φ by φ = L k=1 ĉk B k so that we can also estimate the derivative φ. An estimator of θ can be obtained by looking for the parameter that makes φ satisfy the dierential equation (1.1) in the best possible manner.

Two dierent methods have been proposed, based on a L 2 distance between φ and f (t, φ, θ): The rst one, called the two-step method, was originally proposed by [START_REF] Varah | A spline least squares method for numerical parameter estimation in dierential equations[END_REF], and has been particularly developed in (bio)chemical engineering [START_REF] Madar | Incorporating prior knowledge in cubic spline approximation -application to the identication of reaction kinetic models[END_REF][START_REF] Voit | Decoupling dynamical systems for pathway identication from metabolic proles[END_REF][START_REF] Poyton | Parameter estimation in continuous-time dynamic models using principal dierential analysis[END_REF]. It avoids the numerical integration of the ODE and usually gives rise to simple optimization program and fast procedures that usually performs well in practice.

The statistical properties of this two stage estimator (and several variants) have been studied in order to understand the inuence of nonparametric techniques to estimate a nite dimensional parameter [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF][START_REF] Liang | Parameter estimation for dierential equation models using a framework of measurement error in regression models[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF]. While keeping the same kind of numerical approximation of the solution, [START_REF] Ramsay | Parameter estimation for dierential equations: A generalized smoothing approach[END_REF] proposed a second method based on the generalized smoothing approach for determining at the same time the parameter θ and the nonparametric estimation φ. The essential dierence between these two approaches is that the nonparametric estimator in the generalized smoothing approach is computed adaptively with respect to the parametric model, whereas two-step estimators are model-free smoothing.

We introduce here a new estimator that can be seen as an improvement and a generalization of the previous two-step estimators. It uses also a nonparametric proxy φ, but we modify the criterion used to identify the ODE parameter (i.e. the second step). The initial motivations are

• to get a closed-form expression for the asymptotic variance and condence sets,

• to reduce sensitivity to the estimation of the derivative in Gradient Matching approaches,

• to take into account explicitly time-dependent vector eld, with potential discontinuities in time.

The most notable feature of the proposed method is the use of a variational formulation of the dierential equations instead of the classical point-wise one, in order to generate conditions to satisfy. This formulation is rather general and can cover a greater number of situations: we come up with a generic class of estimator of Dierential Equations (e.g Ordinary, Delay, Partial, Dierential-Algebraic), that can incorporate relatively easily prior knowledge about the true solution. In addition to the versatility of the method, the criterion is built in order to oer computational tractability, that implies that we can give a precise description of the asymptotics and give the bias and variance of the estimator. We also give a way to ameliorate adaptively our estimator and to compute asymptotic condence intervals.

First, we introduce the statistical ODE-based model and main assumptions, we motivate and describe our estimator, and show its consistency. Then, we provide a detailed description of the asymptotics, by proving its root-n consistency and asymptotic normality. Based on the asymptotic approximation, we give a closed-form expression of the asymptotic variance, and we address the problem of obtaining the best variance through the choice of an appropriate weighting matrix. Finally, we provide some insights into the practical behavior of the estimator through simulations and by considering two realdata examples. The objective of the experiments parts is to show the interest of OC with respect to the nonlinear least squares and classical gradient matching estimators.

Examples

We motivate our work in detail by presenting two models that are relatively common and simple but that nevertheless causes diculties for estimation.

Ricatti ODE

The (scalar) Ricatti equation is dened by a quadratic vector eld f (t, x) = a(t) ). The explosions have to be handled in estimation algorithms and this slows down the exploration of the parameter space (which can be dicult for high-dimensional state or parameter spaces). Nevertheless, we show in the experiment part that NLS or Gradient Matching can do well for parameter estimation, but some additional diculties does appear when the time-dependent function c(•) has abrupt changes. We consider the case where c(t

) = c √ t -d 1 [Tr,T ] , T r is a change-point time,
with d > 0. This situation is classical (e.g in engineering) where some input variables t → u(t) modify the evolution of the system φ = f (t, φ(t)) + u(t) (typically it can be the introduction of a new chemical species in a reactor at time T r ), see gure 1.1. The Cauchy-Lipschitz theory for existence and uniqueness of solutions to time-discontinuous ODE is extended straightforwardly with measure theoretic arguments [START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF]. The (generalized) solution is dened almost everywhere and belongs to a Sobolev space. For sake of completeness, we provide a generalized version of the Cauchy-Lipschitz theorem for IVP in algorithms due its intrinsic complexity. Nevertheless, the following DDE is commonly acknowledged as a correct model [START_REF] Ellner | Dynamic Models in Biology[END_REF]23]:

Ṅ = P N (t -τ ) exp (-N (t -τ )/N 0 ) -δN (t) (1.3)
whose parameter tting (of P, N 0 , δ) remains delicate. In particular classical NLS are dicult to use in this setting as the initial condition, which is a function dened on [-τ, 0], is unknown. Alternative solutions, such as Gradient Matching or Bayesian Methods (based on ABC, [START_REF] Wood | Statistical inference for noisy nonlinear ecological dynamic systems[END_REF]) give reliable estimates that reproduce the observed dynamics without estimation of the initial condition. These aforementioned methods use particular statistics or functions of the model that provides high-level information on the parameters. The Orthogonal Conditions estimator has a similar approach for dealing with the estimation of Dierential Equations.

2 Dierential Equation Model and Gradient Matching

ODE models and Gradient Matching

For ease of readability, we focus on a two-dimensional system of ODEs. In our case, as there is no computational and theoretical dierences between the situation d = 2 and d >

Y i = φ * (t i ) + i (2.1)
where 1 , . . . , n are i.i.d with E( i ) = 0 and V ( i ) = σ 2 I 2 . We suppose that the regression function φ * belongs to the Sobolev space G (φ * (0), φ * (1)) = 0 with G a given function, is one of the most common and important one, as it arises in numerous applications (physics, control theory,. . . ). We emphasize that a convenient way to deal theoretically and computationally with BVP, in particular linear second order dierential ODEs, is not based on an adaptation of the IVP theory, but it rather involves elaborated concepts from functional analysis such as weak derivative, variational formulation and Sobolev spaces [START_REF] Conway | A course in functional analysis[END_REF]. If we denote the inner product of L 2 as ∀ϕ, ψ ∈ L 2 ([0, 1]) , ϕ, ψ = ´1 0 ϕ(t)ψ(t)dt, the weak derivative of the function g in H 1

H 1 = {u ∈ L 2 ([0, 1]) | u ∈ L 2 ([0, 1])},
where f = (f 1 , f 2 ) is a vector eld from [0, 1] × X × Θ to R
is not dened point-wise but as the function ġ ∈ L 2 satisfying ġ, ϕ = -g, φ , for all function ϕ in

C 1 with support included in ]0, 1[ (denoted C 1 C (]0, 1[)). Of course, if t → φ (t, x 0 , θ) is a C 1 function on
]0, 1[, the classical derivative φ is also the weak derivative. We introduce then the (weak) variational formulation of the ODE (1.1): a weak solution g to (1.1) is a function in

H 1 such that ∀ϕ ∈ C 1 C (]0, 1[) ˆ1 0 f (t, g(t), θ)ϕ(t)dt + ˆ1 0 g(t) φ(t)dt = 0 (2.3)
This variational formulation is the key of the Finite Elements Method which is the reference approach for solving Boundary Value Problems and Partial Dierential Equations, [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. In the case of ODEs, this formulation is not well used for computing solutions, because the geometry of the (1-D) interval ]0, 1[ is simple, and it is easy to build a spline approximation by collocation that solves approximately the ODE.

Nevertheless, the characterization (2.3) is useful for the statistical inference task, as it enables to give necessary conditions for a good estimate. In particular, we emphasize that we do not solve the ODE, but we want to identify a parameter θ indexing the vector eld f . Hence, we develop a new algorithmic approach, dierent from the one used for solving the direct problem.

Denition

We dene a new gradient matching estimator based on (2.3): starting from a nonparametric estimator φ, computed from the observations (t i , y i ), i = 1, . . . , n, we want to nd the parameter θ that minimizes the discrepancy between the parametric derivative t → f t, φ(t), θ and a nonparametric estimate of the derivative, e.g. φ. A classical discrepancy measure is the L 2 distance, that gives rise to the two-step estimator θTS dened as θTS = arg min θ∈Θ R n,w (θ) where

R n,w (θ) = ˆ1 0 | φ(t) -f t, φ(t), θ | 2 w(t)dt. (2.4)
This estimator is consistent for several usual nonparametric estimators [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF][START_REF] Liang | Parameter estimation for dierential equation models using a framework of measurement error in regression models[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF], but the use of a positive weight function w vanishing at the boundaries (w(0) = w(1) = 0) is needed to get the classical parametric root-n rate. The importance of the weight function w for the asymptotics of θTS is assessed by theorem 3.1 in [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF]. Indeed, if w does not vanish at the boundaries, then θTS does not have a root-n rate, because the asymptotics is then dominated by the nonparametric estimates φ(0) and φ(1).

The usefulness of such weighting function is well acknowledged in nonparametric or semiparametric estimation. For instance, the so-called weighted average derivative is based on a similar weight function in order to get estimators with parametric rate in partial index models [START_REF] Newey | Eciency of weighted average derivative estimators and index models[END_REF].

The use of a nonparametric proxy (instead of a solution to be computed) gives the opportunity to consider parameter estimation in f 1 and in f 2 separately. For this reason and ease of readability, we consider only the estimation of the parameter θ 1 when f can be written

f (t, x, θ) = (f 1 (t, x, θ 1 ), f 2 (t, x, θ 2 )) and θ = (θ 1 , θ 2 ) (θ i ∈ R p i and p 1 + p 2 = p).
The joint estimation of θ = (θ 1 , θ 2 ) can be done by stacking the observations into a single column: there is no consequence on the asymptotics, but the estimator covariance matrix has to be slightly modied in order to take into account the correlations between the two equations f 1 and f 2 . Having said that, we write simply f = f 1 and θ = θ 1 and we consider only one

equation ẋ1 = f (t, x, θ). We use a nonparametric estimator φ = ( φ1 , φ2 ) of φ * : [0, 1] → R 2 .
Starting from (2.3), a reasonable estimator θ should satisfy the weak formulation

∀ϕ ∈ C 1 C (]0, 1[) , ˆ1 0 f t, φ(t), θ ϕ(t)dt + ˆ1 0 φ1 (t) φ(t)dt = 0.
(2.5)

The vector space C 1 C (]0, 1[) is not tractable for variational formulation, and one prefers Hilbert space with a structure related to L 2 . In our case, we use

H 1 0 = {h ∈ H 1 |h(0) = h(1)
= 0} which has a simple description within L 2 : an orthonormal basis is given by the sine functions t → √ 2 sin( πt), ≥ 1 and we have

H 1 0 = ∞ =1 a √ 2 sin ( πt) ∞ =1 2 a 2 < ∞ (2.6)
Hence, it suces to consider a countable number of orthogonal conditions (2.5) dened, for instance, with the test functions ϕ = √ 2 sin( πt), ∀ ≥ 1:

C (θ) : ˆ1 0 f t, φ(t), θ ϕ (t)dt + ˆ1 0 φ(t) φ (t)dt = 0. (2.7)
More generally, we consider a family of orthonormal functions ϕ ∈ H 1 0 , with ≥ 1, and we introduce the vector space F = span{ϕ , ≥ 1}. The vector space F may not be necessarily dense in H 1 0 , as the functions ϕ could be chosen for computational tractability or because of a natural interpretation (for instance B-splines, polynomials, wavelets, ad-hoc functions, . . . ). For this reason, we introduce the orthogonal decomposition of H 1 0 = F ⊕ F ⊥ , where F ⊥ = {g ∈ H 1 0 | g, ϕ = 0, ϕ ∈ F}, and we can have

F = H 1 0 .
In general, an estimator θ satisfying C ( θ) for ≥ 1 also approximately satises (2.5). However in practice, we will use a nite set of orthogonal constraints dened by L test functions (L > p).

In order to discuss the inuence of the choice of F and of nite dimensional subspace spanned by ϕ 1 , . . . , ϕ L we introduce the nonlinear operator E :

(g, θ) → E (g, θ), such that t → E (g, θ) (t) = f (t, g(t), θ).
For all θ in Θ and g in H 1 , the Fourier coecients of E(g, θ) -ġ in the basis (ϕ ) ≥1 are e (g, θ) = E(g, θ) -ġ, ϕ = E(g, θ), ϕ + g, φ , and we introduce the vectors in R L e L (g, θ) = (e (g, θ)) =1..L and e * L (θ) = (e (φ * , θ)) =1..L . Finally, our estimator is dened by minimizing the quadratic form

Q n,L (θ) = e L ( φ, θ) 2 : θn,L = arg min θ∈Θ Q n,L (θ).
(

θn,L is the parameter that almost vanishes the rst L Fourier coecients in the orthogonal decompo-

sition of H 1 0 = F ⊕ F ⊥ : E(g, θ) -ġ = E L (g, θ) + R L (g, θ) + E ⊥ F (g, θ) with E L (g, θ) = L =1 e (g, θ) ϕ , R L (g, θ) = >L e (g, θ) ϕ and E ⊥ F (g, θ) ∈ F ⊥ .
The function E ⊥ F (φ * , θ) represents the behavior of E(g, θ) -ġ at the boundaries of the interval

[0, 1]. As φ approaches φ * asymptotically in supremum norm, the objective function Q n,L (θ) is close to Q * L (θ) = E L (φ * , θ) 2 L 2 .
The discriminative power of Q * L (θ) can be analyzed locally around its global minimum θ * L , as it behaves approximately as the quadratic form

Q * L (θ) ≈ (θ -θ * L ) J * θ,L J * θ,L (θ -θ * L )
where J * θ,L is the matrix in R L×p with entries ´1 0 f θ j (t, φ * (t), θ * L )ϕ (t)dt, for j = 1, . . . , p, = 1, . . . , L.

Boundary Conditions and Construction of Orthogonal Conditions

The construction of the orthogonal conditions e (θ) exposed in the previous section is generic and can be proposed for numerous types of Dierential Equations, in particular for Ordinary and Delay Dierential Equations. Moreover, similar orthogonal conditions could be also derived for solutions of PDEs with a relevant set of test functions ϕ, but this extension is beyond the scope of the present paper. A process for deriving "regular" orthogonal conditions, (i.e that gives rise to root-n consistent estimator, as it is shown in section 4) is to use conditions C (θ) with an integral expression ´1 0 h t, φ(t), θ dt. The function h : (t, x, θ) -→ R must be smooth and must satisfy the remarkable identity ´1 0 h (t, φ * (t), θ * ) = 0. The variational formulation generates functions h (t, x, θ) = (f (t, x, θ) ϕ (t) -φ (t)x) whereas the classical Gradient Matching considers a single function h(t, x, y, θ) = f (t, x, θ) -y 2 ϕ(t), and the variable y is evaluated along the derivative φ(t). The asymptotic analysis shows that the dependency in y can be removed and that h behaves in fact as a function h(t, x, θ).

The OC framework then generalizes the classical TS estimator and gives ways to ameliorate it. Among other, the use of the boundary vanishing function ϕ implies an information loss close to the boundaries.

This loss can be sensible in terms of estimation quality, and should be avoided when the boundary values are known. For instance, for an IVP with known initial condition φ(0) = φ 0 , we can derive an orthogonal condition that takes into account the knowledge of φ 0 . By direct computation, we have

ˆ1 0 h(t, φ(t), θ)dt = ˆ1 0 f (t, φ(t), θ)ϕ(t)dt -[φ(1)ϕ(1) -φ(0)ϕ(0)] + ˆ1 0 φ(t) φ(t)dt.
If φ(1) is unknown, but φ(0) is known, it suces to take ϕ such that ϕ(1) = 0 and ϕ(0) = 0. The orthogonal condition still have the same expression h(t, x, θ). The same adaptation can be done when boundary values of the derivative are known (called Neumann's condition), for instance φ(1) = φ 1 is known. Indeed, the ODE gives a relationship between the second order derivative φ and the state φ, as

φ(t) = ∂ x f (t, φ, θ)f (t, φ, θ).
By choosing ϕ such that ϕ(0) = 0 and by Integration By Part, the following identity

ϕ(1)φ 1 = ˆ1 0 ∂ x f (t, φ, θ)f (t, φ, θ)ϕ(t)dt + ˆ1 0 f (t, φ, θ) φ(t)dt
gives a new condition that exploits the behavior of the solution at the boundary. Obviously, these conditions can be successfully used if the nonparametric proxy satises the boundary conditions of interest. At the contrary, it seems rather dicult to integrate such information about the boundary within the criterion R n,w (θ). The orthogonal conditions introduced in the previous section are a direct exploitation of the ODE model, and the introduction of the space F is a way to deal with the problem of the choice of the number of conditions and their type. Nevertheless, it would be useful to introduce model specic conditions h(t, φ(t), θ) which are known to have a vanishing integral for θ = θ * . Our estimator can be thought as a Generalized Method of Moments estimator, but where Moments do characterize curves and not probability distributions. A similar idea has been developed recently in the context of functional data analysis [START_REF] James | Curve alignment by moments[END_REF].

Consistency of the Orthogonal Conditions estimator

In order to obtain precise results with closed-form expression for the bias and variance estimators, we consider series estimators, i.e. estimators expressed as φj = K k=1 ĉk,j p kK = ĉj p K , where p K = (p 1K , . . . , p kK ) is a vector of approximating functions and the coecients ĉj = (ĉ k,j ) k=1..K are computed by least squares. For notational simplicity, we use the same functions (and the same number K) for estimating φ * 1 and φ * 2 . We denote P K = (p kK (t i )) 1≤i,k≤n,K the design matrix and Y j = (y i,j ) i=1..n the vectors of observations. Hence, the estimated coecients ĉj = P K P K † P K Y j (where † denotes a generalized inverse) gives rise to the so-called hat matrix H = P K P K P K † P K and the vector of smoothed observations is φj = HY j , j = 1, 2. One can typically think of regression splines, [START_REF] Ruppert | Semiparametric regression. Cambridge series on statistical and probabilistic mathematics[END_REF]. We introduce now the conditions required for the denition and consistency of our estimator.

Condition C1: (a) Θ is a compact set of R p and θ * is an interior point of Θ, X is an open subset of

R 2 ; (b) (t, x) → f (t, x, θ * ) is L 2 -Lipschitz
(t) ≤ ζ 0 (K) and K = K(n) such that ζ 0 (K) 2 K/n -→ 0 as n -→ ∞ ; (c) There are α, c 1,K , c 2,K such that φ * j -p K c j,K ∞ = sup [0,1] φ * j (t) -p K (t) c j,K = O(K -α ).
Condition C3: There exists D > 0, such that the D-neighborhood of the solution range

D = {x ∈ R 2 | ∃t ∈ [0, 1], |x -φ * (t)| < D} is included in X and f is C 2 in (x, θ) on D × Θ for t in [0, 1] a.e.
Moreover, the derivatives of f w.r.t x and θ (with obvious notations) f x , f θ , f xx , f xθ and f θθ are L 2 uniformly bounded on D × Θ by L 2 functions hx , hθ , hxθ , hxx and hθθ (respectively).

Condition C4: Let (ϕ ) ≥1 be an orthonormal sequence of C 1 functions in H 1 0 .

Condition C5:

θ * is the unique global minimizer of Q * F and inf |θ-θ * |> Q * F (θ) > 0.
Condition C6: There exists L 0 such that for L ≥ L 0 , J θ,L (g, θ) is full rank in a neighborhood of (φ * , θ * ).

Condition C1 gives the existence and uniqueness of a solution φ * in H 1 to the IVP for θ = θ * and x(0) = φ * (0). If f is continuous in t and x, then the derivative φ * (t) = f (t, φ * (t), θ * ) can be dened on ]0, 1[ and is also continuous, see appendix A. More generally, C1 does apply when there is a discontinuous input variable, such as in the Ricatti example described in section 1.2.1.

Under condition C2 (satised among others by regression splines with ζ 0 (K) = √ K), it is known that the series estimator φj are consistent estimators of φ * j for usual norms, in particular φj -φ *

j ∞ = O P ζ 0 (K) K /n + K -α (theorem 1, [ 24]
). If φ * is C s and we use splines then α = s and φ -φ

* ∞ = O P K / √ n + K 1/2-s .
Condition C3 is here to control the continuity and regularity of the function E involved in the inverse problem. Moreover, it provides uniform control needed for stochastic convergence.

Condition C4 is a sucient condition for deriving independent conditions C (θ), and normalization is useful only to avoid giving implicitly more weight to a condition w.r.t. the other conditions.

Condition C5 means that θ * is a global and isolated minima of Q * F (θ), which is standard in Mestimation [START_REF] Van Der | Asymptotic Statistics. Cambridge Series in Statistical and Probabilities Mathematics[END_REF], but can be hard to check in practice. Indeed, the parametric identiability of ODE models can be hard to show, even for small systems. No general and practical results do exist for assessing the identiability of an ODE model [START_REF] Miao | On identiability of nonlinear ode models and applications in viral dynamics[END_REF]: it is useful to discriminate between ODE identiability, statistical identiability and practical identiability. The latter being the most useful but almost impossible to check a priori. The essential meaning of condition C5 is that the addition of more and more orthogonal conditions should lead to a perfect and univocal estimation of the true parameter. From our experience and by numerical computations, we can check that Q * L (θ) has a unique minima in θ * in a region of interest, for L big enough (usually L 2 × p). The natural criterion for estimating θ and for identiability analysis is

Q * (θ) = E (φ * , θ) -E (φ * , θ * ) 2 L 2 but E ⊥ F (φ * , θ) 2
L 2 is withdrawn and we use the quadratic form Q * F (θ) in order to avoid boundary eects. This is needed in order to get a parametric rate of convergence, as in the original two-step criterion (2.4).

As a consequence, we lose a piece of information brought by the trajectory t → φ * (t) and we have to be sure that the parameter θ has a low inuence on E ⊥ F (φ * , θ)

2 L 2 .
A favorable case is that it is almost constant on Θ, so that Q * and Q * F are essentially the same functions, with the same global minimum and the same discriminating power. In practice, we can check that C5 is approximately satised by computing numerically the criterion E L φ(•, θn,L ), θ , in a neighborhood of θn,L , for L ≥ L.

Finally, Condition C6 is about the inuence of the number of test functions used. We use only the rst L Fourier coecients of E(g, θ) -ġ to identify the parameter θ, but this might not be sucient to discriminate between two parameters θ and θ . In a way, we perform dimension reduction but we need to be sure that we have an exact recovery when L goes to innity: we expect that the global minimum

θ * L of |e * L (θ)| 2 is close to the global minimum θ * of Q * F (θ) = E F (φ * , θ) 2
L 2 (found under condition C5).

We introduce the Jacobian matrices J θ,L (g, θ) in R L×p with entries ´1 0 f θ j (t, g(t), θ)ϕ (t)dt and J x,L (g, θ)

in R L×d with entries ´1 0 f x i (t, g(t), θ)ϕ (t)dt. For this reason, we suppose that J * θ,L is full rank, so that Q * L (θ) is locally strictly convex, with a unique local minimum θ * L .

The Jacobian matrix introduced in condition C6 is classical in sensitivity analysis (in ODE models). Usually, the sensitivity matrix used is the Jacobian of the least squares criterion (similar to J θ,L φ(•, θ), θ ); it enables to check a posteriori the identiability of the parameter θ. Conversely, local non-identiable parameter (sloppy parameters, [START_REF] Ryan N Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF]) can be detected in that case.

Theorem 3.1. If conditions C1 to C6 are satised, then

θn,L -θ * L = O P (1)
and the bias

B L = θ * L -θ * tends to zero as L → ∞.
In particular, if we use the sine basis and if

E (φ * , θ) is in H 1 for all θ, then B L = o 1 L .
Remark 3.1. The convergence rate of the bias B L can be rened according to the test functions ϕ : if we use B-splines, the bias is controlled by the meshsize ∆ = max j>1 (τ j -τ j-1 ) of the sequence of knots τ j , j = 1, . . . , L dening the spline spaces, see section 6 in [START_REF] Schumaker | Spline Functions: Basic Theory[END_REF].

Remark 3.2. In practice, we have B L = 0 for medium-size L, around 2 × d × p.

Asymptotics

We give a precise description of the asymptotics of θn,L (rate, variance and normality), by exploiting the well-known properties of series estimators. We consider the linear case, then we extend the obtained results to general nonlinear ODEs. We show in a preliminary step that the asymptotics of θn,L -θ * L are directly related to the behavior of e L ( φ, θ * ), which is a classical feature of Moment Estimators.

Asymptotic representation for θn -θ *

L

From the denition (2.8) of θn,L and dierentiability of f , the rst order optimality condition is J θ,L φ, θn,L e L φ, θn,L = 0

(4.1)
from which we derive an asymptotic representation for θn,L , by linearizing e L φ, θn,L around θ * L . We need to introduce the matrix-valued function dened on D×θ such that M L (g, θ) = J θ,L (g, θ) J θ,L (g, θ)

-1

J θ,L (g, θ) , and proposition 4.1 shows that M L ( φ, θn,L ) is also a consistent estimator of M * L .

Proposition 4.1. If conditions C1-C6 are satised, then

J θ,L φ, θn,L J L -1 J θ,L φ, θn,L P -→ M * L = J * θ,L J * θ,L -1 J * θ,L (4.2)
where the matrix J L is the Jacobian J θ,L evaluated at a point θ between θ * and θn,L . Moreover, we have

θn,L -θ * L = -M * L e L ( φ, θ * L ) + o P (1). (4.3)

Linear dierential equations

We consider the parametrized linear ODE dened as where A(θ) and B(θ) are matrices in R L×K with entries A ,k (θ) = ´1 0 (a(t, θ)ϕ (t) + φ (t)) p kK (t)dt

       ẋ1 = a(t, θ 1 )x 1 + b(t, θ 1 )x 2 ẋ2 = c(t, θ 2 )x 1 + d(t, θ 2 )x 2 (4.
and B ,k (θ) = ´1 0 (b(t, θ)ϕ (t)) p kK (t)dt. The gradient of e L ( φ, θ) is J θ,L φ, θ = ∂ θ A(θ)ĉ 1 + ∂ θ B(θ)ĉ 2
where ∂ θ A(θ) and ∂ θ B(θ) are straightforwardly computed by permuting dierentiation and integration.

Although e L ( φ, θ) depends linearly on the observations, we have to take care of the asymptotics as we are in a nonparametric framework and K grows with n. The behavior of linear functionals T ρ ( φ) for several nonparametric estimators (kernel regression, series estimators, orthogonal series) is well known [START_REF] Andrews | Asymptotic normality of series estimators for nonparametric and semiparametric regression models[END_REF][START_REF] Bickel | Nonparametric estimators which can be plugged-in[END_REF][START_REF] Goldstein | Optimal plug-in estimators for nonparametric functional estimation[END_REF][START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF], and in generality it can be shown that such linear forms can be estimated with the classical root-n rate and that they are asymptotically normal under quite general conditions. In the particular case of series estimators, we rely on theorem 3 of [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF] that ensures the root-n consistency and the asymptotic normality of the plugged-in estimators T ρ ( φj ), j = 1, 2 under almost minimal conditions.

We will give in the next section the precise assumptions required for root-n consistency of linear and nonlinear functional of the series estimator. Moreover, the variance of θn,L has a remarkable expression

V e,L (θ) = V e L ( φ, θ) = A(θ)V (ĉ 1 ) A(θ) + B(θ)V (ĉ 2 )B(θ) . (4.6)
We remark that there is no covariance term between ĉ1 and ĉ2 since we assume that V (Y |T = t) is diagonal (assumption C2), but in all generality, we should add 2A(θ)cov(ĉ 1 , ĉ2 )B(θ) . We can use the classical estimates of the variance of ĉ1 and ĉ2 to compute an estimate of V e,L (θ)

V e,L (θ) = A(θ) V (ĉ 1 )A(θ) + B(θ) V (ĉ 2 )B(θ) (4.7)
Thanks to proposition 4.1, we can estimate the asymptotic variance of the estimator θn,L with the con- sistent estimator M L = M L ( φ, θn,L ) and we estimate V θn,L by V θn,L = M L V e L ( φ, θn,L ) M L .

From the asymptotic normality of the plug-in estimate, we can derive condence balls with level 1 -α.

For instance, for each parameter θ i , i = 1, . . . , p:

IC(θ i ; 1 -α) = θn,L i ± q 1-α 2 V θn,L 1/2
ii where q 1-α/2 is the quantile of order 1 -α 2 of a standard Gaussian distribution. Nevertheless, we recall that these condence intervals might be aected by the bias of θn,L depending on L.

Nonlinear dierential equations

We give here general results for the asymptotics of e ( φ, θ) when the functional is linear or not in φ.

In [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF], the root-n consistency and asymptotic normality is obtained if the functional g → e (g, θ) has a continuous Fréchet derivative De (g, θ) with respect to the norm • ∞ . If x → f (t, x, θ) is twice continuously dierentiable for t ∈ [0, 1] a.e. in x and θ in Θ, then we can compute easily its Fréchet

derivative for g ∈ H 1 in the uniform ball g -φ * ∞ ≤ D . For all h ∈ H 1 such g + h -φ * ∞ ≤ D, we have e (g + h, θ) -e (g, θ) = f x (•, g, θ) h, ϕ + h, φ + h f xx (•, g, θ) h, ϕ (4.8)
by a Taylor expansion around g. As in the linear case, we introduce the tangent linear operator

A g (θ) : u → u -a g (t, θ)u with a g (t, θ) = f x 1 (t, g(t), θ) and the function b g (t, θ) = f x 2 (t, g(t), θ).
For all θ, the Fréchet derivative of e (g, θ) (w.r.t to the uniform norm) is the linear operator

h = (h 1 , h 2 ) →De (g, θ).h = h 1 , φ + a g (t, θ)ϕ + h 2 , b g (•, θ)ϕ and satises for all θ ∈ Θ |e (g + h, θ) -e (g, θ) -De (g, θ).h| ≤ C h 2 ∞ because f xx is uniformly dominated on D × Θ.
Moreover, for all (with 0 < < D), for all g, g such

that g -φ * ∞ , g -φ * ∞ ≤ , we have |De (g, θ).h -De (g , θ).h| ≤ ˆ1 0 h(t) f xx (t, g(t), θ) (g(t) -g (t)) ϕ (t)dt ≤ C h ∞ g -g ∞
with C, a constant independent of θ, and g, g (because f xx is uniformly dominated).

As in the linear case, we need to evaluate De (g, θ) on the basis p K . We denote A(g, θ) and B(g, θ)

the matrices in R L×K with entries ´1 0 a g (t, θ)ϕ (t)p kK dt and ´1 0 b g (t, θ)ϕ (t)p kK dt (respectively) and we have the approximation

e L ( φ, θ) = e L (φ * , θ) + A(φ * , θ)ĉ 1 + B(φ * , θ)ĉ 2 + O h 2 ∞ .
(4.9)

We can derive the asymptotic variance of e L ( φ, θ) from (4.9)

V e,L (θ) = A(φ * , θ)V (ĉ 1 ) A(φ * , θ) + B(φ * , θ)V (ĉ 2 )B(φ * , θ) (4.10)
and we can get an estimate V e,L (θ) from the data as in the linear case.

In order to assess the previous discussion and for deriving the root-n rate of our estimator, we introduce the following two conditions:

Condition C7: (a) The times T 1 , . . . , T n have a density π w.r.t. Lebesgue measure such 0

< c < π < C < ∞; (b) E [ 4 ] < ∞.
Condition C8: For = 1, . . . , L, θ ∈ Θ, there exists

βK in R K with fx(•,φ * ,θ)ϕ + φ π -β K p K L 2 -→ 0.
Conditions C7 and C8 are similar to the assumptions given in [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF]. Condition C8 is here to ensure that the Fréchet derivative De (φ * , θ) that drives the asymptotic rate of e (g, θ) (see equation 4.8) can be well approximated in the basis p K as the nonparametric proxy. Then the linearized nonlinear functional of the nonparametric estimator is well approximated by a linear combination of the regression coecients. When we use B-splines with uniform knot sequence, condition C8 can be replaced by the simpler condition C9:

Condition C9: (a) The series estimator is a regression spline with a uniform knot sequence (τ 1,K , . . . , τ N K ,K )

dening the spline basis p K satises max i |τ i+1,K -τ i,K | -→ 0 as K -→ ∞ ; (b) For all θ ∈ Θ, for = 1 . . . L, v : t → fx(t,φ * (t),θ)ϕ (t)+ φ (t) π(t) is C 1 .
Theorem 4.1. If either the following conditions are satised:

p K is a general series estimators Under conditions C1-C8 and if f is a linear vector eld or, f is a nonlinear vector eld and K is chosen such that ζ 0 (K) 4 K 2 n -→ 0 p K is a uniform knot splines Under conditions C1-C2(a),C3-C7,C9 and if f is a linear vector eld and K 2 n -→ 0, or f is a nonlinear vector eld and

K 4 n -→ 0
Then θn,L is such that

√ n θn,L -θ * L N (0, V * L ) (4.11)
with

V * L = M * L V * e,L M * L . (4.12)
where V * e,L = V e,L (θ * L ). The asymptotic variance can be estimated by

M L V e,L ( θn,L ) M L P -→ V * L .
In particular, if we use regression splines and

t → f (t, φ * (t), θ) is C s on [0, 1] with s ≥ 3, then (4.11) holds with K such that √ nK -s → 0 and n -1 K 4 → 0. Moreover, if L = L(n) -→ ∞, n -→ ∞ is chosen such that the bias B L(n) = O(n -1/2
), then we have

θn,L(n) -θ * = O P (n -1/2 ). (4.13)
In particular, this is the case when the test functions ϕ are the sine basis, and

L(n) = O(n α ) with α > 1/2.
This theorem is a direct application of theorem 3 in [START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF] that claims the root-n consistency and asymptotic normality of general (nonlinear) plug-in estimators. The main steps of the proof are given in Supplementary Material I.

5 Experiments

Description of the setting

We compare the NLS estimator θNLS , the Two-Step Estimator (TS) θTS and the OC estimator θOC for varying sample sizes (n = 400, 200, 50) and varying noise levels (high and small). We consider 3 dierent ODEs with dierent mathematical structure: the α-pinene ODE (linear in state and in parameter), the Ricatti ODE (nonlinear in state, linear in parameter) and the FitzHugh-Nagumo ODE (nonlinear in state and in parameter). Experiments on these three dierent models provide a good idea of the behavior of the dierent estimators in terms of the robustness, consistency and eciency. It helps also in assessing the quality of the linear approximation for the asymptotics (in particular for the computation of the covariance matrix).

In the simulations, the noise is homoscedastic and Gaussian, so that the NLS are asymptotically ecient. Hence, the settings n = 200 or n = 400 indicates the eciency loss of the Gradient Matching estimators whereas the small size setting (n = 50) gives some information on the small sample case, where the asymptotic approximations cannot be assessed.

As the standard reference method, the Sum of Squared Errors (SSE) is minimized by a Levenberg-Marquardt algorithm using 20 starting points centered around the true parameter value θ * , and we retain the best minimum. The solution of the ODE is computed by a Runge-Kutta algorithm of order 4, implemented in the Matlab function ode45. Hence, we expect that we obtain the true NLS estimator, and that the estimated variance is the true best one.

The Gradient Matching estimators (TS and OC) use the same regression spline, decomposed on a B-spline basis with a uniform knots sequence ξ k , k = 1, . . . , K. For each dataset (and each dimension), the number of knots is selected by minimizing the GCV criterion, [START_REF] Ruppert | Semiparametric regression. Cambridge series on statistical and probabilistic mathematics[END_REF]. For the plain TS estimator, we use a piecewise ane weight function with w(0) = w(1) = 0, as in [START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF].

The Orthogonal Conditions are dened with the sine basis or B-Splines basis. We have to face with the practical problem of nding the best number of conditions L, that depends on the model and on φ. In each setting, we have xed a minimum and a maximum number of conditions L min and L max ≤ 2 × d × p and we select the OC estimator θn,L that gives the smallest prediction error (i.e that minimizes the SSE): θOC = arg min

L min ≤L≤Lmax n i=1 y i -φ(t i , φ0 , θn,L ) 2
where φ0 = φ(0) is the nonparametric estimate of the initial condition.

We use Monte Carlo simulations, based on N M C = 500 independent draws for comparing the estimators. We compute their Mean Squared Errors θ -θ *

2

. The accuracy of the estimator is roughly estimated by the trace of the covariance matrices of the estimators, denoted T r V ( θ) . Moreover, the reliability of the estimates (and asymptotic approximation) is evaluated with the coverage probabilities of the 95% condence ellipse (except in the case of TS because there is no closed-form for asymptotic variance). For the NLS, the asymptotic variance is computed via the Matlab function nlint. A more detailed analysis of the experiments (including coverage probabilities of condence sets) are given in

×10 -2 M SE T r V ( θ) (n, σ)

Unknown initial condition

In this case, the NLS needs to estimate the initial condition as well, whereas it is not needed for Gradient

Matching estimators and we have the same estimates (for θTS and θOC ) as in the previous section. In this setting, we consider another OC estimator that uses information about the other boundary T = 100.

Indeed, we know that the α-pinene network converges to a stationary point, that is almost reached at time T = 100. Hence the boundary condition φ * (100) = 0 can be used for estimation (section 2.3): if ϕ 1 is a test function with ϕ 1 (100) = 0, we have A 2 < φ * , ϕ 1 > +A < φ * , φ1 >= 0. This gives an additional condition to be satised for the OC estimator, which is denoted as θOC,1 , see section 2.3.

×10 -2 M SE T r V ( θ)
(n, σ) 14) -φ(T r )) = 0 where φ is the antiderivative of ϕ.

When T r is known, we use a cubic B-splines basis with 3 knots at T r , meaning that φ can have a discontinuous derivative at time T r (hence the curve estimation from noisy data is pretty correct at T r ).

The curve is mainly at for t ∈ [0, T r ] and after T r , one can observe a linear behavior: 3 knots are used to estimate the curve, and their positions are selected manually.

When T r is unknown, it is required to estimate θ = (a, c, d , T r ).

           d dt X 1 = ρ m D m (t -τ ) -γ ms d dt X 2 = ρ s D m (t -τ ) -γ sl + γ ms e (X 1 -X 2 ) d dt X 3 = γ sl e (X 2 -X 3 ) -δ l (6.1)
is credible for representing the dynamics of the observations. Model (6.1) is written in log-scale (i.e with X 1 = log(T m E ), X 2 = log(T s E ) and X 3 = log(T l E )), and the parameter θ = (ρ m , ρ s , δ l , γ ms , γ sl ) T has to be estimated from the data. The function D and the delay are known (estimated from the data).

The available data are the variables T m E , T s E and T l E for six dierent subjects and are measured at times T = [0, [START_REF] Bellman | Orbit determination as a multi-point boundary-value problem and quasilinearization[END_REF][START_REF] Bickel | Nonparametric estimators which can be plugged-in[END_REF][START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Brigo | Interest Rate Models -Theory and Practice[END_REF][START_REF] Brunel | Parameter estimation of ode's via nonparametric estimators[END_REF][START_REF] Cappé | Inference in Hidden Markov Models[END_REF][START_REF] Conway | A course in functional analysis[END_REF][START_REF] Ellner | Dynamic Models in Biology[END_REF][START_REF] Hein W Engl | Inverse problems in systems biology[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF][START_REF] Newey | Convergence rates and asymptotic normality for series estimators[END_REF]. Following Wu et al., we stabilize the variance by a log transformation, hence we consider directly the variables X i , i = 1, 2, 3. We assume that each subject share the same true parameter θ * and the same initial conditions: at each time point, we compute the mean of the log-measurement (over the subjects) as pseudo-observations.

We estimate D m with a spline smoother computed with cubic B-Splines and GCV selection for the knots. As in Wu et al, the nonparametric proxy is a regression spline X = X 1 , X 2 , X 3 dened on [START_REF] Bickel | Nonparametric estimators which can be plugged-in[END_REF][START_REF] Gugushvili | Root-n-consistent parameter estimation for systems of ordinary dierential equations: bypassing numerical integration via smoothing[END_REF]; we do not consider earlier times since the inuenza specic CD8+ T cells are not produced before. Since we have a small number of observations, the choice of the knots for the cubic splines is done manually.

Nevertheless for the parameter estimation, we have tested several estimates X (with dierent knots locations), and dierent number of tests functions L: we selected L = 3 or L = 4. The corresponding estimators are denoted θOC 3 and θOC 4 . Moreover, in order to improve the accuracy , we have used a weighted version of the OC estimator, similar to the classical "Generalized Methods of Moments" (this procedure is detailed in section 5 of Supplementary Material I ). The quality of the estimator is evaluated by the SSE:

SSE = 6 s=1 3 d=1 N i=1 y i,d,s -φ d (t i , θ, X(0)) 2
where y i,d,s is the observation at time t i for the s-th subject for the transformed variable X d . As suggested in Wu et al, we use the OC estimates as initial values for NLS estimation. For both estimates, we obtain the same estimator which is then simply denoted as θNLS . We provide three dierent estimates θOC 3 , θOC

and θNLS ; we mention also θref , which is the estimate obtained in Wu et al [START_REF] Wu | Modeling of inuenza-specic cd8+ t cells during the primary response indicates that the spleen is a major source of eectors[END_REF].

Blowy model

The Delay Dierential Equation (1.3) was proposed by Gurney et al [START_REF] Bkythe | Nicholson's blowies revisited[END_REF] to model the dynamics of a population of blowies, from the Nicholson's blowy data [START_REF] Nicholson | The self-adjustement of population to change[END_REF]. These data consists of 350 counts taken every two days during between day 40 = T 0 and day 315 = T 1 . As Gurney did, we take τ = 14.8 days and our aim is to estimate θ = (P, N 0 , δ). The orthogonal conditions derived from the weak form is we have used only the rst 180 observations, see [START_REF] Seifu | Fitting Population Dynamic Models To Time-Series Data By Gradient Matching[END_REF]. For the nonparametric estimation, we have used

Discussion

Among the simulated models we considered (α-pinene, Ricatti), the NLS estimator is often the best estimator in the asymptotic case (and small noise case) in terms of MSE for the parameters. Nevertheless, in some complex case such as unknown initial conditions for α-pinene (with small sample size or high noise level), or Ricatti equation (with known or unknown change point T r ), then TS and OC can oer better statistical performances. The α-pinene model shows the interest of using information on the boundaries in OC (as introduced in section 2.1). Moreover, simulations show that OC can improve on classical TS although it uses only (partial information) about (weak) derivatives. The fact that the NLS can be caught up, even in the very favorable case of a closed-form solution and starting values (for NLS optimization) close to the true parameter indicates that the introduction of Functions Moments oers a competitive estimator to the direct classical for complex case. In the latter case of Ricatti, the TS approaches is uniformly better than NLS, whereas OC is not systematically better than NLS. Ricatti Equation is striking, as it shows that good proxies φ gives a lot of information: when T r is known, the reconstruction of the solution and its derivative is excellent, which gives a clear advantage to the plain TS. Nevertheless, when T r is unknown the derivative estimation is of poor quality around T r , and the TS estimator is unstable and cannot be computed. The same situation occurs for NLS, because of some lack of identiability and dramatic changes in derivative estimation which makes the optimization algorithms inecient. For the inuenza dataset analysis, the two OC estimators give correct parameter estimates from real and sparse data (the simulated ODE have a correct qualitative behavior). When used as starting for NLS, both estimates give the same NLS estimator, which improves (obviously) the we can assess the self-consistency of our estimate. Moreover, the posterior mean is always in the 95% condence set computed for OC.

x 2 + 14 1. 1 .

 2141 b(t)x + c(t) where a(•), b(•), c(•) are time-varying functions. This equation arises naturally in control theory for solving linear-quadratic control problem[START_REF] Sontag | Mathematical Control Theory: Deterministic nite-dimensional systems[END_REF], or in mathematical nance, in the analysis of stochastic interest rate models[START_REF] Brigo | Interest Rate Models -Theory and Practice[END_REF]. We consider one of the simplest case where a is constant, b = 0 and c(t) = c √ t. The objective is to estimate parameters a, c from the noisy observations y i = φ(t i ) + i for t i ∈ [0, 14]. Here the true parameters are a = 0.11, c = 0.09 and φ 0 = -1, and one can see the solution and simulated observations in gure 1.1. Although the solution φ is smooth in the parameters, there exists no closed form and simulations are required for implementing NLS and classical approaches. The hard part in this equation is due to the extreme sensitivity of the squared term in the vector eld: for small dierences in the parameters or initial condition, the solution can explode before reaching the nal time T = Explosions are not due to numerical instability but to the failure of (theoretical) existence of a global solution on the entire interval (e.g the tangent function is solution of φ = φ 2 + 1, φ(0) = 0 and explodes at t = π 2

Figure 1 . 2 :

 12 Figure 1.2: Blowy Data, collected by Nicholson

4 ) 2 ( 4

 424 where a(•, θ), b(•, θ), c(•, θ), d(•, θ) are in L 2 . We focus only on the estimation of the parameter θ = θ 1 involved in the rst equation ẋ1 = a(t, θ)x 1 + b(t, θ)x 2 and we suppose that we have two series estimators φ1 = p K ĉ1 and φ2 = p K ĉ2 satisfying condition C2. The orthogonal conditions are simple linear functionals of the estimators e ( φ, θ) = φ1 , φ + a(•, θ)ϕ + φ2 , b(•, θ)ϕ . Hence the asymptotic behavior of the empirical orthogonal conditions relies on the plug-in properties of φ1 and φ2 into the linear forms T ρ : x → ´1 0 ρ(t)x(t)dt where ρ is a smooth function. Moreover, the linearity of series estimator makes the orthogonal conditions e L ( φ, θ) easy to compute as e L ( φ, θ) = A(θ)ĉ 1 + B(θ)ĉ

Figure 6 . 1 :θOC 1 ,

 611 Figure 6.1: Inuenza model, Estimated curves for X 1 (red), X 2 (green), X 3 (blue); ×: observations, : solution for θOC 1 , •: solution for θOC 2 , solid line: solution with θNLS .

Figure 6 . 2 :

 62 Figure 6.2: Inuenza model, Estimated curves for X 1 (red), X 2 (green), X 3 (blue); solution obtained with OC+NLS, • solution obtained with θref .

∀ϕ ∈ C 1 c

 1 (]a, b[), ˆb a N (u) φ(u)du + P ˆb-τ a-τ N (u)e -N (u) N 0 ϕ(u + τ )du -δ ˆb a N (u)ϕ(u)du = 0 where [a, b] has to be chosen such that: [a, b] ,[a -τ, b -τ ] ⊂ [T 0 , T 1 ]. Due to a change in the dynamics,

Figure 6 . 3 :

 63 Figure 6.3: Solution N of the Nicholson's DDE simulated with the OC estimator (computed with L = 11 conditions -continuous red line). The NLS solution is given by the dashed green curve. The initial function is estimated between day 40 and 55 and the simulation starts after day 55. Drift between data and simulations comes from a chaotic behavior and uncertainty in initial condition (and parameters)

  2, there is no lack of generality by this assumption. We consider noisy observations Y 1 , . . . , Y n ∈ R 2 of the function φ * measured at random times t 1 < • • • < t n ∈ [0, 1]:

  2 , where X ⊂ R 2 .

	The statistical problem can be seen as a noisy version of a parametrized Multipoint Boundary-Value
	Problem (MBVP, [4]). MBVP deals with the existence, uniqueness and computation of a solution φ * to
	equation (1.1), with general boundary conditions φ * (t 1 ) = y 1 , . . . , φ * (t n ) = y n , n ≥ 2. Obviously, MBVP
	is a much more dicult problem than the classical Initial Value Problem although some theoretical
	results do exist in some restricted cases ([3, 27] and references therein). On the computational side,
	numerous algorithms such as collocation, multiple shooting,... have been proposed to solve general
	Boundary Value Problems, [2]. Among them, the 2 points Boundary Value Problem (BVP) where

  (Y i , t i ) are i.i.d. with variance V (Y |T = t) = Σ = σ 2 I 2 ; (b) For every K, there is a nonsingular constant matrix B such that for P K = B K p (t); (i) the smallest eigenvalue of E P K (T )P K (T ) is bounded away from zero uniformly in K and (ii) there is a sequence of constants ζ 0 (K) satisfying sup t P K

and L 2 -Caratheodory (see Supplementary Material I, section 1).

Condition C2: (a)

Table 5 .

 5 1: MSE, Asymptotic Variance for α-pinene model with known Initial Condition

	TS	OC	OC,0	NLS	OC	OC,0	NLS
	(400, 3) 0.72 0.05 0.04 0.02	0.04	0.04	0.02
	(400, 8) 2.28 0.22 0.25 0.10	0.95	1.20	0.12
	(200, 3) 1.19 0.27 0.30 0.03	0.09	0.13	0.03
	(200, 8) 2.95 0.44 0.37 0.18	2.66	2.68	0.27
	(50, 3) 2.39 0.27 0.26 0.16	1.37	1.58	0.16
	(50, 8) 4.54 1.03 0.93 0.68	7.96	7.27	1.68

Table 5 .

 5 2: MSE, Asymptotic Variance for α-pinene model with unknown initial conditions 5.3 Ricatti Equation The true ODE is φ = aφ 2 + c √ t -d 1 [Tr;14] , with a * = 0.11, c * = 0.09, d * = 2 and φ 0 = -1, for t ∈ [0, 14]. For all ϕ in C 1 with ϕ(0) = ϕ(14) = 0, we have < φ, φ > +a < φ 2 , ϕ > +c <

	TS	OC	OC,1	NLS	OC	OC,1	NLS
	(400, 3) 0.25 0.11 0.11 0.07	0.10	0.10	0.06
	(400, 8) 1.07 0.85 0.56 0.50	1.06	0.82	0.61
	(200, 3) 0.6 0.37 0.23 0.14	0.25	0.20	0.14
	(200, 8) 1.64 1.42 0.83 1.34	2.36	1.64	1.54
	(50, 3) 1.33 1.31 0.80 0.69	1.63	1.02	0.76
	(50, 8) 3.64 2.11 1.79 1.96	5.34	2.20	4.38
							√ t, ϕ >

Table 5 .

 5 The OC is no more linear in parameters, but θOC can be computed by solving the general nonlinear program. The Two-Step estimator fails to estimate T r because the derivative of the solution is badly estimated when T r is unknown. OC estimators still give reliable estimates as it uses only φ in the criterion. Some care has to be taken for the knots selection because of unknown T r : when n = 200, 400 we use a uniform grid of 15 knots on [0, 14]. For n = 50, we have used 8 knots uniformly located on [0, 14]. Nevertheless, the nonparametric estimates are too rough to obtaining any correct estimate θTS .Concerning NLS, we were not able to solve the optimization problem and we cannot give Monte Carlo statistics for the evaluation of NLS. NLS collapses in practice because the optimization problem is hard (severely ill-posed problem). Indeed, the Levenberg-Marquardt algorithm becomes very sensitive to initial conditions and gives dierent solutions for very close starting values, even in the neighborhood of the true value θ * . Moreover, we have to face with the problem of explosion of the solutions during the optimization process. In particular, this problem is very delicate because we have to chose (a, c) so that the (potential) explosion of the solution can be balanced by a proper choice of d and T r . Probably, NLS would benet from a specic optimization algorithm that could exploit the particular properties of the ODE, but this is out of the scope of the paper. 4: MSE, Sum Empirical Variance for Parameter estimation for Ricatti with unknown T r

	×10 -2		M SE		T r V θ
	(n, σ)	TS	OC	NLS	OC	NLS
	(400, 0.2) 0.18 0.27 0.58	1.76	0.10
	(400, 0.4) 0.78 1.21 0.94	2.56	0.38
	(200, 0.2) 0.33 0.87 0.57	2.85	0.25
	(200, 0.4) 1.12 2.69 1.12	5.64	0.98
	(50, 0.2) 1.03 1.30 1.54	4.70	1.00
	(50, 0.4) 3.80 4.43 3.94	8.89	4.08
	Table 5.3: MSE , T r V θ	for Parameter estimation for Ricatti Equation with known T r
	×10 -2	M SE( a) M SE( c) M SE( d ) M SE( T r )
	(n, σ)		OC	OC		OC	OC
	(400, 0.2)		0.09	0.00		2.54	1.39
	(400, 0.4)		0.29	0.01		4.27	3.54
	(200, 0.2)		0.21	0.00		4.08	3.18
	(200, 0.4)		0.61	0.01		11.96	6.93
	(50, 0.4)		0.64	0.02		11.20	14.25
	(50, 0.4)		0.77	0.01		17.18	19.40
			×10 -2	M SE T r V ( θ)
			(n, σ)	OC		OC
			(400, 0.2) 4.01		3.97
			(400, 0.4) 8.11		8.02
			(200, 0.2) 7.47		7.35
			(200, 0.4) 19.51	18.94
			(50, 0.2) 26.10	5.14
			(50, 0.4) 37.36	9.49
	6 Real data analysis
	6.1 Inuenza virus growth and migration model
	We consider the ODE model introduced in Wu et. al [42] for the growth and migration of inuenza
	virus-specic eector CD8+ T cells, among lymph node (T	m E ), spleen (T	s E ), and lung (T	l E ) of mice. After
	a model selection process, it turns out that the following model

Table 6 .

 6

				θOC 3	θOC 4	θNLS	θref
			ρ m ρ s δ l γ ms γ sl RMSE 13.5 2.9e-5 2.7e-5 1.5e-5 1.6e -5 4.1e-5 4.7e-5 4.1e-5 4.5e -5 2.0 3.4 3.96 3.7 0.39 0.35 0.157 0.15 0.72 0.81 0.49 0.47 13.9 9.0 9.5
		θOC 3			θOC 4			θNLS
		Low. Bound	Up. Bound	Low. Bound	Up. Bound	Low. Bound	Up. Bound
	ρ m ρ s δ l γ ms γ sl	2.1e-5 0.7e-5 -1.11 0.27 -0.10	3.7e-5 7.4e-5 5.21 0.50 1.55		1.9e-5 0.9e-5 -0.28 0.24 -0.14	3.4e-5 8.4e-5 7.21 0.46 1.76		0.7e-0.5 3.4e-0.5 2.59 0.03 0.39	2.4e-0.5 4.8e-0.5 4.93 0.26 0.55

1: Estimates, RMSE and the 95% condence intervals for dierent L and estimators.

Table 6 .

 6 42 knots located between t = 40 and t = 220. Preliminary tests and comparisons suggests to use the sine basis for the test function ϕ , and we use 2 ≤ L ≤ 15. A simulation is given in gure 6.3 .C Low. Bound Up. Bound Low. Bound Up. Bound Low. Bound Up. Bound 2: Estimates, RSSE and 95% condence intervals for dierent L

				L = 11	L = 9	L = 12
			P	7.81	7.52	7.91
			N 0	381.8	385.9	377.7
			δ	0.154	0.153	0.154
			RSSE 1.7136e+03 1.7557e+03 1.7990e+03
		L = 11			L = 9	L = 12
	P N 0 δ	5.80 303.62 0.10	9.81 459.94 0.20	5.64 306.59 0.11	9.40 465.38 0.19	5.0416 289.36 0.10	10.77 465.98 0.20

O

Acknowledgments

This work was funded by two ANR projects GD2S (ANR-05-MMSA-0013-01) and ODESSA (ANR-09-SYSC-009-01) and received a partial support from the Analysis of Object Oriented Data program (2010-2011) in Statistical and Applied Mathematical Sciences Institute (SAMSI), USA.

Supplementary Material I. This abrupt change causes some diculties in estimating non-parametrically the solution and its derivative, which can make Gradient Matching less precise. We consider then the estimation of the two additional parameters d and T r . Hence, the parameter estimation problem can be seen as a change-point detection problem, where the solution φ still depends smoothly in the parameters. Nevertheless, in the case of the joint estimation of a, c, d and T r , the particular inuence of the parameter T r makes the problem much more dicult to deal with for classical approaches as it is suggested by the objective functions in Supplementary Material II. The variational formulation for model estimation gives a seamless approach for estimating models which possess time discontinuities.

Dynamics of Blowy populations

The modeling of the dynamics of population is a classical topic in ecology an more generally in biology.

Dierential Equations can describe very precisely the mechanics of evolution, with birth, death and migration eects. The case of single-species models is the easiest case to consider, as interactions with rest of the world can be limited, and the acquisition of reliable data is easier. In the 50s, Nicholson measured quite precisely the dynamics of a blowy population, known as Nicholson's experiments [START_REF] Nicholson | The self-adjustement of population to change[END_REF].

The data are relatively hard to model, and it is common to use Delay Dierential Equation (DDE) whose general form is Ṅ (t) = f (N (t), N (t -τ ), θ), in order to account for the almost chaotic behavior Supplementary Materials II: Experiments, Tables and Figures.

α-pinene

A linear ODE with constant coecients is written ẋ = Ax, where A = (A 1 | . . . |A d ). For i = 1, . . . , d, the weak formulation gives the identity Y ϕ i = X ϕ A i to be satised, where X ϕ is a d × L matrix with entries x k , ϕ andY ϕ i is a vector in R L with entries equal to -x i , φ . For illustration, we consider the α-pinene ODE used in [START_REF] Moles | Parameter estimation in biochemical pathways: a comparison of global optimization methods[END_REF] for the comparison of several global optimization algorithms:

The true parameter to be estimated from a completely observed trajectory on [0, 100

As this ODE is linear and time-invariant, we have a closed-form for the solution φ * (t, θ, φ 0 ) = e tA φ 0 that can be directly used for the computation of the NLS estimator.

The test functions used for the OC estimators are B-Splines (with uniform knots sequence) ϕ , = 1, . . . , L with compact support included in ]0, 20[. We consider a varying number of conditions L, i.e 2 ≤ L ≤ 15. Finally, we have two settings for the estimation of θ: when the initial condition φ 0 is known (and equal to (100, 0, 0, 0, 0) as in [START_REF] Rodriguez-Fernandez | Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems[END_REF]), and when φ 0 is unknown and needs to be estimated (for NLS).

Known initial condition

For the OC and TS estimator, we constrain the spline estimator φ to satisfy the condition φ(0) = φ 0 (by adding a linear constraint to the classical least-squares minimization). Moreover, following section 2.3, we integrate the knowledge of the initial condition by adding a test function ϕ 0 which is a B-spline with ϕ 0 (0) = 0. Hence, we dene 2 dierents OC estimators, respectively, θOC,0 and θOC that uses or not (resp.) the knowledge of the initial condition.


Damien Dosimont

Robin Lamarche-Perrin

Lucas Mello Schnorr

Guillaume Huard

Jean-Marc Vincent

Combining Data and Visual Aggregation Techniques to Build a Coherent Spatiotemporal Overview

published or not. The documents may come

L'archive ouverte

Combining Data and Visual Aggregation Techniques to Build a Coherent Spatiotemporal Overview

Damien Dosimont * † ‡ § , Robin Lamarche-Perrin ¶ , Lucas Mello Schnorr , Guillaume Huard † ‡ * § and Jean-Marc Vincent † ‡ * § * Inria † Univ. Grenoble Alpes, LIG, F-38000 Grenoble, France ‡ CNRS, LIG, F-38000 Grenoble, France § Email: firstname.lastname@imag.fr ¶ MPI for Mathematics in the Sciences, 04103 Leipzig, Germany -Email: robin.lamarche-perrin@mis.mpg.de Informatics Institute, UFRGS, Porto Alegre -Email: schnorr@inf.ufrgs.br This article has been rejected from VPA'14 proceedings at the last minute, because of a copyright conflict. Please contact the authors for further information.


Célia Da Costa

Andrea G B Pereira

Tettamanzi

Célia Da

Costa Pereira

Andrea G B Tettamanzi

Syntactic Possibilistic Goal Generation

Introduction

Researchers in the field of AI are growing more and more aware of the importance of integrating the goal generation process when representing the agent's capabilities [START_REF] Hanheide | A framework for goal generation and management[END_REF]. Indeed, the idea of fixing a priori the goals that an agent must achieve is acceptable only for few domains where all possible goals can be determined in advance. On the other hand, a suitable framework for representing the agent's capabilities should consider the fact that an agent must also choose which goals to pursue from the collection of goals which have been generated previously [START_REF] Pokahr | A goal deliberation strategy for bdi agent systems[END_REF][START_REF] Thangarajah | A soft COP model for goal deliberation in a BDI agent[END_REF][START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF].

In most real world situations, information available for an agent is incomplete. This means that the set of the agent's beliefs represents an imprecise description of the real world. Possibility theory [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF] is well-suited for modeling uncertain or vague information by means of a possibility distribution. In [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF], the authors used such a distribution to propose a semantic representation of beliefs allowing thus to express that some worlds (interpretations) are more plausible for a BDI agent than others. As pointed out by Dubois et al. [START_REF] Dubois | Knowledge-driven versus datadriven logics[END_REF], this semantic view of beliefs is also the one developed in the theory of belief change by Gärdenfors [START_REF] Gärdenfors | Knowledge in Flux: Modeling the Dynamics of Epistemic States[END_REF]. In that theory, the logical consequences of the belief base [START_REF] Nebel | A knowledge level analysis of belief revision[END_REF] represents the full set of the agent's beliefsits belief set. In this semantic setting, adding new information to a belief set comes down to some worlds which were previously possible becoming impossible. This means that the more information is available, the smaller the set of worlds considered possible by the agent and the more precise the information held by the agent. It is important to observe that, while this kind of semantic representation is well-suited to a theoretical treatment of the matter, it is not at all adapted to the implementation of an agent framework. This is why, to use a syntactic possibilistic representation of beliefs, in [START_REF] Da | A syntactic possibilistic belief change operator for cognitive agents[END_REF] the authors propose an equivalent syntactic possibilistic belief-change operator.

An algorithm for generating the goals of an agent based on the semantic representation of beliefs and desires was proposed in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF]. However, to the best of our knowledge, no one has yet proposed an algorithm for goal generation based on the syntactic counterparts of the possibilistic representations of beliefs and desires. This paper aims at bridging this gap and at making an efficient implementation of the agent's goal generation and adoption possible. We provide the syntactic version of the possibilistic deliberation algorithm and of the goal-election algorithm proposed in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF].

The rest of the paper is organized as follows. Section 2 presents some basic definitions in possibility theory. Section 3 presents the representation of graded beliefs both in the syntactic and in the semantic settings. Section 4 presents the components of the possibilistic model. Section 5 presents a syntactic representation of desires while in Section 6 an algorithm for generating such kind of desires is proposed. Section 7 presents the two algorithms used for selecting as goals the maximally justified desires among the maximally possible desires. Finally, Section 8 concludes the paper.

Background

In this section, we present some basic definitions in possibility theory that will be used throughout the paper.

Language and Interpretations

We adopt a classical propositional language to develop the theoretical framework used to represent information manipulated by an agent.

Definition 1 (Language) Let A be a finite set of atomic propositions and let L be the propositional language such that A∪{⊤, ⊥} ⊆ L, and, ∀φ,

ψ ∈ L, ¬φ ∈ L, φ ∧ ψ ∈ L, φ ∨ ψ ∈ L.
Additional connectives can be defined as useful shorthands for combination of connectives of L, e.g., φ ⊃ ψ ≡ ¬φ ∨ ψ.

We will denote by Ω = {0, 1} A the set of all interpretations on A. An interpretation ω ∈ Ω is a function ω : A → {0, 1} assigning a truth value p ω to every atomic proposition p ∈ A and, by extension, a truth value φ ω to all formulas φ ∈ L; ω |= φ means that

φ ω = 1 (ω is a model of φ); if S ⊆ L is a set of formulas, ω |= S means ω |= φ for all φ ∈ S; S |= φ means that ∀ω |= S, ω |= φ.
The notation [φ] denotes the set of all models of formula φ ∈ L:

[φ] = {ω ∈ Ω : ω |= φ}. Likewise, if S ⊆ L is a set of formulas, [S] = {ω ∈ Ω : ∀φ ∈ S, ω |= φ} = φ∈S [φ].

Possibility Theory

Fuzzy sets [START_REF] Zadeh | Fuzzy sets[END_REF] are sets whose elements have degrees of membership in [0, 1]. Possibility theory is a mathematical theory of uncertainty that relies upon fuzzy set theory, in that the (fuzzy) set of possible values for a variable of interest is used to describe the uncertainty as to its precise value. At the semantic level, the membership function of such set, π, is called a possibility distribution and its range is [0, 1]. By convention, π(ω) = 1 means that it is totally possible for ω to be the real world, 1 > π(ω) > 0 means that ω is only somehow possible, while π(ω) = 0 means that ω is certainly not the real world.

A possibility distribution π is said to be normalized if there exists at least one interpretation ω0 s.t. π(ω0) = 1, i.e., there exists at least one possible situation which is consistent with the available knowledge.

Definition 2 (Measures) A possibility distribution π induces a possibility measure Π, its dual necessity measure N , and a guaranteed possibility measure ∆. They all apply to a classical set A ⊆ Ω and are defined as follows:

Π(A) = max ω∈A π(ω); (1)
N (A) = 1 -Π( Ā) = min ω∈ Ā{1 -π(ω)}; (2)
∆(A) = min ω∈A π(ω). (3)

Syntactic and Semantic Representations

A possibilistic belief base is a finite set of weighted formulas B = {(φi, αi), i = 1, . . . , n}, where αi is understood as a lower bound of the degree of necessity N ([φi]) (i.e., N ([φi]) ≥ αi). Here, B(φi) = αi means that the degree to which formula φi belongs to the set B is αi. Formulas with αi = 0 are not explicitly represented in the belief base, i.e., only a belief which is somehow believed/accepted by the agent is explicitly represented. The higher the weight, the more certain the formula. From belief base B the degree of belief B(φ) of any arbitrary formula φ ∈ L may be computed as follows [START_REF] Benferhat | A practical approach to revising prioritized knowledge bases[END_REF]:

B(φ) = max{α : Bα |= φ}, (4)
where Bα = {φ : B(φ) ≥ α}, with α ∈ [0, 1], is called an α-cut of B. The meaning of Equation 4 is that the degree to which an agent believes φ is given by the maximal degree α such that φ is entailed only by the formulas whose degree of membership in the base is at least α. This is the syntactic representation of graded beliefs [START_REF] Benferhat | A practical approach to revising prioritized knowledge bases[END_REF] that we will use in this paper. Alternatively, one may regard a belief as a necessity degree induced by a normalized2 possibility distribution π on the possible worlds ω ∈ Ω [2]: π : Ω → [0, 1], where π(ω) is the possibility degree of interpretation ω. In this case, the degree to which a given formula φ ∈ L is believed can be calculated as

B(φ) = N ([φ]) = 1 -max ω |=φ π(ω),
where N is the necessity measure induced by π. This is the semantic representation of graded beliefs proposed in [START_REF] Benferhat | A practical approach to revising prioritized knowledge bases[END_REF] which has been used in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF].

The syntactic and the semantic representations of graded beliefs are equivalent [START_REF] Dubois | Possibilistic logic[END_REF]. Therefore, they may be used interchangeably as convenience demands. This means that, given a belief base B such that, for all α, Bα is consistent, one can construct a possibility distribution π such that, for all φ ∈ L, N ([φ]) = max{α : Bα |= φ}, where Bα is the α-cut of base B. In particular, π may be defined as follows: for all ω ∈ Ω,

π(ω) = 1 -max{α : Bα |= ¬φω}, (5)
where φω denotes the minterm of ω, i.e., the formula satisfied by ω only. Notice that π is normalized. Indeed, since, by hypothesis, for all α, Bα is consistent, there exists an interpretation ω * ∈ Ω, such that, for all α ∈ (0, 1], ω * |= Bα; therefore, π(ω * ) = 1, because no formula φ exists such that ω * |= φ and B(φ) > 0.

A Possibilistic BDI Model

The possibilistic BDI model of agency we adopt is an adaptation of the one used in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF]. The main difference is that we replace the semantic representation of beliefs and desires with syntactic representations in the form of a belief base B and a desire base D. The agent receives information from the outside world and thus forms its beliefs. The "program" of the agent consists of a number of desire-generation rules, contained in a rule-base RJ . These, together with the beliefs, determine a set of goals, which are then fed into a planner module to compute the actions to be performed by the agent.

The internal mental state of the agent is completely described by a possibilistic belief base B and by a set of desire-generation rules RJ . The set B of the formulas believed by the agent is computed from B according to Equation 4.

The set J of the agent's justified desires is generated dynamically through a deliberation process which applies the rules in RJ to the current beliefs and justified desires to produce a desire base D, from which the fuzzy set of justified desires J is computed according to Equation 11.

Finally, the agent rationally elects its goals G * from the justified desires J as the most desirable of the possible sets of justified desires, according to its beliefs. The agent then plans its actions to achieve the elected goals G * by means of a planner module, whose discussion lies outside of the scope of this paper.

In the following sections, we give technical details about all the components of the model and about syntactic deliberation and goal election algorithms. To help the reader, we will use an adaptation of the running example of [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF], which follows.

Example Dr. A. Gent has submitted a paper to ECAI 2014 he has written with his co-author I. M. Flaky, who has promised to go to Prague to present it if it is accepted. Dr. Gent knows that, if the paper is accepted, publishing it (which is his great desire), means to pay the conference registration (for his co-author or for himself) and then be ready to go to Prague to present it, in case I. M. is unavailable.

If the paper is accepted (a), Dr. Gent is willing to pay the registration (r); furthermore, if the paper is accepted and Dr. Flaky turns out to be unavailable (q), he is willing to go to Prague to present it (p). Finally, if he knows the paper is accepted and wishes to present it, he will desire to have a hotel room (h) and a plane ticket reserved (t).

Then, one fine day, Dr. Gent receives the notification of acceptance of his paper: the source is fully trustworthy-the program chair of ECAI 2014. However, soon after learning that the paper has been accepted, Dr. Flaky rushes into Dr. Gent's office to inform him that he is no more available to go to Prague; as always, Dr. Gent does not completely trust what Dr. Flaky tells him, as he is well accustomed to his changing mind. A few weeks later, rumors get to Dr. Gent's ear that all the hotels in Prague are full (f ); although Dr. Gent considers this news as yet unverified, he takes notice of it. Let's assume that, combined with his a priori beliefs, like that if the planes are all booked out (b), he might not succeed in reserving a flight, this yields Dr. Gent's beliefs, represented by the following base: B = {(b ⊃ ¬t, 0.9), (f, 0.2), (q, 0.75), (a, 1), (f ⊃ ¬h, 1), (p ⊃ (a ∧ r), 1)}.

Representing Desires

We may regard desires as expression of preference for some states of affairs over some others. Therefore, from a semantic point of view, such preference may be encoded as an assignment of a qualitative utility u(ω) ∈ [0, 1] to every world ω ∈ Ω: u(ω) = 0 means that ω is among the least preferred worlds, u(ω) = 1 means that ω is among the most preferred worlds, and u(ω) > u(ω ′ ) means that ω is preferred to ω ′ . Such qualitative utility u is thus, formally, a possibility distribution. Guaranteed possibility measure, ∆, can be used to associate a preference degree to arbitrary formulas [START_REF] Benferhat | Logical representation and fusion of prioritized information based on guaranteed possibility measures: application to the distance-based merging of classical bases[END_REF]. The set of the agent's justified desires, J, a fuzzy set of formulas in the language of choice, is thus defined based on possibility distribution u, which, unlike π, needs not be normalized, since desires may very well be inconsistent, as follows: for all formulas φ ∈ L,

J(φ) = ∆([φ]) = min ω|=φ u(ω). (7)
J may be extended to sets of formulas in the following way. Let S ⊆ L be a set of formulas, then J(S) = min φ∈S J(φ).

The basic mechanism which determines how desires arise, i.e., which desires are justified and to which extent, is rule-based and may be described in terms of desire-generation rules.

Definition 3 (Desire-Generation Rule) A desire-generation rule r is an expression of the form βr, ψr ⇒ + D φ, 3 where βr, ψr, φ ∈ L.

The unconditional counterpart of this rule is

α ⇒ + D φ, with α ∈ (0, 1].
The intended meaning of a conditional desire-generation rule is: "an agent desires every world in which φ is true at least as much as it believes βr and desires ψr". The intended meaning of an unconditional rule is straightforward: the degree to which the agent desires φ is α.

Given a desire-generation rule r, we shall denote rhs(r) the formula on the right-hand side of r.

Example (continued) Dr. Gent's RJ may be described by the following desire-generation rules:

R1 : a, p ⇒ + D t ∧ h, R2 : a ∧ q, ⊥ ⇒ + D p, R3 : a, ⊥ ⇒ + D r.
The degree of activation of a desire-generation rule depends on the degree to which its left-hand side is satisfied, i.e., the degree to which βr is believed and ψr is desired.

Definition 4 (Rule Activation) Let r = βr, ψr ⇒ + D φ be a desiregeneration rule. The degree af activation of r, Deg(r), is given by Deg(r) = min{B(βr), J(ψr)}.

For an unconditional rule r = αr ⇒ + D φ, Deg(r) = αr. A semantic deliberation algorithm was proposed in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF] which, given a belief base B induced by a possibility distribution π and a set of desire-generation rules RJ , computes the corresponding possibility distribution u.

In order to replace such algorithm with a syntactic deliberation algorithm, the first step is to replace the qualitative utility u with a desire base D, which, by analogy with the belief base B, will be represented as a fuzzy set of formulas. However, whereas the membership degrees of B are interpreted and treated as necessity degrees, the membership degrees of formulas of D are to be interpreted as minimum guaranteed possibilities. In D, each piece of information [φi, αi] expresses that any world satisfying φi is considered satisfactory for the agent to at least a degree αi.

Let supp(D) = {φ : D(φ) > 0}. One important property that a desire base should obey, which is a consequence of the membership degrees representing guaranteed possibilities is that, for all formulas φ, ψ ∈ supp(D),

if φ |= ψ then D(φ) ≥ D(ψ). (8)
Exactly as a belief base B induces a corresponding possibility distribution πB, a desire base D induces a corresponding qualitative utility uD as follows (cf. Definition 12 in [START_REF] Benferhat | Logical representation and fusion of prioritized information based on guaranteed possibility measures: Application to the distance-based merging of classical bases[END_REF]): for all ω ∈ Ω,

uD(ω) = max φ:ω|=φ D(φ). (9)
If ω does not satisfy any formula φ in D then ω is not satisfactory at all for the agent. Formally, if there is no [φi, ai] ∈ D such that ω |= φi then uD(ω) = 0.

The above definition may be understood as follows: every formula φ occurring in the desire base may be regarded as the representative of [φ], the set of its models. It has been proven in [START_REF] Dubois | Knowledge-driven versus datadriven logics[END_REF] that uD, defined as per Equation 9, is the most specific possibility distribution satisfying D(φ) = ∆([φ]) for all formulas φ ∈ supp(D). Now, let JD be the (fuzzy) set of justified desires in D. By definition, the degree of justification of all formulas occurring in the base must be identical to their degree of membership in the base, i.e., for all formulas φ, if D(φ) > 0 then JD(φ) = D(φ).

If D(φ) = α, it means that min ω∈[φ] uD(ω) = α, or, in other terms, that, ∀ω ∈ [φ], uD(ω) ≥ α. Therefore, we get Equation 9.

The next step is to show how, given a desire base D, the degree of justification of any arbitrary desire formula ψ may be calculated. Based on the definition of J, we may write, for all formulas ψ,

JD(ψ) = min ω|=ψ uD(ω) = min ω|=ψ max φ:ω|=φ D(φ). (10)
We now have to eliminate all references to the models ω from the above formula in order to make it "syntactic". We do not want references to the models because we do not want to be obliged to enumerate all the interpretations explicitly in order to compute JD(ψ). This is indeed the reason why the syntactic view for goal generation proposed in this paper should lead to a more efficient implementation than the semantic view already present in the literature. We will show this assertion later in the paper. One way to obtain the elimination of all references to the models is to construct a set of formulas P(D), the "partition" of Ω according to D, containing all the 2 supp(D) conjunctions of positive or negated formulas occurring in D. Then, with the convention that max ∅ = 0,

JD(ψ) = min ξ∈P(D) ψ∧ξ =⊥ max φ∈supp(D) ξ|=φ D(φ) = = min i=0,...,2 n -1 ψ∧ξ i =⊥ max j=1,...,n ξ i |=φ j D(φj). ( 11
)
Proof: The models of the formulas in P(D) form a partition of Ω:

2 n -1 i=0 [ξi] = Ω, ∀i, j, i = j, ξi ∧ ξj = ⊥.
The qualitative utility uD is constant over each

Therefore, instead of minimizing over all ω ∈ Ω, it is sufficient to minimize over all ξi ∈ P(D). Moreover, since we are minimizing uD(ω) over all models of ψ, we should only consider those ξi such that [ξi] ∩ [ψ] = ∅, i.e., such that ξi ∧ ψ = ⊥. ✷ Example (continued) Let's assume Dr. Gent's desire base is D = {(t ∧ h, 0.75), (p, 0.75), (r, 1)}.

We will see later how such base may be derived from the desiregeneration rules. To be able to compute J(ψ) for any arbitrary formula ψ, we may pre-compute P(D) and, for all ξi ∈ P(D), the corresponding term αi = max φ∈supp(D) ξ i |=φ D(φ) in Equation 11:

ξ0 = ¬(t ∧ h) ∧ ¬p ∧ ¬r, α0 = max ∅ = 0, ξ1 = (t ∧ h) ∧ ¬p ∧ ¬r, α1 = max{0.75} = 0.75, ξ2 = ¬(t ∧ h) ∧ p ∧ ¬r, α2 = max{0.75} = 0.75, ξ3 = (t ∧ h) ∧ p ∧ ¬r, α3 = max{0.75, 0.75} = 0.75, ξ4 = ¬(t ∧ h) ∧ ¬p ∧ r, α4 = max{1} = 1, ξ5 = (t ∧ h) ∧ ¬p ∧ r, α5 = max{0.75, 1} = 1, ξ6 = ¬(t ∧ h) ∧ p ∧ r, α6 = max{0.75, 1} = 1, ξ7 = (t ∧ h) ∧ p ∧ r, α7 = max{0.75, 0.75, 1} = 1.
Now if, for instance, we want to compute J(ψ), we will have to compute the minimum, for i such that ψ ∧ ξi = ⊥, of the αi. Therefore, we have, for example, J(t) = 0 and J(r) = 1.

Based on Equation 11, computing the degree of justification of formula ψ given the desire base D requires O(n2 n ) entailment checks, where n is the size of the desire base. Checking whether a formula entails another formula is a logical reasoning problem which may be reduced to the satisfiablility problem, whose computational complexity varies depending on the specific logic considered, but does not depend on n. For instance, satisfiability in propositional logic (also known as Boolean satisfiability) is NP-complete [START_REF] Garey | Computers and Intractability: A guide to the Theory of NP-Completeness[END_REF]; concept satisfiability in description logics goes from polynomial to NEXPTIMEcomplete [START_REF]The Description Logic Handbook: Theory, implementation and applications[END_REF].

Generating Desires

We are now ready to present a syntactic deliberation algorithm, which calculates the desire base given a set of desire-generation rules RJ and a belief set B.

We first recall the semantic deliberation algorithm presented in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF]. Let R ω J = {r ∈ RJ : ω |= rhs(r)} denote the subset of RJ containing just the rules whose right-hand side would be true in world ω and Deg µ (r) the degree of activation of rule r calculated using µ as the qualitative utility assignment. Given a mental state S = π, RJ , the following algorithm computes the corresponding qualitative utility assignment, u.

Algorithm 1 (Semantic Deliberation) INPUT: π, RJ . OUTPUT: u. 1. i ← 0; for all ω ∈ Ω, u0(ω) ← 0; 2. i ← i + 1; 3. For all ω ∈ Ω, ui(ω) ← max r∈R ω J Deg u i-1 (r), if R ω J = ∅, 0, otherwise;
4. if maxω |ui(ω)ui-1(ω)| > 0, i.e., if a fixpoint has not been reached yet, go back to Step 2; 5. For all ω ∈ Ω, u(ω) ← ui(ω); u is the qualitative utility assignment corrisponding to mental state S.

Let Deg X (r) be the degree of activation of rule r based on the desire base X. When X = ∅, the degree of justification of all desire formulas is zero. The belief base B does not change during the deliberation process.

Algorithm 2 (Syntactic Deliberation)

INPUT: B, RJ . OUTPUT: D. 1. i ← 0; D0 ← ∅; 2. i ← i + 1; 3. Di ← {(rhs(r), Deg D i-1 (r)) : r ∈ RJ }; 4. if Di = Di-1, i.e.

, if a fixpoint has not been reached yet, go back

to Step 2; 5. D ← Di and the deliberation has finished.

We now have to prove that this syntactic deliberation algorithm is equivalent to the semantic deliberation algorithm presented in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF]. The hypothesis is summarized in the following diagram:

RJ ✲ Algorithm 1 u ✲ Algorithm 2 D ✲ Equation 9 uD
Proposition 2 We can prove that uD = u.

Proof: We may proceed by induction: we will consider the sequence {ui}i=0,1,... of the possibility distributions constructed by Algorithm 1 and the sequence {uD i }i=0,1,..., whose elements uD i are the possibility distributions induced by the desire bases Di constructed by Algorithm 2, and we will prove that, if ui-1 = uD i-1 , then ui = uD i . Now, ui-1 = uD i-1 means that Deg u i-1 (r) = Deg D i-1 (r) for all rules r. By Equation 9, for all ω ∈ Ω, we may write

uD i (ω) = max φ:ω|=φ Di(φ) = max r∈R J :ω|=rhs(r) Deg D i-1 (r) = max r∈R ω J Deg u i-1 (r) = ui(ω).
This proves the induction step. Finally, it is straightforward to verify that u0 = uD 0 , therefore the two sequences will be identical and so their limits, u and uD, and this concludes the proof. ✷

Example (continued) Let us apply Algorithm 2 to Dr. Gent's mental state: we obtain D0 = ∅, D1 = {(p, 0.75), (r, 1)}, D2 = {(t ∧ h, 0.75), (p, 0.75), (r, 1)}, D3 = {(t ∧ h, 0.75), (p, 0.75), (r, 1)} = D2.

Therefore,

D = {(t ∧ h, 0 .75), (p, 0.75), (r, 1)}. (13)
7 Generating Goals

In [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF], the assumption was made that a rational agent would select as goals the maximally justified desires among the maximally possible desires. In other words, a rational agent should first restrict attention only to those desires that it would be most normal (i.e., unsurprising, likely, . . . ) to expect they might come true and then decide to actively pursue those, among them, that have the highest qualitative utility.

In order to write a goal election algorithm according to such assumption, we need to define the set of desires possible to a given degree.

Definition 5 Given γ ∈ (0, 1], Jγ = {φ ∈ supp(J) : Π([φ]) ≥ γ} is the (classical) subset of supp(J) containing only those desires whose overall possibility is at least γ. We recall that Π([φ]) = 1 -B(¬φ).

We now define a goal set for a given level of possibility γ, as the set of the maximally justified γ-possible desires. Definition 6 (Goal set) The γ-possible goal set is

Gγ = arg max S⊆Jγ J(S) if Jγ = ∅, ∅ otherwise.
We denote by γ * the maximum possibility level such that Gγ = ∅. Then, the goal set elected by a rational agent will be

G * = Gγ * , γ * = max Gγ =∅ γ. ( 14
)
Let Img(π) be the level set 4 of possibility distribution π and Img(u) be the level set of qualitative distribution u. Notice that Img(u) and Img(π) are both finite, independently of Ω being finite, as proven in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF].

The following two algorithms, adapted from [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF], allow an agent to compute Gγ for a given possibility lower bound γ, and the optimal goal set G * , based on a semantic representation of beliefs and desires as two possibility distributions, π and u. We will call them semantic, to distinguish them from the two algorithms that we are going to propose to replace them, which will assume a syntactic representation of both beliefs and desires. This also suggests that we define γ = 1γ as the "impossibility of the goals", and reformulate the goal election algorithm as a search for the least impossible set of maximally justified goals. ✷ Since (i) Algorithm 6 iterates over the level set of the belief base and (ii) Algorithm 5, which is called as a subroutine at each iteration of Algorithm 6, loops over the level set of the desire base, and (iii) the most complex task performed at each iteration of Algorithm 5 is computing the degree of belief of each negated desire formula ψ in a δ-cut of the desire base, we may conclude that the computational cost (in number of entailment checks) of syntactic goal election is

O( Img(B) • Img(D) ) • C,
where C is the number of entailment checks needed for computing B(φ) for an arbitrary formula φ, given the belief base B, which is done using Equation 4, thus giving C = O( Img(B) ). Furthermore, we may observe that Img(B) ≤ B = m and Img(D) ≤ D = n. Therefore, we may conclude that carrying out the syntactic goal election requires O(m 2 n) entailment checks.

The termination of Algorithms 1 and 4 is proved in [START_REF] Da | Belief-goal relationships in possibilistic goal generation[END_REF]; the termination of Algorithms 2 and 6 is a direct consequence of their termination.

Example (continued) We may now apply Algorithm 6 to elect the goals of Dr. Gent, given that his belief base is the one given in Equation 6 and his desire base is the one given in Equation 13; therefore, Img(B) = {0.2, 0.75, 0.9, 1} and Img(D) = {0.75, 1}.

We begin by calling Algorithm 5 with γ = 0: δ is set to max Img(D) = 1, and the corresponding δ-cut of D is in fact the core of D, D1 = {(r, 1)}. Now, B(¬r) = 0 ≤ γ; therefore G0 = {r} and Algorithm 6 terminates immediately with γ * = 1γ = 1, G * = G0 = {r}, i.e., Dr. Gent will elect as his goal just to register to ECAI 2014.

Conclusions

We have proposed a syntactic representation for desires within a possibilistic BDI model of agency; we have shown its equivalence to the semantic representation based on a qualitative utility; we have provided the syntactic equivalent of the deliberation algorithm, which generates the set of justified desires given a set of desire-generation rules and a belief base. We have then provided the syntactic equivalent of the goal election algorithm, which generates the goals as the maximally justified desires among the maximally possible desires.

The cost of computing the degree of justification of a formula and of electing the goals has been given in terms of the basic operation of checking whether a formula entails another formula. Even though the cost of the former task grows exponentially with to the size of the desire base, in practice it is expected to be feasible, given that the size of the desire base depends on the number of desire-generation rules and it is hard to think of applications that would call for a large number of such rules. On the other hand, goal generation is polynomial in the size of the belief and desire bases. Compare this to what happens when using the semantic representation, where 2 L interpretations have to be explicitly represented and iterated over to compute the degree of justification of a formula. In our opinion, a syntactic computation of beliefs, desires, and goals is the only viable alternative to implement possibilitic agents based on languages whose semantics involve an infinite number of possible worlds. This is why we believe the results here presented are a first and important step towards the practical implementation of a possibilistic BDI framework.

Figure 1 .

 1 Figure 1. A schematic illustration of the proposed BDI model. The meaning of the symbols is explained in the text.

Proposition 1

 1 Let supp(D) = {φ1, φ2, . . . , φn} and let P(D) = {ξ0, ξ1, ξ2, . . . , ξ2n-1} with ξ0 = ¬φ1 ∧ ¬φ2 ∧ . . . ∧ ¬φn, ξ1 = φ1 ∧ ¬φ2 ∧ . . . ∧ ¬φn, . . . . . . ξ2n-1 = φ1 ∧ φ2 ∧ . . . ∧ φn.

  [ξi]: for all ω |= ξi, uD(ω) = max j=1,...,n ξ i |=φ j D(φj).

Proof:

  We begin by observing that Img(u) = Img(D), i.e., the level set of the desire base is the same as the level set of the corresponding qualitative utility. Furthermore, since we now have a desire base, we may replace the construction of φ based on the minterms in Step 2 of Algorithm 3 with a more straightforward the disjunction of all the formulas in the δ-cut of the desire base. This also suggests that, instead of returning Gγ = {φ}, it is equivalent, but more intuitive, to return Gγ = D δ . Finally, instead of testing the condition Π([φ]) ≥ γ in Step 3, it is equivalent to test the condition B(¬φ) ≤ 1γ. Now, by the DeMorgan laws and by the properties of necessity, 5 B(¬φ) = B ψ∈D δ ¬ψ = min ψ∈D δ B(¬ψ), which allows us to avoid constructing φ explicitly and to test directly the condition min ψ∈D δ B(¬ψ) ≤ 1γ.

5 N

 5 (A ∩ B) = min{N (A), N (B)}.

  ≥ γ, terminate with Gγ = {φ}; otherwise, 4. δ ← max{α ∈ Img(u) : α < δ}, 0 if no such α exists; 5. if δ > 0, go back to Step 2; 6. terminate with Gγ = ∅.Algorithm 4 (Semantic Goal Election) INPUT: π, u. OUTPUT: G * . 1. γ ← max Img(π) = 1, since π is normalized; 2. compute Gγ by Algorithm 3; 3. if Gγ = ∅, terminate with γ * = γ, G * = Gγ; otherwise, 4. γ ← max{α ∈ Img(π) : α < γ}, 0 if nosuch α exists; 5. if γ > 0, go back to Step 2; 6. terminate with G * = ∅: no goal may be elected. Gγ = ∅, terminate with γ * = 1γ, G * = Gγ; otherwise, 4. γ ← min{α ∈ Img(B) : α > γ}, 1 if no such α exists; 5. if γ < 1, go back to Step 2; 6. terminate with G * = ∅: no goal may be elected.

	Algorithm 3 (Semantic Computation of Gγ)
	INPUT: π, u. OUTPUT: Gγ.
	1. δ ← max Img(u);
	2. determine the least specific formula φ such that J(φ) ≥ δ as fol-
	lows:
	φ ←	φω,
	u(ω)≥δ
	where φω denotes the minterm of ω, i.e., the formula satisfied by
	ω only;
	3. if Π([φ]) Proposition 3 The syntactic versions of Algorithms 3 and 4 are Al-
	gorithms 5 and 6 given below.
	Algorithm 5 (Syntactic Computation of Gγ)
	INPUT: B, D. OUTPUT: Gγ.
	1. δ ← max Img(D);
	2. if min ψ∈D δ B(¬ψ) ≤ γ, terminate with Gγ = D δ ;
	otherwise,
	3. δ ← max{α ∈ Img(D) : α < δ}, 0 if no such α exists;
	4. if δ > 0, go back to Step 2;
	5. terminate with Gγ = ∅.
	Algorithm 6 (Syntactic Goal Election)
	INPUT: B, D. OUTPUT: G * .
	1. γ ← 0;
	2. compute Gγ by Algorithm 5;
	3. if

Univ. Nice Sophia Antipolis, CNRS, I3S, UMR 7271, 06900 Sophia Antipolis, France, email: {celia.pereira,andrea.tettamanzi}@unice.fr

Normalization of a possibility distribution corresponds to consistency of the beliefs.

Note that the implication used to define a desire-generation rule is not the material implication.

The level set of a possibility distribution π is the set of α ∈ [0, 1] : ∃ω such that π(ω) = α.