Laurent Gornet

R Desmorat

Gilles Marckmann

Pierre Charrier

A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis

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

INTRODUCTION

The aim of the present paper is to develop new constitutive models for the multiaxial loadings of elastomers which are both physically motivated, well adapted for numerical problems, and accurate for rubber materials. Hyperelastic models dedicated to rubber materials can be classified into two types of strain energy formulations. The first kind of models is issued from mathematical developments such the well-known Rivlin series or as the Ogden model [START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF]. The second kind of models is the one developed from physical motivations. Such models are based on both physics of polymer chains networks and statistical methods. In this work, we propose a strain energy function expressed in terms of independent strain invariants. The proposed model (GD) is dedicated to multiaxial loadings. This model is successfully identified on both Treloar and Kawabata experiments (Treloar 1944, Kawabata et al. 1981). We present a bridge between the phenomenological strain energy formulation and the physical motivation of the proposed (GD) model. In the Finite Element context, the proposed model can easily be implemented because of its strain invariants formulation. This has been done in computer codes like Cast3M-CEA (French Atomic Agency) and AB-AQUS. Finite Element predictions of a roll restrictor developed by Trelleborg automotive have been carried out with ABAQUS. The continuous damage approach is applied to soften by damage the material parameters of this initial model in order to describe Mullins effect. A new model for Mullins effect (GDM model) is therefore derived. The identification of the GDM model is made thanks to experimental data from [START_REF] Marckmann | Comparison of hyperelastic models for rubberlike materials[END_REF].

STRAIN ENERGY FUNCTIONS

In this section, the strain energy functions that define the proposed model (GD) are briefly recalled. Links between this model and the eight-chains model are highlighted. A phenomenological energy function expressed in term of the second invariant I 2 is also introduced based on physical motivations.

Incompressible GD strain energy density

Assuming that rubber materials are both isotropic and incompressible, the proposed strain energy function W only depends on the two first invariants of the left Cauchy-Green stretch tensor B:

( ) ( ) ( ) ( ) 3 1 1 2 1 1 2 2 2 3 1 1 2 2 2 , 1 3 GD h I W I I W I W I h e dI h dI I - = + = + ∫ ∫ (1)
where h 1 , h 2 , h 3 are the material parameters. In this strain energy, the I 1 part of W describes the global response of the material and is equivalent to the Hart-Smith model [START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF]. The second term that involves I 2 improves the accuracy of the model for multiaxial loading conditions (a general power I 2 -term has been introduced by Lambert-Diani & Rey, 1999). The true stress tensor is defined by the differentiation of the proposed strain energy with respect to B:

A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis

P. Charrier

Modyn Trelleborg, Zone ind. de Carquefou,BP 419,44474 Carquefou Cedex -France ABSTRACT: The present paper focuses on static stiffness modelling of rubber materials for multiaxial loadings with a minimal number of material parameters in order to ensure robustness of both identification and Finite Element analysis. A physically motivated isotropic constitutive model is proposed for multiaxial large stretch rubber deformation. This model is expressed in terms of classical independent strain invariants. Only three parameters are needed to successfully represent both Treloar and Kawabata experiments (Kawabata et al. 1981) with a response equivalent to the one obtained by Ogden six parameters model [START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF]. Mullins effect modelling is finally derived according to a damage mechanics approach (Chagnon et al. 2004).

1 1 2 2 2 2 W W W pI I B B I I I σ   ∂ ∂ ∂ = -+ + -   ∂ ∂ ∂   (2)
( )

1 I Tr B = , ( ) ( ) 2 2 2 1 2 I Tr B Tr B   = -   (3)

Links between GD and eight-chains models

The present part is devoted to the comparison of the GD model to the eight-chains one (Arruda and Boyce 1993) in order to highlight the physical motivations of this proposed phenomenological model. As these constitutive equations are supposed to be qualitatively efficient for the entire range of strains, both small and large strain responses are compared.

Comparison of these models is established using their polynomial expansion in terms of I 1 respectively (equations 4 and 5) for small strain. High order terms are then eliminated:

( ) ( ) ( ) 2 1 3 1 1 1 1 0 3 2 1 ! i i i C W I C I i i ∞ + = = - + ∑ (4)
The first five terms of the development of the eightchains model are:

  = -+ -       + - + -       + - +     (5)
The small strain stiffness of the models are respectively defined by first terms of Equations ( 4) and ( 5). Then, both the strain energy functions reduce to the neo-Hookean expression. The both models are characterized by their ability to describe the strainhardening of the material that takes place under large strains. This strain-hardening phenomenon is mainly due to the extensibility limit of polymer chains. Chagnon et al. (2004) established that the first part of the Hart-Smith strain energy W 1 (I 1 ), which is identical to the first part of the GD model, is equivalent to the eight-chains and Gent models for the entire range of strains.

( )

( ) 2 3 1 3 1 1 1 1 h I W I h e I - ∂ = ∂ (6)
We proposed now to highlight the physical motivation of the second part of strain energy W 2 (I 2 ) of the GD model. As mentioned by Treloar [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], term as function of I 2 in constitutive equations can be seen as corrections of the phantom network theory (terms as function I 1 ). We propose to constrain the eight-chains model by a new network of chains on the surface of the cube (figure 1). The confinement of the eight-chains model is governed by a strain energy potential. This potential constrains the eight-chains cube surface. Let us recall that the surface of the eight-chains model (a cube) is I 2 1/2 and that its increase under deformation is I 2 1/2 -3 1/2 . We therefore define a pressure constrain of the eightchains rubber network. This phenomenon is modeled by the second invariant energy part :

( ) 2 2 2 2 2 1 3 W I h dI I = ∫ or ( ) 2 2 2 2 2 3 W I h I I ∂ = ∂ ( 7
)
where 3h 2 stands for the pressure constrain of entanglement of the eight-chains rubber network.

Compressible GD strain energy density

In this part we introduce the proposed strain energy function for compressible isotropic hyperelastic materials in terms of strain invariants. Based on kinematic assumption, we use a decoupled representation of the strain energy function originally proposed by [START_REF] Flory | Thermodynamic relations for highly elastic materials[END_REF][START_REF] Ogden | Recent advances in the phenomenological theory of rubber elasticity[END_REF].

( ) ( ) ( ) 1 2 1 2 , , , iso vol W I I J W I I W J = + (8) ( ) det J F = , 1 3 F J F - = , T B F F = (9) ( ) 1 I Tr B = , ( ) ( ) 2 2 2 1 2 I Tr B Tr B   = -     (10)
Where W iso (I 1 ,I 2 ) and W vol (J) stand for the isochoric and volumetric elastic strain energy of the material, respectively. The isochoric part is equivalent to the equation ( 1) with modified invariants. The volumetric part is chosen to enforce a nearly quasi incompressible behaviour.

( ) ( )

1 3 2 3 1 2 1 2 1 2 2 1 , h I iso W I I h e d I h d I I - = + ∫ ∫ (11) ( ) 2 1 1 vol W J D = - (12)
The volumetric properties part is presented in [START_REF] Doll | On the developpement of volumetric strain energy function[END_REF]. The true stress tensor is defined by the differentiation of the proposed strain energy with respect to B:

5 7 2 3 3 1 1 2 2 2 2 iso iso iso iso W W W J I B J B I I I σ - - ∂ ∂ ∂   = + -   ∂ ∂ ∂   (13) vol vol W I J σ ∂ = ∂ , vol iso σ σ σ = + (14)
3 EXPERIMENTAL DATA AND IDENTIFICATION OF PARAMETERS In order to compare the efficiency of the models, we choose two complementary data sets issued from classical references [START_REF] Marckmann | Comparison of hyperelastic models for rubberlike materials[END_REF]. The first set is due to Treloar (1944). In the current study, data from Treloar for unfilled natural rubber (cross-linked with 8 parts of S phr) were used. This material exhibits highly reversible elastic response and no stretch-induced crystallization up to 400%. Thus it is well-modeled by hyperelastic constitutive equations. Experimental measures were performed for four different loading conditions: equibiaxial extension of a sheet (EQE), uniaxial tensile extension (UE), pure shear (PS) and biaxial extension (BE). The second data set is due to Kawabata et al. (1981). It was obtained using an experimental apparatus for general biaxial extension testing. In terms of stretch ratios, unfilled polyisoprene specimens were stretched from 1.04 to 3.7 in the first direction (λ 1 ) and from 0.52 to 3.1 in the perpendicular direction (λ 2 ). These values correspond to moderate strain but lead to deformation conditions from uniaxial extension to equibiaxial extension. Here, both experimental data sets are simultaneously considered to compare the models because the two materials are quite similar. Thus, for a given model, a unique set of material parameters must be able to reproduce these data with a good agreement. The parameter identification is performed using genetic algorithms as presented in [START_REF] Marckmann | Comparison of hyperelastic models for rubberlike materials[END_REF] The GD parameters are: h 1 =0.142 MPa, h 2 =1.585 10 -2 MPa and h 3 =3.495 10 -4 . The proposed model GD is able to accurately reproduce the whole "S" shaped response of the material. The model behaves satisfactory under all the presented loadings (figures 2-4).

The Ogden six parameter strain energy density is classically

MULLINS EFFECT

Elastomers present a loss of stiffness after the first loading cycle of a fatigue experiment [START_REF] Mullins | Softening of rubber by deformation[END_REF]. It has been observed that this phenomenon is only dependent on the maximum deformation previously reached in the history of the material. It is quite important to model it because the mechanical behaviour of rubber products is highly modified by this softening phenomenon. Moreover, as the Mullins effect depends on the maximum deformation endured previously, material REV of the product are not identically affected. As a consequence, it is not acceptable to determine experimentally an accommodated hyperelastic constitutive equation for the material. The stress-softening should be explicitly included in the model. The continuum damage me-chanics has often been used to model the Mullins effect even if phenomenon undergoing Mullins effect is not a strictly speaking an irreversible damage phenomenon and can even be modelled without damage [START_REF] Cantournet | Müllins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model[END_REF]. For example, it can be recovered with time and annealing accelerates this recovery. A thermodynamic variable d is introduced to represent the loss of stiffness and the corresponding stress-softening. The general theory of Continuum Damage Mechanics is detailed in [START_REF] Lemaitre | Mechanics of solid materials[END_REF] book. When applied to GD hyperelastic density, it yields to GDM model. This new strain energy function for hyperelastic model with damage variables can be written (eq 17) :

( )

( ) ( ) ( ) ( ) 3 1 2 3 1 2 1 1 2 2 2 1 1 1 2 2 2 3 3 3 1 , 3 1 , 1 , 1 h I GDM

W I I h e dI h dI I h h d h h d h h d

- = + = - = - = - ∫ ∫ % % % % % % (17)
The incompressible state laws (eq. 2) associated with this model are classically obtained by equation 18.

1 1

GDM I W B σ = ∂ = ∂ (18)
We define the thermodynamic forces

i d Y associated
with damage internals variables i d by:

i GDM d i W Y d ∂ -= ∂ (19)
A non standard damage model is build here as the damage thermodynamics forces are not used to describe damage evolution. We consider instead that the quantity governing the damage evolution laws is the maximum of the first invariant. This assumption is achieved according to the physical motivations of maximum strain state endured during the history of the deformation [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF]. In order to do this we can introduce damage criterion functions

f i = I 1 -k i (d i ) (20)
such as f i <0 implies no damage evolution and as damage evolves at f i =0 with then d i -k i -1 (I 1 ). Choosing particular expressions for k i -functions allows to derive the damage evolution laws (equations 19-20) and make the model complete. The proposed GDM model is based on an improved method already proposed in (Chagnon et al. 2004). The evolution equation of the damage variable is expressed thanks to the first strain invariant and presents an exponential form (equations 21). This model is able to represented unloadings Mullins effect responses (figures 11 and 12).

max 1 1 1 1 1 exp I d d η ∞     = - -         (21) max 1 2 2 2 1 exp I d d η ∞     = - -        
where d 1∞, η 1, d 2∞, η 2 and b are material parameters.

I 1 max represents the maximum value of the first strain invariant obtained during loadings. The coupling with damages d 1 and d 2 is similar to the discontinuous damage part of the constitutive equations proposed by [START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden type large strain elastic materials[END_REF]. Considering the mass conservation of the polymer network, which implies that the number of monomer segments per unit volume N n must remain constant [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF], where N is the mean number of monomer per chain and n the number of chain per unit of volume. According to relationship between parameters derived in (Chagnon et al. 2004), this leads to link the mechanical properties 1 h % and 3 h % and shows that d 3 is not an independent thermodynamics damage variable, ( )

3 1 1 d F d = - , ( ) ( ) ( ) 1 2 1 1 3 1 F d b h = - % (22)
The GDM is identified on idealized experimental data as presented in (Chagnon et al. 2004). Behaviour is considered time-independent.

FINITE ELEMENT PREDICTIONS

In the finite element context, the proposed GD and GDM models can easily be implemented because of their strain invariants formulations. This has been done in computer codes Cast3M-CEA (French Atomic Agency) and ABAQUS. The two previous constitutive models were implemented in the finite element context, thanks to the UMAT and UHYPER facilities.

Figure 1 .

 1 Figure 1. Entanglement of the eight-chains rubber network. Eight-chains model completed by surrounding chains (in bold) is the physical motivations of the GD model.

  . The GD model includes entanglement of the eight-chains rubber network and identified responses are depicted on figures 2-4. The GD model responses on figures 2 to 4 are almost equivalent to Ogden six parameters results for bi-axial extensions presented on figures 5 to 7. Harts-Smith identified responses are depicted on figures 8-10.

Figure 2 .

 2 Figure 2. Experimental data (--) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.

Figure 3 .

 3 Figure 3. Experimental data (--) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.

Figure 4 .

 4 Figure 4. Experimental data (.) and GD model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible GD model response.

Figure 5 .

 5 Figure 5. Experimental data (--) and Ogden model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Ogden six model response.

Figure 6 .

 6 Figure 6. Experimental data (--) and Ogden model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Ogden six model response.

Figure 7 .

 7 Figure 7. Experimental data (.) and Ogden model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible Ogden six model response.

  µ 1 =0.63 MPa, µ 2 =1.2 10 -3 MPa, µ 3 =-1 10 -2 MPa, α 1 =1.3, α 2 =5 and α 3 =-2. The conditions µ i α i >0 ensures the positive definite character of the strain energy. In order to recall the influence of W 2 (I 2 ), let us recall the classical strain energy of Hart-Smith (equation 16).

  parameters are: a 1 =0.140 MPa, a 2 =5.25410 -4 and a 3 =1.290.

Figure 8 .

 8 Figure 8. Experimental data (--) and Hart-Smith model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Hart-Smith model response.

Figure 9 .

 9 Figure 9. Experimental data (--) and Hart-Smith model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Hart-Smith model response.

Figure 10 .

 10 Figure 10. Experimental data (.) and Hart-Smith model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible Hart-Smith model response.

Figure 11 .

 11 Figure 11. Tensile behaviour of a hyperelastic incompressible GDM model with stress-softening. Identification is performed for a material used on a roll restrictor.

Figure 12 .

 12 Figure 12. Pure Shear behaviour of a hyperelastic incompressible GDM model with stress-softening. Identification is performed for a material used on a roll restrictor.

Figure 13 .

 13 Figure 13. Engine roll restrictor made of steel and rubber parts.

Finite

  element predictions of a roll restrictor (figure13) developed by Trelleborg automotive have been performed with ABAQUS and the GDM model. This behavior law is able to describe the local loss of stiffness of the material and the non-homogeneity of the structure after the first "first gear full torque" engine acceleration. The material properties used for finite element simulations are not the ones used for real structures. Predictions are performed with the GDM model discussed in previous sections. The 3D model is composed of C3D8H elements. It corresponds to an incompressible finite element hybrid formulation. It appears that the stress softening level is heterogeneous in the roll restrictor rubber parts (figures 14-16). Damage is located on the surface for most part of the structure excepted in the both snubbers.

Figure 14 .

 14 Figure 14. White parts represent the damage level d 1 (Mullins effect) pattern in the roll restrictor during the first "first gear full torque" engine's acceleration.

Figure 15 .

 15 Figure 15. White parts represent the damage level d 2 (Mullins effect) pattern in the roll restrictor during the first "first gear full torque" engine's acceleration.

Figure 17 Figure 16 .

 1716 Figure 17 illustrates the above Mullins observations in respect of a GDM model with isotropic damage evolution law (equation 23). Mullins effect is located on the surface for most part of the structure excepted in the both snubbers.

Figure 17 .

 17 Figure 17. Isotropic damage model (d = d 1 = d 2 , d 3 =0). White parts represent the damage area "d" (Mullins effect) during the first "first gear full torque" engine's acceleration.

  GeM, UMR CNRS 6183,Ecole Centrale Nantes, 1 rue de laNoë, BP 92101-44321 Nantes, Cedex 3, France.

	R. Desmorat
	LMT-Cachan ENS Cachan/CNRS/UPMC/PRES UniverSud Paris, 61 avenue du Président Wilson, F-94230
	Cachan Cedex, France.

L. Gornet & G. Marckmann


Laurent Gornet

Gilles Marckmann

R Desmorat

Pierre Charrier

Unified model for Mullins effect and high cycle fatigue life prediction of rubber materials

INTRODUCTION

The aim of the present paper is to develop new constitutive model for lifetime prediction of elastomers which are physically motivated, well adapted for numerical problems, and accurate for rubber materials under multiaxial high cycle fatigue loadings. We propose a strain energy function expressed in terms of independent strain invariants restricted, in the present paper, to uniaxial high cyclic loadings. A way to predict the fatigue life is to calculate the damage accumulation cycle by cycle using the continuum damage mechanics approach [START_REF] Lemaitre | Mechanics of solid materials[END_REF]. For rubber-like materials, CDM approaches are mainly focused on the Mullins effect (Chagnon et al. 2004, 2006, Cantournet 2002[START_REF] Cantournet | Mullins effect and cyclic stress softening of filled elastomer by internal sliding and friction thermodynamics model[END_REF]. Few models are developed for fatigue phenomena [START_REF] Wang | A damage model for fatigue life prediction under multiaxial loading[END_REF], Ayoub et al. 2012). In the present study, the damage model [START_REF] Gornet | A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis[END_REF] initially developed for the multiaxial Mullins effect (few cycles) is extended to high cycle fatigue loadings. In the Finite Element context, the proposed model can easily be implemented thanks to its strain invariants formulation. Finite Element predictions of samples have been carried out with the Finite Element code Abaqus.

STRAIN ENERGY FUNCTION

In this section, the hyperelastic strain energy function that defines the foundation of the proposed fa-tigue model is briefly recalled. For incompressible materials, a phenomenological energy function is expressed in term of the first I 1 and the second invariant I 2 [START_REF] Gornet | A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis[END_REF]).

Incompressible GD strain energy density

Assuming that rubber materials are both isotropic and incompressible, the proposed strain energy function W only depends on the two first invariants of the left Cauchy-Green stretch tensor B:

( ) ( ) ( ) ( ) 3 1 1 2 1 1 2 2 2 3 1 1 2 2 2 , 1 3 GD h I W I I W I W I h e dI h dI I - = + = + ∫ ∫ (1)
where h 1 , h 2 , h 3 are the material parameters. In this strain energy, the term depending on I 1 of W describes the global response of the material and is similar to the first part of the Hart-Smith model [START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF]. The second term that involves I 2 improves the accuracy of the model for multiaxial loading conditions (a general power I 2 -term has been introduced by Lambert-Diani & Rey, 1999). The true stress tensor is defined by the differentiation of the proposed strain energy with respect to B:

2 1 1 2 2 2 2 W W W pI I B B I I I σ   ∂ ∂ ∂ = -+ + -   ∂ ∂ ∂   (2)

P. Charrier

TrelleborgVibracoustic, Zone ind. de Carquefou,BP 419, ABSTRACT: The study describes the basic principles of a general damage model (GDMF) for Mullins effect and high cycle fatigue loadings of rubber materials and demonstrates its prediction possibilities for simulating the complete fatigue failure phenomenon. The present paper focuses on stiffness modelling of rubber materials for uniaxial and multiaxial static and fatigue loadings with a minimal number of material parameters in order to ensure robustness of the identification. The proposed hyperplastic model is expressed in terms of classical independent strain invariants. Mullins effect and high cycle fatigue loadings are both modelled according to a continuum damage mechanics approach.

( )

1 I Tr B = , ( ) ( ) 2 2 2 1 2 I Tr B Tr B   = -   (3)

Links between GD and eight-chains models

The hyperelastic model is characterized by its ability to describe the strain-hardening of the material that takes place under large strains. This strain-hardening phenomenon is mainly due to the extensibility limit of polymer chains. The first part of the GD model is equivalent to Hart-Smith model and then to the eight-chains [START_REF] Arruda | A three dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF] and Gent models for the entire range of strains (Chagnon G. et al. 2004).

( )

( ) 2 3 1 3 1 1 1 1 h I W I h e I - ∂ = ∂ (4)
We proposed to remind the physical motivation of the second part of strain energy W 2 (I 2 ) of the GD model. As mentioned by [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], term that depends on I 2 in constitutive equations can be seen as corrections of the phantom network theory (I 1 function). We have proposed to constrain the eightchains model by a new network of chains on the surface of the cube (figure 1). The confinement of the eight-chains model is governed by a strain energy potential. This potential constrains the eight-chains cube surface. Let us recall that the surface of the eight-chains model (a cube) is I 2 1/2 and that its increase under deformation is I 2 1/2 -3 1/2 . We therefore define a pressure constrain of the eight-chains rubber network. This phenomenon is modeled by the second invariant energy part:

( ) 2 2 2 2 2 1 3 W I h dI I = ∫ or ( ) 2 2 2 2 2 3 W I h I I ∂ = ∂ (5)
where 3h 2 stands for the pressure constrain of entanglement of the eight-chains rubber network. 3 MULTIAXIAL EXPERIMENTAL DATA AND GD PARAMETERS IDENTIFICATION In order to compare the efficiency of the models, we choose two complementary data sets issued from classical references (Marckmann & Verron 2006). The first set is due to Treloar (1944). In the current study, data from Treloar for unfilled natural rubber (cross-linked with 8 parts of S ) were used. This material exhibits highly reversible elastic response and no stretch-induced crystallization up to 400%. Thus it is well-modeled by hyperelastic constitutive equations. Experimental measurements were performed for four different loading conditions: equibiaxial extension of a sheet (EQE), uniaxial tensile extension (UE), pure shear (PS) and biaxial extension (BE). The second data set is due to Kawabata et al. (1981). It was obtained using an experimental apparatus for general biaxial extension testing. In terms of stretch ratios, unfilled polyisoprene specimens were stretched from 1.04 to 3.7 in the first direction (λ 1 ) and from 0.52 to 3.1 in the perpendicular direction (λ 2 ). These values correspond to moderate strain but lead to deformation conditions from uniaxial extension to equibiaxial extension. Here, both experimental data sets are simultaneously considered to compare the models because the two materials are quite similar. Thus, for a given model, a unique set of material parameters must be able to reproduce these data with a good agreement. The parameter identification is performed using genetic algorithms as presented in [START_REF] Marckmann | Contribution à l'étude des elastomers et des membranes souflées[END_REF]). The GD model assumes entanglement of the eight-chains rubber network and identified responses are depicted on figures 2-4. The GD model responses on figures 2 to 4 are almost equivalent to Ogden six parameters [START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF] results. The GD parameters are: h 1 =0.142 MPa, h 2 =1.585 10 -2 MPa and h 3 =3.495 10 -4 . The proposed GD model is able to accurately reproduce the whole "S" shaped response of the material. The model behaves satisfactory under all the presented loadings (figures 2-4).

A UNIFIED DAMAGE MODEL FOR MULLINS AND FATIGUE MODELLING

Material

The material studied here is a carbon black-filled natural rubber. Its recipe and mechanical characteristics are given in Table 1. Samples are obtained by injection molding in order to reproduce industrial conditions of mass produced parts. The compound has been cured for 7 min, with a mold temperature set to 165°C. All samples were made using the same material batch to limit properties scattering due for example to mixing.

Damage modeling

Elastomers present a loss of stiffness after the first loading cycle of a fatigue experiment [START_REF] Mullins | Softening of rubber by deformation[END_REF]. It has been observed that this phenomenon is only dependent on the maximum deformation previously reached in the loading history of the material.

It is quite important to model it because the mechanical behaviour of rubber products is highly modified by this softening phenomenon. The stress-softening for fatigue loadings should be explicitly included in the model. The continuum damage mechanics has often been used to model the Mullins effect even if phenomenon undergoing Mullins effect is not a strictly speaking an irreversible damage phenomenon. For example, it can be recovered with time and annealing accelerates this recovery. The general theory of Continuum Damage Mechanics is detailed in [START_REF] Lemaitre | Mechanics of solid materials[END_REF] book. When applied to GD hyperelastic density, it yields to GDM model. This new strain energy function for hyperelastic model with damage variables can be written (eq 6) :

( )

( ) ( ) ( ) ( ) 3 1 2 3 1 2 1 1 2 2 2 1 1 1 2 2 2 3 3 3 1 , 3 1 , 1 , 1 h I GDM

W I I h e dI h dI I h h d h h d h h d

- = + = - = - = - ∫ ∫ ɶ ɶ ɶ ɶ ɶ ɶ (6)
The incompressible state laws (eq. 2) associated with this model are classically obtained by equation 7.

1 1

GDM I W B σ = ∂ = ∂ (7)
We define the thermodynamic forces 1, 2,3

iM iF i d d d i • • • = + = (9)
Where the term iM d

• corresponds to Mullins Effect under few cyclic static loading and iF d

• is related to high cyclic fatigue one. Variables i d

• represent the total damage rate for each mechanical constant. The details of Mullins and fatigue damage variables will be discussed hereafter.

Mullins effect damage evolution law

For the Mullins part, a non-standard damage model is built here, the damage thermodynamics forces are not used to describe damage evolution laws. Instead of that, we consider that the quantity governing the damage evolution laws is the maximum value of the first invariant. This assumption is achieved according to the physical motivations of maximum strain state endured during the history of the deformation [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF]. In order to do this we can introduce damage criterion functions

f i = I 1 -k i (d iM )
(10) such as f i <0 implies no damage evolution and as damage evolves at f i =0 with then d iM =k i -1 (I 1 ). Choosing particular expressions for k i -functions allows to derive the damage evolution laws (equations 8-13) and makes the model complete for Mullins effect. This model is able to represent unloading Mullins effect responses (figure 12). Evolution equations of the damage variables are expressed thanks to the first strain invariant.

max 1 1 1 1 1 exp M I d d η ∞     = - -         (11) max 1 2 2 2 1 exp M I d d η ∞     = - -        
where d 1∞, η 1, d 2∞, η 2 and b are material parameters.

I 1 max (or 1 I in the following) represents the maximum value of the first strain invariant obtained during loadings. The coupling with damages d 1M and d 2M is similar to the discontinuous damage part of the constitutive equations proposed by [START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden type large strain elastic materials[END_REF]. Let us consider the mass conservation of the polymer network, which implies that the number of monomer segments per unit volume must remain constant [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF]. For the 8 chains model this assumption leads to: R C N cst = . This equation is derived for the GDM model. Chagnon et al. 2004, suggest a relationship between parameters of the 8 chains model, HartSmith model and the Gent model. These relations lead to equation:

1 2 1 3 m R J C N h   = +     (12)
Where the maximum value of 1 3 Iis m J in the Gent models.

According to relationship between parameters this leads to link the mechanical properties 1 h ɶ and 3 h ɶ and shows that d 3 is not an independent thermodynamics damage variables.

( )

3 1 1 M M d F d = - , ( ) ( ) ( ) 1 2 3 1 3 / 1 M b F d h a h = - ɶ (13)

High cycle fatigue damage evolution law

The proposed high cycle fatigue part of the damage model is based on an improved method already proposed in [START_REF] Gornet | High cycle fatigue damage model for delamination crack growth in CF/epoxy composite laminates[END_REF] for composite materials. Hence, in case of fatigue loadings, we assume that the damage variables of the rubber part will be also governed by the maximum value of the first invariant during the cyclic loadings. In this article, a fatigue model based on damage maximum value of the first invariant release rates is proposed as follows:

1 1 1min 1 1 1 1 1 ( ) , , 0 0 0 0 0 i i C C i F i F F F I I t I g d R I I I d d if I and f t if I or f • • • • •       =         ∂  = =  ≥ ≥ ∂     < <  (14)
Where, F f is a damage loading function and defines the threshold of fatigue damages growth. This function can be written in terms of the first strain invariant as * -=

1 1 I I f F
, where * 1 I is the first invariant threshold value and damage will grow only and only if 0

≥ F f
. This threshold value * 1

I is assumed zero in this paper. Here, i g is a dimensionless function and depends on maximum value of the first invariant on the cycle 1 I . Since the damage growth defined by Equation ( 14) is in rate format, it should be integrated over each time increment in the numerical analysis in or-der to obtain the damage at the end of increment. The damage variable at the end of a time increment t ∆ can be written as:

( ) ( ) ( ) N N iF iF iF N iF iF d N N d N d dt d d N N • +∆ + ∆ = + ∂ = + ∂ ∫ 1 ( , ) ( ) ( )
( ) i t t iF iF iF t N N d I iF n N d t t d t d dt d t • +∆ +∆ Ρ = + ∆ = + = + ∫ ∑ ( 15
)
Here t and t t ∆ + are the times corresponding to end of cycles N and N N ∆ + respectively. Here 1 ( , )

i d I Ρ
represents the small variation of damage within one cycle compared to total number of cycles and is expressed through following form: The effect of load ratio "R", for loads varying between maximum and minimum values of 1 I during a cycle is not taken into account in the present study. The sum over the cycle numbers in Equation ( 15) can be approximated by using numerical integration schemes like trapezoidal rule or Simpson's rule for definite integrals. In the present paper we assume that fatigue damage evolution rate is equivalent for each variable 1F d and 2 F d .

( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 3 1 3 1 2 1 1 1 3 1 2 1 2 1 ( , ) exp iF i I I i i d d I N I d d I I I I α α α α α α α - ∂ Ρ = ∂         =             (
1 2 F F F d d d • • • = = (18)
According to relationship between parameters this leads to link the mechanical properties 1 h ɶ and 3 h ɶ and shows that for the unified Mullins and fatigue d 3 is not an independent thermodynamics damage variable.

( )

3 1 1 F d F d = - , ( ) ( ) ( ) (
)

1 2 3 1 3 / 1 F F F b F d h a N h = - ɶ (19) ( ) N F F a N a e γ - =
The last relationship suggests here that the number of structural monomers in the chain network can decrease during cyclic fatigue loading with the increase of dangling chains. Indeed, the mass conservation written in term of monomer in the rubber matrix (Eq. 12) must not be taken strictly. Damage can injure the chain structural network.

In order to verify the possibilities of the proposed fatigue damage evolution law. The global fatigue response of the cylindrical samples under cyclic uniaxial loadings (equations 20-21) are depicted on

Figures 567for several values of ( )

1 i I α with 1 6.15, 38.76 I λ = = . 2 1 2 3 2 2 1 , 2 , 1 I I I λ λ λ λ = + = + = (20) 2 0 1 2 1 1 2 F W W S I I π λ λ λ   ∂ ∂   = = - +     ∂ ∂     (21)

Sample

A classical test sample was used for this experimental and simulation campaign. The sample is called AE2 and the mesh is presented in Figure 9 (half of the specimen). It is an axisymmetric notched sample, the radius of the notch being equal to 2 mm. It is used in the case of uniaxial and multiaxial loading conditions, because its geometry induces a multiaxial stress state in the notch, even under uniaxial loading conditions.

Monitoring and data acquisition

Both displacement (resp. angle) and force (resp. torque) can be monitored during fatigue test. A test is a series of cycle sequences. In this way, a large variety of tests can be conducted. Experimental data are recorded, such as the minimum and maximum displacements and forces, and the shape of hysteresis cycles. In this paper, the fatigue damage model is identified on idealized fatigue experimental data (Figure 10). The loading is an imposed axial displacement with constant sine amplitude. As example, the damage accumulation in AE2 specimen under tensile fatigue loadings is presented Figure 9. Experimental evolution of the hysteresis loop and GDM predictions obtained for non-relaxing uniaxial tensile experiments are presented on Figure 11. In order to overcome identification problems of the parameters overs the 985 experimental cyclic responses an equivalent cylindrical sample is used to evaluate the local material response in the notch of the specimen. Piola-Kirchhoff stress responses are then normalized. The offset of the cyclic loading, generally attributed to anisotropic and viscous effects, is also vanished since the present model cannot reproduce it.

Identification of the parameters

It must be highlighted that uniqueness of the parameters values is one of the fundamental underlying questions in inverse analysis. Genetic algorithms are used to explore the influence of parameters on the behavior responses. Figure 11

Figure 1 .

 1 Figure 1. Entanglement of the eight-chains rubber network. Eight-chains model completed by surrounding chains (in bold) is the physical motivations of the GD model.

Figure 2 .

 2 Figure 2. Experimental data (--) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.

Figure 3 .

 3 Figure 3. Experimental data (--) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.

Figure 4 .

 4 Figure 4. Experimental data (.) and GD model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible GD model response.

  rubber part under fatigue loadings is taken into account by three damage variables F d 1 , F d 2 and F d 3 . The micro damage crack growth under high cycle fatigue can be considered as a combination of damages due to the Mullins Effect and due to the cyclic variation of the loading, hence total dam-age evolution for three different material constants i h ɶ can be expressed as follows:

  16)Where i α are three functions of the maximum value of the first strain invariant 1 I obtained during the cyclic loading.

Figure 5 .Figure 6 .

 56 Figure 5. Variation of the adimensional Piola-Kirchhoff tensile stress with respect to ultimate Failure Cycles over the wide range of 1 α for 1 6.15, 38.76 I λ =

Figure 7 .

 7 Figure 7. Variation of the adimensional Piola-Kirchhoff tensile stress with respect to ultimate Failure Cycles over the wide range of 3 α for 1 6.15, 38.76 I λ =

Figure 8 .

 8 Figure 8. Tensile-torsion fatigue experiments (climatic chamber is not shown).

Figure 9 .

 9 Figure 9. AE2 sample for static and fatigue loading. The maximum value of 1M d (Mullins effect) is in the center of the sample.

Figure 10 .

 10 Figure 10. Experimental evolution of the hysteresis loop for an enforced displacement relaxing experiment on a AE2 sample. N represents the number of fatigue cycles.

  represents the tensile response of the damage model corresponding to experimental condition of Figure 10. This figure shows how the model is able to accurately reproduce the response of fatigue loadings for a maximal stretch.

Figure 11 .

 11 Figure 11. Simulations of the evolution of the averaging response of the hysteresis loop during fatigue loadings (AE2 under relaxing experiments. N represents the number of fatigue cycles.

Figure 12 Figure 12 .

 1212 Figure 12 represents the response of the identified model for fatigue cyclic extension to λ=10 and for several cycles. Even if the model cannot reproduce the hysteresis response of the material, the global pseudo-hyperelastic response can be accurately reproduced. The envelope of the maximum value of the stress for different extension are given in Figure 13 as function of the number of fatigue cycles for different values of extension. The corresponding curves directly derived from of Eq (20) and (21). The Mullins and cyclic fatigue effect of the damage model are presented (Figure 14). The above results showed the efficiency of our damage model to capture the rubber response under constant amplitude loading conditions. However, the model can also be used with variable amplitude cyclic loadings. The identified parameters for GDMF are given below: h 1 =0.02457 MPa, h 2 =1.49 10 -6 MPa and h 3 =1.1883 10 -11 , a 1 =14.01, b 1 =-1.564, a 2 =0.61209, b 2 =0.07659, a 3 =7.3392, b 3 =-0.05347,a F =2.83,

Figure 13 .

 13 Figure 13. Piola-Kirchhoff functions of number of fatigue cycles (N) for several values of the maximum first invariant values of λmax.

Figure 14 .

 14 Figure14. Model responses for tensile tests : load 1, 2 and 3 respectively from λ=1 to 3.5 (first time), 3.5 (second time) then 6.5 (first time) and 6.5 (second time). Fatigue effect of the response for cycle 50 and 500.

  Unified model for Mullins effect and high cycle fatigue life prediction of rubber materials L.Gornet & G. Marckmann GeM, UMR CNRS 6183, Ecole Centrale Nantes, 1 rue de laNoë, BP 92101-44321 Nantes, Cedex 3, France.

	R. Desmorat
	LMT-Cachan ENS Cachan/CNRS/UPMC/PRES UniverSud Paris, 61 avenue du Président Wilson, F-94230
	Cachan Cedex, France.

Table 1 :

 1 Recipe and mechanical characteristics of the material

	______________________________________________
	Formulation ______________________________________________
	NR	100.00
	Zinc Oxide	9.85
	Plastificant	3.00
	Carbon black	43.0


Ruben Herzog
email: rubenherzog@ug.uchile.cl

Maria-Jose Escobar

Adrian Palacios

Bruno Cessac

Dimensionality Reduction in spatio-temporal MaxEnt models and analysis of Retinal Ganglion Cell Spiking Activity in experiments

INTRODUCTION

Retinal spike response to stimuli is constrained, on one hand by short range correlations (receptive field overlap) and on the other hand by lateral connectivity (cells connectivity). This last effect is difficult to handle from statistics because it requires to consider spatio-temporal correlations with a time delay long enough to take into account the time of propagation along synapses. Although MaxEnt model are useful to fit optimal model (maximizing entropy) under the constraints of reproducing observed correlations, they do address spatio-temporal correlations in their classical form (Ising or higher order interactions but without time delay). Binning in such models somewhat integrates propagation effects, but in an implicit form, and increasing binning severely bias data [START_REF] Cessac | Neural Computation[END_REF]. To resolve this issue we have considered spatio-temporal MaxEnt model formerly developed e.g. by Vasquez et al. [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF]. The price to pay, however is a huge set of parameters that must be fitted to experimental data to explain the observed spiking patterns statistics. There is no a priori knowledge of which parameters are relevant and which ones are contributing to overfitting. We propose here a method of dimension reduction, i.e. a projection on a relevant subset of parameters, relying on the so-called Susceptibility matrix closely related to the Fisher information. In contrast to standard methods in information geometry though, this matrix handle space and time correlations.

We have applied this method for retina data obtained in a diurnal rodent (Octodon degus, having 30% of cones photoreceptors) and a 252-MEA system. Three types of stimuli were used: spatio-temporal uniform light, white noise and a natural movie. We show the role played by time-delayed pairwise interactions in the neural response to stimuli both for close and distant cells. Our conclusion is that, to explain the population spiking statistics we need both short-distance interactions as well as long-distance interactions, meaning that the relevant functional correlations are mediated not only by common input (i.e. receptive field overlap, electrical coupling; spillover) but also by long range connections. (b) Volume V(k, ) obtained for PSA, WN and NM stimuli. In the three cases, k c minimizing V(k, ) is represented with a square. The number of relevant dimensions k c , exceeds the firing rates meaning that spatio-temporal pairwise interactions are needed for the optimal value of V(k, ).

Remark:

The number of optimal dimensions k c required to fit the empirical distribution increases with the stimuli high-order correlations.

METHODS

Recordings: Extracellular recording of the electrical activity of retinal patches from 4 healthy, young individuals of the rodent species Octodon degus were performed ex vivo using a Multielectrode Array (USB-MEA256 from Multichannel Systems, MCS GmbH). The images were projected using a conventional DLP LED projector at 60fps and registered at 20KHz. Sorting was done using Offline Sorter R by Plexon Inc, and the spatio-temporal receptive fields were computed using STA in WN stimulus.

Analysis: A binary raster ω of N = 50 neurons with a bin size of 5ms were used, where ω n (t) = 1 if neuron n spikes at time t, 0 otherwise. The goal is to fit the empirical distribution of ω with a MaxEnt model (eq. 1) and compute its corresponding Susceptibility matrix χ (eq. 2):

H(ω) = L l=0 h l m l (ω) (1), χ ll = ∂µ [ m l ] ∂h l (2)
where m l s are the observables, h l the free parameters to fit, l is parameter label. t ranges from 0 to R -1, where the integer R ≥ 1 is called the range (i.e. memory) of the potential H, in this case, R = 2.

The Susceptibility matrix (eq. 2) tells us how much a slight change on one parameter (h l ) changes the estimation of the observables probabilities (m l ). χ depends on stimuli, increasing their entries as the stimuli complexity increases.

DIMRED: GEOMETRIC INTUITION

The set of indistinguishable distributions around the optimal set of parameters defines an ellipsoid on the parameter space (see a) with a volume V satisfying:

log V(k) = 1 2 S H (k) -log Γ k 2 + 1 + k 2 log 2π T , (3)
S H (k) = - k i=1 log λ i ( 4
)
where λ i s are decreasing χ eigenvalues; log Γ( k 2 + 1) is purely combinatoric; log 2π T characterizes the effect of precision accuracy ( ) and finite sampling. Minimizing eq. ( 1) is yields the optimal set of parameters given a finite sample size (T ) with accuracy .

(a) illustrates model composed by the Rates of two neurons (Rate Dim. 1 and 2) and the spatio-temporal correlations between both (z axis). The ellipsoid is more elongated on the sloppy eigendirections (correlations) and less elongated on the stiff eigendirections (rates). (b) shows the % of the remaining pairwise spatial (solid) and temporal (dashed) parameters as the C.I. threshold increases. Temporal interactions are more abundant than spatial ones for all stimuli at all thresholds. PSA shows less remaining pairwise parameters than dynamic stimuli for all thresholds. Then, the number of spatio-temporal parameters required to optimally fit the empirical distribution increases with the stimuli high-order correlations.

(c) shows the Receptive Field (RF) Spacing (< 1 overlap; = 1 adjacency, horizontal dashed line; > 1 no overlap) as the C.I. threshold increases for spatial (blue) and temporal (black) parameters. Lines are medians and shaded area the interquartile range. We see the same behaviour for all cases: the RF spacing slightly decreases with C.I. threshold, showing that the optimal set of parameter includes both short-range (i.e. RF overlap) and long-range (no overlap) pairwise interactions for all stimuli.

SYNTHETIC DATA: FIRING RATES ON N FIRST DIMENSIONS, CORRELATIONS ON FIRST k c Three synthetic rasters (T = 10 6 , N = 10) were generated using 3 different underlying statistics: (i) Independent, with firing rates, L = N ; (ii) Ising, firing rates and spatial correlations, L = N (N -1)/2; (iii) Pairwise R=2 (PWR2), with firing rates and spatio-temporal correlations, L = N (N -1)/2 + N 2 . For each raster a PWR2 MaxEnt model was fit. The respective χ matrices were filtered with k c = 10, 15 and 127, respectively (see eq. 1 and 2). Then, the C.I. was computed for each parameter in order to link relevant dimensions with relevant parameters.

(a) By construction the Indep. raster shows no pairwise parameters at thresholds as small as 0.03, while Ising raster, at the same threshold, shows only firing rate and spatial correlations parameters. PWR2 shows a small reduction on the % of remaining spatio-temporal parameters for k c = 127 but keeping most of them, whereas using k c = N removes all the pairwise parameters at very low thresholds. The first N dimensions contains mainly information about firing rates, while the next k c -N concentrates information about spatio-temporal correlations.

CONCLUSIONS (i)

The first cut-off of χ spectrum contains mainly information about individual neuron firing rates, while the second cut-off contains information about spatio-temporal correlations. If k c = N , there is no relevant information about spatio-temporal correlations on the data and it can be optimally explained by a model which takes into account only neuron firing rates. (ii) All tested stimuli requires spatio-temporal interactions to optimally fit the evoked retinal population activity. So, as other works have shown on different species (salamander [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF][START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Tkacik | Searching for collective behavior in a large network of sensory neurons[END_REF]; monkey [START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF]; guinea pig [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF]), retinal population activity presents significant spatio-temporal interactions. Moreover, the stimuli high-order correlations increases the number of relevant spatio-temporal interactions on the network. (iii) The optimal set of pairwise parameters, which implies a pair of neurons, includes neurons with overlapped receptive fields as well as distant neurons for all tested stimuli, with small differences between stimuli. Short range correlations can be explained via stimulus correlation and receptive field overlap whether long range more presumably reflect lateral connectivity effects.

χ

  SPECTRAL PROPERTIES UNDER DIFFERENT STIMULI (a) The χ eigenvalue spectrum obtained for PSA, WN and NM stimuli. The graphs shows two cut-offs: The first cut-off, close to the network size N represented as a vertical dashed line, is mainly to the effect of individual firing rates. The second cut-off is associated to the information contained on spatio-temporal pairwise interactions. Minimizing eq. (3), we obtain different values of k c which are represented as squares on each curve.

  Stimuli: (i) Photopic spontaneous activity (PSA): A uniform spatiotemporal invariant full field (15 mins); (ii) White noise (WN): Binary checkerboard pattern (20 mins); and, (iii) Natural Movie (NM): recorded in the animal's natural environment (30s movie, 40 repetitions).PSA WN NM

  (b) shows the value of the volume, V(k, ) as the degrees of freedom (k) increases. The volume reaches a minimum (dots) for certain k named as k c . Increasing k upon k c makes the volume explode. FROM χ TO NETWORK PROPERTIES The effect on the estimation of slight changes in the lth parameter is computed as k=K k=1 λ k v 2 lk , where v lk is the lth entry of the eigenvector v k and λ k is corresponding eigenvalue. The Change Index (C.I. equation on a) is the ratio of the effect before (K = L) and after filtering (K = k c ). C.I.=1 means no filtering effect and C.I.=0 means full filtering effect. C.I. summarizes the effect of filtering on the direct influence of h l on m l estimation and all the indirect influences (h l acts on m l which modifies m l ) (a) shows the distribution of C.I. values for the pairwise parameters, where PSA shows two modes, one close to 0 and other to 1. WN and NM shows one big mode close to 1. Thus, static stimulus (PSA) shows much more pairwise parameters affected by filtering than dynamic ones.

( b )

 b The monomials probabilities show a non-linear relationship with C.I.. Monomials with low probabilities are likely to show low C.I. values. The Indep. raster shows only firing rates parameters with C.I. ∼ 1; Ising raster shows firing rates and spatial correlations parameters with C.I. ∼ 1, but none related to temporal correlations shows high C.I.; Almost all PWR2 parameters show C.I.∼ 1 when filtering with k c , but filtering with k c = N , reduces the C.I. of most of the spatio-temporal parameters at values ∼ 0.2, while most firing rates parameters remain unchanged.

ACKNOWLEDGEMENTS

This work has been supported by CONICYT-Fondecyt 1140403 and 1150638, CONICYT-Basal Project FB0008, Millennium Institute ICM-P09-022-F, ECOS-Conicyt C13E06. Special thanks to Dr. Patricio Orio for the access to the computational cluster of the CNSLab at the CINV. ENAS software (https://enas.inria.fr/) was used for the computation of χ.


Ivan Gentil
email: gentil@ceremade.dauphine.fr

Bruno R Émillard
email: bruno.remillard@hec.ca

Carlo Solutions

O F Feynman-Kac Equations

)

Keywords: Feynman-Kac formulae, sequential Monte Carlo, genetic algorithms, systematic sampling, Markov chains. 2000 Mathematics Subject Classification: Primary 65C35, 60G35 Secondary 60J10

While convergence properties of many sampling selection methods can be proven to hold in a context of approximation of Feynman-Kac solutions using sequential Monte Carlo simulations, there is one particular sampling selection method introduced by [2], closely related with "systematic sampling" in statistics, that has been exclusively treated on an empirical basis. The main motivation of the paper is to start to study formally its convergence properties, since in practice it is by far the fastest selection method available. One will show that convergence results for the systematic sampling selection method are related to properties of peculiar Markov chains.

Introduction

Let (X k ) k≥0 be a non-homogeneous Markov chain on a locally compact metric space E, with transition kernels (K n ) n≥1 and initial law η 0 defined on the Borel σ-field B(E).

Further let B b (E) be the set of bounded B(E)-measurable functions.

Given a sequence (g n ) n≥1 of positive functions in B b (E), suppose that one wants to calculate recursively the following Feynman-Kac formulae (η n ) n≥1 :

η n (f ) = γ n (f ) γ n (1) , f ∈ B b (E), (1)
where

γ n (f ) = E f (X n ) n k=1 g k (X k-1 ) . (2)
Note that most nonlinear filtering problems are particular cases of Feynman-Kac formulae.

Following [START_REF] Crisan | Discrete filtering using branching and interacting particle systems[END_REF] and [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF], let M 1 (E) denotes the set of probability measures on (E, B(E)).

If µ ∈ M 1 (E) and n ≥ 0, let µK n be the probability measure defined on B b (E) by

µK n (f ) = µ(K n f ) = E E f (z)K n (x, dz)µ(dx).
In order to understand the relation between the η n s, for any n ≥ 1, let ψ n : M 1 (E) → M 1 (E) be defined by

ψ n (η)f = η(g n f ) η(g n ) , η ∈ M (E), f ∈ B b (E),
and let Φ n denote the mapping from M 1 (E) to M 1 (E) defined by

Φ n (η) = ψ n (η)K n .
Then it is easy to check that for any n ≥ 1,

η n = Φ n (η n-1 ). (3)
Note also that for any n ≥ 1, the mapping Φ n can be decomposed into ηn = ψ n+1 (η n ),

η n+1 = ηn K n+1 , n ≥ 0, η 0 ∈ M 1 (E). ( 4
)
Further remark that the first transformation, η n → ηn , is non-linear, while the second one, ηn → η n+1 , is linear.

Even if the forward system of equations ( 3) looks simple, it can rarely be solved analytically, and even if this is the case, it would require extensive calculations. This is why algorithms for approximating (η n ) n∈N , starting from η 0 , are so important.

One such method, presented in the remarkable surveys [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF][START_REF] Crisan | A survey of convergence results on particle filtering methods for practitioners[END_REF] and the book of [START_REF] Del Moral | Feynman-Kac formulae[END_REF],

is to build approximations of measures (η n ) n∈N using interacting particle systems. The algorithm uses decomposition [START_REF] Cappé | Comparison of resampling schemes for particle filtering[END_REF], and by analogy with genetics, the first step, which is related to a sampling selection method, is often referred to as the selection step, and the second one is termed the mutation step, while in reality it is a Markovian evolution of the particles. The speed of any algorithm depends on the two steps. In [START_REF] Pitt | Filtering via simulation: auxiliary particle filters[END_REF] the authors focus on the mutation step, but it is clear that it also depends on the sampling selection process.

In this paper, one discusses properties of a particular algorithm that is called "systematic sampling" selection herein, while in the genetic algorithms literature, it has been strangely called "Stochastic universal sampling" selection. It seems to have appeared first in [START_REF] Baker | Reducing bias and inefficiency in the selection algorithm[END_REF]. It has been reintroduced in the filtering literature in [START_REF] Carpenter | Building robust simulation-based filters for evolving data sets[END_REF], see also [START_REF] Künsch | Recursive Monte Carlo Filters: Algorithms and Theoretical Analysis[END_REF].

In what follows, a description of the general algorithm is given in Section 2, with a few examples of sampling selection methods, together with some tools for studying its convergence. In Section 3, one focuses on the systematic sampling selection method, giving some properties, and stating some convergence results and a conjecture, based on results from Markov chains proved in the appendix. Finally, in Section 4, numerical comparisons between sampling selection methods are made through a simple model of nonlinear filtering for noisy black-and-white images.

Algorithm and sampling selection methods

The general algorithm for approximating the solution of ( 3) is first given, following the exposition in [START_REF] Crisan | Discrete filtering using branching and interacting particle systems[END_REF][START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF], while particular sampling selection methods are presented next. Throughout the rest of the paper, it is assumed that for any n ≥ 1,

0 < inf x∈E g n (x) ≤ sup x∈E g n (x) < +∞.

General algorithm

Let N be a integer, representing the number of particles and for any n ≥ 0, let

ξ n = ξ 1 n , • • • , ξ N n
denotes the particles at time n and set

η N n = 1 N N i=1 δ ξ i n .
• At time n = 0, the initial particle system

ξ 0 = ξ 1 0 , • • • , ξ N 0 consist of N
independent and identically distributed particles with common law η 0 .

• For each n ≥ 1, the particle system

ξ n = ξ 1 n , • • • , ξ N n consists of N particles, is
obtained in the following way:

(Sampling/Selection) First calculate the weights vector W n ∈ (0, 1) N , where

W i n = g n (ξ i n-1 ) N i=1 g n (ξ i n-1 ) , i = 1, . . . , N. (5)
Then, select, according to a given sampling selection method, a sample ξn-1 = ξ1 n-1 , . . . , ξN n-1 of size N from ξ n-1 and with weights (W i n ).

(Evolution/Mutation) Given ξn-1 , the new particle system ξ n consists of particles ξ i n chosen independently from law K n ( ξi n-1 , dx), 1 ≤ i ≤ N . In other words, for any z = (z 1 , . . . , z N ) ∈ E N ,

P ξ n ∈ dx| ξn-1 = z = N i=1 K n (z i , dx i ).
Note that in order to describe a sampling selection method, it suffices to define how the numbers M 1 n , . . . , M N n ∈ {0, 1, . . . , N } are randomly selected, with M i n representing the number of times particle ξ i n-1 is appears in the new sample. Therefore, one can write

ηN n-1 = 1 N N i=1 δ ξi n-1 = 1 N N i=1 M i n δ ξ i n-1 .
A sampling selection method will be said to be conditionally unbiased, if for any

i ∈ {1, . . . , N } and any k ≥ 1, E(M i k |ξ k-1 ) = N W i k .
Remark 2.1. Conditional unbiasedness yields the following property:

E η N k f |ξ k-1 = Φ k η N k-1 (f ), f ∈ B b (E). (6)
For, in that case,

E η N n f |ξ n-1 = E E η N n f |ξ n-1 , ξn-1 |ξ n-1 = E 1 N N i=1 M i n K n f (ξ i n-1 ) ξ n-1 = N i=1 W i n K n f (ξ i n-1 ) = Φ n η N n-1 (f ).
The mean square error of a particular sampling selection method can be obtained using the following useful result, obtained by [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF]Theorem 2.36]. Before stating the result, define, for any measurable η with values in M (E),

η 2 2 = sup f ∈B b (E), f ∞≤1 E (ηf ) 2 .
Theorem 2.1. Assume that the sampling selection method is conditionally unbiased and that the following condition is verified for all 1 ≤ k ≤ n: there exists a constant

C k such that for all N -dimensional vectors q 1 , • • • , q N ∈ R N , E   1 N N i=1 (M i k -N W i k )q i 2 ξ k-1   ≤ 1 N C k max 1≤i≤N q i 2 . (7)
Then, for all 1 ≤ k ≤ n, there exists a constant C k such that

η N k -η k 2 2 ≤ C k /N.
In what follows, only conditionally unbiased sampling selection methods are considered. As shown in Remark 3.2, one can see that in general, the systematic sampling selection method defined below does not satisfy condition [START_REF] Crisan | Discrete filtering using branching and interacting particle systems[END_REF] of the previous theorem, while classical sampling selection methods, like the ones listed in Section 2.3, do satisfy it. Therefore, weaker conditions must be imposed in order to obtain mean square convergence. In fact, one has the following result.

Theorem 2.2. Let (a N ) be a sequence such that a N /N → 0, as N → ∞. Assume that the sampling selection method is conditionally unbiased. Then

lim N →∞ a N max 1≤k≤n η N k -η k 2 2 = 0 if and only if, for any f ∈ B b (E), lim N →∞ a N max 1≤k≤n E   1 N N i=1 (M i k -N W i k )f (ξ i k-1 ) 2   = 0. ( 8
)
Moreover, sup N ≥1 a N max 1≤k≤n η N k -η k 2 2 is finite if and only if sup N ≥1 a N max 1≤k≤n E   1 N N i=1 (M i k -N W i k )f (ξ i k-1 ) 2   < ∞.
Proof. Suppose that a N /N → 0 and let f ∈ B b (E) be given. First, note that using the unbiasedness condition, together with (6), one has, for any k ∈ {1, . . . , n},

E (η N k f -η k f ) 2 = E (η N k f -Φ k (η N k-1 )f ) 2 + E (Φ k (η N k-1 )f -η k f ) 2 . ( 9
)
Since g k ≥ c k > 0 by hypothesis, for some positive constant c k , k ≥ 1, it follows that lim

N →∞ a N max 1≤k≤n η N k -η k 2 2 = 0
if and only if for any k = 1, . . . , n, lim

N →∞ a N E (η N k f -Φ k (η N k-1 )f ) 2 = 0. Next, it can be shown easily that for any k = 1, . . . , n, E η N k f -Φ n (η k-1 )f 2 ξ k-1 can be written as E   1 N N i=1 (M i k -N W i k )K k f (ξ i k-1 ) 2 ξ k-1   + 1 N Ψ k (η N k-1 )(K k f 2 -(K k f ) 2 ). Since K k f ∈ B b (E), 0 ≤ 1 N Φ k (η N k-1 )(K k f 2 -(K k f ) 2 ) ≤ 1 N f 2 ∞
, and a N /N → 0, it follows from the calculations above that lim

N →∞ a N max 1≤k≤n E (η N k f -Φ k (η N k-1 )f ) 2 = 0
if and only if (8) holds true. The rest of the proof is similar, so it is omitted.

Systematic sampling

By obvious analogy with systematic sampling in Statistics, the first sampling selection method that is described is simply called "systematic sampling". It appears that this method was first proposed by [START_REF] Baker | Reducing bias and inefficiency in the selection algorithm[END_REF] under the strange name "Stochastic Universal Sampling", in a context of unbiased sampling selection for genetic algorithms. However, nobody formally studied its convergence properties.

As opposed to the definition of [START_REF] Baker | Reducing bias and inefficiency in the selection algorithm[END_REF], the sampling selection method can simply be defined in the following way: For n ≥ 1, let U n a uniform random variable on [0, 1) and note for w ∈ [0, 1], M (w, U n ) := N w + U n , where x denotes the integer part of x. Then

M 1 n := M (W 1 n , U n ), M k n := M (W 1 n + • • • + W k n , U n ) -M (W 1 n + • • • + W k-1 n , U n ), k = 2, . . . , N. Since M (1, U n ) = N , one gets that N i=1 M i n = N .
Therefore the number of particles is always N . Properties of that sampling selection method are examined in Section 3.

Other sampling methods

One can grossly classify the various sampling selection methods into two categories, according as the number of particles is constant or random. The following list is by no means exhaustive. One concentrates on "classical" selection methods. For the first two methods, N is constant, while N n fluctuates in the last two methods. For other sampling selection methods, one may consult [START_REF] Crisan | Discrete filtering using branching and interacting particle systems[END_REF][START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF][START_REF] Del Moral | Feynman-Kac formulae[END_REF] and references therein. See also [START_REF] Cappé | Comparison of resampling schemes for particle filtering[END_REF] for a recent review of comparison of various selections. Note that the last two methods are particular cases of what is known as "Branching selection methods" in the filtering literature. However, there are many others sampling methods, e.g., the tree-sampling method proposed in [START_REF] Crisan | Discrete filtering using branching and interacting particle systems[END_REF]. The latter is more computationally demanding, and since it is not as simple to implement as the other methods, it is not considered here.

Simple random sampling

This selection method is based on simple random sampling without rejection. It follows that

M 1 n , • • • , M N n ∼ Multinomial N, W 1 n , • • • , W N n ,
where W i n 1≤i≤N are given by ( 5). This sampling selection method is the classical selection method and has many interesting properties that have been studied mainly by Del Moral and co-authors, e.g. see [START_REF] Del Moral | Feynman-Kac formulae[END_REF]. In particular, conditions (i)-(ii) of Theorem 2.1 are met; also one can prove a Central Limit Theorem and Large Deviations Properties.

The remainder stochastic sampling

This algorithm was first introduced by [START_REF] Brindle | Genetic algorithms for function optimization[END_REF] in a context of unbiased sampling selection for genetic algorithms; see also [START_REF] Baker | Adaptative selection methods for genetic algorithms[END_REF][START_REF] Baker | Reducing bias and inefficiency in the selection algorithm[END_REF] for comparisons between sampling selection methods in the latter context. It is also defined as "Residual sampling" by [START_REF] Liu | Sequential Monte Carlo methods for dynamic systems[END_REF]. It is a little much faster to implement than the simple random sampling selection method, it satisfies conditions (i)-(ii) of Theorem 2.1, and recently, [START_REF] Douc | Limit theorems for weighted samples with applications to sequential monte carlo methods[END_REF] investigated some of its convergence properties. See also [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF] and the references therein. To describe the selection method, first define Ñ = N -

N i=1 N W i n = N i=1 N W i n ,
where {x} stands for the fractional part of x, i.e. {x} = xx . Next, allocate the (possibly) remaining Ñ particles via simple random sampling, i.e.

M 1 n -W 1 n , • • • , M N n -W N n ∼ Multinomial Ñ, W 1 n , • • • , W N n , with W i n = N W i n / N j=1 N W j n , 1 ≤ i ≤ N .
In fact Ñ is expected to be the half of N , then the speed of the selection method is almost the same as for the simple random sampling.

Binomial sampling

As stated before, for this sampling selection method and the next one, the number of particles at time n is random and it is denoted by N n , n ≥ 0.

Of course, N 0 is fixed. For n ≥ 1, and given ξ n-1 and

N n-1 , M 1 n , . . . , M Nn-1 n are independent and M i n ∼ Bin(N n-1 , W i n ), for 1 = 1, . . . , N n-1 . It follows that N n = Nn-1 i=1 M i n .
This sampling selection method is a little bit faster than the simple random sampling selection method, but a major drawback is that there is no control on the number of particles. Moreover, P (N n = 0) > 0.

Bernoulli sampling

The Bernoulli sampling selection selection method was introduced in [START_REF] Crisan | Convergence of a branching particle method to the solution of the Zakai equation[END_REF]. See also [START_REF] Crisan | A particle approximation of the solution of the Kushner-Stratonovitch equation[END_REF][START_REF] Crisan | Exact rates of convergence for a branching particle approximation to the solution of the Zakai equation[END_REF] for additional properties of the sampling selection selection. It is worth noting that M i n takes the same values as in the systematic sampling selection method, provided N n-1 = N . In fact, for n ≥ 1, and given ξ n-1

and N n-1 , M 1 n , . . . , M
Nn-1 n are independent, where M i n is defined by

M i n = N n-1 W i n + ε i n , ε i n ∼ Ber( N n W i n ), 1 ≤ i ≤ N n-1 .
Note that N n ≥ 1 and that the following alternative representation also holds:

M i n = N (W 1 n + • • • + W i n ) + U i n -N (W 1 n + • • • + W i-1 n ) + U i n ,
where U 1 n , . . . , U

Nn-1 n are independent and U i n ∼ Unif([0, 1)), given ξ n-1 , N n-1 .

Some properties and results for systematic sampling selection

Throughout the rest of paper, the selection method is the one defined in Section 2.2. Let's start first with some elementary properties of systematic sampling selection.

Lemma 3.1. Suppose that U n is uniformly distributed over [0, 1). Then, conditionally

on ξ n-1 , one has, for any i ∈ {1, • • • , N }, M i n -N W i n ∼ Ber( N W i n ). ( 10
)
In particular, for any

i ∈ {1, • • • , N }, E M i n |ξ n-1 = N W i n .
Proof. It suffices to show that whenever U ∼ Unif([0, 1)) and x, y ≥ 0, then

U + x + y -U + x -y is a Bernoulli random variable with parameter p = {y}.
To this end, first note that V = {U + x} is also uniformly distributed on [0, 1). Next,

U + x + y -U + x -y = {U + x} + {y} = V + {y} ∼ Ber({y}),
hence the result holds.

Remark 3.1. Using the same proof as in Lemma 3.1, then, conditionally on ξ n-1 , one

obtains M i n + • • • + M j n -N (W i n + • • • + W j n ) ∼ Ber N (W i n + • • • + W j n ) , for any i ≤ j ∈ {1, • • • , N }.
Note also that since the sampling selection method is unbiased, i.e. condition (i) of Theorem 2.1 is satisfied, then for any n ≥ 1, one has formula [START_REF] Crisan | Exact rates of convergence for a branching particle approximation to the solution of the Zakai equation[END_REF].

To obtain L 2 convergence of the algorithm based on the systematic sampling selection method, one would like to apply Theorem 2.1 of [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF]. All sampling selections presented in Section 2.3 satisfies property [START_REF] Crisan | A survey of convergence results on particle filtering methods for practitioners[END_REF]. If N n is random, there is an similar condition to [START_REF] Crisan | A survey of convergence results on particle filtering methods for practitioners[END_REF]. But as shown next, systematic sampling behaves differently. Remark 3.2. Inequality ( 7) is not verified in general for the systematic sampling selection method. Here is an illustration. Suppose that N = 2m and let, for any

i ∈ {1, • • • , N/2}, W 2i n = 1/(2N ), and W 2i-1 n = 3/(2N ).
Then one can check that for

any 1 ≤ i ≤ N/2, M 2i-1 n = 1 if U n ∈ [0, 1/2), M 2i-1 n = 2 if U n ∈ [1/2, 1), M 2i n = 1 if U n ∈ [0, 1/2), M 2i n = 0 if U n ∈ [1/2, 1). Next, if 1 ≤ i ≤ N/2, set q 2i = 1 and q 2i-1 = -1. It follows that E   1 N N i=1 (M i k -N W i k )q i 2 ξ k-1   = 1 4 ,
showing that inequality ( 7) is false.

The description of all selections in Section 2.2 and 2.3 prove that the systematic sampling selection method if the less computationally demanding. At the light of Remark 3.2 one can see that the L 2 -convergence of the particle system with the systematic sampling selection method to the Feynman-Kac formulae is not easy. However one believes that the following holds true.

Conjecture 3.1. Suppose that η 0 and (K n ) n≥1 are absolutely continuous laws and consider that M 1 n , . . . , M N n are obtained using the the systematic sampling selection method. Then, for all f ∈ B b (E) and n ≥ 1, (8) holds with a N ≡ 1, i.e.

lim N →∞ E   1 N N i=1 (M i n -N W i n )f (ξ i n-1 ) 2   = 0. (11)
Note that it follows from Theorem 2.2 that the above conjecture is equivalent to

η N n -η n 2 → 0, as N → ∞, for any n ≥ 0.
In what follows, one tries to motivate why Conjecture 3.1 might be true. To this end, first note that any 1 ≤ k ≤ N ,

M k n -N W k n = N (W 1 n + • • • + W k-1 n ) + U n -N (W 1 n + • • • + W k n ) + U n . Now, set F 0 n = 0 and F k n = k j=1 g n (ξ j n-1 ), 1 ≤ k ≤ N .
For any α > 0 and any

f ∈ B b (E), further define Z N n (f, α) = 1 √ N N k=1 f ξ k n-1 F k-1 n α + U n - F k n α + U n = 1 √ N N k=1 f ξ k n-1 S k-1 n -S k n , where S k n = F k n /α + U n and S 0 n = U n . Then, setting ḡn = 1 N N k=1 g n (ξ k n ) and defining Y N n (f ) = Z N n (f, ḡn ), one has Y N n (f ) = 1 √ N N i=1 (M i n -N W i n )f (ξ i n-1
), so one can rewrite [START_REF] Del Moral | Feynman-Kac formulae[END_REF] in the form lim

N →∞ 1 N E Y N n (f ) 2 = 0.
Unfortunately, working with Y N n appears to be impossibly difficult, so one could work instead with a more tractable quantity, namely Z N n . In the case n = 1, one has a least the following result, which is a first step in proving Conjecture 3.1. Before stating it, recall that D([0, 1]) is the space of càdlàg functions with Skorohod's topology. Theorem 3.2. Assume that the law of g 1 (ξ 1 0 ) is absolutely continuous. Then, for any α > 0 and any f ∈ B b (E), the sequence of processes B N ∈ D([0, 1]) defined by

B N f,α (t) = Z N t 1 (f, α), t ∈ [0, 1],
converges in D([0, 1]) to σB f,α , where B f,α is a Brownian motion and

lim N →∞ E Z N 1 (f, α) 2 = σ 2 .
The proof on Theorem 3.2 is an easy consequence of Theorem A.2 applied with

X k = f (ξ k 0 ), Y k = g 1 (ξ k 0
) and f (x, y, s) = x(y -s). In addition, there is an "explicit" expression for σ 2 . More details can be found in Appendix A.

Remark 3.3. Theorem 3.2 does not prove Conjecture 3.1 in the case n = 1. However, if one is willing to deal with a random number of particles at step n = 1, one obtains the following interesting result: Set N 1 = U 1 + N j=1 g 1 (ξ k 0 )/α , and define

ηN1 0 = 1 N 1 N k=1 M k 1 δ ξ k 0 .
Then, as N → ∞, N 1 /N tends to η 0 (g 1 )/α and lim sup

N →∞ N η N1 1 -η 1 2
2 < ∞. To prove convergence for higher orders, i.e n > 1, one would need results from nonhomogeneous Markov chains. The approach will be examined in a near future, using for example the results of [START_REF] Sethuraman | A martingale proof of Dobrushin's theorem for non-homogeneous Markov chains[END_REF].

Remark 3.4. In order to keep N 1 fixed, one could try to control the term

Z N 1 (g 1 ) - Z N 1 (α). Since √ N (g 1 -α) η 0 (g 2 1 )
Z, where Z ∼ N (0, 1), it follows that

√ N 1 g 1 - 1 α = - η 0 (g 2 1 ) η 2 0 (g 1 ) Z + o P (1).
If one could differentiate term by term, one would then obtain

Z N 1 (g 1 ) -Z N 1 (α) η 0 (K 1 f g 1 ) η 0 (g 2 1 ) η 0 (g 1 ) 2 Z = η 1 (f ) η 0 (g 2 1 ) η 0 (g 1 ) Z, so one could guess that Y N η 1 (f ) η 0 (g 2 1 ) η 0 (g 1 ) Z + B f,α (1)
. On the other hand, if the sequence Z N 1 (α) was tight for α in a closed interval not containing zero, then one would get Z N 1 (g 1 ) -Z N 1 (α) → 0 in probability. There is no indication so far in favor of one of these two approaches.

Numerical comparisons

The numerical comparisons will be done through a simple model of filtering for tracking a moving target using noisy black-and-white images, where the exact filter can be calculated explicitly, that is, η n is known for any n ≥ 1, e.g. [START_REF] Gentil | Filtering of images for detecting multiple targets trajectories[END_REF].

Description of the model

One will assume that the target moves on Z 2 according to a Markov chain. Observations consist in black-and-white noisy images of a finite fixed region R ⊂ Z 2 . More precisely, let (X n ) n≥0 be a homogeneous Markovian chain with values in X = {ω ∈ {0, 1} Z : x∈Z ω(x) = 1}. Of course, the position of the target at step n is x 0 if and

only if X n (x 0 ) = 1. Set M (a, b) = P {X n+1 (a) = 1|X n (b) = 1} , a, b ∈ Z 2 . ( 12
)
Note that M describes exactly the movement of the target.

The model for observations Y k ∈ {0, 1} R , k = 1, . . . , n is the following: Given X 0 , . . . , X n , assume that {Y n (x)} x∈R are independent and for any x ∈ R,

P (Y n (x) = 0|X n (x) = 0) = p 0 , P (Y n (x) = 1|X n (x) = 1) = p 1 , (13)
where 0 < p 0 , p 1 < 1. For a more realistic result one can assume that 1/2 < p 0 , p 1 < 1. One wants to compute the distribution of X k conditionally to Y n , where Y n is the sigma-algebra generated by observations Y 1 , . . . , Y n , and Y 0 is the trivial sigmaalgebra. As in Section 2 of [START_REF] Gentil | Filtering of images for detecting multiple targets trajectories[END_REF], note that for any (ω, ω ) ∈ {0, 1} R ⊗ X, the conditional probability

P (Y k = ω|X k = ω ) = Λ(ω, ω ) satisfies Λ(ω, ω ) = p |R|-1 0 (1 -p 1 ) 1 -p 0 p 0 <ω> p 0 p 1 (1 -p 0 )(1 -p 1 ) <ωω >
, where < ω >= x∈R ω(x) and < ωω >= x∈R ω(x)ω (x).

Let P be the joint law of the Markovian targets with initial distribution ν, and the observations, and let Q be the joint law of the Markovian targets with initial distribution ν, and independent Bernoulli observations with mean 1/2. Further let G n be the sigma-algebra generated by Y 1 , . . . , Y n , X 0 , . . . , X n . Then it is easy to check that with respect to G n , P is equivalent to Q and

dP dQ Gn = n j=1 2 |R| Λ(Y j , X j ). Further define L n = n j=1 Λ(Y j , X j ).
Denoting by E P (resp. E Q ) expectation with respect to P (resp. Q), observe that for any f ∈ B b (X), one has

ηn (f ) = E P (f (X n )|Y n ) = E Q (f (X n )L n |Y n ) E Q (L n |Y n ) . ( 14
)
This formula is a consequence of the properties of conditional expectations, and in the context of filtering, ( 14) is known as the Kallianpur-Stribel formula, e.g. [START_REF] Kallianpur | Estimation of stochastic systems: Arbitrary system process with additive white noise observation errors[END_REF].

Denote by K the Markov kernel associated with the Markov chain (X n ) n≥0 defined by M , as in [START_REF] Del Moral | Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF]. One can check that η n and ηn satisfy (4) with g n (x) = Λ(y n , x) and K n = K. Note also that in that case, g n takes only two values which can be assumed to belong to Q because of rounding errors. It follows from Remark A.2 that

sup N0≥1 E N 1 η N1 1 -η 1 2 < ∞.
The results proved in Section 2 of [START_REF] Gentil | Filtering of images for detecting multiple targets trajectories[END_REF] provide an algorithm for computing recursively the exact filter, i.e. the law of X n given Y n . In the next section, one will compare the results from the exact filter with those obtained by the Monte Carlo algorithm described in Section 2 with various sampling methods.

Simulation results

In what follows, R is chosen to be the window of size 100 × 100 defined by R = {0, . . . 99} 2 . To makes things simple, the target starts at (50, 50) and it moves according to a simple symmetric random walk, i.e its goes up, down, right or left to the nearest neighbor with probability 1/4. The estimation of the position of the target is taken to where

S i = S i-1 + Y i , i ≥ 1.
For n ∈ Z, set e n (s) = e 2πins , s ∈ [0, 2), and let ζ n = E (e n (Y )). Further set

N = {n ∈ Z; ζ n = 1}.
Recall that (e n ) ∈Z is a complete orthonormal basis of the Hilbert space H = L 2 ([0, 1), λ) with scalar product (f, g) = 1 0 f (s)ḡ(s)ds and norm

f 2 = (f, f ).
It is easy to check that (Z i ) i≥0 is a Markov chain on Z with kernel K defined by

Kf (x, s) := Z f (x , {s + y})P (dx , dy), f ∈ B b (Z), (15)
and stationary distribution µ = P X ⊗ λ. Note that for any f ∈ L 2 (µ), by Tonelli's theorem, Kf is well defined, it depends only on s ∈ [0, 1), and it belongs to H since

Kf 2 2 ≤ 1 0 Z f 2 (x, {s + y})P (dx, dy)ds = Z 1 0 f 2 (x, u)duP (dx, dy) = Z f 2 (z)µ(dz) = f 2 L 2 (µ) .
Finally, let L and A be the linear bounded operators from L 2 (µ) to H defined by

Lf (s) = n∈N (Kf, e n )e n (s), Af (s) = R f (x, s)P X (dx), s ∈ [0, 1). Theorem A.1. Let f ∈ L 2 (µ) be given and set W N = 1 N N k=1 f (Z k ). Then: (i) If the initial distribution of Z 0 = (X 0 , S 0 ) is µ, then W N converges almost
surely and in mean square to W given by

W = Lf (S 0 ) = n∈N (Kf, e n )e n (S 0 ). (16)
If, in addition,

n∈Z\N |(Kf, e n )| |(Af, e n )| |1 -ζ n | < ∞, (17)
then

N E (W N -W) 2 converges, as N → ∞, to f 2 L 2 (µ) -Lf 2 2 + 2 n∈Z\N (Kf, e n )(Af, e n ) 1 -ζ n . ( 18
) (ii) If the initial distribution of Z 0 = (X 0 , S 0 ) is µ, if N = {0} and n∈Z\{0} |(Kf, e n )| 2 |1 -ζ n | 2 < ∞, (19)
then the sequence of processes B N , defined by ]) to σB, where B is a Brownian motion and σ 2 is given by [START_REF] Liu | Sequential Monte Carlo methods for dynamic systems[END_REF].

B N (t) = √ N W N t -µ(f ) , t ∈ [0, 1], converges in D([0, 1
(iii) If P Y admits a square integrable density h, then the Markov chain is geometrically ergodic, that is, there exists ρ ∈ (0, 1) such that for any f ∈ L 2 (µ),

|K n f (Z 0 ) -µ(f )| ≤ h 2 ρ n-2 f L 2 (µ) , n ≥ 2.
Proof. For simplicity, set ψ = Kf ∈ H. To prove (i), start the Markov chain from µ and denote the law of the chain by Q. Then the sequence (Z n ) n≥0 is stationary, and Birkhoff's ergodic theorem, e.g. [14, Section 6.2], can be invoked to claim that W N converges almost surely and in mean square to some random variable W. To show that W is indeed given by ( 16), it suffices to show that E (W N -W) 2 tends to 0, as

N → ∞. First, note that E(W 2 ) = Lf 2 2 = n∈N |(ψ, e n )| 2 . Next, set ϕ(s) = Af (s).
If n ∈ N , then e n (Y ) = 1 P-a.s., and it follows, by Fubini's theorem, that (ψ, e n ) =

e k ({s + y})P Y (dy) = ζ k e k (s), s ∈ [0, 1).
Hence, for any k ≥ 1 and any n ∈ Z, one obtains

K k e n = ζ k n e n . (20)
Now, using the Markov property of the chain together with [START_REF] Sethuraman | A martingale proof of Dobrushin's theorem for non-homogeneous Markov chains[END_REF], one has

E(W N W) = E(W N W) = 1 N N k=1 n∈N (ψ, e n )E f (Z k )e n (S 0 ) = 1 N N k=1 n∈N (ψ, e n )E K k f (S 0 )e n (S 0 ) = 1 N N k=1 n∈N (ψ, e n ) K k-1 ψ, e n = 1 N j∈Z N k=1 n∈N (ψ, e j )(ψ, e n ) K k-1 e j , e n = 1 N j∈Z N k=1 n∈N (ψ, e j )(ψ, e n )ζ k-1 j (e j , e n ) = 1 N N k=1 n∈N |(ψ, e n )| 2 ζ k-1 n = E(W 2 ),
since, by definition, ζ n = 1, for any n ∈ N .

Next, using stationarity, the Markov property, [START_REF] Sethuraman | A martingale proof of Dobrushin's theorem for non-homogeneous Markov chains[END_REF], and also using identity

2 N 2 N -1 k=1 k j=1 z j-1 = N -1 1 -z - z -z N (1 -z) 2 , z ∈ C, z = 1, it follows that E(W 2 N ) = 1 N E f 2 (Z 0 ) + 2 N 2 N -1 k=1 k j=1 E K j-1 ψ(S 0 )f (Z 0 ) = 1 N f 2 L 2 (µ) + 2 N 2 N -1 k=1 k j=1 K j-1 ψ, ϕ = 1 N f 2 L 2 (µ) + n∈Z (ψ, e n )(ϕ, e n )   2 N 2 N -1 k=1 k j=1 ζ j-1 n   = 1 N f 2 L 2 (µ) + N -1 N n∈N (ψ, e n )(ϕ, e n ) + 2 N 2 n∈Z\N (ψ, e n )(ϕ, e n ) N -1 1 -ζ n - ζ n -ζ N n (1 -ζ n ) 2 = 1 N f 2 L 2 (µ) + N -1 N E(W 2 ) + 2 N 2 n∈Z\N (ψ, e n )(ϕ, e n ) N -1 1 -ζ n - ζ n -ζ N n (1 -ζ n ) 2 .
Collecting the expressions obtained for E(W 2 N ) and E(W N W), one gets

E (W N -W) 2 = 1 N f 2 L 2 (µ) - 1 N E(W 2 ) (21) + 2 N 2 n∈Z\N (ψ, e n )(ϕ, e n ) N -1 1 -ζ n - ζ n -ζ N n (1 -ζ n ) 2 . Since n∈Z\N |(ψ, e n )||(ϕ, e n )| is finite, sup nZ\N N -1 1 -ζ n - ζ n -ζ N n (1 -ζ n ) 2 = N -1 k=1 k j=1 ζ j-1 n ≤ N 2 2 ,
it follows from (21) and the Dominated Convergence Theorem that lim

N →∞ E (W N -W) 2 = 0
and under the additional condition [START_REF] Künsch | Recursive Monte Carlo Filters: Algorithms and Theoretical Analysis[END_REF], one also obtains lim

N →∞ N E (W N -W) 2 = f 2 L 2 (µ) -E(W 2 ) + 2 n∈Z\N (ψ, e n )(ϕ, e n ) 1 -ζ n ,
completing the proof of (i).

The proof of (ii) is inspired by [START_REF] Durrett | Probability: theory and examples second ed[END_REF]. First, note that since N = {0}, Lf = µ(f ) for any f ∈ L 2 (µ) and it follows from (i) that 1 N N k=1 f (Z k ) converges almost surely and in L p to µ(f ), for any 1 ≤ p ≤ 2. Moreover given any f ∈ L 1 (µ), one can find

f n ∈ L 2 (µ) such that f -f n L 1 (µ) < 1 n . It follows that for any n ≥ 1, lim sup N →∞ E 1 N N k=1 f (Z k ) -µ(f ) ≤ 2 n + lim sup N →∞ E 1 N N k=1 f n (Z k ) -µ(f n ) = 2 n .
Since the latter is true for any n ≥ 1, one may conclude that 1

N N k=1 f (Z k ) converges in L 1 to µ(f ). By Birkhoff's ergodic theorem, 1 N N k=1 f (Z k ) converges almost surely to µ(f ). Next, let D be the subset of H defined by D =    h ∈ H; n∈Z\{0} |(h, e n )| 2 |1 -ζ n | 2 < ∞    ,
(I -K) -1 (I -L) on D. Next, if f ∈ D, set g = Ξf . Since Lf = µ(f ), it follows that √ N(W N -µ(f )) = 1 √ N N k=1 [g(Z k ) -Kg({S k-1 )}] + 1 √ N Kg(S 0 ) - 1 √ N Kg({S N }).
Now, setting F k = σ{Z j ; j ≤ k}, the terms ξ k = g(Z k ) -Kg({S k-1 )} are square integrable martingale differences with respect to (F j ) j≥0 , i.e. E(ξ k |F k-1 ) = 0, and because g 2 and (Kg) 2 both belong to L 1 (µ), it follows from (i), as shown above, that

1 N N k=1 E ξ 2 k |F k-1 = 1 N N k=1 Kg 2 ({S k-1 }) -(Kg) 2 ({S k-1 })
converges almost surely to µ(g 2 )µ((Kg) 2 ). Note that since Kg = ΞKf , one has (Kg, Lf ) = 0 and expression ( 18) can be written as

σ 2 = (I -L)f 2 L 2 (µ) + 2(ΞKf, Af ) = (I -K)g 2 L 2 (µ) + 2(Kg, Af ) = (I -K)g 2 L 2 (µ) + 2(Kg, A(I -L)f ) = (I -K)g 2 L 2 (µ) + 2(Kg, A(I -K)g) = µ(g 2 -2gKg + (Kg) 2 ) + 2µ(gKg) -2µ((Kg) 2 ) = µ(g 2 ) -µ((Kg) 2 ).
Finally, because of the stationarity of (ξ k ) k≥1 , it follows that for any > 0,

1 N N k=1 E ξ 2 k I(|ξ k | > √ N ) = E ξ 2 1 I(|ξ 1 | > √ N ) → 0,
as N → ∞. The conditions of Theorem 7.4 in [START_REF] Durrett | Probability: theory and examples second ed[END_REF] are all met, so one may safely conclude that defining the process

B N (t) = √ N W N t -µ(f ) , t ∈ [0, 1], then B N converges in D[0, 1] to σB,
where B is a Brownian motion.

To prove part (iii), note first that since the density h of Y is square integrable, then

N = {0}, sup n≥1 |ζ n | = ρ < 1, ζ n = (e n , h), and h 2 2 = n∈Z |ζ n | 2 . Therefore, for any g ∈ H, n∈Z |(g, e n )| |ζ n | ≤ g 2 h 2 < ∞. It follows that for any k ≥ 2, K k f = K k-1 ψ = n∈Z (ψ, e n )ζ k-1 n e n ,
the latter series converging absolutely. Thus

sup z0∈Z K k f (s) -µ(f ) = sup s∈[0,1) K k-1 ψ(s) -λ(ψ) ≤ n∈Z\{0} |(ψ, e n )||ζ n |ρ k-2 ≤ h 2 f L 2 (µ) ρ k-2 .
This completes the proof of the theorem. The next result is a straightforward extension of the previous theorem. Before stating it, denote by ν the joint law of (Z 1 , S 0 ), where S 0 ∼ Unif([0, 1)).

Remark A.1. Note that if ζ n = 1 for some n > 0, then k → ζ k is n-periodic, so {ζ k ; n ∈ Z \ N } is finite. Therefore sup k∈Z\N |ζ k | = ρ <
Theorem A.2. Suppose that f ∈ L 2 (ν) and set W N = 1 N N k=1 f (Z k , {S k-1 }). Then: (i) If the initial distribution of Z 0 = (X 0 , S 0 ) is µ, then W N converges almost
n∈Z\N |(Kf, e n )| |(Af, e n )| |1 -ζ n | < ∞, (22)
then

N E (W N -W) 2 converges, as N → ∞, to f 2 L 2 (ν) -Lf 2 2 + 2 n∈Z\N (Kf, e n )(Af, e n ) 1 -ζ n . ( 23
) (ii) If N = {0}, if the initial distribution of Z 0 is µ and n∈Z\{0} |(Kf, e n )| 2 |1 -ζ n | 2 < ∞, (24)
then the sequence of processes B N , defined by B N (t) = √ N W N tµ(f ) , t ∈ [0, 1], converges in D([0, 1]) to σB, where B is a Brownian motion and σ 2 is given by (23).

(iii) If P Y admits a square integrable density h, then the Markov chain is geometrically ergodic, that is, there exists ρ ∈ (0, 1) such that for any f ∈ L 2 (µ),

|K n f (Z 1 , S 0 ) -µ(f )| ≤ h 2 ρ n-2 f L 2 (µ) , n ≥ 2.
Remark A.2. For example, suppose that X k is bounded and set f (x, y, s) = x(ys).

Then it is easy to check that for any n ∈ N , (Kf, e n ) = Z×[0,1)

x({y + s}s)P (dx, dy)e n (s)ds

= Z×[0,1)
xu(e n (y) -1)e n (u)P (dx, dy)du = 0, since P (e n (Y ) = 1) = 1. It follows from Theorem A.2 that

W N = 1 N N k=1 X k ({S k } -{S k-1 })
converges to 0 almost surely and in mean square.

Furthermore, if card(N ) > 1 then condition (22) holds and sup

N ≥1 N E(W 2 N ) < ∞,
while if P Y is absolutely continuous, then condition (24) holds true and √ NT N converges in law to a centered Gaussian random variable with variance σ 2 given by (23).

f 1 0

 1 (x, {s + y})e n (s)P (dx, dy)ds = Z f (x, u)e n (u)P (dx, dy)du = (ϕ, e n ). As a result, E(W 2 ) = n∈N |(ϕ, e n )| 2 . Next, using the fact that for any k ∈ Z, and any s, y ∈ [0, 1), one has e k ({s + y}) = e k (s + y) = e k (s)e k (y), and it follows that Ke k (s) = Z e k ({s + y})P (dx, dy) = [0,1)

  and let Ξ be the operator from D to H that satisfies Ξh = n∈Z\{0} (h, e n ) 1ζ n e n . Note that since (I -K)Ξh = (I -L)h, then Ξ = (I -K) -1 (I -L) on D. Let D be the set of all f ∈ L 2 (µ) such that f satisfies (19), i.e. Kf ∈ D. Then Ξ can be extended to a mapping from D to L 2 (µ) viz. Ξf = (I -L)f + ΞKf . Using KL = LK = L, one obtains that Ξ =

  surely and in mean square to W given by[START_REF] Kallianpur | Estimation of stochastic systems: Arbitrary system process with additive white noise observation errors[END_REF] whereKf (s) = Z f (x, {s + y}, s)P (dx, dy).If Af (s) = Z f (x, s, {s -y})P (dx, dy), s ∈ [0, 1) and if in addition,

  1 and condition (17) is satisfied. Also, if P Y has a non degenerate absolutely continuous part, then N = {0} and sup n≥1 |ζ n | = ρ < 1, so condition (17) holds true.

Table 1 :

 1 Differences in L 2 , between the particle filters and the optimal filter (OF) for one target performing a simple symmetric random walk in images of size 100 × 100 with 10% of errors.

	n	1	2	3	4	5	100
	N 0 = 1000 SRS 0.27	0.78	0.79 0.81 1.0	1.1
	RSS 0.12	0.37	0.64 0.65 0.68 1.3
	SyS	0.07	0.28	0.31 0.85 0.54 1.1
	BiS	0.04 0.058 0.087 0.15 0.17 0.93
	BeS	0.07	0.16	0.21 0.65 0.69 1.3
	N 0 = 10000 SRS 0.031 0.11	0.26 0.40 0.36 1.2
	RSS 0.039 0.12	0.40 0.50 0.91 1.2
	SyS 0.025 0.13	0.22 0.38 0.40 1.3
	BiS 0.028 0.049 0.082 0.15 0.17 0.96
	BeS 0.028 0.07	0.11 0.23 0.25 1.1
	N 0 = 30000 SRS 0.018 0.10	0.23 0.22 0.50 1.4
	RSS 0.17	0.10	0.23 0.21 0.56 1.4
	SyS 0.018 0.098 0.17 0.20 0.25 1.0
	BiS 0.027 0.048	0.8	0.15 0.17 0.95
	BeS 0.018 0.056 0.10 0.17 0.22 0.97
	N 0 = 50000 SRS 0.013 0.094 0.26 0.22 0.37 1.2
	RSS 0.021 0.068 0.18 0.26 0.39 1.3
	SyS 0.012 0.095 0.26 0.24 0.45 1.2
	BiS	0.02	0.04	0.08 0.15 0.17 0.70
	BeS 0.017 0.052 0.086 0.16 0.20 1.2

Table 2 :

 2 Mean absolute error for one target performing a simple symmetric random walk in images of size 100 × 100 with 10% of errors.

	n	[2, 100] [10, 100] [30, 100]
	OF	4.1	2.9	0.8
	N 0 = 1000 SRS	57.4	60.3	56.2
	RSS	51.8	53.6	43.8
	SyS	42.7	43.7	36.9
	BiS	54.0	56.5	45.5
	BeS	13.8	12.1	6.9
	N 0 = 10000 SRS	76.0	81.0	64.1
	RSS	1.9	0.8	0.5
	SyS	2.4	0.7	0.5
	BiS	6.4	6.7	0.5
	BeS	77.5	82.5	85.4
	N 0 = 30000 SRS	8.7	3.0	0.5
	RSS	2.4	0.6	0.4
	SyS	3.9	1.5	0.4
	BiS	4.0	2.1	0.5
	BeS	8.1	6.2	0.4
	N 0 = 50000 SRS	3.9	3.4	0.9
	RSS	10.2	5.2	0.7
	SyS	5.0	2.5	0.6
	BiS	4.8	2.1	0.8
	BeS	3.6	1.5	0.3

Acknowledgments:

The authors thank the referees for helpful remarks and suggestions that lead to an improved version of the article.

be the mean of the various measures. The simulations were performed with p 0 = 0.9 and p 1 = 0.9, that there are 10% of errors in pixels.

We use two methods to compare the sampling methods described in section 2, i.e. simple random sampling (SRS), remainder stochastic sampling (RSS), systematic sampling (SyS), binomial sampling (BiS), and Bernoulli sampling (BeS).

The first method is a comparison between the sampling method and the Optimal Filter (OF). For the systematic sampling we illustrate Conjecture 3.1 and compute η N nη n 2 for different values of n (1, 2, 3, 4, 5 and 100) and N (1000, 10000, 30000 and 50000). One find that the conjecture seems to be true. We also illustrate with the others sampling methods and observe that the convergence of all samplings methods are quite similar. These results are reported in Table 1.

The second method illustrates how the sampling methods are efficient to estimate the exact position of the target. In that setting, we compare also with the result obtained with the optimal filter. We compute the mean absolute error between the estimated position and the true one was calculated over several time intervals, namely [START_REF] Baker | Reducing bias and inefficiency in the selection algorithm[END_REF]100], [START_REF] Crisan | A particle approximation of the solution of the Kushner-Stratonovitch equation[END_REF]100] and [30,100]. The number of particles N takes values 1000, 10000, 30000 and 50000. The results are reported in Table 2.

According to these results, one may conclude that the algorithm based on the systematic sampling selection method performs quite well, provided the number of particles is large enough. Surprisingly, the Monte Carlo based approximate filters seem to perform better than the optimal filter. However the difference may not be statistically significant. Next, based on the results of Table 2 for the time interval [30, 100], note that when the target is precisely detected, the error seems to stabilize near zero, indicating that the distance between η N n and η n might be uniform on n. Finally, other simulations performed with several moving targets seem to indicate that the algorithm based on systematic sampling also give impressive results.

Appendix A. Convergence results for a Markov chain

Suppose that (X i , Y i ) i≥1 are independent observations of (X, Y ) ∈ Z := R × [0, 1) of law P , with marginal distributions P X and P Y respectively. Further let λ denotes Lebesgue's measure on [0, 1). Given Z 0 = (X 0 , S 0 ) ∈ R × [0, 1), set Z i = (X i , {S i }),