Nathalie Verdière
email: verdiern@univ-lehavre.fr

Carine Jauberthie

Parameter estimation procedure based on input-output integro-differential polynomials. Application to the Hindmarsh-Rose model

This paper deals with a parameter estimation method based on input-output integro-differential polynomials. From the Rosenfeld-Groebner algorithm, some differential relations depending only on the inputs, the outputs and the parameters of the model are obtained. A pretreatment consisting in some integrations of these relations permits to obtain new ones. The latter contain essentially integrals depending only on the model inputs/outputs and their higher order derivatives are lower than the initial relations. Therefore, they are less sensitive to the noise on the measurements compared to the initial ones. Integrating permits also to attenuate the effect of the noise improving by the same the parameter estimates. However, even if the numerical estimation algorithm provides a very good value of the parameters, the latter can lead to an incorrect behavior of the model output. Indeed, in biological or physical applications, a little change of some parameter values can lead to a radical change of the model behavior as for the Hopf bifurcation. A Hopf bifurcation refers to a radical change of the model output dynamic due to a parameter crossing a reference value. Therefore, an algorithm is proposed in this paper to test the reliability of any parameter estimation procedure with respect to the dynamic of the system. More precisely, from a given noise on the output(s), it consists in calculating the probability that the result of a parameter estimation algorithm will permit to reproduce the correct behavior of the model output. Finally, this algorithm is applied on the estimation procedure based on the input-output integrodifferential polynomials and on the Hindmarsh-Rose model, a slow-fast model able to reproduce the main behaviors of a neuron and presenting a Hopf bifurcation.

I. INTRODUCTION

Parameter estimation of dynamical systems is an important step in the modeling to make a model usable. Several global or local methods can be found in the literature [START_REF] Kiranyaz | Optimization Techniques: An Overview[END_REF], [START_REF] Walter | Identification of parametric models from experimental data[END_REF]. They are mainly focused on the optimization of some criterion function over the parameter space and required a first initial guess or at least a definition interval. This first initial guess can be obtained by an industrial process or by a numerical procedure and may determine the convergence of a local parameter estimation algorithm. This paper proposes to focus on a parameter estimation method based on differential polynomials in order to obtain a first initial guess of the parameters without, a priori, any knowledge of them.

Afterwards, the general dynamical systems have the following form:

Γ θ ẋ(t, θ) = f (x(t, θ), θ) + u(t)g(x(t, θ), θ), y(t, θ) = h(x(t, θ), θ).

Here x(t, θ) ∈ IR n denotes the state variables. The model outputs are given by y(t, θ) ∈ IR m . The set of model parameters to be estimated are given by θ ∈ U p the unknown parameter vector (U p is an open subset in IR s ).

The functions f (., θ), g(., θ) and h(., θ) are real, rational and analytic for every θ ∈ U p on M (a connected open subset of IR n such that x(t, θ) ∈ M for every θ ∈ U p and every t ∈ [0, T ]). u(t) ∈ R l denotes the input variables. In the case of uncontrolled system, u is equal to 0. We also let x 0 = (x 0,i ) i=1,...,n the vector of initial conditions for x(t, θ) whose some components can depend on the parameters to be estimated.

From system Γ θ , some differential relations can be obtained in using the Rosenfeld-Groebner elimination algorithm and are called input-output (IO) polynomials [START_REF] Boulier | Computing representation for radicals of finitely generated differential ideals[END_REF]. From the latter, an identifiability study of Γ θ can be done and is recommended before putting in place a parameter estimation procedure. Indeed, this study ensures that the parameters can be uniquely inferred from the output measurements of the system. However, IO polynomials may contain sometimes high-order derivatives which is a factor of instability in numerical procedures. Indeed, estimate derivatives from noisy measurements is an ill-posed problem. The main challenge is also to obtain new relations with lower order derivatives. Lots of authors have proposed different ways to obtain integro-differential input-output (ID-IO) relations since they are less sensitive to the noise compared to the initial IO polynomials (for example [START_REF] Boulier | An algorithm for converting nonlinear differential equations to integral equations with an application to parameter estimation from noisy data[END_REF], [START_REF] Loeb | More about process identification[END_REF], [START_REF] Sira-Ramirez | Algebraic Identification and Estimation Methods[END_REF], [START_REF] Verdière | A distribution input-output polynomial approach for estimating parameters in nonlinear models, application to a chikungunya model[END_REF]). Indeed, ID-IO polynomials may contain lower order derivatives and integrating permits to attenuate the structured noise whose amplitude is unknown [START_REF] Fliess | Questioning some paradigms of signal processing via concret examples[END_REF]. In order to obtain ID-IO polynomials, [START_REF] Loeb | More about process identification[END_REF] multiplied them by modulating functions and then integrate the new relation whereas [START_REF] Sira-Ramirez | Algebraic Identification and Estimation Methods[END_REF] used iterative integrals in the frequency domain. In [START_REF] Verdière | A distribution input-output polynomial approach for estimating parameters in nonlinear models, application to a chikungunya model[END_REF], a distribution approach was proposed with a complete study on the numerical error due to this method. In this paper, for the parameter estimation, we focus on the simple use of iterative integrals leading to relations with no derivative in our case. Even in the case of a good parameter estimate, measurement noise is an important factor of uncertainty and error in parameter estimation procedures. To have the probability that the parameter estimate captures the dynamic of the system with respect to the noise on the output(s), an algorithm is proposed in this paper. To apply these results, the Hindmarsh-Rose (HR) model resulting from a simplification and a generalization of the Hodgkin-Huxley model is considered (see [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF], [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF] for more details). The aim of such models is to better understand neuron activity as presented in [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], [START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF], [START_REF] Tokuda | Parameter estimation for neuron models[END_REF]. Its particularity is to present a Hopf bifurcation with respect to its slow-fast parameter. To estimate parameters of the HR neuron model, methods were already proposed. In [START_REF] Tokuda | Parameter estimation for neuron models[END_REF], the authors give two approaches.

The first one concerns the synchronization based parameter estimation and consists in solving a least-square problem subjects to constraints. The second approach is based on the adaptive observers as in [START_REF] Steur | Parameter estimation in Hindmarsh-Rose neurons[END_REF]. This method aims to find a dynamical system so that it synchronizes with the measured voltages from a real neuron. In [START_REF] Schumann-Bishoff | State and parameter estimation using unconstrained optimization[END_REF], authors use nonlinear optimization exploiting the particular structure of the relevant cost function. However, for main of these methods, a first initial guess is needed. That is why, we propose to use the approach based on ID-IO polynomials to obtain this first initial guess.

The paper is organized as follows. In Section II, the parameter estimation procedure based on the ID-IO polynomials is presented as well as the procedure to calculate the probability that a parameter estimate will permit to reproduce the correct behavior of the model output with respect to the noise. Section III concerns the HR model. Its first part presents this slow-fast model and its bifurcation. In the second part, an identifiability study and a parameter estimation are achieved on it. The end of Section III ends with an analysis of the effect of the noise on the estimation of its bifurcation parameter. Finally, Section IV concludes the paper.

II. PARAMETER ESTIMATION METHOD

A. Differential Algebra and identifiability

An identifiability study ensures that two different parameter vectors will lead to two different model outputs. The definition can be formalized as follows:

Definition 1: The model Γ θ is said globally identifiable at θ if a time t 1 ∈ [0, T ] exists such that for all t ∈ [0, t 1 ], for all θ ∈ U, θ = θ then y(t, θ) and y(t, θ) are different. The identifiability study can be analyzed through relations linking input, output and parameters of model Γ θ . They are obtained using the Rosenfeld-Groebner algorithm implemented in the package DifferentialAlgebra of Maple [START_REF] Boulier | Computing representation for radicals of finitely generated differential ideals[END_REF]. From the elimination order consisting in eliminating unobservable variables, the algorithm provides IO polynomials having the following form (See [START_REF] Denis-Vidal | An algorithm to test identifiability of non-linear systems[END_REF]):

P i (y, u, θ) = m i 0 (y, u) + qi j=1 γ i j (θ)m j,i (y, u) = 0 (1)
where (γ i j ) 1≤j≤qi are rational in θ, γ i u = γ i v (u = v), (m j,i ) 1≤j≤qi are differential polynomials with respect to y and u. m i 0 = 0 and i from 1 to m. According to [START_REF] Denis-Vidal | An algorithm to test identifiability of non-linear systems[END_REF], the number of relations is the number of observations. Afterwards, to lighten the notations, only one output is considered so that the index i is omitted. Let us denote s the time derivative of highest order of y in P .

From the coefficients of the IO polynomial P and the initial conditions (if any), we can define the real-valued function φ by φ : θ = (θ 1 , . . . , θ p ) → (γ 1 , . . . , γ q , y(0, θ), . . . . . . , y (s-1) (0, θ)). Under some technical conditions omitted in this paper (for more details, the interested reader can refer to [START_REF] Denis-Vidal | An algorithm to test identifiability of non-linear systems[END_REF]), the identifiability result is the following:

Proposition 1: Assume that the functional determinant P (y, u, θ) = det(m j (y, u), j = 1, . . . , q) ≡ 0 1 and the coefficient of y (s) is not equal to zero at t = 0. The model Γ θ is globally identifiable if and only if the function φ is injective on U p .

B. Estimation procedure

The measurement noise is taken into account by assuming that the output y is disturbed by a centered Gaussian noise η so that the measurement z of y has the following form: z(t) = y(t, θ) + η(t) where θ represents the "true" parameter vector value.

If (t k ) 1≤k≤N is a sequence of discrete times, we denote {z k = z(t k ), k = 1, . . . , N } the set of measure- ments at (t k ) 1≤k≤N and u k = u(t k ) the associated inputs. Let f a real-valued function and I ν (f ), ν ∈ IN, the integrated function such that I ν (f ) = t t-τ τ1 τ1-τ . . . τν τν-1-τ f (τ ν )dτ ν . . . dτ 1 .
Its approximated value by a numerical procedure from the measurements z k , k = 1, . . . , N will be denoted I est ν (f ). Due to the linearity of the integral, a new relation is obtained from P and can be rewritten:

I ν (P ) = I ν ( m0 (y, u)) + q j=1 γj (θ)I ν ( mj (y, u)).
Afterwards, I ν (P ) is called the ID-IO polynomial. Evaluating the expression I ν (P ) at each t k leads to a linear system whose unknown are the coefficients γj (θ), j = 1, . . . , q. Indeed, I ν (P ) is linear with respect to the parameter blocks γ1 (θ), . . . , γq (θ). The following system can also be deduced:

A Θ = b

1 To verify this assumption, it is sufficient to verify the linear independence of the m j (y, u), j = 1, . . . , q by checking the functional determinant given by the Wronskian [START_REF] Denis-Vidal | Some effective approaches to check identifiability of uncontrolled nonlinear systems[END_REF]

     m 1 (y, u)
. . . mq(y, u) m 1 (y, u) (1) . . . mq(y, u) (1) . . .

m 1 (y, u) (q-1) . . . mq(y, u) (q-1)      (2)
is not identically equal to zero. In fact, the monomials are linearly independent if there exists a time point at which the Wronskian is non-zero.

where Θ = (γ 1 (θ), . . . , γq (θ)) T . The kth line of A and b are respectively given by (A) k = (I est ν ( mj (y, u))) j=1,...,q and b k = -I est ν ( m0 (y, u)). Remark 1: The linear system that can be deduced directly from P may contain derivatives of high order. Indeed, the elimination algorithm of Rosenfeld-Groebner is based on addition, differentiation in time and multiplication by any polynomials in x, u, y and θ. Or, integrating permits to decrease the order of some derivatives but also to attenuate the structured noise whose amplitude is unknown. By definition, a structured noise is solution of a given homogeneous linear differential equation (See [START_REF] Fliess | Questioning some paradigms of signal processing via concret examples[END_REF] for more detailed). Finally, whatever the integral method applied on P , the new relation has always given better results in the case of the parameter estimation (See [START_REF] Loeb | More about process identification[END_REF], [START_REF] Sira-Ramirez | Algebraic Identification and Estimation Methods[END_REF], [START_REF] Verdière | A distribution input-output polynomial approach for estimating parameters in nonlinear models, application to a chikungunya model[END_REF] for some examples). Remark that the estimated parameters obtained by this way may be used as an initial guess for a local algorithm to improve the results. Even in the case of a good relative error between the exact and the estimated parameters, the behavior of the output and the simulated output may not be the same. For example, if the model has a Hopf bifurcation, near the critical point, the model output can either tend to a stable equilibrium point or be periodic. Thus, the estimated parameter vector can lead to a deterioration of the dynamic of the model output during the time as it will be seen at section III-C. Remark too that we can not guarantee the result obtained by a local algorithm with a first initial guess for which the system does not present the good dynamic (see Remark 2). In order to test the reliability of any parameter estimation procedure, we propose, in the following section, an algorithm to calculate the probability to have an estimated parameter providing the correct behavior of the model output.

C. Probability of obtaining the correct behavior

Let θ the parameter vector of Model Γ θ to estimate. Suppose that its kth component, θ k , is a bifurcation parameter. With no restriction, we can suppose that k = 1. Let θ c 1 its bifurcation value that is the value for which a little change leads to a complete different behavior of Model Γ θ . For a given error level, to determine the probability for the estimated value of θ 1 not being in the accepted behavioral interval, 1000 simulated data are generated from the output model. For each simulated data, a parameter estimation procedure is used so that the parameter vector θ containing the estimated value θ1 of θ 1 is calculated. The comparison between θ c 1 and θ1 permits to conclude if θ1 is in the accepted behavioral interval. The procedure is summed up below.

1) Use the nominal parameter values to numerically solve the ODE model to get the solution of the output or measurement variables at the experimental design time points. 2) Generate M sets (e.g., 1000) of simulated data from the output or measurement model with a given measurement error level. 3) For each of the M simulated data sets, estimate parameter vector θ using the method developed at section II for example. 4) Count the number of estimates θ that are in the accepted behavioral interval.

III. HINDMARSH-ROSE MODEL

The model of Hindmarsh-Rose (HR) results from a simplification and a generalization of the Hodgkin-Huxley model [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF], [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]. This model is a slow-fast model and is able to reproduce the rich dynamics of a neuron, such as spiking, bursting and chaotic behaviors. The aim of such models is to better understand neuron activity [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], [START_REF] Izhikevich | Dynamical systems in neuroscience[END_REF], [START_REF] Tokuda | Parameter estimation for neuron models[END_REF]. The HR model [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF] reads as follows (where θ = (a, b, d, ε) T ):

       ẋ1 (t, θ) = x 2 (t, θ) + ax 1 (t, θ) 2 -x 1 (t, θ) 3 -x 3 (t, θ) + I(t), ẋ2 (t, θ) = 1 -dx 1 (t, θ) 2 -x 2 (t, θ), ẋ3 (t, θ) = ε(b(x 1 (t, θ) -c x1 ) -x 3 (t, θ)).
(3)

x 1 represents the membrane potential and is the output y := x 1 of the HR model, x 2 the recovery variable associated with the fast current due to the passage of the Na + or K + ions and x 3 the adaptation current associated with the slow current due to the passage of the Ca + ions. Parameters a, b, d and are experimentally determined from measurements of membrane potentials, while c x1 is the x 1 -coordinate of the leftmost equilibrium of the two-dimensional system given by the first two equations of (3) when I(t) = 0 and x 3 (t, θ) = 0. I corresponds to the applied current (in Ampere) and afterwards it is supposed constant. Finally, parameter ε represents the ratio of time scales between fast and slow fluxes across the membrane of a neuron.

The authors in [START_REF] Corson | Hopf bifurcations in a chain of coupled hindmarsh-rose system[END_REF] show that the parameter ε presents a Hopf bifurcation. A Hopf bifurcation refers to the local birth or death of a periodic solution near an equilibrium point due to a change of one parameter value. Therefore, the HR model can present different dynamics with respect to the value of ε. Let ε c its bifurcation value. In [START_REF] Corson | Hopf bifurcations in a chain of coupled hindmarsh-rose system[END_REF], authors give an explicit formula of the bifurcation parameter ε c in function of a, b and d. For example, for a = 3 , b = 4 , d = 5 , I(t) = 3.25A, System (3) admits one equilibrium point and the constant bifurcation parameter ε c is approximatively equal to 0.125912. Fig. 1 and 2 represent the (x 1 , x 2 , x 3 ) view of the phase portrait for two different values of ε, before and after the Hopf bifurcation. For ε greater than ε c the solution is periodic and disappears when ε increases such that the system has a stable equilibrium point. The time series of system (3) for the two previous parameters of ε are represented at Fig. 3 and4 respectively.

A. Identifiability study of the HR model

To obtain the IO polynomial of the HR model, the package DifferentialAlgebra of Maple is used. I being an input of the system, the elimination order

[θ] ≺ [y, I] ≺ [x 1 , x 2 , x 3 ] is chosen. It consists in eliminating x 1 ,
x 2 , x 3 in order to obtain a polynomial depending only on y, θ and I. The Rosenfeld-Groebner algorithm provides a polynomial with a derivative of order 3 and with 24 expressions. To obtain a simpler polynomial, let u 1 (t) = e -t , v 1 (t) = e -t t 0

x 2 1 (s, θ)e s ds. In integrating the second Equation of (3), one gets x 2 (t, θ) = (x 2 (0, θ) -1)u 1 (t) + 1 -dv 1 (t).

Consider now the following system:

               ẋ1 (t, θ) = x 2 (t, θ) + ax 2 1 (t, θ) -x 3 1 (t, θ) -x 3 (t, θ) + I(t), x 2 (t, θ) = (x 2 (0, θ) -1)u 1 (t) + 1 -dv 1 (t), ẋ3 (t, θ) = ε(b(x 1 (t, θ) -c x1 ) -x 3 (t, θ)), u1 (t) = -u 1 (t), v1 (t) = -v 1 (t) + x 2 1 (t). (4)
System (3) completed with the initial conditions (x 1 (0, θ), x 2 (0, θ), x 3 (0, θ)) is equivalent to System (4) completed with (x 1 (0, θ), x 2 (0, θ), x 3 (0, θ), 1, 0). Considering now the elimination order

[θ] < [y, u 1 , v 1 , I] ≺ [x 1 , x 2 , x 3 ],
we obtain (the time variable t is omitted):

P := ÿ + 3 y 2 ẏ + (-ε y 0 + ε + y 0 -1) u 1 +(d ε -d) v 1 + ε (x 3 + ẏ) + (-a ε + d) y 2 -2 a y ẏ + b ε y -b c y ε -ε(1 + I) = 0. (5)
The function φ is defined by: φ : θ → (-ε

y 0 + ε + y 0 - 1, d ε -d, ε, -a ε + d, -2 a, b ε, -b c x ε -ε(1 + I)) is clearly injective so the model (3) with y = x 1 is identifiable.
Remark that y 0 can be estimated too using the polynomial P .

B. ID-IO polynomial

In order to decrease the derivative order of polynomial P defined by ( 5), P is integrated twice given the following relation (the variable θ is omitted):

I2(P ) := t t-τ τ 1 τ 1 -τ ÿ(τ2)dτ2dτ1 +3 t t-τ τ 1 τ 1 -τ y(τ2) 2 ẏ(τ2)dτ2dτ1 +(-ε y0 + ε + y0 -1) t t-τ τ 1 τ 1 -τ u1(τ2)dτ2dτ1 +(d ε -d) t t-τ τ 1 τ 1 -τ v1(τ2)dτ2dτ1 -ε t t-τ τ 1 τ 1 -τ (-y(τ2) 3 -ẏ(τ2))dτ2dτ1 +(-a ε + d) t t-τ τ 1 τ 1 -τ y(τ2) 2 dτ2dτ1 -2 a t t-τ τ 1 τ 1 -τ y(τ2) ẏ(τ2)dτ2dτ1 +b ε t t-τ τ 1 τ 1 -τ y(τ2)dτ2dτ1 -(I ε + b cy ε + ε) t t-τ τ 1 τ 1 -τ 1dτ2dτ1. (6)
Some double integrals can be simplified:

t t-τ τ 1 τ 1 -τ ẏ(τ2)dτ2dτ1 = t t-τ (y(τ1) -y(τ1 -τ ))dτ1, t t-τ τ 1 τ 1 -τ ÿ(τ2)dτ2dτ1 = y(t) -2 y(t -τ ) + y(t -2τ ), t t-τ τ 1 τ 1 -τ y(τ2) 2 ẏ(τ2)dτ2dτ1 = t t-τ (y(τ1) 3 -y(τ1 -τ ) 3 )dτ1, 2 t t-τ τ 1 τ 1 -τ ẏ(τ2)y(τ2)dτ2dτ1 = t t-τ (y(τ1) 2 -y(τ1 -τ ) 2 )dτ1, t t-τ τ 1 τ 1 -τ dτ2dτ1 = τ 2 .
In substituting these relations in the ID-IO polynomial (6), the following relation which does not contain no derivative is obtained:

I2(P ) = y(t) -2 y(t -τ ) + y(t -2τ ) + t t-τ (y(τ1) 3 -y(τ1 -τ ) 3 )dτ1 +(-ε y0 + ε + y0 -1) t t-τ τ 1 τ 1 -τ u1(τ2)dτ2dτ1 +(d ε -d) t t-τ τ 1 τ 1 -τ v1(τ2)dτ2dτ1 -ε - t t-τ τ 1 τ 1 -τ y(τ2) 3 dτ2dτ1 + t t-τ (y(τ1) -y(τ1 -τ ))dτ1 +(-a ε + d) t t-τ τ 1 τ 1 -τ y(τ2) 2 dτ2dτ1 -a t t-τ (y(τ1) 2 -y(τ1 -τ ) 2 )dτ1 +b ε t t-τ τ 1 τ 1 -τ y(τ2)dτ2dτ1 -(I ε + b cy ε -d + ε)τ 2 .
Evaluating this ID-IO polynomial at each t k , k = 1, . . . , N provides a linear system of the form A Θ = b where Θ = (-(-ε

y 0 + ε + y 0 -1), -(d ε -d), ε, -(-a ε + d), a, -b ε, I ε + b c y ε -d + ε) and b k = z k -2 z k-τ + z k-2τ + I est 1 (y(.) 3 -y(. -τ ) 3 ), A k = I est 2 (u 1 ), I est 2 (v 1 ), -I est 2 (y 3 ) +I est 1 (y(.) -y(. -τ )), I est 2 (y 2 ), I est 1 (y(.) 2 -y(. -τ ) 2 ), I est 2 (y), τ 2 ,
A k and b k represent the k th line of A and b respectively (k from 1 to 7). System A Θ = b is solved using the QR factorization that does not necessitate an initial guess of the unknown parameters.

C. Parameter estimation

The trapezoidal formulas are used to estimate the previous integrals. Let (t k ) k=1,...,N the discretization points of the time interval [0, T ], h = T /(N -1). Then, for a simple integral, one gets:

T 0 f (s)ds = h 1 2 f (t 1 ) + 1 2 f (t N ) + N -1 i=2 f (t i )
and for a double integral f (s)ds [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF] For the simulations, the following values are taken: a = 3, b = 4, d = 5, I = 3.25A, ε = 0.12, (x 1 (0, θ), x 2 (0, θ), x 3 (0, θ)) = (0.2; 0.7; 4). The time interval is [0, 100] with a step size h = 0.01s. The integrals are evaluated in using 29 points. This number gave the best results after several numerical tests. The outputs are perturbed by an additive white centered Gaussian noise with a standard deviation σ given by the three successive values: 0.0001, 0.0005, 0.001. For each of this value, the estimate of â, b, d, ε are given in Table 1. Table 2 provides first, the relative errors between the exact and the estimated parameters, that is between θ = (ε, a, b, c) and θ = (ε, â, b, ĉ) then, the relative errors between the measured output and the simulated output obtained with θ. The standard deviation of the Gaussian noise is given by a factor of the model output. Model outputs obtained with the different estimated parameter values given in Table 1 are represented at Figure 5.

For the three deviations, the relative error on the estimated parameters is good contrary to the relative error on the outputs. For the two standard deviations 0.0005 and 0.001, the behavior of the system with the approximated value ε does not correspond to the one with the initial value ε (for the first one, a longer time interval is required). Indeed, while initial system presents a periodic behavior, the simulated outputs obtained with ε tend to a stable equilibrium point.

Remark 2: The use of the local Levenberg-Marquard algorithm implemented in scilab did not permit to obtain a value of ε in the accepted behavioral interval in the case of σ = 0.001. Indeed, we obtain the approximated value ε = 0.128 and ε > ε c > ε.

In order to calculate the probability that the approximated value gives the correct behavior of the system, the algorithm given at Section II-C is applied. The results are gathered in Figure 6. For example, for ε = 0.10, 100% of estimates are in the accepted behavioral interval when σ ≤ 0.0002. If Fig. 5. Outputs of the model obtained with the parameter values given at Table 1. The black solid line represents the measured output. The light blue dash line, dark blue line and red dash line represent respectively the outputs obtained with σ = 0.0001, σ = 0.005 and σ = 0.001. Fig. 6. Percentage of estimated values ε out of the accepted behavioral interval. In blue points, 0% of estimated values are out of the accepted behavioral interval. In blue plus, between 0 (excluded) and 5 percents of estimated values are not in the accepted behavioral interval. In blue cross, between 5 and 10 percents of estimated values are not in the accepted behavioral interval. In red asterisk, between 10 and 30 percents of estimated values are not in the accepted behavioral interval. In red circle, between 30 and 50 percents of estimated values are not in the accepted behavioral interval. Finally, in red square, more than 50% are not in the accepted behavioral interval.

σ ∈ [0.0003; 0.0004[, less than 5% of estimated parameters are not accepted behavioral values. Then, the percentage of accepted behavioral values decreases until σ = 0.001 for which more that 50% of estimated parameters are not accepted behavioral values. For a value ε equal to 0.12, as soon as there is some noise, more than 30% of estimated values are not in the accepted behavioral interval.

IV. CONCLUSION

The numerical method proposes in this paper permits to obtain a very good first initial guess without any knowledge of the parameter values. However, even if the relative error on the parameters is small, the relative error between the measured outputs and the simulated ones can be important due to a change of the output behavior as seen as for the Hindmarsh Rose model. Indeed, it appears that the effect of the standard deviation influences the expected behavior. More generally, it appears that the consideration of the noise is essential for some particular parametric models since a little error on the parameter values can lead to a completely different behavior of the model. Future works concern the study of i) a similar analysis using other forms of noise, ii) the theoretical evaluations of upper bounds of bias and variance of parameter estimates with respect to intensity of the noise.

Fig. 1 .

 1 Fig. 1. (x 1 , x 2 , x 3 ) view of the phase portrait of system (3) when ε = 0.12 < εc.

Fig. 2 .

 2 Fig. 2. (x 1 , x 2 , x 3 ) view of the phase portrait of system (3) when ε = 0.13 > εc.

Fig. 3 .

 3 Fig.3. Time series of system (3) when ε = 0.12 < εc.

Fig. 4 .

 4 Fig. 4. Time series of system (3) when ε = 0.13 > εc.

Table 1 :

 1 Parameter values obtained with the method presented at Section II.

	σ	ε	â	b	d
	0.0001	0.1204	3.0013 3.9712 4.9757
	0.0005	0.1390	2.9371 3.9790 5.0537
	0.001	0.15437 2.8076 4.4480 5.1488
	σ	Relative errors	Relative errors
		on the parameters on the outputs
	0.0001	0.005	0.0161
	0.0005	0.012	0.1949
	0.001	0.072	0.4923

Table 2 :

 2 Relative errors between θ and θ (that is θ -θ 2 / θ 2) obtained with the method presented at Section II.


Guillaume Chauvet

Properties of Chromy's sampling procedure

proposed a unequal probability sampling algorithm, which enables to select a sample in one pass of the sampling frame only. This is the default sequential method used in the SURVEYSELECT procedure of the SAS software. In this article, we study the properties of Chromy sampling. We prove that the Horvitz-Thompson is asymptotically normally distributed, and give an explicit expression for the second-order inclusion probabilities. This makes it possible to estimate the variance unbiasedly for the randomized version of the method programmed in the SURVEYSELECT procedure.

1 Introduction [START_REF] Chromy | Sequential sample selection methods[END_REF] proposed a fixed-size unequal probability sampling design which is strictly sequential, in the sense that a sample is selected in one pass of the sampling frame only. This algorithm benefits from a stratification effect, in the sense that the selected units are well spread over the population like with systematic sampling. The drawback is that many second-order inclusion probabilities are zero, making unbiased variance estimation not possible. [START_REF] Chromy | Sequential sample selection methods[END_REF] therefore proposed to partially randomize the order of the units in the population before applying the sampling algorithm. This randomization is sufficient to guarantee that the second-order inclusion probabilities are positive. The randomized Chromy algorithm is the default sequential method currently available in the SURVEYSELECT procedure of the SAS software. The method has been extensively used for sample surveys, see for example [START_REF] Mills | Disparities in retail marketing for menthol cigarettes in the united states, 2015[END_REF], [START_REF] Radwin | national postsecondary student aid study (npsas: 16): Student financial aid estimates for 2015-16[END_REF], [START_REF] Schmitt | Peer reviewed: Relationships between theoretically derived short-term outcomes and support for policy among the public and decision-makers[END_REF] and [START_REF] Rust | Tobacco marketing at snap-and wic-authorized retail food stores in the united states[END_REF] for recent examples.

So far, the properties of this sampling algorithm have not been fully investigated, and this is the purpose of the current paper. We prove that Chromy sampling is equivalent to ordered pivotal sampling [START_REF] Deville | Unequal probability sampling without replacement through a splitting method[END_REF][START_REF] Chauvet | On a characterization of ordered pivotal sampling[END_REF], in the sense that both algorithms lead to the same sampling design. This leads to the Horvitz-Thompson being consistent and asymptotically normally distributed, under weak assumptions. This also leads to an explicit formula for the second-order inclusion probabilities, making unbiased variance estimation possible for randomized Chromy sampling.

The paper is organized as follows. In Section 2, the notation and the assumptions are given. Chromy sampling and ordered pivotal sampling are introduced in Sections 3 and 4. The equivalence between both sampling designs is proved in Section 5, and the properties of Chromy sampling are studied. The results of a small simulation study are given in Section 6. We conclude in Section 7. The proofs are gathered in the Appendix.

Notation and assumptions

We consider a finite population U of N sampling units that may be represented by integers k = 1, . . . , N. Denote by π = (π 1 , . . . , π N ) ⊤ a vector of probabilities, with 0 < π k < 1 for any unit k ∈ U, and with n = k∈U π k the expected sample size. A random sample S is selected in U by means of a sampling design p(•) with parameter π, in the sense that the expected number of draws for unit k is π k . We let I k denote the number of times that unit k is selected in the sample, and we note I = (I 1 , . . . , I N ) ⊤ .

The set of probabilities π may be defined proportionally on some positive auxiliary variable known for any unit in the population, which leads to unequal probability sampling with probabilities proportional to size (π-ps) (Särndal et al., 1992, Section 3.6.2). The sampling algorithm proposed by [START_REF] Chromy | Sequential sample selection methods[END_REF] may handle π k 's greater than 1, in which case a same unit may be selected several times in the sample. In this paper, we focus on the fairly usual situation when all π k 's lie between 0 and 1, which means that the sampling design is without replacement. In this case, we may interpret π k as the probability for unit k to be included in the sample, and I k is the sample membership indicator for unit k. We are interested in the total t y = k∈U y k of some variable of interest taking the value y k for unit k ∈ U. The Horvitz-Thompson estimator is

tyπ = k∈U y k π k I k . (2.1)
In order to study Chromy sampling, some additional notation is needed. We let V k = k l=1 π l denote the cumulated inclusion probabilities up to unit k, with V 0 = 0. The integer part of V k is the largest integer smaller than V k , and is denoted as V I k . The difference between V k and its integer part V I k is the fractional part, and is denoted as

V F k = V k -V I k . For example, if V k = 3.6 we have V I k = 3 and V F k = 0.6, and if V k = 4.0 we have V I k = 4 and V F k = 0. A unit k is a cross-border if V k-1 ≤ i and V k > i for some integer i = 1, . . . , n -1.
The cross-border units are denoted as k i , i = 1, . . . , n -1, and we note

a i = i -V k i -1 and b i = V k i -i.
The cross-border units define a partition of the population into microstrata U i , i = 1, . . . , n, which are defined as

U i = {k ∈ U; k i-1 ≤ k ≤ k i } with k 0 = 0 and k n = N + 1. (2.2)
The quantities are presented in Figure 1 for illustration.

We consider the following assumptions, which are the same than in Chauvet and Le Gleut (2019):

H1: There exists some constants 0 < f 0 and f 1 < 1 such that for any k ∈ U:

f 0 n N ≤ π k ≤ f 1 . (2.

3)

H2: There exists some constant C 1 such that:

k∈U π k y k π k - t y n 4 ≤ C 1 N 4 n -3 . (2.4) i -1 i V k i-1 -1 V k i-1 V k i-1 +1 . . . V k i -2 V k i -1 V k i V k i +1 a i-1 b i-1 a i b i π k i-1 π k i-1 +1 . . . π k i -1 π k i π k i +1 U i-1 U i U i+1
Figure 1: Inclusion probabilities and cross-border units in microstratum U i , for population U H3: There exists some constant C 2 > 0 such that:

n i=1 k∈U i α ik y k π k - l∈U i α il y l π l 2 ≥ C 2 N 2 n -1 , (2.5)
where for any unit k ∈ U i we take:

α ik =    b i-1 if k = k i-1 , π k if k i-1 < k < k i , a i if k = k i ,
with the convention that b 0 = a n = 0.

It is assumed in (H1) that the first-order inclusion probabilities are bounded away from 1. This is not a severe restriction in practice, since a unit with an inclusion probability close to 1 is usually placed into a take-all stratum, i.e. the probability is rounded to 1 and the unit is not involved in the selection process. It is also assumed in (H1) that the first-order inclusion probabilities have a lower bound of order n/N, which ensures that no design-weight d k = 1/π k is disproportionately larger than the others.

Under the condition (H1), the condition (H2) holds in particular if the variable y has a finite moment of order 4. This seems a fair assumption in practice, unless the variable of interest is heavily skewed like in wealth sur-veys, for example.

Assumption (H3) is somewhat technical, and is used to ensure that the Horvitz-Thompson estimator is not close to being degenerate. This assumption is only needed to prove a central-limit theorem for the Horvitz-Thompson estimator. It requires that the dispersion within the microstrata does not vanish. For example, it does not hold if y k is proportional to π k .

3 Chromy sampling [START_REF] Chromy | Sequential sample selection methods[END_REF] proposed a sampling algorithm which is strictly sequential, in the sense that the units in the population are successively considered for possible selection, and the decision for the unit is made at once. The method is presented in Algorithm 1. Let us denote

S c ∼ Chr(π; U) (3.1)
for a sample selected by means of Chromy sampling with parameter π in the population U. The method was originally proposed for π-ps sampling where the inclusion probabilities are defined proportionally on some auxiliary variable, but it is applicable to any set of inclusion probabilities.

The method proceeds by considering at each step k = 1, . . . , N the unit k for possible selection, and by computing its probability of selection conditionally on the number of units already selected. This algorithm defines a fixed-size sampling design, and [START_REF] Chromy | Sequential sample selection methods[END_REF] proves that the parameter π defining the inclusion probabilities is exactly matched. Note that the case

V F k > V F k-1
(Step 2.a) corresponds to the treatment of a non cross-border unit, while the case

V F k-1 ≥ V F k (
Step 2.b) corresponds to the treatment of a cross-border unit. For illustration, the complete probability tree for Chromy sampling on a small population is given in Appendix A.

This algorithm allocates the sample regularly in the population, as stated in Proposition 1. The proof is given in Appendix B. At any step k of the procedure, the number of units selected is equal to the sum of inclusion probabilities up to rounding, a property which is sometimes coined as spatial balancing [START_REF] Grafström | Spatially balanced sampling through the pivotal method[END_REF].

Algorithm 1 Chromy sampling with parameter π in the population U 1. At step k = 1, take I 1 = 1 with probability π 1 .

2. At step k = 2, . . . , N:

(a) If V F k > V F k-1
, then we take the transition probabilities

P r I k = 1 k-1 l=1 I l = V I k-1 = V F k -V F k-1 1 -V F k-1 , P r I k = 1 k-1 l=1 I l = V I k-1 + 1 = 0. (b) If V F k-1 ≥ V F
k , then we take the transition probabilities

P r I k = 1 k-1 l=1 I l = V I k-1 = 1, P r I k = 1 k-1 l=1 I l = V I k-1 + 1 = V F k V F k-1
.

Proposition 1. For any k = 1, . . . , N, we have

V I k ≤ k l=1 I l ≤ V I k + 1. (3.2)
A drawback of the method is that, by construction, two non cross-border units inside the same microstratum U i may not be selected jointly in the sample. Therefore, many second-order inclusion probabilities are equal to 0 and the variance of the Horvitz-Thompson estimator may not be unbiasedly estimated.

For this reason, [START_REF] Chromy | Sequential sample selection methods[END_REF] proposed to use a randomized procedure, which is as follows. The population U is viewed as a closed loop, and we consider the set Σ c of the N possible circular permutations, each of which using a different unit as the first one. The first permutation is σ 1 , with the natural order 1, . . . , N. For k = 2, . . . , N the k-th permutation is σ k where the units are in the order k, . . . , N, 1, . . . , k -1. The randomized Chromy sample S rc is selected as

S rc ∼ Chr(π σ k ; U σ k ), (3.3)
with σ k a random permutation selected in Σ c with probability π k /n, with U σ k the population ordered with respect to σ k , and π σ k the vector of probabilities ordered accordingly. This is the algorithm currently implemented in the SURVEYSELECT procedure of the SAS software.

Ordered Pivotal sampling

Ordered pivotal sampling [START_REF] Fuller | Sampling with random stratum boundaries[END_REF][START_REF] Deville | Unequal probability sampling without replacement through a splitting method[END_REF][START_REF] Chauvet | On a characterization of ordered pivotal sampling[END_REF] is presented in Algorithm 2. This is a succession of duels between units, and at each step the two first units remaining in the population are considered. If the sum of their probabilities is lower than 1 (rejection step), one of the unit is randomly discarded while the other gets the sum of their probabilities. If the sum of their probabilities is greater than 1 (selection step), one of the unit is randomly selected while the other goes on with the residual probability. For illustration, the complete probability tree for ordered pivotal sampling on a small population is given in Appendix A.

The pivotal sample is selected in at most N -1 steps. Pivotal sampling is a

Algorithm 2 Ordered pivotal sampling with parameter π in the population U 1. Initialize with π(0) = π N .

2. At step t = 1, . . . , T :

(a) Initialize with π(t) = π(t -1).

(b) Take k < l the two first units in the population such that

π k (t -1) / ∈ {0, 1} and π l (t -1) / ∈ {0, 1} (c) If π k (t -1) + π l (t -1) ≤ 1 (rejection step), then do: {π k (t), π l (t)} = {π k (t -1) + π l (t -1), 0} with prob. p(t) {0, π k (t -1) + π l (t -1)} with prob. 1 -p(t),
where

p(t) = π k (t -1) π k (t -1) + π l (t -1) . (d) If π k (t -1) + π l (t -1) > 1 (selection step), then do: {π k (t), π l (t)} = {1, π k (t -1) + π l (t -1) -1} with prob. p(t) {π k (t -1) + π l (t -1) -1, 1} with prob. 1 -p(t)
where

p(t) = 1 -π l (t -1) 2 -π k (t -1) -π l (t -1)
.

3. The algorithm stops at step T when all the components of π(T ) are 0 or 1. Take I = π(T ).

particular case of the cube method [START_REF] Deville | Efficient balanced sampling: the cube method[END_REF], which enables to perform balanced sampling, i.e. to select samples such that the Horvitz-Thompson estimator exactly matches the known totals for some auxiliary variables. Pivotal sampling has found uses in spatial sampling, since it enables to spread well the sample over space: see for example [START_REF] Grafström | Spatially balanced sampling through the pivotal method[END_REF] for the so-called local pivotal method, [START_REF] Chauvet | Asymptotic properties of pivotal sampling with application to spatial sampling[END_REF] for the so-called pivotal tesselation method, or [START_REF] Benedetti | Spatially balanced sampling: a review and a reappraisal[END_REF] for a recent review on spatial sampling methods. Pivotal sampling is also of use in Monte Carlo methods [START_REF] Gerber | Negative association, ordering and convergence of resampling methods[END_REF].

Properties of Chromy sampling

We first prove in Theorem 1 that Chromy sampling and ordered pivotal sampling are equivalent, which is the main result of the paper. The proof is lengthy, and given in Appendix C.

Theorem 1. Ordered pivotal sampling and Chromy sampling with the same parameter π induce the same sampling design.

By using the characterization of Chromy sampling given in Theorem 1, the mean-square consistency of the Horvitz-Thompson estimator stated in equation (5.1) of Theorem 2 is a direct consequence of Theorem 2 in Chauvet (2017). The central-limit theorem stated in equation (5.2) is a direct consequence of Theorem 1 in [START_REF] Chauvet | Asymptotic properties of pivotal sampling with application to spatial sampling[END_REF].

Theorem 2. Suppose that the sample S c is selected by means of Chr(π; U).

If assumption (H2) holds, then

E N -1 ( tyπ -t y ) = O(n -1 ). (5.1)
If in addition assumptions (H1) and (H3) hold, then

tyπ -t y V ( tyπ ) -→ L N (0, 1), (5.2)
where -→ L stands for the convergence in distribution.

It also follows from Theorem 1 that Chromy sampling is a negatively associated sampling design [START_REF] Joag-Dev | Negative association of random variables with applications[END_REF]. This implies that the Sen-Yates-Grundy conditions are satisfied. From Theorem 2 in [START_REF] Bertail | Bernstein-type exponential inequalities in survey sampling: Conditional poisson sampling schemes[END_REF], the Horvitz-Thompson also satisfies a Bennett/Bernsteintype exponential inequality.

From Theorem 5.1 in [START_REF] Chauvet | On a characterization of ordered pivotal sampling[END_REF], and from the computation given in [START_REF] Deville | Une nouvelle (encore une!) méthode de tirage à probabilités inégales[END_REF], it is possible to give an explicit expression for the secondorder inclusion probabilities under Chromy sampling. This is the purpose of Theorem 3.

Theorem 3. Let k and l be two distinct units in U. If k and l are two non cross-border units that belong to the same microstratum U i , then

π kl = 0,
if k and l are two non cross-border units that belong to distinct microstrata U i and U j , respectively, where i < j, then

π kl = π k π l {1 -c(i, j)} , if k = k i-1
and l is a non cross-border unit that belongs to the microstratum U j where i ≤ j, then

π kl = π k π l 1 -b i-1 (1 -π k ) {π k (1 -b i-1 )} -1 c(i, j) ,
if l = k j-1 and k is a non cross-border unit that belongs to the microstratum U i where i < j, then

π kl = π k π l 1 -(1 -π l )(1 -b j-1 )(π l b j-1 ) -1 c(i, j) ,
if k = p i-1 and l = p j-1 , where i < j, then

π kl = π k π l 1 -b i-1 (1 -b j-1 )(1 -π k )(1 -π l ) {π k π l b j-1 (1 -b i-1 )} -1 c(i, j) , where c(i, j) = j-1 l=i c l , c l = a l b l {(1 -a l )(1 -b l )} -1 and with c(i, i) = 1.
It is clear from Theorem 3 that many second-order inclusion probabilities are equal to zero for Chromy sampling. It also makes possible to compute the second-order inclusion probabilities for randomized Chromy sampling. For any units k = l ∈ U, let us denote π rc kl their second-order inclusion probability under randomized Chromy sampling. Then:

π rc kl = i∈U π i n π σ i kl , (5.3)
with π σ i kl the joint selection probabilities of units k and l with the permutation σ i , i.e. when Chromy sampling is applied to the population U σ i with parameter π σ i . A SAS IML subroutine to compute the second-order inclusion probabilities in (5.3) is available as Supplementary Material.

We evaluate in the simulation study performed in Section 6 a variance estimator making use of second-order inclusion probabilities computed from equation (5.3). For illustration, we give in this Section a small example. We consider a population U of size N = 8, with the parameter π = (0.2, 0.4, 0.7, 0.4, 0.6, 0.6, 0.3, 0.8) ⊤ , which leads to a sample of size n = 4. From equation (5.3), we obtain the following matrix of second-order inclusion probabilities (rounded to three decimal places):

(π rc kl ) =             0.
            .
We also selected 10 6 samples by means of the SURVEYSELECT procedure with the option METHOD=PPS_SEQ, which leads to randomized Chromy sampling. These 10 6 samples are used to obtain a simulation-based approximation of the matrix of second-order inclusion probabilities, which is given below (rounded to three decimal places):

(π rc kl,sim ) =             0.
           
.

It is clear that both matrices are almost identical.

Simulation study

We conducted a simulation study in order to evaluate variance estimation and interval estimation for randomized Chromy sampling. The set-up is inspired from [START_REF] Chauvet | Examining some aspects of balanced sampling in surveys[END_REF]. We generate 2 populations of size N = 500, each consisting of an auxiliary variable x and 4 variables of interest y 1 , . . . , y 4 . In the first population, the x-values are generated according to a Gamma distribution with shape and scale parameters 2 and 2; in the second population, the x-values are generated from a log-normal distribution with parameters 0 and 1.7. The x-values are then shaped and scaled to lie between 1 and 10.

Given the x-values, the values of the variables of interest are generated according to the following models:

linear : y 1k = α 10 + α 11 (x k -µ x ) + σ 1 ǫ k , quadratic : y 2k = α 20 + α 21 (x k -µ x ) 2 + σ 2 ǫ k , (6.1) exponential : y 3k = exp{α 30 + α 31 (x k -µ x )} + σ 3 ǫ k , bump : y 4k = α 40 + α 41 (x k -µ x ) 2 -α 42 exp -α 43 (x k -µ x ) 2 + σ 4 ǫ k ,
where µ x is the population mean of x, and where ǫ k follows a standard normal distribution N (0, 1). The population mean µ y and the population dispersion S 2 y for the two populations and the four variables of interest are given in Table 1.

In each population, we computed inclusion probabilities proportional to the x-values, according to the formula

π k = n x k l∈U x l , (6.2)
with n = 50, 100 or 200. The range of inclusion probabilities is given in Table 2. In some cases, equation (6.2) leads to inclusion probabilities greater than 1 for some units. In such case, the corresponding units are selected with certainty (π k = 1), and the other probabilities are recomputed. For the first population, 5 units are selected with certainty with n = 200. For the second population, 2 units are selected with certainty with n = 100 and 9 units are selected with certainty with n = 100. We select B = 1, 000 samples by means of randomized Chromy sampling, using the SURVEYSELECT procedure. For each sample and each variable of interest, we compute the Horvitz-Thompson estimator tyπ , and the Sen-Yates-Grundy variance estimator

V ( tyπ ) = 1 2 k =l∈Src π k π l -π rc kl π rc kl y k π k - y l π l 2 , (6.3)
where the second-order inclusion probabilities are given by equation ( 5.3). To evaluate the properties of this variance estimator, we compute the relative bias

RB{ V ( tyπ )} = 100 × B -1 B b=1 Vb ( tyπb ) -V ( tyπ ) V ( tyπ ) , (6.4)
where Vb ( tyπb ) denotes the variance estimator in the b-th sample, and where V ( tyπ ) is the exact variance, computed by using the exact second-order inclusion probabilities given in (5.3). As a measure of stability, we use the Relative Root Mean Square Error

RRMSE{ V ( tyπ )} = 100 × B -1 B b=1 Vb ( tyπb ) -V ( tyπ ) 2 1/2 V ( tyπ ) .
Finally, we compute the error rate of the normality-based confidence intervals with nominal one-tailed error rate of 2.5 % in each tail.

The simulation results are given in Table 3. The Sen-Yates-Grundy variance estimator is almost unbiased in all cases considered, except for the first population with n = 50 where the variance estimator is slightly positively biased. The Relative Root Mean Square Error diminishes as n increases, as expected.

We note that the coverage rates are not well respected with n = 50, which is likely due to the small sample size and to the instability of the variance estimator. When the sample size increases, the coverage rates become close to the nominal level.

Conclusion

In this paper, we have studied Chromy's sampling algorithm. We proved that it is equivalent to ordered pivotal sampling, which enables in particular computing the second-order inclusion probabilities for the randomized Chromy algorithm programmed in the SURVEYSELECT procedure. The results in our simulation study confirm that the variance estimator based on the secondorder probabilities computed from Deville's formulas show almost no bias for moderate sample sizes. The number of computations for the second-order inclusion probabilities of randomized Chromy sampling is of order N 3 . Formula (5.3) is therefore tractable in case of a small population, for example when Chromy sampling is used to select a set of Primary Sampling Units in a multistage survey (e.g. [START_REF] Rust | Tobacco marketing at snap-and wic-authorized retail food stores in the united states[END_REF]. Otherwise, we may resort to a simulation-based approximation of equation ( 5.3).

A Chromy sampling and pivotal sampling on an example

We consider the population U = {1, 2, 3, 4, 5} with π = (0.4, 0.8, 0.5, 0.6, 0.7).

The complete probability tree for Chromy sampling with parameter π is given in Figure 2. For example, at the first step the unit 1 is selected with probability 0.4, and discarded with probability 0.6. If unit 1 is selected, then 1 l=1 I l = V I 1 + 1. Since unit 2 is a cross-border unit, we follow Step 2.b of Algorithm 1 and unit 2 is selected at the next step with probability

V F 2 /V F 1 = 0.2/0.4 = 1/2.
The complete sampling design is

S c =                       
{1, 2, 4} with proba. 3/35, {1, 2, 5} with proba. 4/35, {1, 3, 4} with proba. 3/56, {1, 3, 5} with proba. 1/14, {1, 4, 5} with proba. 3/40, {2, 3, 4} with proba. 9/56, {2, 3, 5} with proba. 3/14, {2, 4, 5} with proba. 9/40.

(A.1) Now, we consider ordered pivotal sampling on the same population U with the same parameter π. The complete probability tree for pivotal sampling is given in Figure 3. For example, at the first step, the units 1 and 2 fight with respective probabilities 0.4 and 0.8. With probability (1 -0.8)/(2 -0.4 -0.8) = 1/4, unit 1 is selected and unit 2 gets the residual probability 0.2, and with the complementary probability unit 2 is selected and unit 1 gets the residual probability 0.2. In the first case, unit 2 faces unit 3 with respective probabilities 0.2 and 0.5. With probability 0.2/(0.2 + 0.5) = 2/7, unit 2 gets the sum of the probabilities and unit 3 is discarded, and with the complementary probability unit 3 gets the sum of the probabilities and unit 2 is discarded.

It follows from straightforward computations that the complete sampling design is the same as for Chromy sampling.

(0.4, 0.8, 0.5, 0.6, 0.7) (1, 0.2, 0.5, 0.6, 0.7)

I 1 = 1 I 2 = 1 I 3 = 0 I 4 = 1 & I 5 = 0 I 4 = 0 & I 5 = 1 I 2 = 0 I 3 = 1 I 4 = 1 & I 5 = 0 I 4 = 0 & I 5 = 1 I 3 = 0 I 4 = 1 & I 5 = 1 I 1 = 0 I 2 = 1 I 3 = 1 I 4 = 1 & I 5 = 0 I 4 = 0 & I 5 = 1 I 3 = 0 I 4 = 1 & I 5 =
(1, 0.7, 0, 0.6, 0.7)

(1, 1, 0, 0.3, 0.7)

(1, 1, 0, 1, 0)

(1, 1, 0, 0, 1)

(1, 0.3, 0, 1, 0.7)

(1, 1, 0, 1, 0)

(1, 0, 0, 1, 1)

(1, 0, 0.7, 0.6, 0.7)

(1, 0, 1, 0.3, 0.7)

(1, 0, 1, 1, 0)

(1, 0, 1, 0,

(1, 0, 1, 1, 0)

(1, 0, 0, 1, 1) (0.2, 1, 0.5, 0.6, 0.7) (0.7, 1, 0, 0.6, 0.7)

(1, 1, 0, 0.3, 0.7)

(

(1, 1, 0, 0, 1) The proof is by induction. For k = 1, we have I 1 ∈ {0, 1} and V I 1 = 0, so the property holds. Suppose that the property holds at k -1, namely

(0.3, 1, 0, 1, 0.7) (1, 1, 0, 1, 0) (0, 1, 0, 1, 1) (0, 1, 0.7, 0.6, 0.7) (0, 1, 1, 0.3, 0.7) (0, 1, 1, 1, 0) (0, 1, 1, 0, 1) (0, 1, 0.3, 1, 0.7) (0, 1, 1, 1, 0) (0, 1, 0, 1, 1)
V I k-1 ≤ k-1 l=1 I l ≤ V I k-1 + 1. (B.1)
First note that if k is a cross-border unit, we have

V I k = V I k-1 + 1 and if not, we have V I k = V I k-1 . • Suppose that k-1 l=1 I l = V I k-1 .
If k is a non cross-border unit, we have

V I k-1 ≤ k-1 l=1 I l + I k ≤ V I k-1 + 1 ⇔ V I k ≤ k l=1 I l ≤ V I k + 1.
If k is a cross-border unit, we obtain from Algorithm 1 that I k = 1, and

k l=1 I l = V I k-1 + 1 = V I k . • Suppose that k-1 l=1 I l = V I k-1 + 1. If k is a non cross-border unit, then from Algorithm 1, we have I k = 0 and k l=1 I l = V I k-1 + 1 = V I k .
If k is a cross-border unit, we obtain from I k ∈ {0, 1}

V I k-1 + 1 ≤ k-1 l=1 I l + I k ≤ V I k-1 + 2 ⇔ V I k ≤ k l=1 I l ≤ V I k + 1.

C Proof of Theorem 1

The proof proceeds in two main steps. We first prove that Chromy sampling may be alternatively seen as the result of a two-stage sampling procedure, inside a population U c of clusters which is introduced in Section C.1. We then consider in Section C.2 the ordered sample, which is given by the selected units ranked with respect to the natural order in the population, and we give the transition probabilities between the selected units. These results are used in Section C.3 to prove Theorem 1.

C.1 Clustered population

The N sampling units in the population U are grouped to obtain a population U c = {u 1 , . . . , u 2n-1 } of clusters. There are the clusters associated to the cross-border units (n -1 singletons), denoted as u 2i with associated probability φ 2i = π k i . There are the n clusters of non cross-border units that are between two consecutive integers, denoted as u 2i-1 with associated probability

φ 2i-1 = V k i -1 -V k i-1 .
The vector of inclusion probabilities in the population U c is denoted as φ = (φ 1 , . . . , φ 2n-1 ) ′ . For illustration, useful quantities for population U c are presented in Figure 4.

i - Proof. It is sufficient to prove that the transition probabilities given in Algorithm 1 are the same under Chromy sampling with parameter π and under the two-stage sampling procedure. We use the following notation: for any unit k ∈ U, recall that I k is the sample membership indicator under Chromy sampling with parameter π; for any cluster u i ∈ U c , J c,i is the sample membership indicator under Chromy sampling with parameter φ; for any unit k ∈ U, J k is the sample membership indicator under the two-stage procedure. We first note that, by definition of the two-stage procedure, we have

1 i i + 1 a i-1 b i-1 a i b i a i+1 b i+1 φ 2i-2 φ 2i-1 φ 2i φ 2i+1 φ 2i+2 U i-1 U i U i+1 U i+2
J c,i = k∈u i J k for any u i ∈ U c . (C.1)
We now consider the case when k ∈ U is a cross-border unit, k i say. The corresponding cluster is u 2i , and in such case

J k i = 1 ⇔ J c,2i = 1. (C.2)
We obtain successively

P r J k i = 1 k i -1 l=1 J l = P r J k i = 1 2i-1 j=1 J c,j from equation (C.1) = P r J c,2i = 1 2i-1 j=1 J c,j from equation (C.2) = 1 if 2i-1 j=1 J c,j = i -1, b i 1-a i if 2i-1 j=1 J c,j = i, from Step 2.b in Algorithm 1 = P r I k i = 1 k i -1 l=1 I k , (C.3)
where the last line is obtained again from Algorithm 1.

We now consider the case when k is not a cross-border unit, and belongs to the cluster u 2i-1 , say. We begin by computing the quantities

P r   J k = 1 2i-2 j=1 J c,j , k-1 l=k i-1 +1 J l   .
Note that from Proposition 1, 2i-2 j=1 J c,j may only take the values i -1 and i, and since we select at most one unit l inside u 2i-1 , k-1 l=k i-1 +1 J l may only take the values 0 and 1. If k-1 l=k i-1 +1 J l = 1, we have J k = 0 since we select at most one unit l inside u 2i-1 . Therefore

P r   J k = 1 2i-2 j=1 J c,j = i -1, k-1 l=k i-1 +1 J l = 1   = 0. (C.4)
If 2i-2 j=1 J c,j = i, Algorithm 1 implies that the cluster u 2i-1 and therefore k may not be selected. Therefore

P r   J k = 1 2i-2 j=1 J c,j = i, k-1 l=k i-1 +1 J l = 0   = 0. (C.5)
For the same reason, we may not have simultaneously 2i-2 j=1 J c,j = i and k-1 l=k i-1 +1 J l = 1. Finally, we have

P r   J k = 1 2i-2 j=1 J c,j = i -1, k-1 l=k i-1 +1 J l = 0   = P r J k = 1, k-1 l=k i-1 +1 J l = 0 2i-2 j=1 J c,j = i -1 P r k-1 l=k i-1 +1 J l = 0 2i-2 j=1 J c,j = i -1 = P r J k = 1 2i-2 j=1 J c,j = i -1 P r k-1 l=k i-1 +1 J l = 0 2i-2 j=1 J c,j = i -1 , (C.6)
where the last line in (C.6) follows from the fact that if J k = 1, we necessarily have k-1 l=k i-1 +1 J l = 0. We compute the numerator and the denominator in (C.6) separately. The numerator is

P r J k = 1 2i-2 j=1 J c,j = i -1 = P r J c,2i-1 = 1 2i-2 j=1 J c,j = i -1 P r J k = 1 2i-2 j=1 J c,j = i -1, J c,2i-1 = 1 = P r J c,2i-1 = 1 2i-2 j=1 J c,j = i -1 P r (J k = 1 |J c,2i-1 = 1 ) . From Algorithm 1, we have P r J c,2i-1 = 1 2i-2 j=1 J c,j = i -1 = 1-b i-1 -a i 1-b i-1
, and from the definition of the two-stage procedure P r

(J k = 1 |J c,2i-1 = 1 ) = π k 1-b i-1 -a i . This leads to P r J k = 1 2i-2 j=1 J c,j = i -1 = π k 1 -b i-1 . (C.7)
From the definition of the two-stage procedure, the denominator in (C.6) is

  J k = 1 2i-2 j=1 J c,j = i -1, k-1 l=k i-1 +1 J l = 0   = π k 1 -b i-1 -l<k∈u 2i-1 π l = V F k -V F k-1 1 -V F k-1
.

(C.9) From (C.4), (C.5) and (C.9), we obtain P r

  J k = 1 2i-2 j=1 J c,j + k-1 l=k i-1 +1 J l = i -1   = V F k -V F k-1 1 -V F k-1 , P r   J k = 1 2i-2 j=1 J c,j + k-1 l=k i-1 +1 J l = i   = 0. (C.10)
From equation (C.1), we have (C.11) and from (C.10), this implies that

2i-2 j=1 J c,j + k-1 l=k i-1 +1 J l = k i-1 l=1 J l + k-1 l=k i-1 +1 J l = k-1 l=1 J l ,
P r J k = 1 k-1 l=1 J l = V F k -V F k-1 1-V F k-1 if k-1 l=1 J l = i -1, 0 if k-1 l=1 J l = i, = P r I k = 1 k-1 l=1 I l , (C.12)
where the last line in (C.12) follows from Step 2.b in Algorithm 1. This completes the proof.

C.2 Ordered sample

We use the same notation as in Section C.1, and consider a sample S c selected in U c by means of Chromy sampling with parameter φ. In this case, we let V c,i = i j=1 φ j denote the cumulated inclusion probabilities up to unit u i , with V c,0 = 0. The integer part of V c,i is denoted as V I c,i , and the difference between V c,i and its integer part V I c,i is the fractional part, denoted as V F c,i .

Let us denote by X 1 < . . . < X n the selected units, ranked with respect to the natural order in the population U. The transition probabilities between the ranked selected units are given in Proposition 3.

Proposition 3. Let S c denote a sample selected by Chromy sampling with parameter φ in the population U c . The transition probabilities between the Now, we consider the case when X i = u 2i-1 . We have

X i = u 2i-1 ⇒ 2i-2 j=1 J c,j = i -1 and J c,2i-1 = 1 ⇒ 2i-1 j=1 J c,j = i.
(C.17)

Since at each step of Chromy sampling, the conditional probabilities only depend on the number of units already selected, this leads to P r(X i+1 = u j |X i = u 2i-1 ) = P r X i+1 = u j 2i-1 j=1 J c,j = i , (C.18) which is identical to equation (C.13). Therefore, the second equation in Proposition 3 follows.

Finally, we consider the case when X i = u 2i , which is equivalent to 2i j=1 J c,j = i. If j = 2i + 1:

Figure 2 :

 2 Figure 2: Probability tree for Chromy sampling on a population U of size N = 5

Figure 3 :

 3 Figure 3: Probability tree for ordered pivotal sampling on a clustered population U c 20

Figure 4 :

 4 Figure 4: Inclusion probabilities and cross-border units in microstrata U i and U i+1 for population U c

  P r(X i+1 = u 2i+1 |X i = u 2i ) = P r X i+1 = u i+1b i 1b i . (C.19)

Table 1 :

 1 Population mean and population dispersion for two populations and four variables of interest

	linear	quadratic exponential	bump
	µ y1 S 2 y1	µ y2	S 2 y2	µ y3	S 2 y3	µ y4 S 2 y4
	Population 1 10.1 13.1 11.7 69.8 10.3	18.2	12.3 84.6
	Population 2 10.1 5.3	8.8	8.3	9.9	4.4	5.2 18.0

Table 2 :

 2 Range of inclusion probabilities proportional to x for two populations and three sample sizes

	n = 50	n = 100	n = 200
	Min Max Min Max Min Max
	Population 1 0.030 0.305 0.060 0.609 0.122 1.000
	Population 2 0.079 0.791 0.160 1.000 0.325 1.000

Table 3 :

 3 Relative bias (in % ), Relative Root Mean Square Error (in % ) of the variance estimator and Coverage Rate of the normality-based confidence intervals for three populations and four variables of interest

			n = 50			n = 100			n = 200
		RB RMSE Cov. Rate RB RMSE Cov. Rate RB RMSE Cov. Rate
	linear	-0.8	76	11.20	-0.5	46	7.20	-1.0	28	5.50
	Population 1 quadratic	3.2	118	13.50	-0.6	66	8.60	-2.0	37	7.00
	exponential -0.4	85	11.40	-0.9	50	7.10	-1.1	30	6.50
	bump	3.7	103	11.30	-0.4	59	8.30	-2.1	32	7.00
	linear	-0.7	56	8.20	0.6	31	7.30	0.0	17	4.80
	Population 2 quadratic	-0.7	56	8.60	0.6	30	7.70	0.0	17	4.80
	exponential -0.6	55	8.10	0.5	30	7.50	0.0	17	4.90
	bump	-1.4	61	11.10	0.8	35	6.90	0.9	19	4.30

ordered sampled units X 1 < . . . < X n are:

Proof. We first consider the case when

where the second line in (C.13) follows from the fact that, from Step 2.a of Algorithm 1, 2i-2 j=1 J c,j = i implies that J c,2i-1 = 0. If j = 2i:

where the last but one line in (C.15) follows from Step 2.a of Algorithm 1 and from equation (C.14). If j = 2i + 2, we obtain similarly:

This gives the first equation in Proposition 3.

This completes the proof.

C.3 Proof of Theorem 1

From Lemma 3.1 in [START_REF] Chauvet | On a characterization of ordered pivotal sampling[END_REF] and our Proposition 2, Chromy sampling and ordered pivotal sampling have the same two-stage characterization. It is therefore sufficient to prove that they lead to the same sampling design when sampling with parameter φ in the clustered population U c . However, from equations (4.2)-(4.4) in [START_REF] Chauvet | On a characterization of ordered pivotal sampling[END_REF] and our Proposition 3, both Chromy sampling and ordered pivotal sampling have the same transition probabilities between ordered sampled units, which means that the induced sampling designs are identical.