Sophie Donnet

Adeline Samson

EM algorithm coupled with particle filter for maximum likelihood parameter estimation of stochastic differential mixed-effects models

Keywords: Mixed models, Stochastic Differential Equations, SAEM algorithm, Particle Filter, PMCMC, Stochastic volatility, Time-inhomogeneous SDE

Biological processes measured repeatedly among a series of individuals are standardly analyzed by mixed models. These biological processes can be adequately modeled by parametric Stochastic Differential Equations (SDEs). We focus on the parametric maximum likelihood estimation of this mixed-effects model defined by SDE. As the likelihood is not explicit, we propose a stochastic version of the Expectation-Maximization algorithm combined with the Particle Markov Chain Monte Carlo method. When the transition density of the SDE is explicit, we prove the convergence of the SAEM-PMCMC algorithm towards the maximum likelihood estimator. Two simulated examples are considered: an Ornstein-Uhlenbeck process with two random parameters and a time-inhomogeneous SDE (Gompertz SDE) with a stochastic volatility error model and three random parameters. When the transition density is unknown, we prove the convergence of a different version of the algorithm based on the Euler approximation of the SDE towards the maximum likelihood estimator.

Introduction

Biological processes are usually measured repeatedly among a collection of individuals or experimental animals. The parametric statistical approach commonly used to discriminate between the inter-subjects variability (the variance of the individual regression parameters) and the residual variability (measurement noise) is the mixed-effects model methodology [START_REF] Pinheiro | Mixed-effect models in S and Splus[END_REF]. The choice of the regression function depends on the context. Dynamical processes are frequently described by deterministic models defined as ordinary differential equations. However, these functions may not capture the exact process, since responses for some individuals may display local "random" fluctuations. These phenomena are not due to error measurements but are induced by an underlying biological process that is still unknown or unexplained today. The perturbation of deterministic models by a random component leads to the class of stochastic differential equations (SDEs). In this paper we develop a parametric estimation method for mixed model defined by an SDE process, called SDE mixed model. Note that SDE mixed models can be viewed as an extension of state space models with random parameters. Except in the linear case, the likelihood of standard mixed models is intractable in a closed form, making the estimation a hard task. [START_REF] Davidian | Nonlinear models to repeated measurement data[END_REF], Wolfinger (1993), [START_REF] Pinheiro | Mixed-effect models in S and Splus[END_REF] and [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] develop various methods for determinitic regression function. When the regression term is the realization of a diffusion process, the estimation is made more complicated by the difficulties deriving from the SDEs. Parametric estimation of SDE with random parameters (without measurement noise) has been studied by [START_REF] Ditlevsen | Stochastic vs. deterministic uptake of dodecanedioic acid by isolated rat livers[END_REF], [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] or [START_REF] Oravecz | A hierarchical Ornstein-Uhlenbeck model for continuous repeated measurement data[END_REF] in a Bayesian context. SDE mixed model estimation, including measurement noise modeling, is more complex and has received little attention. [START_REF] Overgaard | Non-linear mixed-effects models with stochastic differential equations: Implementation of an estimation algorithm[END_REF] and Tornøe et al. (2005) build estimators based on an extended Kalman filter but without proving the convergence of the algorithm. [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] propose to use a stochastic version of the Expectation-Maximisation algorithm. However, the method involves Markov Chain Monte Carlo (MCMC) samplers which have proved their slow mixing properties in the context of state space models with random parameters. For state space models with fixed known parameters, the Sequential Monte Carlo (SMC) algorithms have demonstrated their efficiency. However, these techniques are difficult to extent to the case of unknown parameters [START_REF] Casarin | Online data processing: comparison of Bayesian regularized particle filters[END_REF]. Recently, [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] developed a powerful algorithm combining MCMC and SMC, namely the Particle Markov Chain Monte Carlo (PMCMC), taking the advantage of the strength of its two components. Besides, its convergence toward the exact distribution of interest is established. To our knowledge, PMCMC methods have been used for parameter estimation only in a Bayesian context and for models without random individual parameters. We propose to focus on maximum likelihood estimation for SDE mixed models by combining a PMCMC algorithm with the SAEM algorithm. When the transition density of the SDE is explicit, we prove the convergence of the SAEM-PMCMC algorithm towards the maximum likelihood estimator. Two simulated examples are considered, an Ornstein-Uhlenbeck mixed SDE model observed with additive noise and a stochastic volatility mixed model with a time-inhomogeneous SDE, showing the increase in accuracy for the SAEM-PMCMC estimates compared with the SAEM-MCMC algorithm [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF]. When the transition density is unknown, we prove the convergence of a different version of the algorithm towards the maximum likelihood estimator. The paper is organized as follows. In Section 2, we present the SDE mixed model. In Section 3, an estimation method is proposed when the transition density of the SDE is explicit. Section 4 shows numerical results based on simulated data in this case. In Section 5, we extend the approach to general SDEs and prove the convergence of the proposed algorithm. In Section 6 we conclude and discuss the advantages and limitations of our approach.

Mixed model defined by stochastic differential equation

Model and notations

Let y = (y ij ) 1≤i≤n,0≤j≤Ji denote the data, where y ij is the noisy measurement of the observed process for individual i at time t ij ≥ t i0 , for i = 1, . . . , n, j = 0, . . . , J i . Let y i,0:Ji = (y i0 , . . . , y iJi ) denote the data vector of subject i. We consider that the y ij 's are governed by a mixed-effects model based on a stochastic differential equation defined as follows:

y ij = g(X ij , ε ij ), ε ij ∼ i.i.d. N (0, σ 2 ), (1)
dX it = a(X it , t, φ i )dt + b(X it , t, φ i , γ)dB it , X iti0 |φ i ∼ π 0 (•|φ i ) (2) φ i ∼ i.i.d. p(φ; β). (3)
In equation ( 1), (ε ij ) i,j are i.i.d Gaussian random variables of variance σ 2 representing the measurement errors and g is the known error model function. The function g(x, ε) = x + ε leads to an additive error model. Multiplicative error models or stochastic volatility models can be considered with g(x, ε) = x(1 + ε) or g(x, ε) = f (x) exp(ε) with a known function f . For each individual i, (X ij ) 0≤j≤Ji are the values at discrete times (t ij ) 0≤j≤Ji of the continuous process (X it ) t≥ti0 = (X it (φ i )) t≥ti0 defined by equation ( 2). This process is proper to each individual through the individual parameters φ i ∈ R p involved in the drift and volatility functions a and b and an individual Brownian motion (B it ) t≥ti0 . Functions a and b -which may be functions of time, leading to time-inhomogeneous SDEs -are common to the n subjects. They are assumed to be sufficiently regular (with linear growth) to ensure a unique strong solution [START_REF] Oksendal | Stochastic differential equations: an introduction with applications[END_REF]. Besides, b depends on an unknown additional volatility parameter γ. The (B it ) and φ k are assumed to be mutually independent for all 1 ≤ i, k ≤ n. The initial condition X i0 of process (X it ) t≥ti0 , conditionally to the individual parameter φ i is random with distribution π 0 (•|φ i ). In the following, X i,0:Ji = (X i0 , . . . , X iJi ) denotes the vector of the i-th process realization at observation times (t ij ) for i = 1, . . . , n and X = (X 1 , . . . , X n ) denotes the whole vector of processes at observation times, for all individuals. The individual parameters φ i are assumed to be distributed with a density p(φ; β) depending on parameter β. We denote Φ = (φ 1 , . . . , φ n ) the vector of all individual parameter vectors.

For simplicity purpose, we restrict to a same number of observations per subject: J i = J for all i. The extension to the general case is straightforward.

Likelihood function

Let θ = (β, γ, σ) be the parameter vector of interest which belongs to some open subset Θ of the Euclidean space R q with q the number of unknown parameters. Our objective is to propose a maximum likelihood estimation of θ. The likelihood function is well-defined under the assumption of existence of a strong solution to (2). We denote x → p(x, ts, s|x s , φ i ; θ) the density of X it given φ i and X is = x s , s < t with respect to the Lebesgue measure. This allows us to write the likelihood of model (1-3) as

p(y; θ) = n i=1 p(y i,0:J |X i,0:J ; σ)p(X i,0:J |φ i ; γ)p(φ i ; β)dφ i dX i,0:J , (4)
where p(y|x; σ), p(x|φ; γ) and p(φ; β) are density functions of the observations given the diffusion process, the diffusion process given the individual parameters and the individual parameters, respectively. By independence of the measurement errors (ε ij ), we have

p(y i,0:J |X i,0:J ; σ) = J j=0 p(y ij |X ij ; σ) where p(y ij |X ij ; σ)
is the density function associated to the error model (1). The Markovian property of the diffusion process (X it ) implies p(X i,0:

J |φ i ; θ) = π 0 (X i0 |φ i ) J j=1 p(X ij , ∆ ij , t ij-1 |X i j-1 , φ i ; γ),
where ∆ ij = t ijt ij-1 . Except for a mixed Wiener process with drift and additive noise, the integral (4) has no-closed form.

Maximizing (4) with respect to θ yields the corresponding Maximum Likelihood Estimate (MLE) θ of the SDE mixed model (1-3). As the likelihood (4) is not explicit, it has to be numerically evaluated or maximized. This requires eitheir to know the transition densities explicitely (see Section 3) or to approximate it by an Euler scheme (see Section 5).

SDE with explicit transition density: estimation method

In this section, we assume that:

(M0) the SDE defined by ( 2) has an explicit transition density p(x, ts, s|x s , φ i ; θ).

Even with this assumption, the estimation is complex because the n random parameters φ i and random trajectories (X i,0:J ) are unobserved. This statistical problem can be viewed as an incomplete data model. The observable vector y is thus considered as part of a so-called complete vector (y, X, Φ).

Estimation algorithm (SAEM algorithm)

The Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is useful in situations where the direct maximization of the marginal likelihood θ → p(y; θ) is more complex than the maximization of the conditional expectation of the complete likelihood

Q(θ|θ ) = E [log p(y, X, Φ; θ)|y; θ ] ,
where p(y, X, Φ; θ) is the likelihood of the complete data (y, X, Φ). The EM algorithm is an iterative procedure: at the m-th iteration, given the current value θ m-1 of the parameters, the E-step is the evaluation of Q m (θ) = Q(θ | θ m-1 ) while the M-step updates θ m-1 by maximizing Q m (θ). For cases where the E-step has no closed form, [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] propose the Stochastic Approximation EM algorithm (SAEM) and replace the E-step by a stochastic approximation of Q m (θ). The E-step is thus divided into a simulation step (S-step) of the non-observed data (X (m) , Φ (m) ) with the conditional distribution p(X, Φ |y; θ m-1 ) and a stochastic approximation step (SA-step):

Q m (θ) = Q m-1 (θ) + α m log p y, X (m) , Φ (m) ; θ -Q m-1 (θ) ,
where (α m ) m∈N is a sequence of positive numbers decreasing to zero. To fulfill the convergence conditions of the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF], we consider the exponential case. More precisely, we assume:

(M1) The parameter space Θ is an open subset included in a compact set of R q . The complete likelihood p(y, X, Φ) belongs to the exponential family i.e. log p(y, X, Φ; θ) = -ψ(θ) + S(y, X, Φ), ν(θ) ,

where ψ and ν are two functions of θ, S(y, X, Φ) is known as the minimal sufficient statistics of the complete model, taking its value in a subset S of R d and •, • is the scalar product on R d .

In this case, the SA-step of the SAEM algorithm reduces to the approximation of E[S(y, X, Φ)|y; θ ]. At step S, a simulation under the conditional distribution p(X, Φ|y; θ m-1 ) is required. However, this distribution is likely to be a complex distribution, resulting in the impossibility of a direct simulation of the non-observed data (X, Φ). [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] suggest to realize the simulation step via a Markov Chain Monte Carlo (MCMC) scheme, resulting in the SAEM-MCMC algorithm. The S-step consists in constructing a Markov chain with p(X, Φ|y; θ m-1 ) as ergodic distribution at the m-th iteration. This S-step by a Markov kernel is detailed in the following subsection.

Simulation of the latent variables (X, Φ) given y

We know deal with the construction of a Markov Chain with ergodic distribution p(X, Φ|y;

θ m-1 ). First note that conditionally to the observations y, the individuals are independent. As a consequence, the simulation of each (X i,0:J , φ i ) conditionally to y 0:J can be performed separately. To ease the reading, in this subsection, we focus on the simulation of the non observed variables (X i,0:J , φ i ) of a single individual. Therefore, we omit the index i in the following. Moreover, we denote θ = (β, γ, σ) the current values of the parameters and omit the index m of the SAEM algorithm iteration. Standard version of the MCMC alternately simulates X 0:J under the distribution p(X 0:J |y 0:J , φ; γ, σ) and φ under the distribution p(φ|y 0:J , X 0:J ; γ, σ). If the conditional distributions p(X 0:J |y 0:J , φ; γ, σ) and p(φ|y 0:J , X 0:J ; γ, σ) cannot be simulated easily, we resort to Metropolis Hastings algorithms for each component of φ and each time component of X 0:J [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF]. However, standard MCMC algorithms have reached their limits in high dimensional context: they do not exploit the Markovian structure of the data and have proved slow mixing properties. We propose to replace them by the Particle Markov Chain Monte Carlo (PMCMC) algorithm proposed by [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF]. The idea of the PMCMC algorithm has been first proposed by [START_REF] Beaumont | Estimation of population growth or decline in genetically monitored populations[END_REF], formalized by [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF] and developped in the context of state-space models by [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF].

Let us consider an ideal Metropolis-Hastings algorithm updating conjointly φ and X 0:J conditionally to y 0:J . A new candidate (X c 0:J , φ c ) would be generated with a proposal distribution:

q(X c 0:J , φ c |X 0:J , φ; θ) = q(φ c |φ)p(X c 0:J |y 0:J , φ c ; γ, σ). and accepted with probability: ρ(X c 0:J , φ c |X 0:J , φ) = min 1, q(φ|φ c ) q(φ c |φ) p(y 0:J |φ c ; γ, σ)p(φ c ; β) p(y 0:J |φ; γ, σ)p(φ; β) , However, because of the complexity of the model, we are not able to simulate exactly under the conditional distribution p(X 0:J |y 0:J , φ; γ, σ) and the marginal likelihood p(y 0:J |φ; γ, σ) has no closed form. These two points can be tackled by an approximation through a particle filter, also called Sequential Monte Carlo (SMC) algorithm. The SMC algorithm produces a set of K particles (X (k) 0:J ) k=1...K and respective weights (W (k) 0:J ) k=1...K approximating the conditional distribution p(X 0:J |y 0:J , φ; γ, σ) by an empirical measure

Ψ K J = K k=1 W (k) 0:J δ X (k) 0:J
where δ 0 is any probability measure on R J+1 .

Algorithm 1 (SMC algorithm).

• Time j = 0: sample X (k) 0

∼ π 0 (•|φ) ∀ k = 1, . . . , K and compute the weights:

w 0 X (k) 0 = p y 0 , X (k) 0 |φ; γ, σ , W 0 X (k) 0 = w 0 X (k) 0 K k=1 w 0 X (k) 0 • Time j = 1, . . . , J: sample K iid variables X (k) 0:j-1 according to the distribution Ψ K j-1 = K k=1 W (k) 0:j-1 δ X (k) 0:j-1
. Then for each k = 1, . . . , K, sample

X (k) j ∼ q •|y j , X (k) 0:j-1 , φ; γ, σ and set X (k) 0:j = (X (k) 0:j-1 , X (k) j ).
Finally compute and normalize the weights:

wj " X (k) 0:j " = p " yj, X (k)
j |y0:j-1, X (k) 0:j-1 φ; γ, σ " q " X (k) j |yj, X (k) 0:j-1 , φ; γ, σ " , Wj(X (k) 0:j ) = wj " X (k) 0:j " P K k=1 wj " X (k) 0:j "
The simulation of one trajectory X 0:J (called a "run of SMC algorithm" ) under the approximation of p(X 0:J |y 0:J , φ; γ, σ) is directly achieved by randomly choosing one particle among the K particles with weights (W (k) 0:J ) k=1...K . Besides, the marginal distribution p(y 0:J |φ; γ, σ) can be estimated through the weights

p K (y 0:J |φ; γ, σ) = 1 K K k=1 w 0 X (k) 0 J j=1 1 K K k=1 w j X (k) 0:j . (5)
As a consequence, [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] propose the PMCMC algorithm:

Algorithm 2 (PMCMC algorithm).

• Initialization : starting from φ(0), generate X 0:J (0) by a run of SMC algorithmwith K particles-targeting p(X 0:J |y 0:J , φ(0); γ, σ) and estimate p(y 0:J |φ(0); γ, σ) by p K (y 0:J |φ(0); γ, σ)

• At iteration = 1, . . . , N (a) Sample a candidate φ c ∼ q(•|φ( -1)) (b) By a run of SMC algorithm with K particles, generate X c 0:J targeting p(•|y 0:J , φ c ; γ, σ) and compute p K (y 0:J |φ c ; γ, σ) estimating p(y 0:J |φ c ; γ, σ) (c) Set (X 0:J ( ), φ( )) = (X c 0:J , φ c ) and p K (y 0:J |φ( ); γ, σ) = p K (y 0:J |φ c ; γ, σ) with probability b ρ K (X c 0:J , φ c |X0:J ( -1), φ( -1)) = min  1, q(φ( -1)|φ c ) q(φ c |φ( -1)) b p K (y0:J |φ c ; γ, σ)p(φ c ; β) b p K (y0:J |φ( -1); γ, σ)p(φ( -1); β) ff If the candidate is not accepted, then set (X 0:J ( ), φ( )) = (X 0:J ( -1), φ( -1)) and p K (y 0:J |φ( ); γ, σ) = p K (y 0:J |φ( -1); γ, σ)
Remark 1. The proposal distributions q(φ c |φ( )) and q(X j |y j , X j-1 , φ c ; γ, σ) used in the SMC algorithm play a crucial role to ensure good mixing properties of the PMCMC algorithm. They are discussed with more details for the two simulated examples (Section 4).

The most remarkable property of PMCMC is that the distribution of interest p(X 0:J , φ|y 0:J ; γ, σ) is left invariant by the transition kernel, whatever the number of particles K, the ergodicity being reached under weak assumptions. More precisely, let Λ denote the auxiliary variables generated by the SMC algorithm, including the set of generated trajectories, the resampling indices and the indice of the particle randomly picked to obtain "a run of SMC algorithm". [START_REF] Andrieu | The pseudo-marginal approach for efficient Monte Carlo computations[END_REF] prove that PMCMC is an exact MCMC algorithm on (X 0:J , φ). Under general assumptions, the stationary distribution π K (Λ, φ) of the PM-CMC algorithm is such that its marginalized distribution over the auxiliary variables Λ is exactly the distribution p(X 0:J , φ|y 0:J ; γ, σ), independently on K: the PMCMC algorithm generates a sequence (X 0:J ( ), φ( )) whose marginal distribution L K (X 0:J ( ), φ( )|y 0:J ; θ) is such that for all θ ∈ Θ and for all K > 0,

||L K (X 0:J ( ), φ( )|y 0:J ; θ) -p(X 0:J , φ|y 0:J ; θ)|| T V ---→ →∞ 0.
where || • || T V is the total variation distance [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF].

SAEM-PMCMC algorithm and convergence

We now combine algorithms PMCMC and SAEM:

Algorithm 3 (SAEM-PMCMC algorithm).

• Iteration 0: initialization of θ 0 and s 0 = E S(y, X, Φ)|y; θ 0 .

• Iteration m = 1, . . . , M: S-Step: * For each individual i,

• Initialize the PMCMC algorithm by setting φ i (0) as the expectation of p(φ, β m-1 ) • Run N iterations of the PMCMC algorithm with K particles at each iteration, targeting p(X i,0:J , φ i |y i,0:J

; θ m-1 ) * Set X (m) = (X (m) 1,0:J , . . . , X (m)
n,0:J ) and

Φ (m) = (φ (m) 1 , . . . , φ (m)
n ) the simulated non observed variables SA-Step: update of s m-1 using the stochastic approximation scheme:

s m = s m-1 + α m S(y, X (m) , Φ (m) ) -s m-1 (6) M-Step: update of θ m-1 by θ m = arg max θ (-ψ(θ) + s m , ν(θ) ) .
As the PMCMC algorithm can be viewed as a standard MCMC algorithm, the convergence of SAEM-PMCMC can be proved using [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] result. We recall the assumptions of [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF].

(M2) The functions ψ(θ) and ν(θ) are twice continuously differentiable on Θ.

(M3) The function s : Θ -→ S defined as s(θ) = S(y, X, Φ)p(X, Φ|y; θ)dXdΦ is continuously differentiable on Θ.

(M4) The function (θ) = log p(y, θ) is continuously differentiable on Θ and

∂ θ p(y, X, Φ; θ)dXdΦ = ∂ θ p(y, X, Φ; θ)dXdΦ. (M5) Define L : S × Θ -→ R as L(s, θ) = -ψ(θ) + s, ν(θ) . There exists a function θ : S -→ Θ such that ∀θ ∈ Θ, ∀s ∈ S, L(s, θ(s)) ≥ L(s, θ).
(SAEM1) The positive decreasing sequence of the stochastic approximation

(α m ) m≥0 is such that m α m = ∞ and m α 2 m < ∞.
(SAEM2) : Θ → R and θ : S → Θ are d times differentiable, where d is the dimension of S(y, X, Φ).

(SAEM3) For all θ ∈ Θ, ||S(y, X, Φ)|| 2 p(X, Φ|y; θ)dXdΦ < ∞ and the function Γ(θ) = Cov θ (S(X, Φ)) is continuous.
(SAEM4) S is a bounded function.

(SAEM5) The transition probability Π θ of the PMCMC algorithm is Lipschitz in θ and generates a uniformly ergodic chain. The Markov chain (X (m) , Φ (m) ) m≥0 takes its values in a compact subset.

We comment the different assumptions. Assumption (M0) is a strong hypothesis, which is partially relaxed in Section 5. Assumptions (M1-M5), (SAEM1-SAEM3) are standard and not restrictive. Assumption (SAEM4) and the compacity assumption of (SAEM5) are the most restrictive and not really realistic. They could be relaxed using a principle of random boundaries presented in [START_REF] Allassonnière | Construction of bayesian deformable models via stochastic approximation algorithm: A convergence study[END_REF]. Assumption (SAEM5) is verified by the PMCMC algorithm depending on the proposal distributions. Indeed, the Lispchitz property holds if the complete likelihood is continuously derivable which is the case under (M2-M3) and if θ remains in a compact set, which is realistic in practice. About the uniform ergodicity, [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] prove that if the noise density p(y|X; σ) is bounded above

sup y,X p(y|X; σ) < M σ , (7)
PMCMC inherits the convergence properties of the corresponding ideal MCMC algorithm.

For instance, if q(φ c |φ) = p(φ c ; β) then the kernel of the ideal MCMC q(X c 0:J , φ c ) = p(φ c ; β)p(X c 0:J |y 0:J , φ c ; γ, σ) is independent. For this kernel, the ratio p(X0:J ,φ|y0:J ) q(X0:J ,φ)

is bounded if (7) holds, ensuring the uniform ergodicity [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF].

Theorem 1. Assume that (M0-5), (SAEM1-5) hold. Under the assumption that {s m } m≥0 takes its values in a compact subset, the sequence θ m supplied by the SAEM-PMCMC algorithm converges a.s. towards a (local) maximum of the log-likelihood (θ) = log p(y; θ).

Remark that despite the SMC approximation in the PMCMC algorithm, the fact that the marginal stationary distribution of MCMC is the exact conditional distribution p(X, Φ|y; θ) is sufficient to prove the convergence of SAEM-PMCMC.

SDE with explicit transition density: simulation study

The respective performances of the SAEM implemented with a standard MCMC and of the SAEM combined with the PMCMC are compared to two models of various complexity: the Ornstein-Uhlenbeck process and the time-inhomogeneous Gompertz process. In order to validate our results, we perform a large scale simulation study with various sets of parameters and designs (n, J). Moreover, we study the influence of the number of particles and of the proposals involved the SMC algorithm. We compare the results using a biais and RMSE criteria. More precisely, for each condition (parameters, design, number of iterations) 100 datasets are generated. The corresponding estimate is obtained on each data set using both the SAEM-PMCMC algorithm and the SAEM-MCMC algorithm. Let θ r denote the estimate of θ obtained on the r-th simulated dataset, for r = 1, . . . , 100 by the corresponding algorithm. The relative bias 1 100 r ( θ rθ)/θ and relative root mean square error (RMSE) 1 100 r ( θ rθ) 2 /θ 2 for each component of θ are computed for both algorithms. The SAEM algorithm requires initial value θ 0 and the choice of the sequence (α m ) m≥0 . The initial values are chosen arbitrarily as the convergence of the SAEM algorithm little depends on the initialization. The step of the stochastic approximation scheme is chosen as recommended by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]: α m = 1 during the first iterations 1 ≤ m ≤ M 1 , and α m = 1 (m-M1) 0.8 during the subsequent ones. Indeed, the initial guess θ 0 might be far from the maximum likelihood value and the first iterations with α m = 1 allow the sequence of estimates to converge to a neighborhood of the maximum likelihood estimate. Subsequently, smaller step sizes during M -M 1 additional iterations ensure the almost sure convergence of the algorithm to the maximum likelihood estimate. We implement the SAEM algorithm with M 1 = 60 and M = 100 iterations.

Example 1: Ornstein-Uhlenbeck process 4.1.1. Ornstein-Uhlenbeck mixed model

The Ornstein-Uhlenbeck process has been widely used in neuronal modeling, biology, and finance (see e.g. [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Consider an SDE mixed model driven by the Ornstein-Uhlenbeck process and an additive error model

y ij = X ij + ε ij , ε ij ∼ N (0, σ 2 ), dX it = - X it τ i -κ i dt + γdB it , X 0 = 0 where κ i ∈ R, τ i > 0. We set φ i = (log(τ i ), κ i ) the vector of individual random parameters. We assume that log(τ i ) ∼ i.i.d. N (log(τ ), ω 2 τ ), κ i ∼ i.i.d. N (κ, ω 2 κ ).
The parameter vector is θ = (log τ, κ, ω τ , ω κ , γ, σ). This model can be easily discretized, resulting in a state-space model with random individual parameters:

y ij = X ij + ε ij , ε ij ∼ N (0, σ 2 ), X ij = X i ij-1 e -∆ijτi -κ i τ i (1 -e -∆ij τi ) + η ij , η ij ∼ N 0, γ 2 τ i (1 -e -∆ij τi ) .
The vector X i = (x i1 , . . . , x iJ ) conditional on φ i is Gaussian with mean vector m iX and covariance matrix G iX equal to

m iX = τ i κ i 1 -e - t i1 τ i , . . . , τ i κ i 1 -e - t iJ τ i , G iX = τiγ 2 2 1 -e - 2 min(t ij ,t ik ) τ i e - |t ij -t ik | τ i 1≤j,k≤J , (8)
where is the transposed vector. Altough this SDE is linear, the Gaussian transition density

p(X ij , ∆ ij , t ij-1 |X ij-1 , φ i ; θ)
is a nonlinear function of φ i . Thus, the likelihood has no closed form and the exact estimator of θ is unavailable. This example is a toy example. First, we compare the performances of the SAEM algorithm implemented with a simple MCMC to those of the SAEM-PMCMC algorithm. Next, we compare the influence of the number of particles K and the choice of the proposal distribution q(X tij | X ij-1 , y ij , φ i ; θ) in the SAEM-PMCMC.

Simulation step of the SAEM algorithm

Whereas the conditional distribution p(X i,1:J |, y i,0:J , φ i ; θ) is Gaussian, the joint distribution p(X i,1:J , φ i |, y i,0:J , φ i ; θ) is not explicitly and we have to resort to a approximate simulation at the S-step.

A first solution is to implement a standard MCMC algorithm, alternatively simulating under the distributions p(X i,1:J |, y i,0:J , φ i ; θ) and p(φ i |, y i,0:J , X i,1:J ; θ) for each subject. The posterior distribution p(X i,1:J |, y i,0:J , φ i ; θ) is Gaussian with easily computable mean vector and variance matrix derived from (8). Similarly, the posterior distribution of κ i is Gaussian with explicit mean and variance. On the other hand, the posterior distribution of τ i is not explicit and we use a Metropolis-Hastings step with a random walk proposal.

If we consider implementing a PMCMC kernel at the S-step, we have to choose two proposals q(X tij | X ij-1 , y ij , φ i ; θ) and q(•|φ i ). In this particular linear example, the proposal q(X ij | X ij-1 , y ij , φ i ; θ) can be the optimal proposal, the exact posterior density p(X ij |X ij-1 , y ij , φ i ; θ), which minimizes the variance of the particle weights. Indeed, this distribution is explicit for the Ornstein-Uhlenbeck, Gaussian with conditional mean and variance easily computable.

As an alternative, we also consider the transition density p(X ij |X ij-1 , φ i ; θ) as proposal.

The proposal q(•|φ) for the individual parameters φ within the PMCMC algorithm is a classical random walk on each component of vector φ.

Maximization step of the SAEM algorithm

The SAEM algorithm is based on the computation of the sufficient statistics for the maximization step. The statistics for the parameters µ = (log τ, κ), Ω = diag(ω 2 τ , ω 2 κ ) and σ are the three classic ones for mixed models (see e.g. [START_REF] Samson | The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model[END_REF]:

S 1 (y, X, Φ) = n i=1 φ i , S 2 (y, X, Φ) = n i=1 φ i φ i , S 3 (y, X, Φ) = n i=1 (y i -X i ) (y i -X i ).
Let s 1m , s 2m , s 3m denote the corresponding stochastic approximated conditional expectations at iteration m of SAEM. The M step for these parameters reduces to

µ (m) = 1 n s 1m Ω (m) = 1 n s 2m - 1 n 2 s 1m s 1m σ (m) = 1 n s 3m
The sufficient statistic corresponding to the parameter γ depends on the SDE. For the Ornstein-Uhlenbeck, we have

S 4 (y, X, Φ) = n i=1 J j=1 X i itij -X i itij-1 e -∆ijτi -κ i τ i (1 -e -∆ij τi ) 2 .
The M step is thus γ (m) = 1 nJ s 4m . Convergence assumptions (M0-5), (SAEM1-2) are easily checked on this example. (SAEM3) is implied by (SAEM4). ( SAEM4) is not theoretically verified even if in practice, this is not a problem. Finally, given our choice of proposals for PMCMC, (SAEM5) also holds, the compacity being verified in practice.

Simulation design and results

Two different designs are used for the simulations with equally spaced observation times: n = 20, J = 40, ∆ =0 .5 and n = 40, J = 20, ∆ = 1, t 0 = 0. Two sets of parameter values are used. The first set is log(τ ) = 0.6, κ = 1, ω τ = 0.1, ω κ = 0.1, γ = 0.05, σ = 0.05. The second set uses greater variances and is log(τ ) = log(10), κ = 1, ω τ = 0.1, ω κ = 0.1, γ = 1, σ = 1. For each design and each set of parameter values, one hundred datasets are simulated. The SAEM-PMCMC and SAEM-MCMC algorithms are initialized. For the first set of parameters, we set log(τ ) 0 = 1.1, κ 0 = 1.5, ω τ 0 = 0.5, ω κ0 = 0.5, γ 0 = 0.25 and σ 0 = 0.25, i.e. the initial standard deviations are 5 times greater than the true standard deviations. For the second set of parameters, we set log(τ ) 0 = log(10) + 0.5, κ 0 = 1.5, ω τ 0 = 0.5, ω κ0 = 0.5, γ 0 = 5 and σ 0 = 5. The SAEM-MCMC algorithm is implemented with N = 100 MCMC iterations. Several values of N and K are used in the SAEM-PMCMC algorithm. Figure 1 presents the convergence of the SAEM-PMCMC algorithm for one dataset simulated with the second set of parameters, n = 40 and J = 20. This illustrates the low dependence of the initialization of SAEM and the quick convergence in a small neighborhood of the maximum likelihood. Table 1 presents the bias and RMSE (%) of the SAEM-PMCMC and SAEM-MCMC algorithms obtained for the two designs and two sets of parameters. The results are almost identical for both algorithms and very satisfactory (bias less than 3% and RMSE less than 10%). Figure 2 shows the density of the estimators obtained by both algorithms with n = 40, J = 20 and the second set of parameters. The SAEM-PMCMC algorithm has slightly greater bias than SAEM-MCMC for the standard deviations of the random effects (ω τ and ω κ ). On the contrary, the bias for γ and σ are always lower with SAEM-PMCMC than with SAEM-MCMC. In this particular linear model, the standard MCMC on (X, Φ) is "ideal" as the distribution p(X i,0:J |y i,0:J , φ i ) can be sampled exactly. As a consequence, the SAEM-PMCMC can not be hopped to perform better. However, results show that the use of a particle approximation of p(X i,0:J |y i,0:J , φ i ) does not deteriorate the quality of estimation. Besides, note that the design (n, J) of the study affects very little the estimation quality when γ and σ are small (first set of parameters). When γ and σ are larger (second set of parameters), the RMSE are greater with fewer measures per subject (n = 40, J = 20) than with more (n = 20, J = 40). In a second part, we study the influence of the particles number K in the SMC algorithm. Table 2 presents the bias and RMSE obtained on the 100 datasets simulated with the first set of parameters, n = 40, J = 20 and different implementations of the SAEM-PMCMC algorithm. We successively use K = 25, K = 50 and K = 100 particles and either N = 50, N = 100 or N = 200 PMCMC iterations, using the exact posterior distribution q(X ij |X ij-1 , y ij , φ i ; θ) = p(X ij |X ij-1 , y ij , φ i ; θ). The results are almost identical. On the contrary, the choice of the proposal q(X ij |y ij , X ij-1 , φ; γ, σ) in the SMC algorithm may affect the result. Table 2 presents results obtained with the transition density p(X ij |X ij-1 , φ; γ, σ) as proposal q. The algorithm fails to converge due to the degeneracy of the particles. We present the results obtained by increasing significantly the number of par- ticles (K = 1000), but we observe the degeneracy of the algorithm on 90% of the simulated datasets. Morevoer, the bias and RMSE are greater than with the posterior conditional distribution p(X ij |y ij , X ij-1 , φ; γ, σ).

In conclusion, although the proposal q seems to be crucial in these simulations, the number of particles K has less influence on the convergence of the estimation algorithm. This is in concordance with the theoretical results proved by [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF].

However, as emphasized by [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF], increasing K can make arbitrarily small the probability of visiting unfavourable states by the Markov chain, which correspond to large values of the ratio target density to proposal density. But our simulations show that when choosing the "optimal" proposal density q, namely the exact posterior distribution, the ratio of the target density to the proposal density seems to have no larges values and thus the number of particles K has a very low influence. This low influence of K can also be explained by the fact that the SAEM only requires few iterations of PMCMC, without convergence of the Markov chain to the stationary distribution. The convergence is made over the iterations of SAEM.

(Xt ij |Xt ij-1 , yij, φi; θ) = p(Xt ij |Xt ij-1 , yij
(Xt ij |Xt ij-1 , yij, φi; θ) = p(Xt ij |Xt ij-1 , φi; θ), K =

Example 2: stochastic volatility Gompertz model 4.2.1. Stochastic volatility Gompertz mixed model

The Gompertz model is a well-known growth model (see e.g. [START_REF] Jaffrézic | Modelling variances with random effects in non linear mixed models with an example in growth curve analysis[END_REF]. Recently, [START_REF] Donnet | Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[END_REF] proposed a stochastic version of this model to take into account random fluctuations in growth process. The stochastic volatility Gompertz mixed model is the following one:

y ij = X ij (1 + ε ij ), ε ij ∼ N (0, σ 2 ), dX it = B i C i e -Cit X it dt + γX it dB it , X i0 = A i e -Bi (9)
where

A i > 0, B i > 0, C i > 0.
We set φ i = (log A i , log B i , log C i ) the vector of individual random parameters. The initial conditional is a function of the random individual parameters φ i . We assume that log

A i ∼ i.i.d. N (log A, ω 2 A ), log B i ∼ i.i.d. N (log B, ω B ), log C i ∼ i.i.d. N (log C, ω 2 C ) and the parameter of interest is θ = (log A, log B, log C, ω A , ω B , ω C , γ, σ)
. By the Itô formula, the conditional expectation and variance of log(X it )|φ i , for t ≥ 0, are

E(log X it |φ i ) = log A i -B i e -Cit - 1 2 γ 2 t, Var(log X it |φ i ) = γ 2 t.
The transition density of (log(X t )) is explicit, Gaussian and equals to

log p(log X ij , ∆ ij , t ij-1 | log X ij-1 , φ i ; θ) = -1 2 log(2πγ 2 ∆ ij ) -1 2 (log Xij -log Xij-1+Bie -Ct ij-1 (e -C i ∆ ij -1)-1 2 γ 2 ∆ij ) 2 γ 2 ∆ij .
which is a nonlinear function of φ i . As a consequence, the likelihood has no closed form and the exact estimator of θ is unavailable.

Simulation step of the SAEM algorithm

First, the S step of the SAEM algorithm is tackled with a standard MCMC algorithm. Due to the multiplicative structure of the observation noise, the distribution p(X i,0:J |y i,0:J , φ i , θ) is no more explicit and we have to resort to a random walk proposal to simulate the trajectories X i,0:J conditionally to the observations y i,0:J . The components of X i,0:J are updated time by time using a random walk proposal. A random walk proposal is also used to simulate the individual parameters φ i from the conditional distribution. Secondly, we implement the SAEM-PMCMC algorithm and the proposal q(X ij |X ij-1 , y ij , φ i ; θ) has to be chosen. We propose to approximate the ideal proposal p(X ij |X ij-1 , y ij , φ i ; θ) by a Gaussian distribution with mean and variance deduced from the true ones. More precisely, we consider the following proposal on log X tij

q(log X tij | log X ij-1 , y ij , φ i ; θ) = N (m Xij ,post , Γ Xij ,post ), with Γ Xij ,post = σ -2 + (γ 2 ∆ ij ) -1 -1 , m Xij ,post = Γ Xij ,post µ Xij ,post , µ Xij ,post = log y ij σ 2 + 1 γ 2 ∆ ij log X ij-1 -B i e -Citij-1 (e -Ci∆ij -1) - γ 2 ∆ ij 2 .
and then take the exponential to obtain a candidate for X tij . The proposal q(•|φ) for the individual parameters φ within the PMCMC algorithms is a classical random walk on each component of the vector φ.

Maximization step of the SAEM algorithm

Within the SAEM algorithm, for parameters µ = (log A, log B, log C), Ω = diag(ω 2 A , ω 2 B , ω 2 C ) and σ, the sufficient statistics and the M step are the same as in Example 1. As the parameter γ appears both in the expectation and the variance of log(X t )|φ, its estimator is not the same as in Example 1. The sufficient statistic corresponding to γ is

S 4 (y, log X, Φ) = n i=1 J i=1 ∆ ij log X itij -log X itij-1 + B i e -Citij-1 (e -Ci∆ij -1) 2 .
Let s 4m denote the stochastic approximation of this sufficient statistic at iteration m of the SAEM algorithm. For the sake of simplicity, we assume that the step size ∆ ij is a constant ∆. The estimator γ m at iteration m is deduced by maximizing the complete likelihood:

γ m = 2 ∆ -1 + 1 + s 4m nJ
When ∆ ij is not a constant, the estimator is more complex but also explicit.

Simulation results

We propose to simulate two experimental designs with the same observation times for all individuals: n = 20, J = 40, ∆ = 0.01 and n = 40, J = 20, ∆ = 0.02 (t 0 = 0), respectively. As in [START_REF] Donnet | Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[END_REF], we use the following population parameters: log A = log(3000), log B = log(5), log C = log(14), ω A = 0.1, ω B = 0.1, ω C = 0.1. Two levels of variability γ and σ are taken respectively equal to 0.4, 1/ √ 20 and 0.8, 1/ √ 15 and referred to as first and second set of parameters thereafter. For each design and each set of parameters, 100 datasets are simulated and θ is estimated with the two algorithms (SAEM-PMCMC and SAEM-MCMC). The SAEM algorithm is initialized with log A 0 = 8.21, log B 0 = 1.81, log C 0 = 2.84, ω A0 = 0.5, ω B 0 = 0.5, ω C 0 = 0.5, γ 0 = 1.2, σ 0 = 0.66. The PMCMC is implemented with N = 100 iterations and K = 50 particles. The MCMC is implemented with N = 100 iterations. Bias and RMSE (%) are presented in Table 3. Both algorithms give satisfactory results overall (bias and RSME ≤ 20% except for γ). Nevertheless, we clearly see that the SAEM-PMCMC algorithm performs better than the SAEM-MCMC algorithm. Indeed, for the population means (log A, log B, log C), the population standard deviations (ω A , ω B , ω C ) and σ, the bias and RMSE are smaller for SAEM-PMCMC even if SAEM-MCMC results are already small. The difference is obvious for γ, which is estimated with a large bias by SAEM-MCMC (up to a bias of 77%) while this is not the case for SAEM-PMCMC (bias < 7%). This phenomenon is accentuated when the variabilities γ and σ increase (second set of parameters). The better performances of SAEM-PMCMC versus SAEM-MCMC can also be observed in Figure 3 where we plot the densities of the parameters θ among the 100 simulated datasets for the second set of parameters with n = 20 and J = 40. The improved quality of estimation has a computational cost since it takes less than 8 minutes for the SAEM-MCMC and nearly 23 minutes for the SAEM-PMCMC algorithm to supply results. However, estimating precisely the population parameters is a key issue in population studies (such a genetic specification) and a bias can have real consequences on the conclusions.

Theoretical results for general SDE requiring Euler approximation

The class of SDEs considered under assumption (M0) is limited. For a general SDE, we propose to approach the process by an Euler-Maruyama scheme, leading to an approximate model with an explicit (Gaussian) transition density.

Euler approximation

For the sake of simplicity, we now consider that ∀i = 1 . . . n, ∀j = 0 . . . J, t ij = j and that the observations times are regularly spaced with interval ∆. Let r be an integer ≥ 1. On each interval time [t j , t j+1 [, we define an approximate diffusion by the recursive Euler-Maruyama scheme of step size 1/r by : X i0 (r) ∼ π 0 (•|φ i ) and

X ij (0) = X ij-1 (r), X ij (u + 1) = X ij (u) + 1 r a(X ij (u), t j-1 + u/r, φ i ) (10) + 1 r b(X ij (u), t j-1 + u/r, φ i , γ)U u
where the variables U u are iid unit centered Gaussian. Then X ij (r) is a realization of the approximate diffusion at time t j . We consider the error model

y ij = g(X ij (r), ε ij ) (11)
Equations ( 10), ( 11) and ( 3) define the so-called approximate mixed model. Under this model, the likelihood of the observations y is :

p (r) (y; θ) = n i=1
p(y i,0:J |X i,0:J (r); σ)p (r) (X i,0:J (r)|φ i ; γ)p(φ i ; β)dφ i dX i,0:J (r),

where p (r) (X i,0:J (r)|φ i ; θ) is the density of the Euler approximation X ij (r) given φ i .

SAEM-PMCMC on the approximate Euler model

A first approach to estimate θ is to use the algorithm SAEM-PMCMC presented in Section 3 on the approximate model. In that case, Theorem 1 implies the following corollary:

Corollary 1. Assume that (M1-5), (SAEM1-5) hold for the approximate model. Under the assumption that {s m } m≥0 takes its values in a compact subset, the sequence { θ m } m≥0 generated by the algorithm SAEM-PMCMC on the approximate Euler model converges to a (local) maximum of the likelihood p (r) (y; θ) of the approximate model.

The main limit of this result is that the SAEM-PMCMC converges to the maximum of the Euler approximate likelihood p (r) (y; θ) instead of the maximum of the exact likelihood p(y; θ). [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] prove that there exists a constant C such that the total variation distance between p (r) (y; θ) and the exact likelihood p(y; θ) can be bounded by

||p (r) (y; θ) -p(y; θ)|| T V ≤ C r
This results says that the two likelihoods are close but this does not imply that the two maxima are close. A solution would be to consider r → ∞, so that the Euler approximation converges to the exact diffusion. This requires r to increase along the iterations of the SAEM algorithm, also implying the Euler approximate model and the likelihood to change with the iterations. The concept of stationary distribution of PMCMC is no more valid in this context. We could imagine a reversible jump approach for the simulation step of the SAEM algorithm. But even in this case, as the approximate likelihood changes with the iterations, it is not possible to prove the convergence of the SAEM algorithm.

From a practical point of view, the choice of r is sensitive because it adds intermediate latent variables to the model. A large value of r is prefered as the Euler approximate likelihood is closer to the exact likelihood. However, it is well known that a large volume of latent variables is difficult to handle within a MCMC or PMCMC approach. This has been discussed by [START_REF] Roberts | On inference for partially observed nonlinear diffusion models using the Metropolis-Hastings algorithm[END_REF] in a general diffusion context and by [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] for SDE mixed models. Furthermore, as emphasized by [START_REF] Roberts | On inference for partially observed nonlinear diffusion models using the Metropolis-Hastings algorithm[END_REF], when r → ∞, MCMC algorithms with r additional latent variables are unable to correctly estimate γ. Indeed, in that case, the density p (r) (x j |x j-1 , φ; γ) of the Euler approximate diffusion converges to the likelihood of the continuous path (X t ) t∈ [tj-1,tj] given by the Girsanov formula. It is well-known that it is not possible to estimate the volatility from this continuous path likelihood.

Exact theoretical convergence for a SAEM algorithm coupled with a SMC algorithm

To circumvent the problem of the estimation of γ, we restrict to the estimation of

θ = (β, σ)
with γ assumed to be known. To estimate θ, we propose a theoretical algorithm coupling SAEM with a naive SMC. This algorithm is known to have bad numerical properties (this is the reason why we do not implement it) but we prove that it supplies a sequence ( θ m ) m∈N converging a.s. towards a local maximum of the exact likelihood log p(y; θ). This theoretical result yields to a consequent improvement of Corollary 1. Now, we present the SMC that we propose to use in the simulation step of the SAEM algorithm. In the literature, SMC algorithms have originally a filtering purpose to approximate the distribution p(X 0:J |φ, y 0:J ; θ). However, they have been rapidly extended to the case where the parameters φ are also sampled. In the following, we present a naive SMC extension targeting the distribution p (r) (X 0:J , φ|y 0:J ; θ) of the approximate model.

Algorithm 4 (Naive SMC algorithm with parameters sampling).

• Time j = 0:

∀ k = 1, . . . , K, sample φ (k) 0 ∼ p(•; β) and X (k) 0 ∼ π 0 (•|φ (k) 0
) and compute the weights:

W 0 φ (k) 0 , X (k) 0 ∝ p (r) y 0 , X (k) 0 , φ (k) 0 ; θ • Time j = 1, . . . , J: sample K iid variables (φ (k) j-1 , X (k) 0:j-1 ) according to the law Ψ K j-1 = K k=1 W (k) 0:j-1 1 φ (k) j-1 ,X (k) 0:j-1 . Then for each k = 1, . . . , K, sample X (k) j ∼ q (r) •|y j , X (k) 0:j-1 , φ (k) j-1 ; γ, σ and set φ (k) j = φ (k) j-1 and X (k) 0:j = (X (k) 0:j-1 , X (k) j ).
Finally compute and normalize the weights:

W j φ (k) j , X (k) 0:j ∝ p (r) y j , X (k) j |y 0:j-1 , X (k) 0:j-1 φ (k) j-1 ; γ, σ q (r) X (k) j |y j , X (k) 0:j-1 , φ (k) j-1 ; γ, σ ,
We consider a SAEM algorithm in which the simulation step is performed by this naive SMC sampler. To ensure the convergence of this estimation algorithm towards the exact maximum likelihood, we need the number of particles K and the step size 1/r of the Euler approximation to vary along the iterations of the SAEM algorithm. Thus K m and r m denote the number of particles and the step size of the Euler scheme at iteration m of the SAEM algorithm. We denote by ( θ m ) m∈N the sequence of estimates supplied by this algorithm. Technical assumptions (E1-4) are required to prove the convergence of the naive SMC sampler and are presented in Appendix.

Theorem 2. Assume that there exists a constant c > 1 such that K m varies along the iterations:

K m = O (log(m) c ) and that r m = min{r ∈ N | r ≥ √ K m }.
Assume that (M1-5), (SAEM1-4), (E1-4) hold for the exact model. Then, with probability 1, lim m→∞ d( θ m , L) = 0 where L = {θ ∈ Θ, ∂ θ (θ) = 0} is the set of stationary points of the exact log-likelihood (θ) = log p(y; θ).

Theorem 2 is proved in Appendix 7.2. This theoretical result is noteworthy. Indeed, even if the SAEM-SMC algorithm is performed on the approximate Euler model, the algorithm converges to the maximum of the exact likelihood. This is due to the convergence of the Euler approximate SMC towards the exact filter distribution. This powerful result has first been proved by Del [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF]. We propose an extension of their results to our SMC (see Lemma 1 in Appendix 7.1). Then, by generalizing the proof of convergence of the SAEM algorithm to an approximate simulation step, we are able to deduce the convergence of the estimates to the maximum of the exact likelihood. The naive SMC algorithm 4 provides an asymptotically consistent estimate of the target distribution p (r) (X 0:J , φ|y 0:J ; θ) under very weak assumptions but has proved bad properties in practice (in particular a degeneracy of the parameters due to the resampling step). It can be improved by including a MCMC step on the parameters φ [START_REF] Doucet | An introduction to sequential Monte Carlo methods[END_REF]. Although these algorithms remains poorly efficient in practice (explaining the developement of the PMCMC algorithm for instance), they allow to prove this theorem.

Discussion

The stochastic differential mixed-effects models are quite widespread in the applied statistics field. However, in the absence of efficient and computationally reasonable estimation methods, a simplifying assumption is often made: either the observation noise or the volatility term are standardly neglected. In this paper we present an EM algorithm combined with a Particular Monte Carlo Markov Chain method to estimate parameters in stochastic differential mixed-effects models including observation noise. We prove the convergence of the algorithm towards the maximum likelihood estimator when the transition density is explicit. This proof is classical as the PMCMC acts as an exact marginal MCMC. On the contrary, when we consider SDE with unexplicit transition density, we prove the convergence of the SAEM-PMCMC algorithm to the maximum likelihood of an approximate model obtained by a Euler scheme. We are not able to prove the convergence of this algorithm to the exact maximum likelihood. But we propose the convergence of a SAEM-SMC algorithm to the exact maximum likelihood. The suggested method supplies accurate parameter estimation in a really moderate computational time on practical examples. Moreover, it is not restricted to homogeneous time SDE. On various simulated datasets, we illustrate the superiority of this method over the SAEM algorithm combined with a standard MCMC algorithm. This efficiency is due to the fact that the PMCMC algorithm takes advantage of the Markovian properties of the non-observed process. A major advantage of this methodology is its automatic implementation. Indeed the generation of the non-observed process does involves less tuning parameters than a standard MCMC, such as the size of the random move in the random walk Metropolis-Hastings algorithm. For example, the number of particles is not crucial, as illustrated in Example 1. In this paper we present a PMCMC algorithm where all the particles are re-sampled at each iteration. Many other resampling distributions have been proposed in the literature, trying to achieve an optimal procedure. West (1993) developed an effective method of adaptive importance sampling to address this issue. The procedure developed by [START_REF] Pitt | Auxiliary variable based particle filters[END_REF] is similar in spirit and has real computational advantages. A stratified resampling is proposed in [START_REF] Kitawaga | Monte carlo filter and smoother for non-gaussian nonlinear state space models[END_REF]. Each resampling distribution implies a specific update of the weights. Tornøe, C., R. Overgaard, H. Agersø, H. Nielsen, H. Madsen, and E. Jonsson (2005).

Stochastic differential equations in NONMEM: implementation, application, and comparison with ordinary differential equations. Pharm. Res. 22 (8), 1247-58. West, M. (1993). Approximating posterior distributions by mixtures. J. Roy. Statist. Soc. Ser. B 55 (2), 409-422.

Wolfinger, R. (1993). Laplace's approximation for nonlinear mixed models. Biometrika 80 (4), 791-795.

Appendix

Convergence of the Euler approximate naive SMC

Let us introduce some notations. For any bounded Borel function f taking real values, we denote π jθ f = E (f (X 0:j , φ)|y 0:j ; θ) the conditional expectation under the exact smoothing distribution of the exact model p(X 0:j , φ|y 0:j ; θ) and Ψ are continuous on C. E3. f is a function of (φ, X 0:J ) taking its values into R, uniformly bounded E4. There exist κ 1 and κ 2 independent of σ such that p(y|x, σ) ≤ κ 1 and such that p (r) y j , X Then, for any ε > 0, for any j = 1, . . . , J, there exist constants C 1 and C 2 independent of θ such that P Ψ

Ψ (K) j,θ f -Ψ (K) j,θ f = 1 K K k=1 η k -E(η k |G) with η k = f (φ (k) j , X (k
) 0:j ). By deviation inequalities, under assumption (E3), we obtain that there exist two constants such that

P(|Ψ (K) j,θ f -Ψ (K) j,θ f | ≥ δ) ≤ A 1 exp -K A 2 δ 2 f 2
Next, we introduce the kernels H j such that :

H j f (φ, X 0:j-1 ) = p(y j |x j )p(x j |X 0:j-1 , φ)p(φ)f (φ, X 0:j-1 , x j )dφdx j

Thus we have

π j f = π j-1 H j f π j-1 H j 1
and the same definitions on the Euler approximate model, with kernels H K j by H K j f (φ, X 0:j-1 ) = p(y j |x j )p (r) (x j |X 0:j-1 , φ)p(φ)f (φ, X 0:j-1 , x j )dφdx j such that π K j f = E r (f (φ, X 0:j-1 )|y 0:j ; θ) the conditional expectation under the approximate model is

π K j f = π K j-1 H K j f π K j-1 H K j 1
Then, considering the σ-filed G generated by the variables φ Assumption (E4) implies that there exist two constants C 1 , C 2 such that

|η k | ≤ C 1 ||f ||, E(|η k | 2 |G) ≤ C 2 ||f || 2
Deviation inequalities imply that there exists a constant C such that

P W (k) 0:j f -Ψ K j H j f ≥ δ ≤ 2e -K δ 2 C||f || 2
With similar arguments than Del [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] and assumptions (E1-2), we prove by induction that there exist two constants A 1 and A 2 such that

P(|π jθ f -Ψ (K) j,θ f | ≥ δ) ≤ A 1 exp -K δ 2 A 2 f 2
where A 1 = 3 j+2 , A 2 = C(8ρ) j+1 with C = 2 + 2α and ρ = 2κ J+1 /ε J+1 , where α is the constant quantifying the error induced by the Euler scheme and ε J+1 = H J . . . H 1 1. With regularities assumptions and assumptions (E1-2), we can prove that the constants are independent of θ.

Figure 1 .

 1 Figure 1. Convergence of the SAEM-PMCMC estimates log(τ ) m , b κm, c ωτ m , c ωκ m , b γm and b σm along the SAEM iterations for one dataset simulated with the Ornstein-Uhlenbeck mixed model, with the second set of parameters, n = 40 and J = 20. Horizontal lines represent the true values.

  , φi; θ) and K = 25, K = 50 or K = 100 particles, N = 50, N = 100 or N = 200 PMCMC iterations and prior proposal q

Figure 2 .

 2 Figure 2. Ornstein-Uhlenbeck mixed model, n = 40, J = 20 and second set of parameters. Densities of the estimators bθ obtained with the SAEM-PMCMC (plain line) and SAEM-MCMC algorithms (dashed line) on 100 simulated datasets, for each component of (log τ, κ, ω τ , ωκ, γ, σ). The true parameter value is the vertical line.

Figure 3 .

 3 Figure 3. Stochastic volatility Gompertz mixed model, n = 20, J = 40 and second set of parameters. Densities of the estimators b θ obtained with the SAEM-PMCMC (plain line) and SAEM-MCMC algorithms (dashed line) on 100 simulated datasets, for each component of (log A, log B, log C, ωA, ωB, ωC, γ, σ). The true parameter value is th vertical line.

  f (φ(k) , X (k) 0:J ) the conditional expectation under the naive SMC distribution of the approximate model with K particles and a Euler step size r = min{r ∈ N | r ≥ √ K}. Lemma 1. Assume that E1. The functions a and b are two times differentiable with derivatives w.r.t x and φ all orders up to 2 uniformly bounded w.r.t φ E2. Θ is included in a compact set C. Besides, the functions θ → π(φ; β) and σ → p(y|x; σ)

  fπ j,θ f ≥ δ ≤ C 1 exp -K δ 2 C 2Proof. Following Del[START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] proof, let us define the empirical measure Ψ has been defined in algorithm 4. Set the σ-filed G generated by the variables φ

.

  We have to bound |η k | and E(|η k | 2 |G).

Table 1 .

 1 Ornstein-Uhlenbeck mixed model: biais and RMSE (%) for b θ obtained by the SAEM-PMCMC and SAEM-MCMC algorithms on 100 simulated datasets with two designs (n = 20, J = 40 and n = 40, J = 20) and two sets of parameters. PMCMC is implemented with K = 50 particles and the exact posterior proposal q(Xt ij |Xt ij-1 , yij, φi; θ) = p(Xt ij |Xt ij-1 , yij, φi; θ).

	Parameters		log(τ )	κ	ωτ	ωκ	γ	σ
	True value		0.62	1.00	0.10	0.10	0.05	0.05
	n = 20, J = 40
	SAEM-PMCMC	Bias	0.38	0.05 -2.55 -3.29 -0.84	0.07
		RMSE	0.48	0.30	1.22	1.34	0.93	0.44
	SAEM-MCMC	Bias	-0.02	0.29 -0.60 -2.38	3.34 -1.43
		RMSE	0.49	0.30	1.23	1.29	0.96	0.45
	n = 40, J = 20
	SAEM-PMCMC	Bias	0.11	0.01 -2.75 -1.59 -0.73	0.50
		RMSE	0.33	0.20	0.82	0.99	0.94	0.52
	SAEM-MCMC	Bias	-0.24	0.22 -1.46 -1.12	2.96 -1.37
		RMSE	0.34	0.21	0.81	1.02	0.96	0.55
	True value		2.30	1.00	0.10	0.10	1.00	1.00
	n = 20, J = 40
	SAEM-PMCMC	Bias	0.55 -0.63 -1.49 -1.92	0.15 -1.33
		RMSE	0.42	0.97	0.51	0.55	0.64	0.81
	SAEM-MCMC	Bias	-0.41	1.24 -0.62 -1.25	0.44 -1.61
		RMSE	0.35	0.84	0.52	0.54	0.56	0.72
	n = 40, J = 20
	SAEM-PMCMC	Bias	0.28 -0.18 -0.73 -1.04 -0.47 -0.03
		RMSE	3.19	7.03	3.91	4.38	5.45	9.81
	SAEM-MCMC	Bias	0.26 -0.19 -0.99 -1.17 -1.59	2.64
		RMSE	2.32	5.39	3.89	4.05	4.96	9.42

Table 2 .

 2 Ornstein-Uhlenbeck mixed model, n = 20, J = 40: biais and RMSE (%) for b θ obtained by the SAEM-PMCMC algorithm on 100 simulated datasets. The PMCMC algorithm is successively implemented with posterior proposal q

  1000 particles and N = 100 iterations.

	Parameters log(τ )	κ	ωτ	ωκ	γ	σ
	True value	0.62 1.00	0.10	0.10	0.05	0.05
	SAEM-PMCMC with posterior proposal, K = 25, N = 100
	Bias	0.26 0.11	0.32	-4.84	0.64 -0.41
	RMSE	0.49 0.30	1.28	1.30	1.04	0.45
	SAEM-PMCMC with posterior proposal, K = 50, N = 100
	Bias	0.38 0.05 -2.55	-3.29	-0.84	0.07
	RMSE	0.48 0.30	1.22	1.34	0.93	0.44
	SAEM-PMCMC with posterior proposal, K = 100, N = 100
	Bias	0.04 0.24 -0.97	-4.11	6.86 -2.68
	RMSE	0.48 0.29	1.24	1.29	1.15	0.51
	SAEM-PMCMC with posterior proposal, K = 50, N = 50
	Bias	0.28 0.04	4.07 -11.61	10.06 -3.32
	RMSE	0.49 0.30	1.40	1.55	1.38	0.55
	SAEM-PMCMC with posterior proposal, K = 50, N = 200
	Bias	-0.06 0.31 -1.27	-2.68	6.10 -2.36
	RMSE	0.48 0.30	1.16	1.30	1.15	0.49
	SAEM-PMCMC with prior proposal, K = 1000, N = 100
	Bias	-0.63 0.66 -2.78	-4.73	15.42 -5.44
	RMSE	0.53 0.34	1.43	1.77	2.17	0.83

Table 3 .

 3 Stochastic volatility Gompertz mixed model: bias and RMSE (%) of b θ obtained by the SAEM-PMCMC and SAEM-MCMC algorithms, on the 100 datasets with two designs (n = 20, J = 40 and n = 40, J = 20) and two sets of parameters.

	Parameters		log A log B log C	ωA	ωB	ωC	γ	σ
	True value		8.01	1.61	2.64	0.10	0.10	0.10	0.40	0.22
	n = 20, J = 40
	SAEM-PMCMC	Bias	0.01	-0.07	0.03 -0.69	-2.31	-1.79	2.77 -0.04
		RMSE	0.61	1.16	0.91	2.90	9.71	5.98 11.15	3.15
	SAEM-MCMC	Bias	2.04 1.713 -1.99 -8.00	-5.00	-7.55 61.93 -2.56
		RMSE	2.16	2.07	2.19	8.15	10.87	9.14 63.78	4.30
	n = 40, J = 20
	SAEM-PMCMC	Bias	0.05	-0.30	0.12 -1.92	-2.38	-1.41	6.79 -0.81
		RMSE 0.797	1.64	1.36	4.53	14.88	9.10 18.54	3.31
	SAEM-MCMC	Bias	2.37	1.90 -1.98 -6.76	-4.58	-8.17 77.66 -2.00
		RMSE	2.52	2.43	2.33	7.20	14.80	10.73 79.74	4.01
	True value		8.01	1.61	2.64	0.10	0.10	0.10	0.80	0.26
	n = 20, J = 40
	SAEM-PMCMC	Bias	0.56	0.72 -0.44 -1.11	-6.56	-1.97	1.67 -0.15
		RMSE	1.47	1.90	1.67	3.90	13.10	6.33 11.20	3.20
	SAEM-MCMC	Bias	4.76	4.64 -3.67 -7.58 -11.92 -11.49 62.90 -2.60
		RMSE	4.92	4.95	4.01	7.83	16.55	12.08 75.11	4.84
	n = 40, J = 20
	SAEM-PMCMC	Bias	0.36	0.38 -0.47 -0.70	-2.16	-1.31	0.37	0.82
		RMSE	1.07	1.41	1.43	2.21	8.84	4.64	9.55	3.57
	SAEM-MCMC	Biais	4.37	4.11 -3.60 -8.15	-8.05 -10.05 49.70 -4.33
		RMSE	4.46	4.31	3.81	8.26	11.38	10.55 52.36	5.64

Acknowledgments

Sophie Donnet has been financially supported by the french program ANR BigMC.

Proof of Theorem 2

Proof. At iteration m, the simulation step provides (X (m) , Φ (m) ) simulated with an approximate distribution

where Ψ (Km) i,J, θm is the empirical measure obtained by the naive SMC for individual i. The stochastic approximation step of the SAEM-SMC algorithm can be decomposed into:

J, θ(sm) S(y, X, Φ) -E(S(y, X, Φ)|y, θ(s m )) Following Theorem 2 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] on the convergence of Robbins-Monro scheme, the convergence of SAEM-SMC to the set of stationary points of the exact log-likelihood (θ) is ensured if we prove the following assertions:

exists and is finite with probability 1. (c) lim m→∞ R m = 0 with probability 1.

• To prove assertion 1, we first remark that, as γ is known and θ = (β, σ), we have

Thus under (M1-M5), (SAEM2), the functions h(s) = s( θ(s))-s and

With similar arguments than Lemma 2 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF], we can also prove that the set V ({s, F (s) = 0}) is of zero measure.

• Assertion 2 is proved following Theorem 5 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. Note F = {F m } m≥0 the increasing family of σ-algebra generated by all the random variables used by the SMC algorithm until the mth iteration. By construction of the simulation step, the simulated variables (X (1) , Φ (1) ), . . . , (X (m) , Φ (m) ) given θ 0 , . . . , θ m are independent. Consequently, we have that for any positive Borel function f ,

We can then prove that m k=1 α k e k is a martingale, bounded in L 2 under assumptions (M5) and (SAEM1-4).

• We now prove the almost sure convergence of R m (assertion 3). The term R m is of dimension d. It converges a.s. i.f.f. each of its component converges a.s. As a consequence, for the sake of simplicity, we suppose that d equals 1. Applying Borel-Cantelli lemma, it is sufficient to prove that m≥0 P (|R m | ≥ ) < ∞ for all > 0.

We have

J, θ(sm) Sπ J, θ(sm) S where π J,θ S = E(S(X, Φ)|y; θ) First, by independence of the individual i = 1 . . . n, we can decompose S into a sum of n terms S(X, Φ) = n i=1 S i (X i,0:J , φ i ). and write:

J, θ(sm) S iπ i,J, θ(sm) S i ≥ δ n with π i,J, θ(sm) S i = E(S i (X i,0:J , φ i )|y i,0:J ; θ). Now, we bound each term of the sum by Lemma 1 and under assumptions (SAEM4), (E1-E2). This implies the almost sure convergence of (R m ) towards 0.