Keywords: Monte Carlo method, randomization, stratification, variance reduction, random walk, numerical quadrature, simulation, diffusion equation, coagulation equation

particles are first sampled from the initial size distribution. A time step is fixed and, in every step and for each particle, a coalescence partner is chosen and a random number decides if coalescence occurs. If the particles are ordered in every time step by increasing sizes and if the random numbers are replaced by stratified points, a variance reduction is observed, when compared to the results of usual MC algorithm.

Finally, we review the work done and we give some perspectives for future work.

i

Résumé

Les méthodes de Monte Carlo (MC) sont des méthodes de simulation probabilistes basées sur l'utilisation des nombres aléatoires pour résoudre avec des ordinateurs des problèmes des sciences appliquées et des techniques. Leur principe de base est simple et repose sur le théorème de limite central : un paramètre, écrit sous forme de l'espérance d'une variable aléatoire, peut être approché par des évaluations répétées et indépendantes de cette variable. L'erreur de la méthode est évaluée par la variance de l'estimateur. Une façon d'améliorer la convergence est de réduire l'ordre de variance des estimateurs.

Le présent travail analyse des méthodes de réduction de la variance et examine leur efficacité pour l'intégration numérique et la résolution d'équations différentielles et intégrales.

Nous présentons d'abord les méthodes MC stratifiées et les méthodes d'échantillonnage par hypercube latin (LHS : Latin Hypercube Sampling). Parmi les méthodes de stratification, nous privilégions la méthode simple (MCS) : l'hypercube unité I s := [0; 1) s est divisé en N sous-cubes d'égale mesure, et un point aléatoire est choisi dans chacun des sous-cubes. Cette génération peut être vue comme une hasardisation du point milieu. Nous analysons la variance de ces méthodes pour le problème de la quadrature numérique. Nous étudions particulièrment le cas de l'estimation de la mesure d'un sous-ensemble de I s . La variance de la méthode MCS peut être majorée par O(1/N 1+1/s ). Les résultats d'expériences numériques en dimensions 2, 3 et 4 montrent que les majorations obtenues sont précises.

Nous proposons ensuite une méthode hybride entre MCS et LHS, qui possède les propriétés de ces deux techniques, avec un point aléatoire dans chaque sous-cube et les projections des points sur chacun des axes de coordonnées également réparties de manière régulière : une projection dans chacun des N sous-intervalles qui divisent I := [0; 1) uniformément. Cette technique est appelée Stratification Sudoku (SS) à cause des propriétés de ses points rappelant ainsi une grille de Sudoku. Dans le même cadre d'analyse que précédemment, nous montrons que la variance de la méthode SS est majorée par O(1/N 1+1/s ) ; des expériences numériques en dimensions 2, 3 et 4 valident les majorations démontrées.

Nous nous intéressons ensuite à deux problèmes d'équations aux dérivées partielles pour lesquels nous proposons une approche de la méthode de marche aléatoire utilisant les techniques de réduction de variance précédentes.

Nous proposons un algorithme de résolution de l'équation de diffusion, avec un coefficient de diffusion constant ou non-constant en espace. On utilise des particules échantillonnées suivant la distribution initiale, qui effectuent un déplacement gaussien à chaque pas de temps. On ordonne les particules suivant leur position à chaque étape et on rem-Résumé place les nombres aléatoires qui permettent de calculer les déplacements par les points stratifiés utilisés précédemment. On évalue l'amélioration apportée par cette technique sur des exemples numériques.

Nous utilisons finalement une approche analogue pour la résolution numérique de l'équation de coagulation, qui modélise l'évolution de la taille de particules pouvant s'agglomérer. Les particules sont d'abord échantillonnées suivant la distribution initiale des tailles. On fixe un pas de temps et, à chaque étape et pour chaque particule, on choisit au hasard un partenaire de coalescence et un nombre aléatoire qui décide de cette coalescence. Si l'on classe les particules suivant leur taille à chaque pas de temps et si l'on remplace les nombres aléatoires par des points stratifiés, on observe une réduction de variance par rapport à l'algorithme MC usuel.

Finalement, nous faisons un bilan du travail effectué et nous proposons quelques perspectives pour des travaux futurs.

Mots-clés : Méthode de Monte Carlo, hasardisation, stratification, réduction de variance, marche aléatoire, quadrature numérique, simulation, équation de diffusion, équation de coagulation.

iii Abstract Monte Carlo (MC) methods are probabilistic simulation methods using random numbers to solve on computers problems from applied sciences and techniques. Their principle is based on the use of the central limit theorem : One estimates a parameter, written as the expected value of a random variable, by repeated and independent evaluations of the variable. The error of the method is approximated through the variance of the estimator. One way to improve the convergence is to reduce the order of the variance of the estimators.

In the present work, we analyze variance reduction methods and we test their efficiency for numerical integration and for solving differential or integral equations.

First, we present stratified MC methods and Latin Hypercube Sampling (LHS) technique. Among stratification strategies, we focus on the simple approach (MCS) : the unit hypercube I s := [0; 1) s is divided into N subcubes having the same measure, and one random point is chosen in each subcube. This technique can be seen as a randomization of the midpoint rule. We analyze the variance of the method for the problem of numerical quadrature. The case of the evaluation of the measure of a subset of I s is particularly detailed. The variance of the MCS method may be bounded by O(1/N 1+1/s ). The results of numerical experiments in dimensions 2, 3 and 4 show that the upper bounds are tight.

Then we propose an hybrid method between MCS and LHS, that has properties of both approaches, with one random point in each subcube and such that the projections of the points on each coordinate axis are also evenly distributed : one projection in each of the N subintervals that uniformly divide the unit interval I := [0; 1). We call this technique Sudoku Sampling (SS) due to the properties of the points recalling a Sudoku grid. Conducting the same analysis as before, we show that the variance of the SS method is bounded by O(1/N 1+1/s ) ; the order of the bound is validated through the results of numerical experiments in dimensions 2, 3 and 4.

Next, we focus on two problems of partial differential equations for which we propose an approach of the random walk method using the variance reduction techniques previously analyzed.

We propose an algorithm for solving the diffusion equation with a constant or spatiallyvarying diffusion coefficient. One uses particles that are sampled from the initial distribution ; they are subject to a Gaussian move in each time step. The particles are renumbered according to their positions in every step and the random numbers which give the displacements are replaced by the stratified points used above. The improvement brought by this technique is evaluated in numerical experiments.

An analogous approach is finally used for numerically solving the coagulation equation ; this equation models the evolution of the sizes of particles that may agglomerate. The

Introduction

Les méthodes de Monte Carlo (MC) ont connu leur essor lors de la Seconde Guerre Mondiale, notamment dans le cadre du projet américain Manhattan portant sur le développement de l'arme nucléaire. Leur principe consiste à modéliser le problème considéré sous forme de l'intégrale d'une fonction en utilisant les outils du calcul des probabilités puis de la simuler numériquement en l'approchant par la moyenne de la fonction en des points pseudo-aléatoires. Ces points sont engendrés par ordinateur et ont des propriétés semblables à celles des points aléatoires. Il est difficile de savoir exactement qui parmi les chercheurs von Neumann, Ulam, Fermi et Metropolis, leur a donné le nom de Monte Carlo qui fait référence aux jeux de hasard pratiqués dans la principauté de Monaco [START_REF] Metropolis | The Monte Carlo method[END_REF][START_REF] Hammersley | Monte Carlo Methods[END_REF]. L'accessibilité des ordinateurs et la facilité de leur utilisation dans la résolution de problèmes réalistes ont popularisé ces méthodes [START_REF] Dautray | Méthodes probabilistes pour les équations de la physique[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF]. Aujourd'hui, elles comptent sans doute parmi les outils les plus utilisés par les ingénieurs. Elles le sont en physique des particules où des simulations probabilistes permettent d'estimer la forme d'un signal ou la sensibilité d'un détecteur, et pour des problèmes de transport ou de diffusion [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]. En mathématiques financières, elles sont couramment mises en oeuvre pour calculer des produits financiers élaborés [START_REF] Glasserman | Monte Carlo in Financial Engineering[END_REF].

À partir de l'estimateur de Monte Carlo du paramètre étudié, il est possible de construire un intervalle de confiance qui mesure la précision de la méthode. La taille de cet intervalle est de l'ordre de N -1/2 , pour N points de simulation, et elle est indépendante de la dimension de l'espace, ce qui les rend avantageuses pour des problèmes de dimension élevée [START_REF] Davis | Méthods of numerical integration[END_REF][START_REF] Evans | Approximating Integrals via Monte Carlo and Deterministic Methods[END_REF]. Cependant, ces méthodes ont plusieurs inconvénients. En particulier, elles ne permettent pas d'exploiter les éventuelles propriétés de régularité de la fonction à intégrer et leur taux de convergence est faible, puisque, pour réduire l'erreur d'un facteur de 100, il faut multiplier N par 10 000. De plus, les simulations MC donnent souvent des résultats très bruités.

Dans l'objectif de réduire la variance d'un estimateur MC, qui est de l'ordre de N -1 , différentes techniques ont été mises au point, parmi lesquelles nous citons : la méthode des variables antithétiques, des variables de contrôle, d'échantillonnage préférentiel, d'échantillonnage stratifié [START_REF] Fishman | Concepts, Algorithms, and Applications[END_REF][START_REF] Madras | Lectures on Monte Carlo Methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF].

Une autre possibilité pour améliorer la convergence des méthodes MC est de remplacer les nombres pseudo-aléatoires, qui simulent des variables aléatoires indépendantes, par des nombres quasi-aléatoires obtenus à partir de suites réparties le plus uniformément possible dans le domaine d'intégration considéré. On parle alors de méthodes quasi-Monte Carlo [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF][START_REF] Niederreiter | Quasi-Monte Carlo methods and pseudo-random numbers[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF] Drmota | Sequences, Discrepancies and Applications[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. L'ordre de convergence de ces méthodes pour le problème de l'intégration numérique est de N -1+ǫ . En se limitant à des fonctions régulières et à des ensembles de points particuliers, on peut même avoir un ordre en N -α+ǫ avec α > 0 arbitrairement grand [START_REF] L'ecuyer | Quasi-Monte Carlo methods with applications in finance[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. Pour ces méthodes déterministes, nous disposons d'une borne d'erreur due à Koksma et Hlawka qui fait intervenir la variation de la fonction à intégrer et qu'il est en général difficile de caluler précisément. De plus cette borne est pessimiste car l'ordre de convergence obtenu numériquement est meilleur que l'ordre théorique. Pour y remédier, on peut employer des nombres quasi-aléatoires hasardisés, c'est-à-dire faisant intervenir un paramètre aléatoire, grâce auquel on peut calculer la variance du résultat.

L'objectif de ce travail de thèse est la mise au point, le développement et la validation par des tests numériques de méthodes de simulation stratifiées dont les points sont aléatoires et possèdent également des propriétés de répartition uniforme comparables à celles des points QMC dans le but de réduire la variance par rapport à une méthode MC classique.

Au premier chapitre, nous présentons les méthodes MC, QMC et des méthodes quasi-Monte Carlo hasardisées. Nous les comparons dans des calculs numériques d'intégrales.

Au chapitre 2, une méthode de Monte Carlo stratifiée simple (MCS) est introduite et analysée. Nous considérons la stratégie suivante [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF] : le cube unité de dimension s, I s = [0, 1) s , est partitionné en N sous-cubes identiques et dans chacun d'entre eux nous engendrons une variable aléatoire distribuée uniformément. Cela revient à effectuer une hasardisation des points de quadrature de la méthode des points milieux. Cette méthode d'échantillonnage stratifié a été analysée dans [START_REF] Cheng | The problem of dimensionality in stratified sampling[END_REF] pour l'intégration des fonctions régulières. Nous nous intéressons ici à l'intégration de fonctions indicatrices de sous-ensembles de I s à l'aide de l'ensemble stratifié ainsi construit, puisque certains problèmes de simulation peuvent être transformés en des évaluations de telles intégrales. Nous établissons une borne de la variance de l'estimateur dont nous prouvons l'optimalité à travers des expériences numériques. La majoration de variance obtenue est de l'ordre de N -1-1/s , ce qui prouve que cette méthode est bien une technique de réduction de la variance, au moins pour N suffisamment grand. Dans le cas où la frontière du domaine considéré est représentée par une fonction f : I s-1 → I, la borne d'erreur s'exprime en fonction de la variation de f . L'analyse de cette approximation par une méthode quasi-Monte Carlo a été faite dans [START_REF] Niederreiter | Diskrepanz und Distanz von Maßen bezüglich konvexer und Jordanscher Mengen[END_REF][START_REF] Lécot | Error bounds for quasi-Monte Carlo integration with nets[END_REF]. Nous analysons également la méthode d'échantillonnage par hypercube latin (LHS) présentée dans [START_REF] Mckay | A comparaison of three methods for selecting values of imput variables in the analysis of output from a computer code[END_REF]. Nous majorons la variance de l'estimateur du volume d'un rectangle à côtés parallèles aux axes de coordonnées et nous comparons les différentes méthodes dans des expériences numériques.

Dans le chapitre 3, nous développons une version hybride de la méthode stratifiée simple et de la méthode LHS. Il s'agit de construire un ensemble stratifié de points de I s qui possède certaines des propriétés de répartition uniforme des réseaux utilisés dans les méthodes QMC. L'ensemble des points ainsi construit rappelle une grille de Sudoku, c'est pourquoi nous parlons de méthode de Monte Carlo stratifiée Sudoku. Nous démontrons une borne de la variance de l'estimateur de l'intégrale d'une fonction indicatrice. Des expériences numériques illustrent l'avantage de cette nouvelle méthode puisque les vitesses de convergence sont meilleures que celles des méthodes MC, LHS et MCS.

Le chapitre 4 est consacré à la simulation de l'équation de diffusion dans des milieux homogène et hétérogène par des méthodes de marche aléatoire [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF][START_REF] Hald | Convergence of random methods for a reaction-diffusion equation[END_REF][START_REF] Ghoniem | Grid-free simulation of diffusion using random walk methods[END_REF][START_REF] Sherman | A Monte Carlo method for scalar reaction diffusion equations[END_REF][START_REF] Fogelson | Optimal smoothing in function-transport particle methods for diffusion problems[END_REF]. Pour l'équation de diffusion en milieu infini, à coefficient de diffusion constant D 0 , nous considérons un ensemble de particules qui, à chaque pas de temps ∆t, effectuent des déplacements aléatoires. Un ensemble stratifié Sudoku est utilisé pour sélectionner les particules qui vont se déplacer et déterminer les amplitudes de leurs déplacements. Dans le cas d'un coefficient de diffusion qui varie en fonction de l'espace, D(x), il ne suffit pas de remplacer D 0 par D(x) dans les expressions des déplacements [START_REF] Hunter | On the use of random walk models with spatially variable diffusivity[END_REF]. Une correction des pas gaussiens effectués par les particules a été mise en oeuvre par Farnell et Gibson dans [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium : correction to the Gaussian steplength[END_REF] pour une simulation MC. Une version QMC a été mise au point dans [START_REF] El Haddad | Quasi-Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium[END_REF]. Nous proposons ici une version stratifiée de cet algorithme corrigé. Dans les deux cas, une étape de renumérotation des particules est essentielle pour la convergence du schéma numérique et permet de bénéficier des propriétés de répartition de l'ensemble stratifié. L'efficacité des deux algorithmes est mise en évidence à travers des essais numériques dans lesquels nous les comparons à l'algorithme de marche aléatoire classique. Les résultats obtenus indiquent que le nouveau schéma converge plus rapidement qu'un algorithme de Monte Carlo.

Enfin, au chapitre 5, nous nous intéressons à la résolution numérique des équations de coagulation discrète et continue. Pour chaque équation, un algorithme stratifié est proposé. Il est inspiré des algorithmes MC proposés par Babovsky [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF] et QMC de [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF][START_REF] Lécot | A quasi-stochastic simulation of the general dynamics equation for aerosols[END_REF] et consiste à approcher la solution de l'équation par une somme de mesures de Dirac, celles-ci représentent la distribution de masse qui est solution de l'équation. Une première étape consiste à calculer un ensemble de N points dont la distribution approche la loi de probabilité de la condition initiale et qui, dans la suite, vont évoluer selon la dynamique de l'équation de coagulation pour approcher la solution à tout instant. Un schéma d'Euler est utilisé pour discrétiser l'équation et à chaque pas de temps, le problème est formulé comme l'intégration d'une indicatrice d'un sous-ensemble de I 3 . Un ensemble de points stratifiés Sudoku est alors utilisé pour effectuer une quadrature et déterminer l'évolution des particules. Comme au chapitre 4, un tri des particules selon leur masse croissante est effectué à chaque pas de temps. La comparaison de la méthode stratifiée avec la méthode de marche aléatoire met en évidence l'efficacité de la première à travers des essais numériques.

Finalement, dans la conclusion, nous résumons les résultats obtenus et nous donnons des pistes de recherche pour des études futures. Les méthodes MC comptent sans doute parmi les outils les plus utilisés par les ingénieurs d'aujourd'hui. Leur principe est simple : il consiste à échantillonner la quantité à estimer à l'aide de points (pseudo-)aléatoires. Dans le cas de l'intégration numérique sur I s := [0, 1) s , on approche

I = I s f (x)dx
par la moyenne I des valeurs prises par la fonction f sur un ensemble de N variables aléatoires indépendantes X i , i = 1, . . . , N , où chaque X i est de loi uniforme sur [0, 1) s :

X i U [0, 1) s . C'est-à-dire I := 1 N N i=1 f (X i ).
En calculant la variance de l'estimateur (sans biais) ou en utilisant le théorème de limite central, on montre que l'ordre de convergence de la méthode est O(N -1/2 ). Les méthodes MC sont aussi utilisées non seulement pour calculer des intégrales [START_REF] Hammersley | Monte Carlo Methods[END_REF][START_REF] Davis | Méthods of numerical integration[END_REF][START_REF] Fishman | Concepts, Algorithms, and Applications[END_REF][START_REF] Evans | Approximating Integrals via Monte Carlo and Deterministic Methods[END_REF][START_REF] Madras | Lectures on Monte Carlo Methods[END_REF][START_REF] Tuffin | La Simulation de Monte Carlo[END_REF] mais aussi pour simuler des phénomènes de la physique, de la mécanique (des fluides ou des solides), de la chimie, de la biologie [START_REF] Dautray | Méthodes probabilistes pour les équations de la physique[END_REF][START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF] et estimer des produits financiers [START_REF] Glasserman | Monte Carlo in Financial Engineering[END_REF].

Ces méthodes présentent plusieurs inconvénients. En particulier elles ne permettent pas, sans modification, d'exploiter les éventuelles propriétés de régularité de la fonction. Un autre inconvénient est leur faible taux de convergence : pour réduire l'erreur d'un facteur de 100, il faudrait multiplier N par 10 000. Comme les estimateurs sont sans biais, les résultats sont donc très bruités.

Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

Les versions déterministes des méthodes MC sont les méthodes QMC consistant à remplacer les variables aléatoires indépendantes par des éléments d'ensembles de points à faible discrépance. Ceux-ci sont construits de façon à couvrir le domaine d'intégration I s le plus uniformément possible. Sous certaines hypothèses de régularité sur f , l'ordre de convergence des méthodes QMC est O(N -1+ǫ(s) ) [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF][START_REF] Niederreiter | Quasi-Monte Carlo methods and pseudo-random numbers[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF] Tezuka | Uniform Random Numbers : Theory and Practice[END_REF][START_REF] Drmota | Sequences, Discrepancies and Applications[END_REF][START_REF]Random and Quasi-Random Point Sets[END_REF][START_REF] Niederreiter | Constructions of (t, m, s)-nets and (t, s)-sequences[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF].

Un inconvénient des méthodes QMC est que les majorations d'erreur déterministes sont souvent trop larges donc inexploitables, quand on les compare aux intervalles de confiance des méthodes MC. Dans la pratique, les méthodes QMC peuvent concurrencer les méthodes MC en dimension petite (s ≤ 3) ou modérée (4 ≤ s ≤ 10). Mais en grande dimension (s ≥ 11), l'approche QMC pure devient inefficace.

Pour améliorer les caractéristiques et les performances des méthodes QMC, les spécialistes (académiques, industriels ou financiers) ont essayé d'introduire des paramètres aléatoires dans les ensembles de points déterministes. Cette démarche, assez récente, permet de bénéficier de la précision des méthodes QMC et de disposer d'intervalles de confiance pour les estimateurs, comme pour les méthodes MC. Ces méthodes hybrides sont souvent appelées méthodes quasi-Monte Carlo randomisées, que nous préférons nommer hasardisées. Le caractère aléatoire des points de ces méthodes permet l'évaluation de la variance de l'estimateur en faisant un certain nombre de répétitions indépendantes. Parmi les synthèses récentes sur les méthodes quasi-Monte Carlo hasardisées sont les publications de P. L'Ecuyer et C. Lemieux [START_REF] L'ecuyer | Recent advances in randomized quasi-Monte Carlo methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF].

Les méthodes de Monte Carlo

Nous commencons par l'introduction des méthodes de Monte Carlo dans le cadre de l'intégration numérique.

Intégration numérique multi-dimensionnelle

Dans de nombreux problèmes scientifiques, il s'agit de calculer une intégrale

I(f ) = I s f (x)dx,
où s est une dimension. Il n'est pas toujours possible de l'obtenir exactement : on cherche alors à l'approcher.

Les formules de quadrature classiques permettent d'approcher l'intégrale d'une fonction par une somme pondérée de ses valeurs en différents points. Parmi ces méthodes, les plus élémentaires sont les formules des rectangles, des trapèzes, des points milieux, de Simpson.

Le détail de la formule des points milieux composée va permettre de mettre en évidensce une caractéristique faisant intervenir la dimension s du problème. La version la plus simple utilise un entier m et un pas h := 1/m ; alors

I(f ) ≈ S(f ), où S(f ) := h s m-1 j 1 =0 . . . m-1 js=0 f (2j 1 + 1) h 2 , . . . , (2j s + 1) h 2 .
L'erreur est définie par :

E(f ) := I(f ) -S(f ).
Si f est de classe C 2 sur I s , l'erreur est de l'ordre

O 1 m 2 = O 1 N 2/
R s |f (x)|dπ(x) < +∞.

On note

I π (f ) := R s f (x)dπ(x).
Soit (X n ) n≥1 une suite de variables aléatoires s-dimensionnelles sur un espace probabilisé (Ω, T , P), équidistribuées de loi π. Alors chaque f • X n est une variable aléatoire réelle intégrable et l'on a

E[f • X n ] = Ω f • X n (ω)dP(ω) = R s f (x)dπ(x) = I π (f ).
Si (Y n ) n≥1 est une suite de variables aléatoires réelles équidistribuées, intégrables, de moyenne µ, on appelle moyenne empirique de N termes et on note Y N la variable aléatoire réelle :

Y N := 1 N (Y 1 + Y 2 + . . . + Y N ).
Alors

E[Y N ] = 1 N (E[Y 1 ] + E[Y 2 ] + . . . + E[Y N ]) = E[Y 1 ] = µ.
Ce qui prouve que Y N est un estimateur sans biais de µ. Dans le cas de l'intégration MC,

f • X N = 1 N (f • X 1 + f • X 2 + . . . + f • X N )
est un estimateur sans biais de I π (f ). La convergence de la méthode de Monte Carlo est assurée par la loi forte des grands nombres. → I π (f ).

Intervalle de confiance

Sous des hypothèses plus fortes, la méthode de Monte Carlo fournit aussi un intervalle qui contient l'approximation de l'intégrale avec une probabilité (de limite) donnée. C'est la définition de ce qu'on appelle l'intervalle de confiance. On fait donc l'hypothèse supplémentaire :

R s |f (x)| 2 dπ(x) < +∞ et l'on note σ 2 π (f ) := R s |f (x)| 2 dπ(x) - R s f (x)dπ(x) 2 . (1.2.1)
Chaque variable aléatoire réelle f • X n est alors de carré intégrable et

E[(f • X n ) 2 ] = Ω (f • X n (ω)) 2 dP (ω) = R s |f (x)| 2 dπ(x), donc Var(f • X n ) = σ 2 π (f ). Si (Y n ) n≥1 est
une suite indépendante de variables aléatoires réelles équidistribuées de carré intégrable, de moyenne µ et de variance σ 2 , on a :

Var(Y N ) = 1 N 2 N n=1 Var(Y n ) = σ 2 N .
On appelle variance empirique (de N termes) et on note S 2 N la variable aléatoire réelle

S 2 N := 1 N -1 N n=1 (Y n -Y N ) 2 = 1 N -1 N n=1 Y n 2 -N Y N 2 .
On a

E[S 2 N ] = N N -1 (E[Y 2 1 ] -E[Y 2 N ]) = N N -1 (σ 2 + µ 2 -E[Y on a E[S 2 N ] = N N -1 σ 2 + µ 2 - σ 2 N -µ 2 = σ 2 .

Donc S

2 N est un estimateur sans biais de σ 2 . Dans le problème de l'intégration numérique

S 2 π,N (f ) := 1 N -1 N n=1 (f • X n ) 2 -N (f • X n ) 2
est un estimateur sans biais de σ 2 π (f ). L'intérêt de intervalle de confiance est basé sur le théorème de limite central qui permet de contrôler la précision de la méthode de Monte Carlo.

Théorème 1.2 (Théorème de limite central) Si (Y n ) n≥1 est une suite indépendante de variables aléatoires réelles équidistribuées dans L 2 , on a

1 √ N N n=1 Y n -E[ N n=1 Y n ] L → N (0, σ 2 ), où σ 2 est la variance commune des Y n .
Soit (Y n ) n≥1 une suite indépendante de variables aléatoires équidistribuées dans L 2 . On note µ et σ 2 l'espérance et la variance commune des Y n et on pose :

Z N := Y N -µ σ/ √ N .
On remarque que

E[Z N ] = 0 et Var(Z N ) = 1. En utilisant le fait que √ N (Y N -µ) L → N (0, σ 2 ), on obtient : Z N L → N (0, 1).
Donc, pour a > 0 lim

N →∞ P(|Z N | ≤ a) = 2Φ(a) -1, où Φ(x) := 1 √ 2π x -∞
e -y 2 /2 dy, x ∈ R est la fonction de répartition de la loi normale centrée réduite N (0, 1). On a donc lim

N →∞ P |Y N -µ| ≤ σa √ N = 2Φ(a) -1.

Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

(R s , B(R s )) et f est une fonction mesurable sur (R s , B(R s )), telle que R s (f (x)) 2 dπ(x) < +∞.
La suite (X n ) n≥1 est une suite indépendante de variables aléatoires s-dimensionnelles équidistribuées de loi π, donc la famille (f •X n ) n≥1 est une suite indépendante de variables aléatoires réelles équidistribuées, d'espérance I π (f ) et de variance σ 2 π (f ). Le théorème de limite central implique, que pour tout a > 0, lim

N →∞ P 1 N N n=1 f • X n -I π (f ) ≤ σ π (f )a √ N = 2Φ(a) -1. L'intervalle aléatoire 1 N N n=1 f • X n - σ π (f )a √ N , 1 N N n=1 f • X n + σ π (f )a √ N
s'appelle l'intervalle de confiance de niveau de confiance

β := 2Φ(a) -1.
La Table 1.1 donne les valeurs du niveau de confiance en fonction de a. On a la conséquence suivante du théorème de limite central, utile en pratique.

Proposition 1.1 Soit (Y n ) n≥1 une suite indépendante de variables aléatoires réelles équidistribuées et dans L 2 , d'espérance µ et de variance σ 2 . Si

Y N := 1 N N n=1 Y n et S 2 N := 1 N -1 N n=1 Y 2 n -N Y 2 N ,
on a, pour tout a > 0, lim

N →∞ P Y N -µ ≤ S N a √ N = 2Φ(a) -1, où S N := S 2 N .
Preuve. La loi forte des grands nombres (théorème 1.1) implique :

S N p.s. → σ.
D'après le théorème de limite central (théorème 1.2), pour tout a > 0, lim

N →∞ P(|Z N | ≤ σa) = 2Φ(a) -1, où Z N := √ N (Y N -µ).
Soit 0 < α < a. On a

{|Z N | ≤ S N a} = {|Z N | ≤ σ(a + α), |Z N | ≤ S N a} + {|Z N | > σ(a + α), |Z N | ≤ S N a} ⊂ {|Z N | ≤ σ(a + α)} + {σ(a + α) < S N a} donc P(|Z N | ≤ S N a) ≤ P(|Z N | ≤ σ(a + α)) + P σα a < S N -σ .
De même

{|Z N | ≤ σ(a -α)} = {|Z N | ≤ σ(a -α), |Z N | ≤ S N a} + {|Z N | ≤ σ(a -α), |Z N | > S N a} ⊂ {|Z N | ≤ S N a} + {S N a < σ(a -α)} donc P(|Z N | ≤ σ(a -α)) ≤ P(|Z N | ≤ S N a) + P σα a < σ -S N .
Par conséquent

P(|Z N | ≤ σ(a -α)) -P σα a < σ -S N ≤ P(|Z N | ≤ S N a) ≤ P(|Z N | ≤ σ(a + α)) + P σα a < S N -σ . Soit ε > 0. Comme lim α→0 α>0 Φ(a -α) = Φ(a) = lim α→0 α>0 Φ(a + α), il existe 0 < α < a tel que Φ(a) -ε < Φ(a -α) < Φ(a) < Φ(a + α) < Φ(a) + ε. Il existe N 1 ∈ N tel que si N ≥ N 1 , P(|Z N | ≤ σ(a -α)) ≥ 2Φ(a -α) -1 -ε ≥ 2Φ(a) -1 -3ε, P(|Z N | ≤ σ(a + α)) ≤ 2Φ(a + α) -1 + ε ≤ 2Φ(a) -1 + 3ε.
→ σ donc il existe N 2 ∈ N tel que si N ≥ N 2 , P σ -S N > σα a ≤ P σ -S N > σα a < ε, P S N -σ > σα a ≤ P S N -σ > σα a < ε. Finalement pour tout N ≥ max(N 1 , N 2 ) 2Φ(a) -1 -4ε ≤ P(|Z N | ≤ S N a) ≤ 2Φ(a) -1 + 4ε donc |P(|Z N | ≤ S N a) -2Φ(a) + 1| ≤ 4ε.

Par conséquent lim

N →∞ P(|Z N | ≤ S N a) = 2Φ(a) -1,
qui est le résultat. Dans le problème d'intégration, sous les mêmes conditions sur la suite (X n ) n≥1 , on a pour tout a > 0, lim

N →∞ P 1 N N n=1 f • X n -I π (f ) ≤ S π,N (f )a √ N = 2Φ(a) -1, où S π,N (f ) := S π,N (f ) 2 . L'intervalle aléatoire 1 N N n=1 f • X n - S π,N (f )a √ N , 1 N N n=1 f • X n + S π,N (f )a √ N
est l'intervalle de confiance pratique de niveau de confiance β = 2Φ(a) -1. En effet la variance σ 2 π (f ) est tout autant inconnue que I π (f ). La variance contrôle l'erreur de la méthode MC. L'usage d'échantillons stratifiés de points aléatoires peut réduire cette variance [START_REF] Fishman | Concepts, Algorithms, and Applications[END_REF][START_REF] Evans | Approximating Integrals via Monte Carlo and Deterministic Methods[END_REF][START_REF] Madras | Lectures on Monte Carlo Methods[END_REF][START_REF] Glasserman | Monte Carlo in Financial Engineering[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF]. Nous présentons dans la suite cette technique d'échantillonnage.

Méthode d'échantillonnage stratifié

Soit N > 0 un entier et soient X 1 , . . . , X N des variables aléatoires indépendantes de loi π. On note

S := 1 N N n=1 f • X n .
On a

E[S] = I π (f ) et Var(S) = σ 2 π (f ) N .
Alors S est l'estimateur de Monte Carlo de I π (f ).

Soit {D i : 1 ≤ i ≤ p}
une partition de R s formée d'ensembles boréliens :

R s = p i=1 D i et i = j ⇒ D i ∩ D j = ∅.
On note

a i := D i dπ(x) et I i (f ) := D i f (x)dπ(x). Alors p i=1 a i = 1 et p i=1 I i (f ) = I π (f ).
Si

π i := 1 D i a i π,
alors π i est une probabilité sur R s . Soient n 1 , . . . , n p des entiers strictement positifs et pour chaque i :

1 ≤ i ≤ p soient X (i) 1 , . . . , X (i)
n i des variables aléatoires de loi π i . On suppose que la famille de variables aléatoires {X

(i) k : 1 ≤ i ≤ p , 1 ≤ k ≤ n i } est indépendante. Soit T i := 1 n i n i k=1 f • X (i) k , 1 ≤ i ≤ p et T := p i=1 a i T i . Comme E[T i ] = 1 n i n i k=1 E[f • X (i) k ] = E[f • X (i) 1 ] = Ω f • X (i) 1 (ω)dP(ω) = R s f (x)dπ i (x) = 1 a i D i f (x)dπ(x) = I i (f ) a i , on a E[T ] = p i=1 a i I i (f ) a i = I π (f ),
donc T est un estimateur sans biais de I π (f ) : on l'appelle l'estimateur stratifié de I π (f ). On note

σ 2 i (f ) := R s (f (x)) 2 dπ i (x) - R s f (x)dπ i (x) 2 .

Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

Comme

Var(T i ) = 1 n 2 i n i k=1 Var(f • X (i) k ) = 1 n i Var(f • X (i) 1 ) = 1 n i Ω (f • X (i) 1 ) 2 (ω)dP(ω) - Ω f • X (i) 1 (ω)dP(ω) 2 = 1 n i R s (f (x)) 2 dπ i (x) - R s f (x)dπ i (x) 2 = σ 2 i (f ) n i , on a Var(T ) = p i=1 a 2 i Var(T i ) = p i=1 (a i σ i (f )) 2 n i .
On suppose dans la suite

p i=1 n i = N. (1.2.2)
Pour la méthode d'allocation proportionnelle, on choisit

n i = N a i , 1 ≤ i ≤ p.
Alors

Var(T ) = 1 N p i=1 a i σ 2 i (f ) = 1 N p i=1 a i R s (f (x)) 2 dπ i (x) - 1 N p i=1 I i (f ) 2 a i = 1 N p i=1 D i (f (x)) 2 dπ(x) - 1 N p i=1 I i (f ) 2 a i = 1 N R s (f (x)) 2 dπ(x) - 1 N p i=1 I i (f ) 2 a i = Var(S) + I π (f ) 2 N - 1 N p i=1 I i (f ) 2 a i donc Var(S) = Var(T ) + 1 N p i=1 a i I i (f ) a i -I π (f ) 2 ≥ Var(T ).
La méthode de Monte Carlo d'échantillonnage stratifié avec allocation proportionnelle ne peut pas être plus mauvaise que la méthode de Monte Carlo simple, car la variance de l'estimateur de la première ne peut pas dépasser celle de la seconde. On note

δ := 1 N p i=1 a i σ i (f ) et V (n 1 , . . . , n p ) := p i=1 (a i σ i (f )) 2 n i . Si (n 1 , . . . , n p ) ∈ N * p vérifie (1.2.2), V a 1 σ 1 (f ) δ , . . . , a p σ p (f ) δ = 1 N p i=1 a i σ i (f ) 2 ≤ 1 N p i=1 (a i σ i (f )) 2 n i p i=1 n i = V (n 1 , . . . , n p ),
en utilisant l'inégalité de Cauchy-Schwarz. Le choix optimal de (n • la borne d'erreur est probabiliste,

• la convergence est lente quand N → ∞,

• la régularité de la fonction à intégrer n'apparaît pas dans l'estimation,

• la génération d'une suite de variables aléatoires est problématique.

Une voie pour éviter ces inconvénients en conservant des avantages des méthodes MC est de remplacer les suites de points aléatoires par des suites déterministes ayant une meilleure répartition uniforme ; ces suites sont appelées suites à discrépance faible. C'est le principe des méthodes quasi-Monte Carlo.

Les méthodes quasi-Monte Carlo

On va de nouveau étudier l'approximation

I s f (x)dx ≈ 1 N N n=1 f (x n ),
où f est une fonction définie sur I s , à valeurs dans R et (x n ) 1≤n≤N est une famille de points de I s . On cherche à obtenir une majoration d'erreur meilleure que O(N -1/2 ). Les ouvrages [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF][START_REF] Niederreiter | Quasi-Monte Carlo methods and pseudo-random numbers[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF] Tezuka | Uniform Random Numbers : Theory and Practice[END_REF][START_REF] Drmota | Sequences, Discrepancies and Applications[END_REF][START_REF]Random and Quasi-Random Point Sets[END_REF][START_REF] Niederreiter | Constructions of (t, m, s)-nets and (t, s)-sequences[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF] sont les références classiques sur les méthodes quasi-Monte Carlo.

D N (B, X) := A(B, X) N -λ s (B).
La discrépance de l'ensemble X est définie par

D N (X) := sup J |D N (J, X)|,
où J parcourt l'ensemble des sous-intervalles semi-ouverts de I s :

J = s i=1 [u i , v i ).
La discrépance à l'origine de X est donnée par

D ⋆ N (X) := sup J ⋆ |D N (J ⋆ , X)|,
où J ⋆ parcourt l'ensemble des sous-intervalles semi-ouverts de I s de la forme

J ⋆ = s i=1 [0, w i ).
Si X := (x n ) n≥1 est une suite infinie de points de I s , on note D N (X) la discrépance de l'ensemble {x 1 , x 2 , . . . , x N } des N premiers termes de la suite, et on note D ⋆ N (X) la discrépance à l'origine de cet ensemble.

Discrépance et discrépance à l'origine sont des erreurs d'intégration des fonctions indicatrices des intervalles. Elles sont liées dans le résultat suivant.

Proposition 1.2 Les discrépances D N (X) et D ⋆ N (X) vérifient D ⋆ N (X) ≤ D N (X) ≤ 2 s D ⋆ N (X).
En dimension s = 1, on a des formules explicites de calcul des discrépances.

Proposition 1.3 Si X := {x 1 , x 2 , . . . , x N } ⊂ I avec 0 ≤ x 1 ≤ x 2 ≤ • • • ≤ x N < 1 on a D ⋆ N (X) = 1 2N + max 1≤n≤N x n - 2n -1 2N , D N (X) = 1 N + max 1≤n≤N n N -x n -min 1≤n≤N n N -x n .
On peut conclure du résultat précédent que l'ensemble

2n -1 2N : 1 ≤ n ≤ N
est celui qui a la plus petite discrépance à l'origine possible. Notre but est de majorer l'erreur d'intégration :

I f (x)dx - 1 N N n=1 f (x n )
pour une classe de fonctions assez générale. On supposera que f est à variation bornée sur I. On note la variation de f

V (f ) := sup Ξ p j=1 |f (ξ j ) -f (ξ j-1 )|
où Ξ = {ξ 0 , ξ 1 , . . . , ξ p } parcourt l'ensemble des familles de points tels que

0 = ξ 0 < ξ 1 < • • • < ξ p = 1.
Toute fonction à variation bornée sur I est intégrable au sens de Riemann sur I.

Le résultat suivant est dû à J.F. Koksma (1942).

Proposition 1.4 Si f a une variation bornée V (f ) sur I et si X := {x 1 , x 2 , . . . , x N } est un ensemble de points de I on a

I f (x)dx - 1 N N n=1 f (x n ) ≤ V (f )D ⋆ N (X).
Pour étendre l'inégalité de Koksma au cas multi-dimensionnel, il faut disposer d'une notion de variation des fonctions de plusieurs variables [START_REF] Hobson | The Theory of Functions of a Real Variable and the Theory of the Fourier's Series[END_REF]. Pour une fonction f :

I s → R et pour w, w ′ ∈ I s et 1 ≤ i ≤ s, on note T i w f la restriction de f à l'hyperplan x i = w i , ∆ i w,w ′ f := T i w ′ f -T i w f.

Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

Si K = {i 1 , . . . , i k } ⊂ {1, . . . , s}, on note

T K w f := T i 1 w • • • T i k w f et ∆ K w,w ′ f := ∆ i 1 w,w ′ • • • ∆ i k w,w ′ f.
On pose T w f := T {1,...,s}

w f et ∆ w,w ′ f := ∆ {1,...,s} w,w ′ f. Si 0 = x 0,i < x 1,i < • • • < x n i ,i = 1, pour 1 ≤ i ≤ s
définit une partition de I s en sous-intervalles, et si a = (a 1 , . . . , a s ) est un indice multidimensionnel (avec des entiers a i , 0 ≤ a i < n i ), on note

x a := (x a 1 ,1 , . . . , x as,s ) et a+ = (a 1 + 1, . . . , a s + 1).

V (k) (T K c 1 f ), où K c := {1, . . . , s} \ K.
L'analogue de l'inégalité de Koksma est le résultat suivant, dû à E. Hlawka [START_REF] Hlawka | Funktionen von beschränkter Variation in der Theorie der Gleichverteilung[END_REF] (voir aussi sa présentation dans [START_REF] Zaremba | Some applications of multidimensional integration by parts, Annales Polonici Mathematici[END_REF]). Théorème 1.3 (Koksma-Hlawka) Si f a une variation bornée V (f ) sur I s au sens de Hardy et Krause, et si X := {x 1 , x 2 , . . . , x N } est un ensemble de N points de I s on a

I s f (x)dx - 1 N N n=1 f (x n ) ≤ V (f )D ⋆ N (X).

Points quasi-aléatoires

L'analyse d'erreur des méthodes quasi-Monte Carlo d'intégration numérique montre que l'utilisation d'ensembles de points à faible discrépance garantit des erreurs petites (pour des fonctions à variation bornée). Dans la suite, il est commode de faire commencer les indices des ensembles ou suites à zéro.

On a vu qu'à une dimension, l'ensemble de N points qui a la plus petite discrépance à l'origine est l'ensemble des points x n , 0 ≤ n < N , où

x n := 2n + 1 2N
(Proposition 1.3). Cette discrépance vaut 1/2N . Concernant les suites infinies, on sait qu'en dimension 1, il n'existe aucune suite dont la discrépance décroisse plus rapidement que O((Ln N )/N ) [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] ; on connaît des suites ayant exactement cette décroissance, comme les suites de Van

φ b (n) = ∞ j=0 a j (n)b -j-1 , n ≥ 0.
La suite de Van der Corput en base b est la suite (x n ) n≥0 où

x n := φ b (n), n ≥ 0.
Une extension en dimension quelconque est donnée dans la définition suivante [START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals[END_REF].

x n := φ b 1 (n), φ b 2 (n), . . . , φ bs (n) , n ≥ 0.
H. Niederreiter a montré dans [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] que si b 1 , . . . , b s sont premiers entre eux deux à deux, la discrépance à l'origine de la suite de Halton S vérifie pour tout N ≥ 2 :

D * N (S) ≤ c s (log N ) s N + O (log N ) s-1 N , où c s := 1 2 s s i=1 (b i -1) log b i .
La majoration de la discrépance à l'origine de ces ensembles a été améliorée par Atanassov qui a établi dans [START_REF] Atanassov | On the discrepancy of the Halton sequence[END_REF] une majoration de la discrépance à l'origine avec une constante

c ′ s := 1 s!2 s s i=1 (b i -1) log b i ,
qui tend vers 0 quand s tend vers l'infini. En se limitant à des ensembles finis, on a l'ensemble de Hammersley de bases b 1 , b 2 , . . . , b s-1 [START_REF] Hammersley | Monte Carlo methods for solving multivariable problems[END_REF] :

(x n ) 0≤n<N où x n := n N , φ b 1 (n), φ b 2 (n), . . . , φ b s-1 (n) , 0 ≤ n < N.
Dans le cas où b 1 , . . . , b s-1 sont premiers entre eux deux à deux, la discrépance à l'origine d'un ensemble de Hammersley à N points est de l'ordre de (log N ) s-1 /N [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF].

Des comparaisons pratiques des suites de points quasi-aléatoires sont données dans [START_REF] Morokoff | Quasi random sequences and their discrepancies[END_REF].

1.3.3 Réseaux (t, m, s) et suites (t, s)
Actuellement, les ensembles de points qui ont la plus faible discrépance sont obtenus à partir de la théorie des réseaux (t, m, s) et des suites (t, s) en base b [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF]Random and Quasi-Random Point Sets[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF].

Définition 1.5 Soit b ≥ 2 un entier.
Un intervalle élémentaire en base b est un intervalle de la forme (1) Si X est un réseau (t, m, s) en base b, alors X est un réseau (u, m, s) en base b, pour tout entier u tel que t ≤ u ≤ m.

J = s i=1 a i b d i , a i + 1 b d i , avec des entiers d i ≥ 0 et 0 ≤ a i < b d i pour 1 ≤ i ≤ s. Soient 0 ≤ t ≤ m
(2) Si X est une suite (t, s) en base b, alors X est une suite (u, s) en base b, pour tout tout entier u tel que t ≤ u. est un réseau (0, m, 1) en base b.

(2) La suite de Van der Corput en base b est une suite (0, 1) en base b.

Les premiers réseaux et suites en base 2 ont été proposés par Sobol [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. La généralisartion en base quelconque est due à Faure, d'abord en dimension 1 [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension un)[END_REF], puis en dimension quelconque [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF] ; voir aussi l'article récent [START_REF] Lemieux | New perspectives on (0, s)-sequences[END_REF]. x n := (x n,1 , x n,2 , . . . , x n,s ),

x n,i := ∞ j=1 y (i) n,j b -j , y (i) n,j := ∞ k=j-1 k j -1 b k-j+1 i a k (n) ∈ Z b et n = ∞ k=0 a k (n)b k , a k (n) ∈ Z b
est le développement en base b de n.

Suites de Niederreiter

H. Niederreiter a proposé dans [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] un procédé général de construction de suites-(t, s) en base q. On se limite ici au cas où q est une puissance première, c'est-à-dire q = p α avec p premier et α ≥ 1. On note F q le corps fini à q éléments. On choisit :

(S1) des bijections ψ r : Z q → F q pour tout r ≥ 0, telles que ψ r (0) = 0 pour r suffisamment grand ;

(S2) des bijections η i,j : F q → Z q pour tout 1 ≤ i ≤ s et j ≥ 1, telles que η i,j (0) = 0 pour tout 1 ≤ i ≤ s et j suffisamment grand ;

(S3) des éléments c

(i) j,r ∈ F q pour tout 1 ≤ i ≤ s, j ≥ 1 et r ≥ 0, tels que c (i) j,r = 0 pour tout 1 ≤ i ≤ s, tout r ≥ 0 et j suffisamment grand.
On note F q ((x -1 )) le corps des séries de Laurent formelles sur F q . Tout élément L ∈ F q ((x -1 )) s'écrit

L = ∞ k=w t k x -k , où w ∈ Z et t k ∈ F q pour tout k ≥ w. On choisit p 1 , . . . , p s ∈ F q [x
] des polynômes irréductibles 1 , unitaires 2 distincts. Des tableaux de polynômes irréductibles sont donnés dans les deux livres [START_REF] Lidl | Introduction to Finite Fields and Their Applications[END_REF][START_REF] Lidl | Finite Fields[END_REF]. On pose e i := deg(p i ) pour tout 1 ≤ i ≤ s. Pour 1 ≤ i ≤ s, j ≥ 1 et 0 ≤ k < e i , on écrit le développement en série de Laurent :

x k p i (x) j = ∞ r=0 a (i) (j, k, r)x -r-1 .
Ce développement permet de définir pour tout

1 ≤ i ≤ s, j ≥ 1 et r ≥ 0 des éléments c (i) j,r de F q : c (i) j,r := a (i) (Q(i, j) + 1, k(i, j), r), (1.3.1) où j -1 = Q(i, j)e i + k(i, j) avec 0 ≤ k(i, j) < e i .
On remarque que pour tout 1 ≤ i ≤ s et r ≥ 0 les éléments c (i) j,r sont nuls pour j suffisamment grand.

On peut maintenant définir la suite de Niederreiter S = {x n : n ≥ 0} dans

I s . Soit n ∈ N et n = ∞ r=0 a r (n)q r
son développement en base q, avec a r (n) ∈ Z q pour r ≥ 0. Pour tout n ≥ 0, 1 ≤ i ≤ s, soit

x (i) n := ∞ j=1 y (i) n,j b -j , où, pour tout j ≥ 1 y (i) n,j := η i,j ∞ r=0 c (i) j,r ψ r (a r (n)) ∈ Z b .
Les éléments x n sont alors définis par :

x n := (x (1) n , . . . , x (s) n ), n ≥ 0. Le résultat fondamental sur ces suites est le suivant. Théorème 1.4 La suite de Niederreiter est une suite (t(q, s), s) en base q, où t(q, s) := s i=1 (e i -1).

1 Un polynôme est irréductible si dans toute factorisation en deux polynômes, l'un est constant. 2 Un polynôme est unitaire si le coefficient du monôme de plus haut degré est égal à 1.

Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

Afin d'obtenir la meilleure majoration de discrépance, il faut minimiser la valeur de t(q, s) pour s et q fixés. Pour cela on range les polynômes irréductibles unitaires de F q [x] par degré croissant : p 1 , p 2 , . . . et on choisit les s premiers.

Cas particulier. Dans le cas où q est une puissance première et s une dimension ≤ q, on peut choisir des polynômes de la forme p i (x) = xb i pour tout 1 ≤ i ≤ s, où b 1 , . . . , b s sont des éléments distincts de F q . Avec ce choix, on a t(q, s) = 0.

La relation (1.3.1) s'écrit alors

c (i) j,r = a (i) (j, 0, r), 1 ≤ i ≤ s, j ≥ 1, r ≥ 0.
Ces éléments s'obtiennent en écrivant le développement en série de Laurent :

1 p i (x) j = 1 x j (1 -b i x -1 ) j = x -j ∞ r=0 r + j -1 j -1 b r i x -r = ∞ r=j-1 r j -1 b r-j+1 i x -r-1 .
Il en résulte que, pour

1 ≤ i ≤ s et j ≥ 1, on a c (i) j,r =    0 pour 0 ≤ r < j -1, r j -1 b r-j+1 i pour r ≥ j -1,
en convenant que 0 0 = 1 ∈ F q . Quand q est premier, on identifie Z q à F q et on choisit les bijections ψ r et η i,j égales à l'identité : on retrouve alors la construction de H. Faure [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF].

Les ensembles digitaux et les suites digitales généralisent les réseaux (t, m, s) et les suites (t, s) [START_REF] Tezuka | Uniform Random Numbers : Theory and Practice[END_REF][START_REF]Random and Quasi-Random Point Sets[END_REF][START_REF] Dick | Digital Nets and Sequences, Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. On présente dans ce qui suit leur principe de construction.

(1) Soient s ≥ 1, b ≥ 2 et k ≥ 1 trois entiers. On considère :

• un anneau commutatif unitaire R de cardinal b ;

• des bijections ψ r :

Z b → R, pour 0 ≤ r ≤ k -1 ;
• des bijections η j,l : R → Z b , pour 1 ≤ j ≤ s et 1 ≤ l ≤ k ;

• des matrices (dites génératrices) C 1 , . . . , C s de dimension k × k d'éléments de R.

Pour i = 0, . . . , b k -1, on écrit la représentation de i en base b :

i = k-1 r=0 a r b r , où a r ∈ Z b . Soit y = (ψ 0 (a 0 ), . . . , ψ k-1 (a k-1 )) t ∈ R k et (b j,1 , b j,2 , . . . , b j,k ) t = C j • y, avec b j,l ∈ R. Pour j = 1, . . . , s soit u i,j = η j,1 (b j,1 ) b + η j,2 (b j,2 ) b 2 + • • • + η j,ℓ (b j,k ) b k .
L'ensemble P = {u i = (u i,1 , . . . , u i,s ), i = 0, . . . , b k -1} est appelé réseau digital sur R en base b.

(2) Soient s ≥ 1 et b ≥ 2 deux entiers. On considère :

• un anneau commutatif unitaire R de cardinal b ;

• des bijections ψ r : Z b → R, pour r ≥ 0, vérifiant ψ r (0) = 0, pour r suffisamment grand ;

• des bijections η j,l : R → Z b , pour 1 ≤ j ≤ s et l ≥ 1 ;

• des matrices (dites génératrices) C 1 , . . . , C s d'éléments de R, d'indices dans N * ×N * .
Pour i ≥ 0, on écrit la représentation de i en base b :

i = ∞ r=0 a r b r , où a r ∈ Z b . Soit y = (ψ 0 (a 0 ), ψ 1 (a 1 ), . . .) t ∈ R N * et (b j,1 , b j,2 , . . .) t = C j • y, avec b j,l ∈ R. Pour j = 1, . . . , s soit u i,j = η j,1 (b j,1 ) b + η j,2 (b j,2 ) b 2 + • • • . La suite P = {u i = (u i,1 , . . . , u i,s ), i ≥ 0} est appelée suite digitale sur R en base b. Exemple 1.2 Soient s ≥ 1 et b ≥ 2,
puissance première supérieure à s. Les suites de Faure généralisées [START_REF] Tezuka | Uniform Random Numbers : Theory and Practice[END_REF] sont obtenues en utilisant la construction précédente dans laquelle les matrices génératrices sont de la forme :

C j = A j P j-1 , j = 1, . . . , s,
où P est la transposée de la matrice de Pascal : La suite de Faure originale est obtenue en prenant b le plus petit entier premier ≥ s et les matrices A j toutes égales à la matrice identité I.

P i,j := j -1 i -1 ∈ F b et A j est
Exemple 1.3 Une construction de suites, élaborée par Faure et Tezuka [START_REF] Faure | Another random scrambling of digital (t, s)-sequences[END_REF][START_REF] Tezuka | I-binomial scrambling of digital nets and sequences[END_REF] consiste à multiplier à droite les matrices génératrices des suites (t, s) par des matrices régulières triangulaires supérieures (NUT : non-singular upper triangular). Les nouvelles matrices génératrices sont :

C ′ j = C j U j , où les U j , j = 1, . . . , s sont des matrices NUT. Comme la multiplication par les matrices U j ne préserve pas les propriétes de répartition dans les intervalles élémentaires, la suite ainsi obtenue n'est généralement pas de type (t, s). Une possibilité pour conserver la propriété d'uniformité des suites est de choisir U j sous la forme :

U j = γ j U,
où U est une matrice NUT fixée et γ j ∈ F b , γ j = 0. Les matrices génératrices deviennent alors :

C ′ j = γ j C j U. On a le résultat suivant [START_REF] Faure | Another random scrambling of digital (t, s)-sequences[END_REF] : si une suite (t, s) en base b est engendrée par les matrices génératrices C j , 1 ≤ j ≤ s, alors la suite engendrée avec les matrices

C ′ j = γ j C j U , 1 ≤ j ≤ s est également une suite (t, s) en base b.

Hasardisation

Le but des techniques d'hasardisation est de construire des ensembles (ou suites) de points de faible discrépance, vérifiant :

(1) chaque point dans la suite hasardisée est uniformément distribué sur [0, 1) s , (2) la régularité de la suite hasardisée est celle de la suite déterministe de départ.

Méthodes de décalage linéaire

Les procédés les plus simples d'hasardisation sont les décalages linéaires.

Décalage aléatoire modulo 1

Cette méthode d'hasardisation est aussi connue sous le nom de rotation de Cranley-Patterson [START_REF] Cranley | Randomization of number theoretic methods for multiple integration[END_REF]. Soit P n := {u i ; i = 0, 1, . . . , n -1} un ensemble de points de [0, 1) s et ∆ un vecteur aléatoire s-dimensionnel uniformément distribué sur [0, 1) s . L'ensemble hasardisé P n := { u i ; i = 0, 1, . . . , n -1} est défini par :

u i = (u i + ∆) mod 1.
L'ensemble ainsi obtenu est uniformément distribué sur [0, 1) s .Cependant, cette technique d'hasardisation ne préserve pas les propriétés d'équirépartition de l'ensemble de départ.

Décalage digital b-adique

Soit P n := {u i ; i = 0, 1, . . . , n -1} un réseau (t, m, s) en base b. Une méthode analogue à la précédente consiste à écrire la représentation b-adique du vecteur ∆ et à additionner ses composantes à celles des points u i , en utilisant les opérations sur F b ( [START_REF] Matoušek | On the L 2 -discrepancy for anchored boxes[END_REF][START_REF] Hong | Algorithm 823 :implementing scrambled digital sequences[END_REF]). Plus précisément, si ∆ = (∆ 1 , . . . , ∆ s ) avec

∆ j = ∞ l=1 d j,l b -l , u i,j = ∞ l=1 u i,j,l b -l , on calcule u i ⊕ ∆ = ( u i,1 , . . . , u i,s ), où u i,j = ∞ l=1 ((u i,j,l + d j,l ) mod b)b -l .
L'ensemble hasardisé est alors P n = { u i ; i = 0, . . . , n -1}. Ce type de décalage digital est le plus approprié aux réseaux (t, m, s) puisqu'il préserve les propriétés d'uniformité et les valeurs du paramètre t.

Ensembles de points à faible discrépance brouillés

Cette méthode a été proposée par Owen dans [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF][START_REF] Owen | Monte Carlo variance of scrambled net quadrature[END_REF]. L'idée est de perturber les digits des ensembles de points à faible discrépance à l'aide de permuations aléatoires de digits, tout en préservant les propriétés d'équirépartition. Soit P = (u i ) i un ensemble de points de [0, 1) s . On considère la représentation b-adique de chacune des composantes du i ème terme u i = (u i,1 , . . . , u i,s ) de P :

u i,j = ∞ k=1 u i,j,k b -k , où 0 ≤ u i,j,k < b, ∀i, j, k.
La version hasardisée de P est l'ensemble P dont le i ème élément u i = ( u i,1 , . . . , u i,s ) est donné par :

u i,j = ∞ k=1 u i,j,k b -k ,
où les u i,j,k sont définis comme suit :

u i,j,1 = π j (u i,j,1 ) u i,j,2 = π j,u i,j,1 (u i,j,2 )
. . .

u i,j,k = π j,u i,j,1 ,...,u i,j,k-1 (u i,j,k ),
1. Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique les π j,u i,j,1 ,...,u i,j,ℓ : Z b → Z b étant des permutations aléatoires indépendantes uniformément distribuées sur l'ensemble des b! permutations possibles de {0, 1, . . . , b -1}. Notons que la permutation appliquée au coefficient u i,j,k dépend des valeurs des coefficients u i,j,1 , . . . , u i,j,k-1 . Ce qui fait que l'implémentation numérique de ces méthodes demande des espaces de stockage et des temps de calcul importants.

Dans [START_REF] Owen | Monte Carlo variance of scrambled net quadrature[END_REF] Owen a montré que si P est un réseau (t, m, s) (respectivement une suite (t, s)) en base b alors P est un réseau (t, m, s) (respectivement une suite (t, s)) en base b presque sûrement et que les points de l'ensemble P sont uniformément distribués sur [0, 1) s .

Des alternatives permettant la réduction de la place mémoire ont été proposées par Matoušek [START_REF] Matoušek | On the L 2 -discrepancy for anchored boxes[END_REF], Tezuka [START_REF] Tezuka | Uniform Random Numbers : Theory and Practice[END_REF], Hong et Hickernell [START_REF] Hong | Algorithm 823 :implementing scrambled digital sequences[END_REF]. L'idée est d'appliquer des permutations affines aux différents coefficients comme suit. Si n = b k , on considère :

• s matrices régulières triangulaires inférieures L 1 , . . . , L s , à indices dans N * × N * , dont les éléments sont choisis aléatoirement et indépendamment dans F b , avec des éléments diagonaux non nuls ;

• s vecteurs (à indice dans N * ) d 1 , . . . , d s dont les composantes sont indépendantes et uniformement distribuées dans F b .
Les coefficients u i,j,1 , u i,j,2 . . . de la j ème composante u i,j de l'élément u i sont donnés par :

   u i,j,1 u i,j,2 . . .    = L j    u i,j,1 u i,j,2 . . .    + d j , (1.4.1)
toutes ces opérations étant effectuées dans F b . Le décalage digital garantit que chaque point obtenu est uniformément distribué dans [0, 1) s [START_REF] Owen | Variance with alternative scrambling of digital nets[END_REF]. Cette technique préserve les propriétés d'équirépartition de l'ensemble original. En effet, dans [START_REF] Matoušek | Geometric Discrepancy : An Illustrated Guide[END_REF] Matousěk a montré que la valeur du paramètre t de la suite ainsi obtenue ne dépasse pas celle de la suite originale, le brouillage peut donc potentiellement améliorer la qualité de la suite. Dans nos essais numériques, nous nous sommes intéressés tout particulièrement à trois types d'hasardisation par permutations affines. Elles diffèrent par le choix des matrices L j :

• La matrice de brouillage linéaire aléatoire de Matoušek (Matousěk scrambling) [START_REF] Matoušek | On the L 2 -discrepancy for anchored boxes[END_REF] :

      h 11 0 0 0 . . . g 21 h 22 0 0 . . . g 31 g 32 h 33 0 . . . g 41 g 42 g 43 h 44 . . . . . . . . . .      
• La matrice de brouillage I-binomial de Tezuka (Tezuka scrambling) [START_REF] Tezuka | I-binomial scrambling of digital nets and sequences[END_REF] :

      h 1 0 0 0 . . . g 2 h 1 0 0 . . . g 3 g 2 h 1 0 . . . g 4 g 3 g 2 h 1 . . . . . . . . . .      
• La matrice de brouillage d'Owen (Owen scrambling) [START_REF] Owen | Variance with alternative scrambling of digital nets[END_REF] :

      h 1 0 0 0 . . . h 1 h 2 0 0 . . . h 1 h 2 h 3 0 . . . h 1 h 2 h 3 h 4 . . . . . . . . . .      

Essais numériques

Il s'agit dans cette partie de valider les approches quasi-Monte Carlo hasardiseés et de les comparer aux schémas de Monte Carlo et quasi-Monte Carlo traditionnels dans des essais numériques.

Nous nous limitons au calcul d'intégrales multi-dimensionnelles dont la valeur exacte est connue par les trois méthodes citées ci-dessus. Nous traçons les courbes d'erreur puis évaluons les valeurs des variances, pour différentes dimensions. La différence entre les trois méthodes est due à la nature des points utilisés dans les simulations.

• Pour les essais du type Monte Carlo, nous utilisons le générateur de nombres pseudoaléatoires MRG32k3a, mis au point par L'Ecuyer [START_REF] L'ecuyer | Good parameters and implementations for combined multiple recursive random number generators[END_REF].

• La suite déterministe utilisée dans les essais QMC est celle de Faure [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF].

• Pour les suites QMC hasardisées, nous utilisons la suite de Faure hasardisée par les méthodes de Matoušek, I-binomiale et d'Owen détaillées au paragraphe précédent.

Le nombre N de points considérés est une puissance de la base b où b est le plus petit entier supérieur à la dimension s.

Les intégrales calculées sont les suivantes :

J s = [0,1] s 1 s s i=1 |4x i -2|dx 1 . . . dx s , I = [0,1] 4 (e x 1 x 2 x 3 x 4 -1)dx 1 dx 2 dx 3 dx 4 .
Les valeurs exactes de ces intégrales sont :

J s = 1, ∀s ≥ 1 et I = 5 4 -2 log 2 + π 2 24 .
On calcule l'erreur, qui est la valeur absolue de la différence entre la valeur exacte et la valeur approchée.

Estimation de J s

On calcule des valeurs approchées de l'intégrale J s pour différentes valeurs de s. Les résultats sont représentés sur les Figures 1.2 Les résultats montrent que les résultats obtenus par les méthodes QMC et QMC hasardisées sont nettement meilleurs que ceux obtenus par la méthode MC.

Pour estimer la variance de l'estimateur des approches MC et QMC hasardisées, nous répétons le calcul de l'intégrale J s de façon indépendante M = 100 fois pour chaque dimension. Nous obtenons ainsi des copies indépendantes : J s,1 , J s,2 , . . . , J s,100 de l'estimateur. La variance est alors estimée par la variance empirique de l'échantillon ainsi obtenu. Les résultats sont donnés dans la Table 1.2 (on a choisi le plus grand nombre N de points).

Il est clair que la substitution des suites aléatoires par des suites hasardisées permet de réduire la variance.

Estimation de I

De même, l'intégrale I calculée par les différentes méthodes montre une nette amélioration due à l'utilisation des suites hasardisées par rapport aux résultats obtenus avec la suite de Faure ou une approximation MC traditionnelle (voir Figure 1.6).

De même, la variance des estimateurs obtenus par les méthodes QMC hasardisées est plus petite que celle de la méthode MC, comme le montre la Table 1.3. On a choisi le plus grand nombre N de points et on a calculé la variance échantillonnale, pour M = 100 répétitions indépendantes.

Conclusion

Dans ce chapitre, on a présenté les méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique. Ensuite, on a décrit des méthodes d'hasardisation des suites déterministes. On a testé l'efficacité de ces dernières méthodes dans des essais numériques d'évaluation d'intégrales multi-dimensionnelles en les comparant aux méthodes de Monte Carlo et quasi-Monte Carlo traditionnelles.

On a remarqué que les versions hasardisées des méthodes QMC permettent d'améliorer les résultats. En fait les erreurs obtenues sont généralement plus petites que celles obtenues par une méthode MC classique, de même que les variances des estimateurs correspondants.

Nous rappelons que le développement des méthodes de Monte Carlo et quasi-Monte Carlo peut être suivi dans deux séries de biennales. La première est celle des conférences internationales Monte Carlo and Quasi-Monte Carlo Methods [START_REF]Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing[END_REF][START_REF] Niederreiter | Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF] Niederreiter | Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo Methods[END_REF][START_REF]Monte Carlo and Quasi-Monte Carlo Methods[END_REF]. La seconde est celle des congrès IMACS Seminar on Monte Carlo Methods [START_REF]Selection of papers presented at the IMACS Seminar on Monte Carlo Methods[END_REF][START_REF] Dimov | IMACS Sponsored Special Issue : The Second IMACS Seminar on Monte Carlo Methods[END_REF][START_REF]Special Issue : 3rd IMACS Seminar on Monte Carlo Methods[END_REF][START_REF] Sabelfeld | Selection of papers presented at the International Conference IV IMACS Seminar on Monte Carlo Methods[END_REF][START_REF] Mascagni | Special Issue : Fifth IMACS Seminar on Monte Carlo Methods[END_REF][START_REF] Alexandrov | Special Issue : The Sixth IMACS Seminar on Monte Carlo Methods[END_REF][START_REF]Selected papers from the Seventh IMACS Seminar on Monte Carlo Methods[END_REF][START_REF]Monte Carlo Methods and Applications, Proceedings of the 8th IMACS Seminar on Monte Carlo Methods[END_REF].

Chapitre 2

Les méthodes de Monte Carlo stratifiées et LHS

Introduction

Dans ce chapitre nous commençons par introduire les méthodes de Monte Carlo stratifiées. Ensuite nous nous intéressons à un cas particulier de stratification pour évaluer des intégrales et plus précisément, les intégrales de fontions indicatrices, car certaines méthodes de simulation peuvent être interprétées comme des successions de quadratures multi-dimensionnelles : voir par exemple [START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF][START_REF] Lécot | Quasi-Monte Carlo methods for estimating transient measures of discrete time Markov chains[END_REF][START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF]. Nous proposons donc un estimateur de Monte Carlo stratifié (MCS) simple et nous analysons l'ordre de convergence de la variance pour le calcul d'un volume multi-dimensionnel. Ensuite nous présentons la méthode d'échantillonnage par hypercube latin (LHS : Latin Hypercube Sampling). Nous majorons la variance de l'estimateur du volume d'un rectangle à côtés parallèles aux axes de coordonnées. Nous présentons les résultats de quelques essais numériques qui permettent d'une part de valider l'étude théorique faite et d'autre part de comparer l'efficacité de ces méthodes à celle du MC classique. Une partie des résultats de ce chapitre a été publiée dans [START_REF] El Haddad | Stratified Monte Carlo integration[END_REF].

Méthode de Monte Carlo stratifiée

Soit s ≥ 1 un entier, λ s la mesure de Lebesgue sur R s et I s := [0, 1) s l'hypercube unité s-dimensionnel. On considère une fonction de carré intégrable f : I s → R dont on veut évaluer l'intégrale

J := I s f (x)dλ s (x).
La méthode de Monte Carlo usuelle consiste à considérer une famille indépendante de variables aléatoires {U 1 , . . . , U N } de loi uniforme sur I s . L'estimateur de Monte Carlo simple de J est alors

X := 1 N N k=1 f • U k .
C'est un estimateur sans biais puisque

E[X] = 1 N N k=1 E[f • U k ] = 1 N N k=1 I s f (x)dλ s (x) = J . Par ailleurs, la famille {f • U 1 , . . . , f • U N } étant indépendante, on a Var(X) = 1 N 2 N k=1 Var(f • U k ), avec
Var(f • U k ) = Ω f • U k (ω) 2 dP(ω) - Ω f • U k (ω)dP(ω) 2 = I s (f (x)) 2 dλ s (x) - I s f (x)dλ s (x) 2 . D'où Var(X) = 1 N I s (f (x)) 2 dλ s (x) - 1 N I s f (x)dλ s (x) 2 = σ 2 (f ) N , (2.2.1)
où σ 2 (f ) est associé à la probabilité 1 I s (x)λ s par (1.2.1).

Stratification simple

Le principe de la stratification consiste à subdiviser le domaine d'intégration I s en des sous-ensembles et à engendrer des variables aléatoires dans chacun d'eux, ce qui garantit une meilleure répartition des points. Dans la méthode que nous allons présenter, nous engendrons une variable aléatoire dans chacun des sous-domaines. Soit {C 1 , . . . , C N } une partition de I s , c'est-à-dire une famille de sous-ensembles de I s vérifiant :

I s = N k=1 C k et ∀k, ℓ ∈ {1, 2, . . . , N } k = ℓ ⇒ C k ∩ C ℓ = ∅.
On suppose de plus que ces sous-ensembles ont la même mesure :

λ s (C 1 ) = • • • = λ s (C N ) = 1 N .
Pour 1 ≤ k ≤ N , soit V k une variable aléatoire de loi uniforme sur C k . On suppose que la famille {V 1 , . . . , V N } est indépendante. L'estimateur de Monte Carlo stratifié simple de J est défini par

Y := 1 N N k=1 f • V k .
Il s'agit d'un estimateur sans biais de J . En effet

E[Y ] = 1 N N k=1 E[f • V k ] = 1 N N k=1 Ω f • V k (ω)dP(ω) = 1 N N k=1 1 λ s (C k ) R s f (x)1 C k (x)dλ s (x) = N k=1 C k f (x)dλ s (x) = J .
Cette technique de stratification a été initialement proposée par Haber [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF] et analysée dans [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF] et [START_REF] Cheng | The problem of dimensionality in stratified sampling[END_REF]. On a d'abord (voir plus loin) : pour toute fonction f de carré intégrable

Var(Y ) ≤ Var(X).
Dans [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF], Haber a introduit la notion de régularité d'intervalle. Soit

H = s i=1 [a i , b i )
un intervalle dans I s ; son module de régularité est défini par

ρ(H) = s i=1 (b i -a i ) max 1≤i≤s (b i -a i ) s .
Une famille (H ℓ ) ℓ d'intervalles est dite régulière s'il existe une valeur ρ 0 telle que

∀ℓ ρ(H ℓ ) ≥ ρ 0 .
Il a alors établi les deux résultats suivants.

Proposition 2.1 Soit f une fonction continue sur I s ; si la famille

{C k : 1 ≤ k ≤ N, N ∈ N * }
des intervalles de la partition de I s , pour tous les N , est régulière, alors

lim N →∞ Var(Y ) Var(X) = 0.
Proposition 2.2 Soit f une fonction dérivable sur I s telle que Df soit continue et bornée sur

I s ; si la famille {C k : 1 ≤ k ≤ N, N ∈ N * } est régulière, alors il existe deux constantes positives K 1 et K 2 (indépendantes de f ) telles que pour tout ǫ > 0 K 1 -ǫ √ 12 I s Df (x) 2 2 dλ s (x) 1/2 1 N 1/2+1/s ≤ Var(Y ) ≤ K 2 + ǫ √ 12 I s Df (x) 2 2 dλ s (x) 1/2 1 N 1/2+1/s ,
pour N suffisamment grand. Dans le cas de partitions avec des hypercubes,

K 1 = K 2 = 1.
Cheng et Davenport ont établi dans [START_REF] Cheng | The problem of dimensionality in stratified sampling[END_REF] un résultat analogue :

Proposition 2.3 Soit f une fonction lipschitzienne sur I s ; si la famille {C k : 1 ≤ k ≤ N } est formée d'hypercubes, alors Var(Y ) = O 1 N 1+2/s .

Nous allons d'abord comparer les variances des estimateurs

X et Y dans un cas gé- néral (en supposant f de carré intégrable). Puisque la famille {f • V 1 , . . . , f • V N } est indépendante, on a Var(Y ) = Var 1 N N k=1 f • V k = 1 N 2 N k=1 Var(f • V k ).
Or

Var(f • V k ) = Ω (f • V k (ω)) 2 dP(ω) - Ω f • V k (ω)dP(ω) 2 = 1 λ s (C k ) R s (f (x)) 2 1 C k (x)dλ s (x) - 1 λ s (C k ) R s f (x)1 C k (x)dλ s (x) 2 = N C k (f (x)) 2 dλ s (x) -N 2 C k f (x)dλ s (x) 2 .
On a, d'après l'inégalité de Cauchy-Schwarz,

I s f (x)dλ s (x) 2 ≤ N N k=1 C k f (x)dλ s (x) 2 , et par suite Var(Y ) = 1 N N k=1 C k (f (x)) 2 dλ s (x) - N k=1 C k f (x)dλ s (x) 2 = 1 N I s (f (x)) 2 dλ s (x) - N k=1 C k f (x)dλ s (x) 2 ≤ 1 N I s (f (x)) 2 dλ s (x) - 1 N I s f (x)dλ s (x) 2 . ( 2
Var(Y ) ≤ Var(X).
La méthode de Monte Carlo stratifiée proposée par Haber [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF] est donc bien une technique de réduction de la variance.

Estimation de l'intégrale d'une indicatrice

Dans ce paragraphe, on cherche à majorer la variance de l'estimateur MCS dans le cas où l'intégrande est une fonction indicatrice. On se limitera à considérer un échantillonnage stratifié particulier dans lequel les sous-domaines sont des pavés identiques : on appelle cette méthode la méthode de Monte Carlo stratifiée simple. Soit A un sous-ensemble mesurable de I s . On restreint l'étude au cas où f := 1 A . Il est clair que dans ce cas J = λ s (A). On note X l'estimateur de Monte Carlo de J :

X := 1 N N k=1 1 A • U k , où N ∈ N * et
Var(X) = 1 N λ s (A)(1 -λ s (A)) ≤ 1 4N .
On veut majorer la variance d'un cas particulier de l'estimateur de Monte Carlo stratifié de Haber [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF]. On suppose qu'il existe un entier n

∈ N * tel que N = n s . Pour k = (k 1 , . . . , k s ) ∈ N s avec 1 ≤ k i ≤ n, on note C k = s i=1 k i -1 n , k i n et on considère une variable aléatoire V k uniformément distribuée dans C k . On suppose que la famille {V k : k = (k 1 , . . . , k s ) ∈ N s , 1 ≤ k i ≤ n} est indépendante : cf.
Y := 1 N k 1 A • V k où k parcourt l'ensemble {k = (k 1 , . . . , k s ) ∈ N s : 1 ≤ k i ≤ n}.
Pour calculer la variance de Y , on impose une condition de régularité sur le domaine A. On introduit les deux ensembles suivants (cf. Figure 2.2) : pour tout réel ε > 0, soit

A ε := {x ∈ I s : ∃y ∈ A ; x -y ∞ < ε} et A -ε := {x ∈ I s : ∀y ∈ I s \A ; x -y ∞ ≥ ε}.
Ces notions ont été introduites par Niederreiter et Wills [START_REF] Niederreiter | Diskrepanz und Distanz von Maßen bezüglich konvexer und Jordanscher Mengen[END_REF] et reprises par Lécot [START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF], dans un contexte quasi-Monte Carlo dans les deux cas. On a évidemment les inclusions

A -ε ⊂ A ⊂ A ε . 0 0.
(A ε \ A), λ s (A \ A -ε )) ≤ γ(ε), alors Var(Y ) ≤ 1 2N γ 1 N 1/s . Preuve. Soit K = {k = (k 1 , . . . , k s ) ∈ N s : 1 ≤ k i ≤ n, C k ∩ A = ∅ et C k ⊂ A}. Pour k ∈ K on a C k ⊂ A 1/n \ A -1/n . En effet, -d'une part, si x ∈ C k , comme C k ∩ A = ∅, il existe y ∈ A tel que x -y ∞ < 1 n et donc x ∈ A 1/n ; -d'autre part, si x ∈ C k , comme C k ⊂ A, il existe y ∈ I s avec y ∈ C k et y / ∈ A donc x -y ∞ < 1 n et y ∈ I s \ A et donc x / ∈ A -1/n . Par conséquent k∈K C k ⊂ A 1/n \ A -1/n .
Il en résulte que

1 N #K ≤ 2γ 1 n , où #K désigne le cardinal de K.
La famille des variables aléatoires V k étant indépendante, on a

Var(Y ) = 1 N 2 k Var(1 A • V k ) = 1 N 2 k N λ s (C k ∩ A)(1 -N λ s (C k ∩ A)).
En remarquant que le terme

λ s (C k ∩ A)(1 -N λ s (C k ∩ A)) est nul si C k ∩ A = ∅ ou si C k ⊂ A, la variance de Y s'écrit Var(Y ) = 1 N 2 k∈K N λ s (C k ∩ A)(1 -N λ s (C k ∩ A)). (2.2.3) D'où, en remarquant que 0 ≤ N λ s (C k ∩ A) ≤ 1 et que la fonction x ∈ [0, 1] → x(1 -x) atteint son maximum de 1/4 en 1/2, Var(Y ) ≤ 1 4N 2 #K, et le résultat suit.
Remarque 2.1 Sous les hypothèses de la proposition précédente, si γ est une fonction linéaire, on obtient :

Var(Y ) = O 1 N 1+1/s ,
ce qui montre un gain de N -1/s par rapport à la méthode MC classique.

Dans [6], Cheng et Davenport ont étudié le cas particulier où

A = A f,y := {x ∈ I s : f (x) ≤ y}, avec f : I s → R et y ∈ R. On note S f,y := {x ∈ I s : f (x) = y}.
Ils ont alors établi le résultat suivant.

Var(Y ) = O 1 N 1+1/s .
Le résultat reste valable si l'on suppose que S f,y est une union finie de sous-ensembles S ℓ ayant les propriétés indiquées dans les hypothèses de la proposition.

Parce qu'ils apparaissent dans les simulations Monte Carlo d'équations, on considère dans cette étude des ensembles dont la frontière peut être paramétrée par une fonction régulière. Soit f : I s-1 → I une fonction borélienne et

A f := {(x ′ , x s ) ∈ I s : x s < f (x ′ )} (cf. Figure 2.
3). On veut estimer l'intégrale

I = I s 1 A f (x)dλ s (x) = I s-1 f (x ′ )dλ s-1 (x ′ ).
L'objectif est d'établir une majoration de la variance de l'estimateur MCS de I semblable à celui de la proposition 2.4 mais avec une condition sur la fonction f et non sur l'ensemble A f . On utilise la notion de variation au sens de Hardy-Krause de f , introduite dans le premier chapitre.

On a le résultat suivant [START_REF] Lécot | Error bounds for quasi-Monte Carlo integration with nets[END_REF].

0 1 1 A f Figure 2.3 -Graphe d'une fonction f et domaine A f correspondant en dimension s = 2.
C k = s i=1 k i -1 n i , k i n i . Pour chaque k, soient y k , z k ∈ C k . Alors k |f (z k ) -f (y k )| ≤ V (f ) s i=1 n i s i=1 1 n i . Preuve. Si y ∈ Īs soit T i y f la restriction de f à l'hyperplan x i = y i . Si y, z ∈ Īs soit ∆ i y,z f := T i z f -T i y f. Si K = {i 1 , . . . , i α } ⊂ [1, s], on pose T K y f = T i 1 y . . . T iα y f et ∆ K y,z f = ∆ i 1 y,z . . . ∆ iα y,z f. On note T y f = T [1,s] y f . Soit 1 := (1, . . . , 1), k+ := (k 1 + 1, . . . , k s + 1), x k := k 1 -1 n 1 , . . . , k s -1 n s .

Les méthodes de Monte Carlo stratifiées et LHS

En combinant les identités suivantes :

f (z k ) -f (y k ) = s i=1 ∆ i y k ,z k T {1,...,i-1} y k T {i+1,...,s} z k f, T w f = s α=0 (-1) α K⊂[1,s] #K=α ∆ K w,x k+ T K c x k+ f, pour w ∈ C k , T x k+ f = s α=0 (-1) α K⊂[1,s] #K=α k i <ℓ i ≤n i i∈K ∆ K x ℓ ,x ℓ+ T K c 1 f.
On obtient le résultat annoncé. Une majoration de la variance de l'estimateur MCS de I est la suivante.

Proposition 2.6 Soit f : I s-1 → I une fonction à variation bornée au sens de Hardy et

Krause et {V k : 1 ≤ k i ≤ n} une famille indépendante de N variables aléatoires, où V k est de loi uniforme sur C k . Soit Y := 1 N k 1 A f • V k . Alors Var(Y ) ≤ s -1 4 V (f ) + 1 2 1 N 1+1/s . Preuve. Pour k = (k 1 , . . . , k s ) ∈ N s avec 1 ≤ k i ≤ n, on note k ′ = (k 1 , . . . , k s-1 ) et C ′ k ′ := s-1 i=1 k i -1 n , k i n .
Un calcul analogue à celui effectué pour établir l'identité (2.2.3) donne

Var(Y ) ≤ 1 4N 2 k ′ #{k s : C (k ′ ,ks) ∩ A f = ∅ et C (k ′ ,ks) ⊂ A f }.
On examine les deux cas suivants :

• si C (k ′ ,ks) ∩ A f = ∅, alors il existe x ′ k ′ ∈ C ′ k ′ tel que k s < nf (x ′ k ′ ) + 1 ; • si C (k ′ ,ks) ⊂ A f , alors il existe y ′ k ′ ∈ C ′ k ′ tel que nf (y ′ k ′ ) < k s . Par conséquent Var(Y ) ≤ 1 4N 2 k ′ (n |f (x ′ k ′ ) -f (y ′ k ′ )| + 2) . Le Lemme 2.1 donne k ′ |f (x ′ k ′ ) -f (y ′ k ′ )| ≤ (s -1)V (f )n s-2 ,
ce qui termine la preuve. La variance de l'estimateur de Monte Carlo stratifié est d'ordre 1/N 1+1/s ce qui représente un gain en précision de facteur 1/N 1/s par rapport à l'estimateur de Monte Carlo simple. Ce gain diminue en dimension élevée, ce qui limite l'application de cette technique de stratification simple à des problèmes de dimension modérée, par exemple pour des simulations nécessitant l'évaluation d'intégrales en dimension s comprise entre 2 et 6 : voir par exemple [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | Low discrepancy sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo scheme using nets for a linear Boltzmann equation[END_REF][START_REF] Lécot | Comparison of quasi-Monte Carlo-based methods for the simulation of Markov chains[END_REF][START_REF] Lécot | Quasi-Monte Carlo methods for estimating transient measures of discrete time Markov chains[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of discrete-time Markov chain on multidimensional state space[END_REF][START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] El Haddad | Quasi-Monte Carlo methods for Markov chains with continuous multi-dimensional state space[END_REF], où des simulations QMC sont proposées et [START_REF] L'ecuyer | Randomized quasi-Monte Carlo simulation of Markov chains with an ordered state space[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF][START_REF] L'ecuyer | On array-RQMC for Markov chains : mapping alternatives and convergence rates[END_REF] pour des simulations QMC hasardisées. C'est aussi le cas des simulations détaillées aux chapitres 4 et 5. Pour des résultats optimaux en dimension plus élevée il est en général intéressant de voir si les intégrandes rencontrés peuvent suggérer une stratification partielle du domaine c'est-à-dire une partition de l'hypercube unité I s le long de quelques coordonnées parmi les s.

Une majoration de la variance du même ordre qu'ici avait été donné dans [START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF]. Le cadre de cette dernière étude était différent à deux points de vue. D'abord l'analyse théorique était limitée à la dimension s = 2 ; ensuite cette analyse était faite dans le cadre de la simulation des chaînes de Markov.

Essais numériques

Dans ce paragraphe, nous validons numériquement les analyses faites précédemment. Nous nous limitons au calcul d'intégrales de fonctions indicatrices en dimensions 2, 3 et 4.

On considère le sous-ensemble de la boule unité s-dimensionnelle

Q := {x ∈ I s : x 2 < 1},
où x 2 désigne la norme euclidienne de x. Dans cet exemple, la valeur exacte de l'intégrale

J = I s 1 Q (x)dλ s (x) = λ s (Q) est connue et vaut λ s (Q) = π s/2 2 s Γ s 2 + 1
, ce qui permet de calculer les erreurs des estimations. On a

Q := {x = (x ′ , x s ) ∈ I s : x s < f (x ′ )}, où f : x ′ = (x 1 , . . . , x s-1 ) → 1 -x 2 1 -• • • -x 2 s-1 .

Courbes d'erreur

On estime la valeur de J par une méthode de Monte Carlo usuelle (MC), la méthode stratifiée simple (MCS) décrite au paragraphe précédent et une méthode quasi-Monte Carlo (QMC) utilisant une suite de Faure [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF]. Les résultats sont représentés sur la Figure 2.4. On remarque qu'ils sont extrêmement bruités. L'erreur de chaque méthode est la valeur absolue de la différence entre la valeur exacte et la valeur approchée. En supposant que les erreurs sont Err = O(N -α ), on détermine les valeurs de α par ajustement linéaire en utilisant la méthode des moindres carrés (les droites sont tracées sur la Figure 2.4). Ces valeurs sont données dans la Table 2.1.

On remarque que les erreurs de la méthode stratifiée simple sont du même ordre que celles de la méthode QMC et sont plus petites que celles de la méthode de Monte Carlo usuelle.

Calcul des variances

On évalue ensuite les variances des estimateurs de J calculés par une méthode de Monte Carlo usuelle ou une méthode de Monte Carlo stratifiée simple. Pour un nombre N de points aléatoires utilisés, on calcule la variance empirique d'un échantillon de M observations de J obtenues en répétant le calcul de façon indépendante.

On a représenté les valeurs des variances des estimateurs MC et MCS en fonction du nombre de points N de la simulation sur la Figure 2.5. On a tracé les courbes des variances pour M = 100 et M = 1 000 répétitions : on voit que la différence des résultats est assez faible.

En supposant que la variance est Var = O(N -β ), on détermine les valeurs de β en effectuant un ajustement linéaire pour M = 1 000. Ces valeurs sont regroupées dans la Table 2.2. On a ajouté la borne théorique donnée par la Proposition 2.4.

Il est clair que l'estimateur de la méthode MCS est meilleur que celui de la méthode MC usuelle, puisque la variance est toujours inférieure et l'ordre β toujours supérieur. De plus, les valeurs de β obtenues dans les simulations sont très voisines des ordres théoriques établis dans les Propositions 2.4 et 2.6, ce qui montre que les bornes démontrées sont optimales (en ce qui concerne l'ordre). Supposant que la variance s'exprime par Var = cN -β (avec c constante) on voit que la constante de la méthode MCS est inférieure à celle de la méthode MC, dans cet exemple.

Méthodes LHS

La méthode d'échantillonnage par hypercube latin (LHS : Latin Hypercube Sampling) est une méthode de stratification qui privilégie les projections uni-dimensionnelles des échantillons. Elle est souvent appliquée dans l'analyse d'incertitude.

Cette technique a été décrite dans [START_REF] Mckay | A comparaison of three methods for selecting values of imput variables in the analysis of output from a computer code[END_REF] et analysée dans [START_REF] Stein | Large sample properties of simulations using Latin hypercube sampling[END_REF][START_REF] Owen | Monte Carlo variance of scrambled net quadrature[END_REF]. Une grille carrée contenant des points d'un échantillon est dite un carré latin s'il n'y a qu'un seul point de l'échantillon dans chaque ligne et chaque colonne. Le principe d'échantillonnage par hypercube latin est une généralisation du concept précédent pour une dimension s quelconque : il y a un point unique de l'échantillon de cardinal N dans chacun des intervalles

I i-1 × ℓ -1 N , ℓ N × I s-i , pour 1 ≤ i ≤ s et 1 ≤ ℓ ≤ N .
La Figure 2.6 représente un tel échantillon en dimension s = 2. 1. Soient σ 1 , . . . , σ s des permutations aléatoires indépendantes de l'ensemble N , de loi uniforme sur S N :

Construction d'ensembles LHS en dimension

1 N ! σ i ∈S N δ σ i
(δ ω est la mesure de Dirac en ω).

2. Soient {(ξ ℓ,1 , . . . , ξ ℓ,s ) : ℓ ∈ N } une famille indépendante de N variables aléatoires uniformément distribuées dans I s .

On suppose de plus que toutes ces variables aléatoires sont indépendantes. Pour ℓ ∈ N , on pose

W ℓ = σ 1 (ℓ) + ξ ℓ,1 -1 N , σ 2 (ℓ) + ξ ℓ,2 -1 N , . . . , σ s (ℓ) + ξ ℓ,s -1 N . Pour 1 ≤ i ≤ s et m i ∈ N , on note
I m i := m i -1 N , m i N et pour m := (m 1 , m 2 , .
. . , m s ) ∈ N s , on note

I m := I m 1 × I m 2 × • • • I ms .
W ℓ,i ∈ m -1 N , m N .
2. La variable W ℓ est de loi uniforme sur I s .

Preuve.

(1) On a

W ℓ,i = σ i (ℓ) + ξ ℓ,i -1 N . Alors, pour 1 ≤ m ≤ N , W ℓ,i ∈ m -1 N , m N ⇔ σ i (ℓ) + ξ ℓ,i -1 N ∈ m -1 N , m N ⇔ m -1 N ≤ σ i (ℓ) + ξ ℓ,i -1 N < m N ⇔ m -1 ≤ σ i (ℓ) + ξ ℓ,i -1 < m ⇔ m ≤ σ i (ℓ) + ξ ℓ,i < m + 1 ⇔ σ i (ℓ) = m, car 0 ≤ ξ ℓ,i < 1. Donc W ℓ,i ∈ m -1 N , m N ⇔ ℓ = σ -1 i (m).
D'où le premier résultat.

(2) Soit A ⊂ I s un ensemble borélien. On a

P({W ℓ ∈ A}) = Ω 1 {ω : W ℓ (ω)∈A} (ω)dP(ω) = Ω 1 ω : σ 1 (ω)(ℓ)+ξ ℓ,1 (ω)-1 N ,..., σs(ω)(ℓ)+ξ ℓ,s (ω)-1 N ∈A (ω)dP(ω) = 1 (N !) s σ 1 ∈S N • • • σs∈S N I s 1 (u 1 ,...,us) : σ 1 (ℓ)+u 1 -1 N ,..., σs(ℓ)+us-1 N ∈A (u)du = 1 (N !) s m 1 ∈N • • • ms∈N σ 1 ∈S N σ 1 (ℓ)=m 1 • • • σs∈S N σs(ℓ)=ms I s 1 (u 1 ,...,us) : m 1 +u 1 -1 N ,..., ms+us-1 N ∈A (u)du = 1 N s m 1 ∈N • • • ms∈N I s 1 (u 1 ,...,us) : m 1 +u 1 -1 N ,..., ms+us-1 N ∈A (u)du, car #{σ i ∈ S N : σ i (ℓ) = m i } = (N -1)!
En faisant les changement de variables :

v i := m i + u i -1 N , 1 ≤ i ≤ s, on obtient P({W ℓ ∈ A}) = 1 N s m 1 ∈N • • • ms∈N N s Im 1 A (v)dv = I s 1 A (v)dv,
ce qui prouve le second résultat. Soit f : I s → R une fonction mesurable de carré intégrable. On veut estimer

J = I s f (x)dλ s (x)
par la méthode LHS.

Proposition 2.8 L'estimateur LHS Z = 1 N ℓ∈N f • W ℓ
est un estimateur sans biais de J .

Preuve. Pour ℓ ∈ N , on a

E[f • W ℓ ] = Ω f σ 1 (ω)(ℓ) + ξ ℓ,1 (ω) -1 N , . . . , σ s (ω)(ℓ) + ξ ℓ,s (ω) -1 N dP(ω) = 1 (N !) s σ 1 ∈S N • • • σs∈S N I s f σ 1 (ℓ) + u 1 -1 N , . . . , σ s (ℓ) + u s -1 N du = 1 (N !) s m 1 ∈N • • • ms∈N σ 1 ∈S N σ 1 (ℓ)=m 1 • • • σs∈S N σs(ℓ)=ms I s f m 1 + u 1 -1 N , . . . , m s + u s -1 N du = 1 N s m 1 ∈N • • • ms∈N I s f m 1 + u 1 -1 N , . . . , m s + u s -1 N du. car #{σ i ∈ S N : σ i (ℓ) = m i } = (N -1)
! En faisant les changements de variables :

v i := m i + u i -1 N , on obtient E[f • W ℓ ] = m 1 ∈N • • • ms∈N Im f (u)du = I s f (u)du = J .
L'espérance de Z s'écrit

E[Z] = 1 N ℓ∈N E[f • W ℓ ] = J ,
donc Z est un estimateur sans biais de J .

Variance de l'estimateur LHS

McKay, Beckman et Conover ont introduit dans [START_REF] Mckay | A comparaison of three methods for selecting values of imput variables in the analysis of output from a computer code[END_REF] l'échantillonnage LHS comme une variante de l'échantillonnage stratifié. Ils ont prouvé le résultat suivant.

Proposition 2.9 Si f est monotone par rapport à chacune de ses variables, alors

Var(Z) ≤ Var(X).
Stein a montré dans [START_REF] Stein | Large sample properties of simulations using Latin hypercube sampling[END_REF] que la variance de l'estimateur LHS est asymptotiquement (N → ∞) inférieure à la variance de l'estimateur MC simple.

Proposition 2.10 Si f est de carré intégrable, alors

Var(Z) ≤ Var(X) + o 1 N .
Une proposition d'Owen [START_REF] Owen | Monte Carlo variance of scrambled net quadrature[END_REF] prouve que la variance d'un échantillon LHS à N points ne peut jamais être supérieure à la variance d'un échantillon MC simple à N -1 points.

Proposition 2.11 Si f est de carré intégrable, alors

Var(Z) ≤ σ 2 (f ) N -1 = N N -1
Var(X).

Nous allons ici améliorer la majoration de la variance de l'estimateur LHS, dans le cas simple où f = 1 A (la différence est négligeable quand N est grand). On supposera N ≥ 2.

Proposition 2.12 Si A ⊂ I s est un ensemble borélien, la variance de l'estimateur LHS

Z = 1 N ℓ∈N 1 A • W ℓ vérifie Var(Z) ≤ 1 N 1 4 + s N -1 N .
Preuve. On a

Var(Z) = 1 N 2 ℓ∈N Var(1 A • W ℓ ) + 1 N 2 ℓ,ℓ ′ ∈N ℓ =ℓ ′ Cov(1 A • W ℓ , 1 A • W ℓ ′ ).
On commence par le premier terme (les variances). On a

Var(1 A • W ℓ ) = E[(1 A • W ℓ ) 2 ] -E[1 A • W ℓ ] 2 = E[1 A • W ℓ ] -E[1 A • W ℓ ] 2 = λ s (A) -λ s (A) 2 = λ s (A) 1 -λ s (A) ≤ 1 4 , en remarquant que 0 ≤ λ s (A) ≤ 1 et que la fonction x ∈ [0, 1] → x(1 -x) atteint son maximum de 1/4 en 1/2.
On passe alors à la majoration du deuxième terme (les covariances). Soient ℓ, ℓ ′ ∈ L tels que ℓ = ℓ ′ .

Cov(1 A • W ℓ , 1 A • W ℓ ′ ) = E[(1 A • W ℓ )(1 A • W ℓ ′ )] -E[1 A • W ℓ ]E[1 A • W ℓ ′ ]. On a E[1 A • W ℓ ] = E[1 A • W ℓ ′ ] = λ s (A). D'un autre côté E[(1 A • W ℓ )(1 A • W ℓ ′ )] = Ω 1 A σ 1 (ω)(ℓ) + ξ ℓ,1 (ω) -1 N , . . . , σ s (ω)(ℓ) + ξ ℓ,s (ω) -1 N 1 A σ 1 (ω)(ℓ ′ ) + ξ ℓ ′ ,1 (ω) -1 N , . . . , σ s (ω)(ℓ ′ ) + ξ ℓ ′ ,s (ω) -1 N dP(ω) = 1 (N !) s σ 1 ∈S N • • • σs∈S N I 2s 1 A σ 1 (ℓ) + u 1 -1 N , . . . , σ s (ℓ) + u s -1 N 1 A σ 1 (ℓ ′ ) + u ′ 1 -1 N , . . . , σ s (ℓ ′ ) + u ′ s -1 N dudu ′ Donc E[(1 A • W ℓ )(1 A • W ℓ ′ )] = 1 (N !) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s σ 1 ∈S N σ 1 (ℓ)=m 1 , σ 1 (ℓ ′ )=m ′ 1 • • • σs∈S N σs(ℓ)=ms, σs(ℓ ′ )=m ′ s I s 1 A m 1 + u 1 -1 N , . . . , m s + u s -1 N du I s 1 A m ′ 1 + u ′ 1 -1 N , . . . , m ′ s + u ′ s -1 N du ′ = 1 (N (N -1)) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s I s 1 A m 1 + u 1 -1 N , . . . , m s + u s -1 N du I s 1 A m ′ 1 + u ′ 1 -1 N , . . . , m ′ s + u ′ s -1 N , du ′ car, si m i , m ′ i ∈ N avec m i = m ′ i , #{σ i ∈ S N : σ i (ℓ) = m i et σ i (ℓ ′ ) = m ′ i } = (N -2)
! En faisant les changements de variables :

v i := m i + u i -1 N , v ′ i := m i + u ′ i -1 N on obtient E[(1 A • W ℓ )(1 A • W ℓ ′ )] = N 2s (N (N -1)) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s Im 1 A (v)dv I m ′ 1 A (v ′ )dv ′ = N s (N -1) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ). D'où Cov[(1 A • W ℓ )(1 A • W ℓ ′ )] = N s (N -1) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) -λ s (A) 2 = N s (N -1) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) - m 1 ,m ′ 1 ∈N • • • ms,m ′ s ∈N λ s (A ∩ I m )λ s (A ∩ I m ′ ).
On en déduit

Cov[(1 A • W ℓ )(1 A • W ℓ ′ )] ≤ N s (N -1) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) - m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) ≤ N s (N -1) s -1 N s (N -1) s 1 N 2s ≤ 1 - N -1 N s ≤ s N .

D'où

Var(Z) ≤ 1 N 2 ℓ∈L 1 4 + 1 N 2 ℓ,ℓ ′ ∈L ℓ =ℓ ′ s N ≤ 1 4N + s(N -1) N 2 ,
qui est le résultat. Si A est un rectangle à côtés parallèles aux axes de coordonnées, on a un résultat de réduction de variance par rapport à la méthode de Monte Carlo. On a vu que si X est l'estimateur de Monte Carlo usuel de J , on a

Var(X) = 1 N λ s (A) 1 -λ s (A) .
Ici encore l'amélioration par rapport à la proposition 2.11 est minime. On supposera dans la suite N ≥ 3.

Proposition 2.13 Si A est un rectangle de I s , à côtés parallèles aux axes de coordonnées, la variance de l'estimateur LHS

Z = 1 N ℓ∈N 1 A • W ℓ vérifie Var(Z) ≤ 1 N λ s (A) 1 -λ s (A) .
Preuve. D'après la démonstration de la proposition précédente, on a

Var(Z) = 1 N 2 ℓ∈N Var(1 A • W ℓ ) + 1 N 2 ℓ,ℓ ′ ∈N ℓ =ℓ ′ Cov(1 A • W ℓ , 1 A • W ℓ ′ ) = 1 N λ s (A) 1 -λ s (A) + 1 N 2 ℓ,ℓ ′ ∈N ℓ =ℓ ′ Cov(1 A • W ℓ , 1 A • W ℓ ′ ).
Il suffit donc de démontrer

∀ℓ, ℓ ′ ∈ N ℓ = ℓ ′ Cov(1 A • W ℓ , 1 A • W ℓ ′ ) ≤ 0.

On a vu

Cov(1

A • W ℓ , 1 A • W ℓ ′ ) = N s (N -1) s m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) -λ s (A) 2 .
Soit

A := s i=1 A i , où A i est un intervalle d'extrémités m i -x -1 i -1 N et m i + n i + x +1 i -1 N , avec m i , n i ∈ N x -1 i , x +1 i ∈ [0, 1). On a alors λ s (A) = 1 N s s i=1 n i + x -1 i + x +1 i . On note [1, s] := {1, 2, . . . , s} et, pour J ⊂ [1, s], on note J c := [1, s] \ J. Alors m 1 ,m ′ 1 ∈N m 1 =m ′ 1 • • • ms,m ′ s ∈N ms =m ′ s λ s (A ∩ I m )λ s (A ∩ I m ′ ) = 1 N 2s J⊂[1,s] ǫ k =±1,k∈J c j∈J n j j∈J n j + x -1 j + x +1 j -1 k∈J c x ǫ k k k∈J c n k + x -ǫ k k .
On montre par récurrence sur s

J⊂[1,s] ǫ k =±1,k∈J c j∈J n j j∈J n j + x -1 j + x +1 j -1 k∈J c x ǫ k k k∈J c n k + x -ǫ k k = s i=1 n i (n i + x -1 i + x +1 i -1) + (n i + x -1 i )x +1 i + (n i + x +1 i )x -1 i .
Il suffit alors de démontrer

N s (N -1) s s i=1 n i (n i + x -1 i + x +1 i -1) + (n i + x -1 i )x +1 i + (n i + x +1 i )x -1 i ≤ s i=1 n i + x -1 i + x +1 i 2 .
Soit, pour tous les n ∈ N et x -1 , x +1 ∈ [0, 1) :

N N -1 n(n + x -1 + x +1 -1) + (n + x -1 )x +1 + (n + x +1 )x -1 ≤ n + x -1 + x +1 2 .
On réduit le polynôme quadratique :

P N,n (x -1 , x +1 ) := (N -1) n + x -1 + x +1 2 -N n(n + x -1 + x +1 -1) + (n + x -1 )x +1 + (n + x +1 )x -1 = N n + (x -1 ) 2 + (x +1 ) 2 -n + x -1 + x +1 2 = 1 2 (N -2) x -1 + x +1 - 2n N -2 2 + N 2 x -1 -x +1 2 + nN (N -2 -n) N -2 .
• Si x -1 = x +1 = 0, on a, pour 1 ≤ n ≤ N :

P N,n (0, 0) = n(N -n) ≥ 0.
• Si x -1 = 0, x +1 ∈ (0, 1), alors 1 ≤ n ≤ N -1 et l'on a :

P N,n (0, x +1 ) = (N -1) x +1 - n N -1 2 + nN (N -1 -n) N -1 ≥ 0.
• Si x -1 ∈ (0, 1), x +1 = 0, alors 1 ≤ n ≤ N -1 et l'on a de manière analogue :

P N,n (x -1 , 0) = (N -1) x -1 - n N -1 2 + nN (N -1 -n) N -1 ≥ 0.
• Si x -1 , x +1 ∈ (0, 1), alors 1 ≤ n ≤ N -2 et l'on a d'après ce qui précède : Ce qui donne le résultat. La méthode des tableaux orthogonaux (orthogonal arrays) permet de généraliser l'échantillonnage LHS. Au lieu de rechercher une stratification sur les projections en dimension 1, on vise une stratification sur les projections en dimension ≤ r (avec r ≤ s). Cette approche a été initiée par Tang [START_REF] Tang | Orthogonal array-based Latin hypercubes[END_REF] et Owen [START_REF] Owen | Orthogonal arrays for computer experiments, integration and visualization[END_REF][START_REF] Owen | Lattice sampling revisited : Monte Carlo variance of means over randomized orthogonal arrays[END_REF]. L'analyse d'Owen de [START_REF] Owen | Monte Carlo variance of scrambled net quadrature[END_REF] indique une variance du même ordre que celle de la méthode LHS (et aussi du schéma de Monte Carlo usuel).

P N,n (x -1 , 0) = 1 2 (N -2) x -1 + x +1 - 2n N -2 2 + N 2 x -1 -x +1 2 + nN (N -2 -n) N -2 ≥ 0.

Essais numériques

Dans ce paragraphe, nous voulons vérifier les analyses précédentes. Comme pour la méthode MC stratifiée simple, nous nous limitons au calcul d'intégrales de fonctions indicatrices en dimensions 2, 3 et 4 et nous considérons le même sous-ensemble de la boule unité s-dimensionnelle :

Q := {x ∈ I s : x 2 < 1}.
On évalue les variances des estimateurs de J = λ s (Q) calculés par une méthode de Monte Carlo usuelle ou la méthode LHS. Pour un nombre N de points aléatoires utilisés, on calcule la variance empirique d'un échantillon de M observations de J obtenues en répétant le calcul de façon indépendante.

On a représenté les valeurs des variances des estimateurs MC et LHS en fonction du nombre de points N de la simulation sur la Figure 2.7. On a tracé les courbes des variances pour M = 100 et M = 1 000 répétitions : on voit que la différence des résultats est assez faible.

En supposant que la variance est de l'ordre de N -γ , on détermine les valeurs de γ en effectuant un ajustement linéaire pour M = 1 000. Ces valeurs sont regroupées dans la Table 2.3.

L'estimateur de la méthode LHS est meilleur que celui de la méthode MC simple, puisque la variance est toujours inférieure. Par contre, l'ordre γ est le même.

Supposant que la variance s'exprime par Var = cN -γ (avec c constante) on voit que la constante de la méthode LHS est inférieure à celle de la méthode MC, dans cet exemple.

Conclusion

Dans ce chapitre nous avons d'abord présenté et analysé une méthode de Monte Carlo stratifiée consistant à utiliser une partition de l'hypercube unité en N sous-domaines

Introduction

Dans ce chapitre nous étudions une méthode d'échantillonnage stratifiée appelée "Sudoku", qui sera utilisée pour la résolution numérique des équations de diffusion et coagulation : voir les chapitres 4 et 5. Cette méthode est intermédiaire entre la méthode de stratification simple (MCS) et la méthode LHS analysées dans le chapitre précédent. Elle possède à la fois les propriétés de stratification multi-dimensionnelle de la méthode MCS et les propriétés de stratification uni-dimensionnelle des méthodes LHS. Nous considérons d'abord le problème de l'intégration numérique des fonctions indicatrices et nous démontrons une borne de la variance de l'estimateur. Dans la situation analysée, les variables aléatoires ne sont pas indépendantes, ce qui rend les calculs plus laborieux que dans le chapitre précédent. L'étude théorique est complétée par des essais numériques. Une partie des résultats de ce chapitre a été publiée dans [START_REF] El Haddad | Extended Latin hypercube sampling for integration and simulation[END_REF].

Stratification Sudoku

Soit s ≥ 1 un entier, f : I s → R une fonction mesurable et λ s la mesure de Lebesgue sur R s . On veut approcher

J := I s f (x)dλ s (x)
par une méthode de Monte Carlo stratifiée. On commence par construire l'ensemble de points aléatoires qui est utilisé dans l'approximation. Soit 2. ξ (1) , . . . , ξ (s) des variables aléatoires uniformément distribuées sur I N ;

3. σ (1) , . . . , σ (s) des bijections aléatoires de {1, . . . , n} s-1 dans {1, . . . , n s-1 }, de loi uniforme sur l'ensemble S de toutes ces bijections.

On suppose que toutes les variables aléatoires définies dans (2) et (3) sont indépendantes.

Pour un vecteur x = (x 1 , . . . , x s ), on note xi := (x 1 , . . . , x i-1 , x i+1 , . . . , x s ). On pose pour tout ℓ ∈ L :

W ℓ = (W ℓ,1 , . . . , W ℓ,s ) avec W ℓ,i := ℓ i -1 n + σ (i) ( li ) -1 N + ξ (i) ℓ N . ( 3
ℓ ′ i -1 n , ℓ ′ i n , 1 ≤ ℓ ′ 1 ≤ n, . . . , 1 ≤ ℓ ′ s ≤ n
il n'y a qu'un seul point.

Dans tout hyper-rectangle de la forme

[0, 1) i-1 × k -1 n s , k n s × [0, 1) s-i , 1 ≤ i ≤ s, 1 ≤ k ≤ n s
il n'y a qu'un seul point.

Preuve.

1. Soient 1 ≤ ℓ ′ 1 ≤ n, . . . , 1 ≤ ℓ ′ s ≤ n. On a les équivalences W ℓ ∈ s i=1 ℓ ′ i -1 n , ℓ ′ i n ⇔ ∀i W ℓ,i ∈ ℓ ′ i -1 n , ℓ ′ i n ⇔ ∀i ℓ ′ i -1 n ≤ ℓ i -1 n + σ (i) ( li ) -1 n s + ξ (i) ℓ n s < ℓ ′ i n ⇔ ∀i ℓ ′ i -1 ≤ ℓ i -1 + σ (i) ( li ) -1 n s-1 + ξ (i) ℓ n s-1 < ℓ ′ i ⇔ ∀i ℓ ′ i ≤ ℓ i + σ (i) ( li ) -1 n s-1 + ξ (i) ℓ n s-1 < ℓ ′ i + 1. Or 0 ≤ σ (i) ( li ) -1 n s-1 + ξ (i) ℓ n s-1 < n s-1 -1 n s-1 + 1 n s-1 = 1. Par conséquent, W ℓ ∈ s i=1 ℓ ′ i -1 n , ℓ ′ i n ⇔ ∀i ℓ i = ℓ ′ i .
Ce qui prouve que dans chaque hypercube on a un point unique.

2. Soit 1 ≤ i ≤ s, 1 ≤ k ≤ n s . On a les équivalences

W ℓ ∈ I i-1 × k -1 n s , k n s × I s-i ⇔ W ℓ,i ∈ k -1 n s , k n s ⇔ k -1 n s ≤ ℓ i -1 n + σ (i) ( li ) -1 n s + ξ (i) ℓ n s < k n s ⇔ k -1 ≤ (ℓ i -1)n s-1 + σ (i) ( li ) -1 + ξ (i) ℓ < k ⇔ k ≤ (ℓ i -1)n s-1 + σ (i) ( li ) + ξ (i) ℓ < k + 1. Puisque 0 ≤ ξ (i) ℓ < 1, on a W ℓ ∈ I i-1 × k -1 n s , k n s × I s-i ⇔ k = (ℓ i -1)n s-1 + σ (i) ( li ) ⇔ k -1 = (ℓ i -1)n s-1 + σ (i) ( li ) -1. Comme 1 ≤ σ (i) ( li ) ≤ n s-1 donc 0 ≤ σ (i) ( li ) -1 < n s-1 et ℓ i -1 est le quotient de la division euclidienne de k -1 par n s-1 . Comme 1 ≤ k ≤ n s donc 0 ≤ k -1 < n s on a 0 ≤ ℓ i -1 < n donc 1 ≤ ℓ i ≤ n.
Par ailleurs σ (i) ( li ) -1 est le reste de la division de k -1 par n s-1 . Cela détermine li . Ainsi il y a un point unique dans tout hyper-rectangle.

D'où le résultat. Cette construction de points est une combinaison des méthodes de Monte-Carlo stratifiée simple (MCS) et d'échantillonnage par hypercube latin (LHS) définies au chapitre précédent. Un exemple d'un ensemble stratifié Sudoku en dimension s = 2, formé de 16 points est représenté sur la Figure 3.1. Si n = 4, il y a un point unique dans chacun des carrés de la forme :

ℓ ′ 1 -1 n , ℓ ′ 1 n × ℓ ′ 2 -1 n , ℓ ′ 2 n , 1 ≤ ℓ ′ 1 ≤ n, 1 ≤ ℓ ′ 2 ≤ n
(en haut), dans chacun des rectangles de la forme :

[0, 1) × k -1 n 2 , k n 2 , 1 ≤ k ≤ n 2
(au milieu) et dans chacun des rectangles de la forme :

k -1 n 2 , k n 2 × [0, 1), 1 ≤ k ≤ n 2
(en bas). Si n 2 = b µ (avec b, µ des entiers), cet ensemble a certaines propriétés des réseaux (0, µ, 2) en base b [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]. On note dans la suite :

• pour ℓ ∈ L et m ∈ M, I ℓ,m := s i=1 ℓ i -1 n + m i -1 n s , ℓ i -1 n + m i n s ;
• pour ℓ ∈ L,

I ℓ := m∈M I ℓ,m = s i=1 ℓ i -1 n , ℓ i n .
L'estimateur stratifié Sudoku de J est défini par

Z := 1 N ℓ∈L f • W ℓ .
Proposition 3.2 L'estimateur stratifié Sudoku Z est un estimateur sans biais de J .

Preuve. Pour ℓ ∈ L, on a

E[f • W ℓ ] = Ω f ℓ 1 -1 n + σ (1) (ω)( l1 ) -1 + ξ (1) ℓ (ω) n s , . . . , ℓ s -1 n + σ (s) (ω)( ls ) -1 + ξ (s) ℓ (ω) n s dP(ω) = 1 (n s-1 !) s σ 1 ,...,σs∈S I s f ℓ 1 -1 n + σ 1 ( l1 ) -1 n s + ξ 1 n s , . . . , ℓ s -1 n + σ s ( ls ) -1 n s + ξ s n s dξ 1 . . . dξ s = 1 (n s-1 !) s m∈M σ 1 ,...,σs∈S ∀i σ i ( li )=m i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s .
Or, pour tout 1 ≤ i ≤ s, on a

#{σ i ∈ S : σ i ( li ) = m i } = (n s-1 -1)! D'où E[f • W ℓ ] = 1 n s(s-1) m∈M I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s = (n s ) s (n s-1 ) s m∈M I ℓm f (u 1 , . . . , u s )du 1 . . . du s ,
en faisant les changements de variables

u 1 := ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , u s := ℓ s -1 n + m s -1 n s + ξ s n s . Donc E[f • W ℓ ] = n s m∈M I ℓm f (u 1 , . . . , u s )du 1 . . . du s = N I ℓ f (u 1 , . . . , u s )du 1 . . . du s .
Par suite

E[Z] = 1 N ℓ∈L E[f • W ℓ ] = ℓ∈L I ℓ f (u 1 , . . . , u s )du 1 . . . du s = J .
C'est le résultat.

Calcul de la variance de l'estimateur stratifié Sudoku

On étudie la variance de l'estimateur Z. À la fin de l'étude, on se limite à f := 1 A , où A est un sous-ensemble mesurable de I s . Les variables aléatoires W ℓ , ℓ ∈ L ne sont pas indépendantes, il faut donc calculer des termes de covariance. On a

Var(Z) = Var 1 N ℓ∈L f • W ℓ = 1 N 2 ℓ∈L Var(f • W ℓ ) + 1 N 2 ℓ,ℓ ′ ∈L ℓ =ℓ ′ Cov(f • W ℓ , f • W ℓ ′ ).
Pour la première somme, on a, d'après la démonstration de la Proposition 3.2 :

Var(f • W ℓ ) = E[(f • W ℓ ) 2 ] -(E[f • W ℓ ]) 2 = N I ℓ (f (u)) 2 du -N 2 I ℓ f (u)du 2 . Donc 1 N 2 ℓ∈L Var(f • W ℓ ) = ℓ∈L 1 N I ℓ (f (u)) 2 du - I ℓ f (u)du 2 .
On analyse les termes de la seconde somme. Soient ℓ, ℓ ′ ∈ L avec ℓ = ℓ ′ . On a

Cov(f • W ℓ , f • W ℓ ′ ) = E[(f • W ℓ )(f • W ℓ ′ )] -E[f • W ℓ ]E[f • W ℓ ′ ].
Ensuite,

E[(f • W ℓ )(f • W ℓ ′ )] = Ω f ℓ 1 -1 n + σ (1) (ω)( l1 ) -1 + ξ (1) ℓ (ω) n s , . . . , ℓ s -1 n + σ (s) (ω)( ls ) -1 + ξ (s) ℓ (ω) n s • f ℓ ′ 1 -1 n + σ (1) (ω)( l′ 1 ) -1 + ξ (1) ℓ ′ (ω) n s , . . . , ℓ ′ s -1 n + σ (s) (ω)( l′ s ) -1 + ξ (s) ℓ ′ (ω) n s dP(ω) = 1 (n s-1 !) s σ 1 ,...,σs∈S I s f ℓ 1 -1 n + σ 1 ( l1 ) -1 n s + ξ 1 n s , . . . , ℓ s -1 n + σ s ( ls ) -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + σ 1 ( l′ 1 ) -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + σ s ( l′ s ) -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s = 1 (n s-1 !) s m,m ′ ∈M m =m ′ σ 1 ,...,σs∈S σ i ( li )=m i , ∀i σ i ( l′ i )=m ′ i , ∀i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + m ′ 1 -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + m ′ s -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s .
On distingue deux cas.

1. S'il existe un unique indice j ∈ {1, . . . , s} tel que ℓ j = ℓ ′ j , alors lj = l′ j . Dans ce cas, on a :

E[(f • W ℓ )(f • W ℓ ′ )] = 1 (n s-1 !) s m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i σ 1 ,...,σs∈S ∀i σ i ( li )=m i ∀i σ i ( l′ i )=m ′ i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + m ′ 1 -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + m ′ s -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s .
On a

#{σ j ∈ S : σ j ( lj ) = m j } = (n s-1 -1)! #{σ i ∈ S : σ i ( li ) = m i , σ i ( l′ i ) = m ′ i } = (n s-1 -2)! Par conséquent E[(f • W ℓ )(f • W ℓ ′ )] = 1 n s(s-1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + m ′ 1 -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + m ′ s -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s .
Donc, en faisant les changements de variables

u 1 := ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , u s := ℓ s -1 n + m s -1 n s + ξ s n s , on obtient E[(f • W ℓ )(f • W ℓ ′ )] = n s(s+1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ . 2. Si ∀i li = l′ i , alors E[(f • W ℓ )(f • W ℓ ′ )] = 1 (n s-1 !) s m,m ′ ∈M ∀i m i =m ′ i σ 1 ,...,σs∈S ∀i σ i ( li )=m i ∀i σ i ( l′ i )=m ′ i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + m ′ 1 -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + m ′ s -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s .
Et comme plus haut

E[(f • W ℓ )(f • W ℓ ′ )] = 1 n s(s-1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i I s f ℓ 1 -1 n + m 1 -1 n s + ξ 1 n s , . . . , ℓ s -1 n + m s -1 n s + ξ s n s dξ 1 . . . dξ s I s f ℓ ′ 1 -1 n + m ′ 1 -1 n s + ξ ′ 1 n s , . . . , ℓ ′ s -1 n + m ′ s -1 n s + ξ ′ s n s dξ ′ 1 . . . dξ ′ s .
Donc, en faisant les mêmes changements de variables que précédemment,

E[(f • W ℓ )(f • W ℓ ′ )] = n s(s+1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ .
Par ailleurs

E[f • W ℓ ]E[f • W ℓ ′ ] = N 2 I ℓ f (u)du I ℓ ′ f (u ′ )du ′ . Donc 1. Si ∃!j ∈ {1, . . . , s} ℓ j = ℓ ′ j alors Cov(f • W ℓ , f • W ℓ ′ ) = n s(s+1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ -n 2s I ℓ f (u)du I ℓ ′ f (u ′ )du ′ . 2. Si ∀i ∈ {1, . . . , s} li = l′ i alors Cov(f • W ℓ , f • W ℓ ′ ) = n s(s+1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ -n 2s I ℓ f (u)du I ℓ ′ f (u ′ )du ′ . Par conséquent 1 N 2 ℓ,ℓ ′ ∈L ℓ =ℓ ′ Cov(f • W ℓ , f • W ℓ ′ ) = 1 n 2s s j=1 ℓ,ℓ ′ ∈L ℓ j =ℓ ′ j , lj = l′ j Cov(f • W ℓ , f • W ℓ ′ ) + 1 n 2s ℓ,ℓ ′ ∈L ∀i li = l′ i Cov(f • W ℓ , f • W ℓ ′ ) = s j=1 ℓ,ℓ ′ ∈L ℓ j =ℓ ′ j , lj = l′ j n s(s-1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ - I ℓ f (u)du I ℓ ′ f (u ′ )du ′ + ℓ,ℓ ′ ∈L ∀i li = l′ i n s(s-1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ - I ℓ f (u)du I ℓ ′ f (u ′ )du ′ . On écrit 1 N 2 ℓ,ℓ ′ ∈L ℓ =ℓ ′ Cov(f • W ℓ , f • W ℓ ′ ) = s j=1 ℓ,ℓ ′ ∈L ℓ j =ℓ ′ j , lj = l′ j V j (ℓ, ℓ ′ ) + ℓ,ℓ ′ ∈L ∀i li = l′ i V s+1 (ℓ, ℓ ′ ), où V j (ℓ, ℓ ′ ) := n s(s-1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ - I ℓ f (u)du I ℓ ′ f (u ′ )du ′ , 1 ≤ j ≤ s, V s+1 (ℓ, ℓ ′ ) := n s(s-1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i I ℓm f (u)du I ℓ ′ m ′ f (u ′ )du ′ - I ℓ f (u)du I ℓ ′ f (u ′ )du ′ . Finalement Var(Z) = 1 N 2 ℓ∈L Var(f • W ℓ ) + 1 N 2 ℓ,ℓ ′ ∈L ℓ =ℓ ′ Cov(f • W ℓ , f • W ℓ ′ ) = ℓ∈L V 0 (ℓ) + s j=1 ℓ,ℓ ′ ∈L ℓ j =ℓ ′ j , lj = l′ j V j (ℓ, ℓ ′ ) + ℓ,ℓ ′ ∈L ∀i li = l′ i V s+1 (ℓ, ℓ ′ ) où V 0 (ℓ) := 1 n s I ℓ (f (u)) 2 du - I ℓ f (u)du 2
On se limite désormais à f := 1 A , où A est un sous-ensemble mesurable de I s . Alors

V 0 (ℓ) = 1 n s I ℓ 1 A (u)du - I ℓ 1 A (u)du 2 = 1 n s λ s (A ∩ I ℓ ) -λ s (A ∩ I ℓ ) 2 . Pour 1 ≤ j ≤ s, V j (ℓ, ℓ ′ ) = n s(s-1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i λ s (A ∩ I ℓm )λ s (A ∩ I ℓ ′ m ′ ) -λ s (A ∩ I ℓ )λ s (A ∩ I ℓ ′ ). Et V s+1 (ℓ, ℓ ′ ) = n s(s-1) (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i λ s (A ∩ I ℓm )λ s (A ∩ I ℓ ′ m ′ ) -λ s (A ∩ I ℓ )λ s (A ∩ I ℓ ′ ). • Soit ℓ ∈ L. Si I ℓ ⊂ A ou si I ℓ ∩ A = ∅, V 0 (ℓ) = 0. Si I ℓ ⊂ A et I ℓ ∩ A = ∅, comme 0 ≤ n s λ s (A ∩ I ℓ ) ≤ 1 et V 0 (ℓ) = 1 n 2s n s λ s (A ∩ I ℓ ) 1 -n s λ s (A ∩ I ℓ ) , on a 0 ≤ V 0 (ℓ) ≤ 1 4n 2s . En notant V 0 := ℓ∈L V 0 (ℓ),
on obtient la majoration

|V 0 | = ℓ∈L I ℓ ⊂A, I ℓ ∩A =∅ V 0 (ℓ) ≤ 1 4n 2s #{ℓ ∈ L : I ℓ ⊂ A, I ℓ ∩ A = ∅}. • Soient 1 ≤ j ≤ s et ℓ, ℓ ′ ∈ L avec ℓ j = ℓ ′ j et lj = l′ j . Si I ℓ ⊂ A, V j (ℓ, ℓ ′ ) = 1 n s (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i λ s (A ∩ I ℓ ′ m ′ ) - 1 n s λ s (A ∩ I ℓ ′ ) = 1 n s m ′ ∈M λ s (A ∩ I ℓ ′ m ′ ) - 1 n s λ s (A ∩ I ℓ ′ ) = 0. De même, si I ′ ℓ ⊂ A, V j (ℓ, ℓ ′ ) = 0.
Et on a aussi, si

I ℓ ∩ A = ∅ ou I ℓ ′ ∩ A = ∅ V j (ℓ, ℓ ′ ) = 0. Si I ℓ ⊂ A, I ℓ ∩ A = ∅, I ℓ ′ ⊂ A et I ℓ ′ ∩ A = ∅, alors 0 ≤ n s(s-1) (n s-1 -1) s-1 m,m ′ ∈M m j =m ′ j ∀i =j m i =m ′ i λ s (A ∩ I ℓm )λ s (A ∩ I ℓ ′ m ′ ) ≤ n s(s-1) (n s-1 -1) s-1 n s-1 n s-1 (n s-1 -1) s-1 1 n s 2 1 n s 2 = 1 n 2s et 0 ≤ λ s (A ∩ I ℓ )λ s (A ∩ I ℓ ′ ) ≤ 1 n 2s . Par suite |V j (ℓ, ℓ ′ )| ≤ 1 n 2s . En notant V j := ℓ,ℓ ′ ∈L ℓ j =ℓ ′ j , lj = l′ j V j (ℓ, ℓ ′ ),
on obtient la majoration

|V j | ≤ 1 n 2s #{(ℓ, ℓ ′ ) ∈ L 2 : ℓ j = ℓ ′ j , lj = l′ j , I ℓ ⊂ A, I ℓ ∩A = ∅, I ′ ℓ ⊂ A, I ′ ℓ ∩A = ∅}. • Soient ℓ, ℓ ′ ∈ L tels que ∀i ∈ {1, . . . , s} li = l′ i . Si I ℓ ⊂ A, V s+1 (ℓ, ℓ ′ ) = 1 n s (n s-1 -1) s m,m ′ ∈M ∀i m i =m ′ i λ s (A ∩ I ℓ ′ m ′ ) - 1 n s λ s (A ∩ I ℓ ′ ) = 1 n s m ′ ∈M λ s (A ∩ I ℓ ′ m ′ ) - 1 n s λ s (A ∩ I ℓ ′ ) = 0.
De même, si

I ′ ℓ ⊂ A, V j (ℓ, ℓ ′ ) = 0.
Et on a aussi, si

I ℓ ∩ A = ∅ ou I ℓ ′ ∩ A = ∅ V j (ℓ, ℓ ′ ) = 0. Si I ℓ ⊂ A, I ℓ ∩ A = ∅, I ℓ ′ ⊂ A et I ℓ ′ ∩ A = ∅, on utilise la décomposition λ s (A ∩ I ℓ )λ s (A ∩ I ℓ ′ ) = m,m ′ ∈M λ s (A ∩ I ℓm )λ s (A ∩ I ℓ ′ m ′ ). Alors V s+1 (ℓ, ℓ ′ ) = m∈M λ s (A ∩ I ℓm )V (ℓ, ℓ ′ , m), où V (ℓ, ℓ ′ , m) := n s(s-1) (n s-1 -1) s m ′ ∈M ∀i m ′ i =m i λ s (A ∩ I ℓ ′ m ′ ) - m ′ ∈M λ s (A ∩ I ℓ ′ m ′ ) = n s(s-1) (n s-1 -1) s -1 m ′ ∈M ∀i m ′ i =m i λ s (A ∩ I ℓ ′ m ′ ) - m ′ ∈M\{m ′ ∈M : ∀i m ′ i =m i } λ s (A ∩ I ℓ ′ m ′ ). On a ∀x ∈ [0, 1] 1 -(1 -x) s ≤ sx donc, en prenant x := 1 n s-1 0 ≤ n s(s-1) (n s-1 -1) s -1 ≤ sn (s-1) 2 (n s-1 -1) s et # M \ {m ′ ∈ M : ∀i m ′ i = m i } = n s(s-1) -(n s-1 -1) s ≤ sn (s-1) 2 . On en déduit 0 ≤ n s(s-1) (n s-1 -1) s -1 m ′ ∈M ∀i m ′ i =m i λ s (A ∩ I ℓ ′ m ′ ) ≤ sn (s-1) 2 (n s-1 -1) s (n s-1 -1) s n s 2 = s n 2s-1 et 0 ≤ m ′ ∈M\{m ′ ∈M : ∀i m ′ i =m i } λ s (A ∩ I ℓ ′ m ′ ) ≤ sn (s-1) 2 n s 2 = s n 2s-1 . Par conséquent |V (ℓ, ℓ ′ , m)| ≤ s n 2s-1 et |V s+1 (ℓ, ℓ ′ )| ≤ n s(s-1) n s 2 s n 2s-1 = s n 3s-1 . En notant V s+1 := ℓ,ℓ ′ ∈L ∀i li = l′ i V s+1 (ℓ, ℓ ′ ), on obtient la majoration |V s+1 | ≤ s n 3s-1 #{(ℓ, ℓ ′ ) ∈ L 2 : ∀i li = l′ i , I ℓ ⊂ A, I ℓ ∩ A = ∅, I ′ ℓ ⊂ A, I ′ ℓ ∩ A = ∅}.
En combinant les majorations précédentes on obtient la majoration de la variance de Z suivante

Var(Z) ≤ |V 0 | + s j=1 |V j | + |V s+1 | ≤ 1 4n 2s #{ℓ ∈ L : I ℓ ⊂ A, I ℓ ∩ A = ∅} + 1 n 2s s j=1 #{(ℓ, ℓ ′ ) ∈ L 2 : ℓ j = ℓ ′ j , lj = l′ j , I ℓ ⊂ A, I ℓ ∩ A = ∅, I ′ ℓ ⊂ A, I ′ ℓ ∩ A = ∅} + s n 3s-1 #{(ℓ, ℓ ′ ) ∈ L 2 : ∀i li = l′ i , I ℓ ⊂ A, I ℓ ∩ A = ∅, I ′ ℓ ⊂ A, I ′ ℓ ∩ A = ∅}.
On obtient un ordre de la variance de Z sous une hypothèse de régularité sur la frontère du domaine A. Proposition 3.3 Soit A un sous-ensemble mesurable de I s . On suppose que, pour tout 1 ≤ i ≤ s, il existe une fonction f i : I s-1 → I lipschitzienne telle que

A = {(u 1 , . . . , u s ) ∈ I s : u i < f i (û i )},
où ûi := (u 1 , . . . , u i-1 , u i+1 , . . . u s ). Soit {W ℓ : ℓ ∈ L} l'ensemble de points défini par (3.2.1) et

Z := 1 N ℓ∈L 1 A • W ℓ l'estimateur stratifié Sudoku de λ s (A). Alors Var(Z) ≤ C + 2 4 + 2s(C + 2) 2 1 N 1+1/s , où C est une constante de Lipschitz (pour la norme du maximum • ∞ ) de chaque fonc- tion f i . Preuve. Pour ℓ ∈ L et 1 ≤ i ≤ s, on note Îℓ,i = 1≤j≤s j =i ℓ j -1 n , ℓ j n .
D'après la géométrie de l'ensemble A, on a les propriétés suivantes :

• si I ℓ ⊂ A, alors il existe ûℓ,i ∈ Îℓ,i tel que nf i (û ℓ,i ) < ℓ i ; • si I ℓ ∩ A = ∅, alors il existe vℓ,i ∈ Îℓ,i tel que ℓ i < nf i (v ℓ,i ) + 1. Par conséquent, on a #{ℓ ∈ L : I ℓ ⊂ A, I ℓ ∩ A = ∅} ≤ n s-1 (C + 2), pour tout j : 1 ≤ j ≤ s #{(ℓ, ℓ ′ ) ∈ L 2 : ℓ j = ℓ ′ j , lj = l′ j , I ℓ ⊂ A, I ℓ ∩ A = ∅, I ′ ℓ ⊂ A, I ′ ℓ ∩ A = ∅} ≤ n s-1 (C + 2) 2 et #{(ℓ, ℓ ′ ) ∈ L 2 : ∀i li = l′ i , I ℓ ⊂ A, I ℓ ∩ A = ∅, I ℓ ′ ⊂ A, I ℓ ′ ∩ A = ∅} ≤ n 2(s-1) (C + 2) 2 .
Le résultat suit en utilisant ces inégalités dans la majoration de la variance de Z obtenue plus haut.

Ce résultat indique un gain de précision par rapport à une méthode MC usuelle d'un facteur de N -1/s . Cet avantage s'amoindrit quand la dimension augmente et semble limiter l'utilisation de cette stratification à des dimensions peu élevées. Faisant la même remarque que pour l'approche de MC stratifié simple, signalons que la technique de stratification Sudoku peut être intéressante pour des simulations nécessitant l'évaluation d'intégrales en dimension s comprise entre 2 et 6 : voir par exemple [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | Low discrepancy sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo scheme using nets for a linear Boltzmann equation[END_REF][START_REF] Lécot | Comparison of quasi-Monte Carlo-based methods for the simulation of Markov chains[END_REF][START_REF] Lécot | Quasi-Monte Carlo methods for estimating transient measures of discrete time Markov chains[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of discrete-time Markov chain on multidimensional state space[END_REF][START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] El Haddad | Quasi-Monte Carlo methods for Markov chains with continuous multi-dimensional state space[END_REF], pour une approche QMC et [START_REF] L'ecuyer | Randomized quasi-Monte Carlo simulation of Markov chains with an ordered state space[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF][START_REF] L'ecuyer | On array-RQMC for Markov chains : mapping alternatives and convergence rates[END_REF] pour une approche QMC hasardisée. C'est aussi le cas des simulations détaillées aux chapitres 4 et 5.

L'ordre théorique est confirmé numériquement dans la section suivante, où l'on remarque que la stratification Sudoku reste efficace même si l'hypothèse sur la frontière de A n'est pas satisfaite.

Q := {x ∈ I s : x 2 ≤ 1}.
On veut calculer λ s (Q) en dimensions 2, 3 et 4 par la méthode stratifiée Sudoku et comparer les résultats obtenus à ceux des méthodes MC, MCS et LHS. Dans un premier temps, le nombre de points N dans les simulations varie comme suit :

• pour s = 2, N = 10 2 , 20 2 , . . . , 400 2 = 160 000 points ;

• pour s = 3, N = 10 3 , 20 3 , . . . , 200 3 = 8 000 000 points ;

• pour s = 4, N = 10 4 , 12 4 , . . . , 40 4 = 2 560 000 points.

Dans chaque cas, on calcule la variance des estimateurs en répétant M fois le calcul de λ s (Q) de façons indépendantes. Pour M = 100, 200, . . . , 1 000 les résultats sont similaires. On a choisi de représenter les courbes de variance pour M = 100 et M = 1 000 dans la Figure 3.2. Sur les résultats obtenus pour M = 1 000, une régression linéaire permet d'évaluer l'ordre de la variance. Les résultats sont indiqués dans la Table 3.1.

L'ordre obtenu dans les essais numériques est voisin de celui établi dans la proposition 3.3. Cela montre que l'hypothèse de régularité imposée sur la frontière du domaine d'intégration est trop restrictive, puisque dans le cas de Q, la fonction

f i : ûi → 1 -u 2 1 -• • • -u 2 i-1 -u 2 i+1 -• • • -u 2 s
n'est pas lipschitzienne sur I s-1 . Nous comparons ensuite les performances, variance et temps de calcul de la méthode SS à celles des méthodes MC, MCS et LHS du chapitre précédent. Le nombre de points N dans les simulations varie comme suit :

• pour s = 2, N = 20 2 , 40 2 , . . . , 400 2 = 160 000 points ;

• pour s = 3, N = 10 3 , 15 3 , . . . , 100 3 = 1 000 000 points ; L'efficacité est définie comme l'inverse du produit de la variance par le temps de calcul [START_REF] Hammersley | Monte Carlo Methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF]. Nous avons comparé les efficacités des méthodes MC, LHS, MCS et SS : voir la Figure 3.6. Celle de la méthode MCS est la meilleure, quand la dimension augmente. En supposant que la variance est O(N -β ), une régression linéaire permet d'évaluer l'ordre de convergence β des variances, en fonction du nombre N de points utilisés. Les résultats sont indiqués dans la Table 3 Les ordres de convergence des variances (définis comme précédemment) obtenus par régression linéaire sont indiqués dans la Table 3.3.

Les résultats obtenus par les méthodes SS et MCO sont très proches. L'avantage de la technique SS réside dans son implémentation simple.

Supposant que la variance s'exprime par Var = cN -β , avec c constante, on voit que les constantes de majoration des différentes méthodes diminuent en passant de la méthode MC aux méthodes LHS, MCS, MCO et SS respectivement. On remarque que les constantes des méthodes MCO et SS sont comparables.

Méthode de Monte Carlo stratifiée Sudoku pour l'intégration numérique

Chapitre 4

Méthode de Monte Carlo stratifiée pour la diffusion

Introduction

L'équation de la chaleur est une équation aux dérivées partielles parabolique, introduite par Fourier en 1811, pour décrire certains phénomènes physiques comme la diffusion. Celle-ci traduit le déplacement de matière dans un système physique, des zones de forte concentration vers celles de faible concentration.

L'équation de diffusion en milieu infini uni-dimensionnel s'écrit L'équation peut se généraliser en dimension quelconque, et la diffusion peut être combinée avec d'autres phénomènes, comme la convection, ou la réaction. Les méthodes numériques classiques permettent d'approcher la solution. Les plus utilisées sont celles des différences finies, des éléments finis ou des volumes finis. Les méthodes de Monte Carlo se sont avérées efficaces pour résoudre ce problème, par exemple en dimension élevée, dans un domaine à frontière irrégulière, quand la diffusion n'est pas prépondérante par rapport à d'autres phénomènes, ou quand on cherche seulement les valeurs de la solution en certains points du domaine d'étude [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF][START_REF] Hald | Convergence of random methods for a reaction-diffusion equation[END_REF][START_REF] Ghoniem | Grid-free simulation of diffusion using random walk methods[END_REF][START_REF] Sherman | A Monte Carlo method for scalar reaction diffusion equations[END_REF][START_REF] Dautray | Méthodes probabilistes pour les équations de la physique[END_REF][START_REF] Fogelson | Optimal smoothing in function-transport particle methods for diffusion problems[END_REF][START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Gobet | A spectral Monte Carlo method for the Poisson equation[END_REF][START_REF] Venkiteswaran | Quasi-Monte Carlo algorithms for diffusion equations in high dimensions[END_REF][START_REF] Venkiteswaran | A QMC approach for high dimensional Fokker-Planck equations modelling polymeric liquids[END_REF][START_REF] Tanré | Some new simulation schemes for the evaluation of Feynman-Kac representations[END_REF]. Elles sont basées sur une interprétation probabiliste de l'équation de diffusion et sur l'utilisation des marches aléatoires de particules. Des variantes quasi-Monte Carlo ont aussi été développées : voir [START_REF] Morokoff | A Quasi-Monte Carlo approach to particle simulation of the heat equation[END_REF][START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF] pour l'équation de la chaleur et [START_REF] Lécot | Grid-free simulation of convection-diffusion[END_REF][START_REF] Lécot | Quasi-random walk methods[END_REF][START_REF] Lécot | Quasirandom walk methods[END_REF][START_REF] Ogawa | A quasi-random walk method for one-dimensional reaction-diffusion equations[END_REF] pour les phénomènes de convection-diffusion ou de réaction-diffusion. L'étude menée dans [START_REF] Morokoff | A Quasi-Monte Carlo approach to particle simulation of the heat equation[END_REF] montre qu'il ne suffit pas de remplacer les points pseudo-aléatoires par des points à discrépance faible pour améliorer les résultats de la simulation, car les points déterministes ne vérifient pas les propriétés d'indépendance requises. Des techniques supplémentaires, notamment celle du tri des particules de simulation permettent de résoudre ce problème et d'accroître l'efficacité des méthodes QMC par rapport aux méthodes MC. Cette technique avait auparavant été utilisée pour la simulation QMC de l'équation de Boltzmann des gaz raréfiés [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | Low discrepancy sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF].

∂c ∂t (x, t) = ∂ ∂x D ∂c ∂x (x, t), x ∈ R, t > 0, ( 4
Dans ce chapitre, une méthode de résolution MC stratifiée est proposée dans le but d'améliorer la convergence de la simulation MC usuelle. Elle suit les étapes principales de [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF] et utilise les ensembles de points de la méthode de Monte Carlo stratifiée Sudoku (SS) développée au chapitre 3. Nous comparons les résultats de la simulation SS avec ceux d'une méthode MC classique et des méthodes MCS et LHS du chapitre 2.

Nous commençons par la description de la méthode dans le cas classique d'un coefficient de diffusion constant. Nous continuons par la mise au point d'une généralisation de la méthode dans le cas d'un milieu hétérogène (coefficient de diffusion variable). Une partie des résultats de ce chapitre a été publiée dans [START_REF] El Haddad | Extended Latin hypercube sampling for integration and simulation[END_REF].

Simulation de l'équation de diffusion

Quand le coefficient de diffusion est constant (milieu homogène) l'équation de diffusion en dimension 1 s'écrit :

∂c ∂t (x, t) = D 0 ∂ 2 c ∂x 2 (x, t), (4.2.1)
avec D = D 0 . À l'instant initial t = 0, on considère une population de particules dont la concentration approche la concentration initiale c 0 . On divise le temps en sous-intervalles de longeur ∆t. À chaque pas de temps, chaque particule effectue un déplacement aléatoire sur l'axe des x, suivant une distribution gaussienne de moyenne 0 et de variance σ 2 (que l'on indique dans la suite). Considérons une particule en x = 0 à l'instant initial. Après le premier pas de temps, on trouve une particule dans l'intervalle de longueur dx autour de x avec une probabilité P 1 (x)dx donnée par

P 1 (x)dx = e -x 2 /2σ 2 √ 2πσ dx.
Durant le pas de temps suivant, la particule se déplace en effectuant un pas de même loi normale. La probabilité de trouver une particule dans l'intervalle de longueur dx après le deuxième pas de temps est :

P 2 (x)dx = dx +∞ -∞ e -y 2 /2σ 2 e -(x-y) 2 /2σ 2 2πσ 2 dy = e -x 2 /4σ 2 √ 4πσ dx,
et ainsi de suite. Après k pas de temps on peut vérifier que

P k (x)dx = e -x 2 /2kσ 2 √ 2kπσ dx.
En choisissant σ 2 = 2D 0 ∆t et en notant t k = k∆t, on a

P k (x) = e -x 2 /4D 0 t k √ 4πD 0 t k .
On retrouve la solution fondamentale de l'opérateur

∂ ∂t -D 0 ∂ 2 ∂t 2 : E(x, t) := e -x 2 /4D 0 t √ 4πD 0 t , x ∈ R, t > 0. (4.2.2)
Elle vérifie (4.2.1) et la condition d'émission initiale en x = 0, puisque

E(x, 0) = δ(x),
où δ est la mesure de Dirac en 0. Si le départ de la particule est en x = x 0 , la solution à a l'instant t est :

E(x, t) := e -(x-x 0 ) 2 /4D 0 t √ 4πD 0 t , x ∈ R, t > 0.
D'après (4.2.2), si une particule est en x = 0 au temps t = 0, la densité de probabilité de trouver une particule à la position x au temps t est

f X (x) = e -x 2 /4D 0 t √ 4πD 0 t , (4.2.3)
qui est la densité de la loi normale centrée et de variance 2D 0 t. Le principe du schéma de Monte Carlo de simulation de l'équation de diffusion est alors le suivant. On approche la concentration initiale par une population de particules. Chaque particule de simulation se déplace pendant un pas de temps ∆t d'une distance X(∆t), où X(t) est une variable aléatoire dont la densité est donnée par (4.2.3). Donc le déplacement ∆x durant le temps ∆t est donné par :

∆x = 2D 0 ∆tZ, (4.2.4)
Z étant une variable aléatoire de loi normale centrée et réduite. Pour apprécier la qualité de cette méthode de simulation, nous avons repris des essais indiqués dans [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium : correction to the Gaussian steplength[END_REF]. Il s'agit de calculer la concentration d'ions en différents points d'un milieu autour d'une source d'ions qui émet des particules (voir [START_REF] Matveev | New and corrected simulations of synaptic facilitation[END_REF] pour le contexte biologique). Dans ce premier exemple, trois émissions de 2 10 particules chacune ont lieu à partir de l'origine x = 0 : la première émission est à l'instant initial t = 0, les deux autres surviennent respectivement aux instants t = 10 et t = 20. La simulation est menée jusqu'au temps final T = 30. Le milieu considéré est homogène avec un coefficient de diffusion constant : D 0 = 4 000 ; on simule l'équation de diffusion avec la méthode de Monte Carlo précédente. La Figure 4.1 représente l'évolution de la concentration dans les intervalles [START_REF] Haber | A modified Monte-Carlo quadrature[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF] (en haut) et [160, 180] (en bas). La courbe obtenue par la méthode MC est comparée aux valeurs exactes calculées à l'aide de la solution fondamentale. Nous voyons que si la simulation suit l'évolution de la concentration, elle introduit des oscillations parasites. La suite du chapitre indique une méthode permettant de réduire ces oscillations, en utilisant dans la simulation Monte Carlo les échantillons de la stratification Sudoku.

Lécot et El Khettabi ont proposé une version quasi-Monte Carlo de cette méthode dans [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF]. Nous rappelons plus loin un résultat de convergence pour cette approche.

Algorithme de simulation Monte Carlo stratifié en milieu homogène

On suppose que le coefficient de diffusion D est constant et que la donnée initiale vérifie

R c 0 (x)dx = 1.
En utilisant la solution fondamentale E de l'opérateur

∂ ∂t -D ∂ 2 ∂x 2 , on obtient que pour tout τ ≥ 0, la solution c de (4.1.1) vérifie c(x, t) = R E(x -w, t -τ )c(w, τ )dw, x ∈ R, t > τ. (4.3.1)
Pour la simulation, on choisit un entier n ∈ N * et on note N = n 2 . On fixe un pas de temps ∆t. On cherche une approximation de la solution c(x, t p ) aux instants t p = p∆t par une somme de masses de Dirac en des positions x p 1 , . . . , x p N à déterminer ; on note cette approximation :

c p (x) := 1 N N k=1 δ(x -x p k ).
La première étape de l'algorithme consiste à échantillonner la distribution initiale c 0 en construisant un ensemble de points {x 0 1 , . . . , x 0 N } pour que c 0 (x) ≈ c 0 (x).

Cela peut s'effectuer par la méthode d'inversion :

x 0 k = C -1 0 2k -1 2N , 1 ≤ k ≤ N,
où C 0 est la fonction de répartition associée à c 0 :

C 0 (x) := x -∞ c 0 (y)dy, x ∈ R.
On suppose qu'à l'instant t p on connait l'approximation de c(x, t p ) :

c p (x) := 1 N N k=1 δ(x -x p k ).
La détermination des positions des particules à l'instant t p+1 se fait en deux étapes :

1. renuméroter les particules de l'étape p suivant leur position ;

2. effectuer une intégration Monte Carlo en utilisant un ensemble de points stratifié.

Nous allons développer chacune de ces étapes.

1. Renumérotation des particules. Les particules sont ordonnées selon leurs positions croissantes

x p 1 ≤ x p 2 ≤ • • • ≤ x p N .
À chaque pas de temps il faut effectuer une intégration numérique, le tri permet de minimiser l'amplitude des sauts de la fonction à intégrer. L'idée de tri a été introduite par C. Lécot pour la résolution de l'équation de Boltzmann, d'abord pour un problème uni-dimensionnel [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | Low discrepancy sequences for solving the Boltzmann equation[END_REF][START_REF] Lécot | A quasi-Monte Carlo method for the Boltzmann equation[END_REF] puis généralisée au cas multi-dimensionnel [START_REF] Lécot | A quasi-Monte Carlo scheme using nets for a linear Boltzmann equation[END_REF].

Elle a aussi été utilisée pour la simulation quasi-Monte Carlo de la diffusion [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF][START_REF] Lécot | Grid-free simulation of convection-diffusion[END_REF][START_REF] Lécot | Quasi-random walk methods[END_REF][START_REF] Lécot | Quasirandom walk methods[END_REF] et des chaînes de Markov [START_REF] Lécot | Comparison of quasi-Monte Carlo-based methods for the simulation of Markov chains[END_REF][START_REF] Lécot | Quasi-Monte Carlo methods for estimating transient measures of discrete time Markov chains[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of discrete-time Markov chain on multidimensional state space[END_REF][START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF][START_REF] El Haddad | Quasi-Monte Carlo methods for Markov chains with continuous multi-dimensional state space[END_REF] ; pour une approche QMC hasardisée de la résolution des chaînes de Markov, voir [START_REF] L'ecuyer | Randomized quasi-Monte Carlo simulation of Markov chains with an ordered state space[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF][START_REF] L'ecuyer | On array-RQMC for Markov chains : mapping alternatives and convergence rates[END_REF].

2. Intégration Monte Carlo stratifiée. En utilisant (4.3.1), on obtient une approximation de la solution à l'instant t p+1 :

cp+1 (x) = 1 N N k=1 E(x -x p k , ∆t).
Cette approximation n'est pas une somme de mesures de Dirac : pour l'obtenir, on effectue une quadrature Monte Carlo. Soit

f (u) := √ 2D∆tΦ -1 (u), u ∈ (0, 1),
où Φ est la fonction de répartition de la loi normale N (0, 1). Pour tout sous-ensemble mesurable A de R, on note 1 A sa fonction indicatrice. On a alors

R cp+1 (x)1 A (x)dx = 1 N N k=1 I 1 A (x p k + f (u))du, (4.3.2)
Pour 1 ≤ k ≤ N , soit 1 I k la fonction indicatrice de l'intervalle

I k := k -1 N , k N .
On transforme le membre de droite de (4.3.2) en une intégrale double :

R cp+1 (x)1 A (x)dx = I 2 C p+1 A (u)du, (4.3.3)
où

C p+1 A (u) := N k=1 1 I k (u 1 )1 A (x p k + f (u 2 )), u = (u 1 , u 2 ) ∈ I × (0, 1).
On considère un ensemble de points stratifié Sudoku : (P1) il y a un point unique de l'ensemble dans chaque intervalle de la forme

W = {W ℓ : ℓ = (ℓ 1 , ℓ 2 ), 1 ≤ ℓ 1 , ℓ 2 ≤ n} défini par W ℓ = (W ℓ,1 , W ℓ,2 ) := ℓ 1 -1 n + σ 1 (ℓ 2 ) -1 n 2 + ξ 1 ℓ n 2 , ℓ 2 -1 n + σ 2 (ℓ 1 ) -1 n 2 + ξ 2 ℓ n 2 , où σ 1 et σ 2 sont
ℓ 1 -1 n , ℓ 1 n × ℓ 2 -1 n , ℓ 2 n , 1 ≤ ℓ 1 ≤ n, 1 ≤ ℓ 2 ≤ n;
(P2) il y a un point unique de l'ensemble dans chaque intervalle de la forme

I × k -1 n 2 , k n 2 et k -1 n 2 , k n 2 × I, 1 ≤ k ≤ n 2 .
On définit l'approximation c p+1 (x) de la solution exacte c(x, t p+1 ) en faisant une intégration MC de l'intégrale du membre de droite de (4.3.3) à l'aide de cet ensemble : pour tout sous-ensemble mesurable

A ⊂ R R 1 A (x)c p+1 (x) = 1 N n ℓ 1 =1 n ℓ 2 =1 C p+1 A (W ℓ ).
Cette équation permet de calculer facilement les nouvelles positions des particules de simulation. On établit une correspondance entre les points de W et les positions ordonnées x p k . Pour u ∈ [0, 1), on note k(u) := ⌊N u⌋, où ⌊•⌋ désigne la partie entière. Les positions des particules à l'instant t p+1 sont alors définies par :

x p+1 k(W ℓ,1 ) = x p k(W ℓ,1 ) + f (W ℓ,2 ), ℓ = (ℓ 1 , ℓ 2 ), 1 ≤ ℓ 1 , ℓ 2 ≤ n. (4.3.4)
Ainsi, la première projection W ℓ,1 de l'ensemble stratifié est utilisée pour sélectionner la particule à déplacer, tandis que la deuxième projection W ℓ,2 détermine le déplacement aléatoire. Remarquons que les propriétés de répartition de l'ensemble stratifié Sudoku assurent que chaque particule est considérée une seule fois pour une diffusion durant chaque pas de temps.

À la place de l'ensemble stratifié Sudoku W, on peut prendre les ensembles stratifiés étudiés au chapitre 2 : l'ensemble stratifié simple de l'approche MCS, ou l'ensemble de la stratégie LHS. Dans ces deux approches, on continue de réordonner les particules par position à chaque pas de temps. Les trois méthodes (SS, MCS et LHS) sont comparées dans des essais numériques plus loin.

La méthode de Monte Carlo classique de simulation de l'équation de diffusion consiste à approcher la distribution initiale puis à faire évoluer les particules selon l'équation

x p+1 k = x p k + f (U k ), 1 ≤ k ≤ N, (4.3.5)
où {U 1 , . . . , U N } est une famille indépendante de variables aléatoires de loi uniforme sur I. Dans cette simulation les particules ne sont pas triées. Un algorithme QMC a été proposé et analysé dans [START_REF] Morokoff | A Quasi-Monte Carlo approach to particle simulation of the heat equation[END_REF][START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF], où une borne d'erreur théorique est prouvée. Pour des entiers b ≥ 2 et m ≥ 1 on pose N := b m le nombre de particules. Soit Y = {y 0 , y 1 , . . .} ⊂ I 2 une suite à faible discrépance. On note :

• Y p = {y q : pN ≤ q < (p + 1)N } ;

• π ′ la projection définie par π ′ (y 1 , y 2 ) := y 1 et π ′′ celle définie par π ′′ (y 1 , y 2 ) := y 2 .

On suppose que l'ensemble des points π ′ (Y p ) est un réseau-(0, m, 1) en base b (voir [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]) et que π ′′ (Y p ) ⊂ (0, 1).

Soit X 0 := {x 0 1 , x 0 2 , . . . , x 0 N } un ensemble de particules tel que la densité c 0 soit approchée par la somme de mesures de Dirac

c 0 (x) := 1 N N k=1 δ(x -x 0 k ) ≃ c 0 (x).
Soit X p l'ensemble de particules à l'instant t p , obtenu en faisant des quadratures quasi-Monte Carlo dans l'algorithme précédent : on remplace les ensembles de points stratifiés Sudoku W par les segments Y p de la suite Y . La notion de discrépance à l'origine par rapport à une densité de probabilité sur R généralise la discrépance à l'origine (par rapport à la mesure de Lebesgue sur I).

Proposition 4.1 La discrépance à l'origine de X p par rapport à c p vérifie :

D ⋆ N (X p ; c p ) ≤ D ⋆ N (X 0 ; c 0 ) + b ⌊m/2⌋ p-1 l=0 D N (Y l ) + p b ⌊m/2⌋ .
Ce résultat indique un ordre (déterministe) de convergence en O(1/N 1/2 ), le même que l'ordre (stochastique) de la méthode de Monte Carlo. Les résultats des expériences numériques repris dans [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF] montrent un ordre en O(1/N 0,71 ). De plus la majoration de la proposition précédente croît linéairement avec le temps. Ce n'est pas le comportement de la discrépance calculée dans les mêmes expériences numériques.

Essais numériques

Dans cette section, nous vérifions la pertinence de l'algorithme défini précédemment en reprenant des essais de [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF] sur l'équation de la diffusion. On considère les équations (4.1.1) et (4.1.2) avec la condition initiale

c 0 (x) := 1 √ π e -x 2 , x ∈ R,
J := a 0 c(x, T )dx, pour a = 4 et T = 1.
La discrétisation temporelle se fait avec un pas de temps ∆t = 1/100. On compare les résultats obtenus par la méthode MC ordinaire à ceux que l'on trouve en utilisant trois stratégies de stratification : MCS, LHS et SS. Pour estimer les variances des estimateurs, on construit un échantillon de valeurs approchées de J de taille M = 5 000 puis on calcule leur variance empirique. Les résultats obtenus par la méthode MC classique, la méthode de stratification simple (MCS), celle de l'échantillonnage par hypercubes latins (LHS) présentées au chapitre 2 et la méthode stratifiée Sudoku (SS) décrite précédemment, sont représentés dans la Figure 4.2. Dans ces simulations, le nombre de particules mises en jeu est de la forme N = (10m) 2 avec 1 ≤ m ≤ 20. Sur la même figure, nous avons comparé les temps de calcul pour les différentes méthodes en fonction du nombre de points de la simulation.

L'efficacité est définie comme l'inverse du produit de la variance par le temps de calcul [START_REF] Hammersley | Monte Carlo Methods[END_REF][START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF]. Nous avons comparé les efficacités des méthodes MC, LHS, MCS et SS : voir la Figure 4.3. Il apparaît que celle de la méthode SS est la meilleure.

On a aussi comparé ces méthodes à la méthode hasardisée proposée par Owen (OS) [START_REF] Owen | Variance with alternative scrambling of digital nets[END_REF] et qui consiste à brouiller une suite à faible discrépance en utilisant des permutations affines (voir chapitre 1). La suite à faible discrépance utilisée est la suite de Faure en base b = 2 et en dimension 2. Les courbes des variances sont tracées dans la Figure 4.4 en fonction du nombre de points N , N = 2 5 , . . . , 2 15 , pour cette méthode. On remarque que la méthode OS fournit la variance la plus petite. Ceci est dû au fait que le brouillage conserve les propriétés de répartition de la suite originale : la suite ainsi obtenue est aussi une suite (0, 2) en base 2.

Si l'on suppose que la variance est de l'ordre de N -β , on peut estimer la valeur de β en effectuant une régression linéaire. Les ordres des variances des différentes méthodes sont donnés dans la Table 4.1.

Pour la méthode SS, nous retrouvons un ordre de variance proche de l'ordre de N -3/2 , qui est celui obtenu théoriquement pour le problème de l'intégration numérique en dimension 2, présenté au chapitre 3. Cela permet de penser que la méthode SS est bien une méthode de réduction de la variance.

Supposant que la variance est exprimée par Var = cN -β , avec c constante, on remarque un décalage entre les constantes des différentes méthodes. Ceci montre que les constantes diminuent en passant de la méthode MC aux méthodes LHS, MCS, SS et MCO respectivement. Nous avons aussi constaté que l'étape de tri est inutile quand on utilise des nombres aléatoires ordinaires (approche MC). Par ailleurs, si l'on utilise une stratégie de stratification sans trier les particules de simulation, on obtient une variance plus grande que celle de la méthode de Monte Carlo ordinaire (avec un ordre de convergence proche).

Dans la suite, nous montrons que l'approche de stratification SS permet d'améliorer les résultats dans le calcul de la diffusion des ions déjà considérée et s'adapte à des problèmes en milieu hétérogène (coefficient de diffusion non constant).

Diffusion dans un milieu hétérogène : correction du pas gaussien

On considère à présent un problème de diffusion dans un milieu hétérogène, où le coefficient de diffusion dépend de la variable d'espace : D = D(x). Dans ce cas, un déplacement de chaque particule d'un pas d'amplitude déterminée par la valeur du coefficient de diffusion au point de départ de la particule : ∆x := 2D(x)∆tZ (où Z est une variable aléatoire normale centrée réduite) conduit à des résultats erronés même pour des petites valeurs de ∆t. On observe une convection artificielle en direction des régions de faible diffusivité : voir [START_REF] Hunter | On the use of random walk models with spatially variable diffusivity[END_REF]. Cette approche intuitive n'est donc pas valable. Pour résoudre un problème biologique de diffusion d'ions [START_REF] Matveev | New and corrected simulations of synaptic facilitation[END_REF], L. Farnell et W.G. Gibson ont proposé dans [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium : correction to the Gaussian steplength[END_REF] une adaptation du schéma de marche aléatoire dans le cas d'un milieu hétérogène qui consiste à introduire un terme correcteur au pas gaussien. Nous reprenons leur méthode avant d'en proposer une variante stratifiée.

Correction du pas gaussien

Soit D(x) > 0 un coefficient de diffusion variable dépendant de la position x. Soit x 0 un réel et f X la densité de la loi gaussienne centrée et de variance 2D(x 0 )∆t :

f X (x) = e -x 2 /4D(x 0 )∆t 4πD(x 0 )∆t .
Soit ∆x = 2D(x 0 )∆tZ le pas gaussien non corrigé, Z étant une variable aléatoire réelle de loi N (0, 1). Si U est une variable aléatoire de loi uniforme sur (0, 1), on peut écrire :

U = ∆x -∞ f X (x)dx. (4.4.1)
De même, si f W (w) est la fonction de densité du pas d'espace exact ∆w, alors

U = ∆w -∞ f W (w)dw. (4.4.2)
L'idée consiste à chercher un terme de correction ǫ tel que :

∆w = ∆x + ǫ. (4.4.3)
En se basant sur les équations (4.4.1), (4.4.2) et (4.4.3), on obtient :

∆x -∞ f X (x)dx = ∆x -∞ f W (w)dw + ∆x+ǫ ∆x f W (w)dw = ∆x -∞ f W (w)dw + ǫf W (∆x) + O(ǫ 2 ).
On a donc une expression du terme de correction du premier ordre :

ǫ = ∆x -∞ f X (x)dx - ∆x -∞ f W (w)dw f W (∆x) . (4.4.4)
Ce terme n'est pas pratiquement utilisable car la distribution f W est généralement inconnue. Pour aboutir à une expression utile, on commence par étudier le cas simple d'un coefficient de diffusion linéaire, avant de passer au cas général, où l'on fait un développement de Taylor.

Cas d'un coefficient de diffusion linéaire

On suppose que

D(x) = D 0 (1 + αx),
où D 0 > 0 et α sont des constantes. Dans ce cas, pour une émission initiale ponctuelle d'intensité 1 au point 0, qui se traduit par c 0 = δ, la solution de l'équation de diffusion s'écrit :

c(x, t) = 1 D 0 |α|t e -(2+αx)/D 0 α 2 t I 0 2 √ 1 + αx D 0 α 2 t , (4.4.5)
I 0 étant la fonction de Bessel modifiée [START_REF] Gradshteyn | Table of integrals, Series, and Products[END_REF]. Par suite l'expression exacte de la fonction de densité est :

f W (w) = 1 D 0 |α|∆t e -(2+αw)/D 0 α 2 ∆t I 0 2 √ 1 + αw D 0 α 2 ∆t .
Pour un petit pas de temps ∆t et pour w = O(∆t), un développement asymptotique de I 0 au voisinage de +∞ permet d'écrire

f W (w) = e -w 2 /4D 0 ∆t √ 4πD 0 ∆t 1 - α 4 w + α 8D 0 ∆t w 3 + O(w 2 ) .
En remplaçant dans l'expression (4.4.4) on obtient

ǫ ≈ ǫ 1 := αD 0 2 ∆t + α 4 (∆x) 2 1 -α 4 ∆x + α 8D 0 ∆t (∆x) 3 . (4.4.6)
Le pas corrigé du déplacement gaussien d'une particule partant de l'origine est alors donné par ∆x + ǫ 1 , où ∆x = √ 2D 0 ∆tZ, Z étant une variable aléatoire de loi N (0, 1).

∂c ∂t (x, t) = ∂ ∂x D ∂c ∂x (x, t), x ∈ R, t > 0, (4.5.1)
c(x, 0) = c 0 (x), x ∈ R, (4.5.2)
où D = D(x). Les instants d'observations sont t p = p∆t ; l'approximation cherchée est de la forme :

c p (x) := 1 N N k=1 δ(x -x p k ),
où x p 1 , . . . , x p N sont des réels à déterminer (les positions des N particules de simulation à l'instant t p ).

L'initialisation se fait comme dans le cas d'un milieu homogène : un ensemble de N réels {x 0 1 , . . . , x 0 N } échantillonne la distribution initiale c 0 . On suppose, qu'à l'instant t p , les particules sont en x p 1 , . . . , x p N . Leur évolution à l'instant suivant se fait en les ordonnant dans un premier temps puis en effectuant une quadrature de Monte Carlo à l'aide d'un ensemble stratifié Sudoku.

1. Renumérotation des particules. Les particules sont ordonnées selon leurs positions croissantes :

x p 1 ≤ x p 2 ≤ • • • ≤ x p N .
2. Intégration Monte Carlo stratifiée. Soit

W = {W ℓ = (W ℓ,1 , W ℓ,2 ) : ℓ = (ℓ 1 , ℓ 2 ), 1 ≤ ℓ 1 , ℓ 2 ≤ n} ⊂ I 2
un ensemble stratifié Sudoku. À tout réel u ∈ (0, 1), on associe l'indice k(u) := ⌊N u⌋ ∈ {0, . . . , N -1}.

On définit, pour x ∈ R et u ∈ (0, 1) :

• le déplacement f (x, u) = 2D(x)∆tΦ -1 (u),
• le terme de correction

ǫ(x, u) = dD dx (x) ∆t 2 + dD dx (x) (f (x, u)) 2 4D(x) 1 - dD dx (x) f (x, u) 4D(x) + dD dx (x) (f (x, u)) 3 8(D(x)) 2 ∆t
.

Les positions des particules à l'instant t p+1 sont alors définies par :

x p+1 k(W ℓ,1 ) = x p k(W ℓ,1 ) + f x p k(W ℓ,1 ) , W ℓ,2 + ǫ x p k(W ℓ,1 ) , W ℓ,2 , pour ℓ = (ℓ 1 , ℓ 2 ), 1 ≤ ℓ 1 , ℓ 2 ≤ n.

Essais numériques

Dans cette partie, nous vérifions la validité de l'algorithme stratifié Sudoku en faisant des expériences numériques. Nous montrons que la technique de stratification approche mieux la solution du problème que le shéma MC classique. Ces essais, repris de [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium : correction to the Gaussian steplength[END_REF], portent sur l'étude de la diffusion des ions de calcium dans un milieu biologique dans lequel la source émet périodiquement des groupes de particules à des instants précis [START_REF] Matveev | New and corrected simulations of synaptic facilitation[END_REF].

Les mêmes problèmes ont été considérés dans [START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium[END_REF]. Nous traitons l'équation de diffusion dans deux cas :

1. le coefficient de diffusion est constant ;

2. le coefficient de diffusion varie en fonction de la position.

Dans le premier cas, nous pouvons comparer les résultats à la solution exacte ; dans le second, la solution de référence est obtenue par un schéma de MC fin (avec un grand nombre de particules). Nous estimons l'évolution de la concentration c des ions dans différentes régions de l'espace. La concentration est calculée en faisant le quotient du nombre de particules dans la région par le nombre total de particules émises. Dans nos essais numériques, nous considérons l'émission de particules par une source ; cette émission est exactement représentée par le schéma (sans approximation numérique).

Coefficient de diffusion constant

Coefficient de diffusion variable

Nous considérons le coefficient de diffusion variable suivant : Dans ce cas, on remarque une déviation considérable des résultats du calcul MC par rapport à la courbe de référence qui devient de plus en plus importante au cours du temps, alors que la méthode SS réduit nettement les oscillations des résultats.

D(x) = D 0 (1 -0.
Une étape importante dans l'algorithme stratifié proposé est le tri des particules de la simulation à chaque pas de temps. Des essais numériques, dont les résultats ne sont pas reportés ici, ont montré que l'utilisation de la technique de stratification sans renuméroter les particules donne des variances supérieures à celles de la méthode MC, avec un ordre de convergence de N -1 . La méthode perd donc son efficacité, pour la même raison que dans le paragraphe précédent : l'amplitude des sauts des fonctions que l'on intègre n'est pas minimisée.

Conclusion

Dans ce chapitre nous avons proposé un schéma de Monte Carlo stratifié pour la simulation de l'équation de diffusion. Ce schéma initialement mis au point pour un coefficient de diffusion constant peut aussi être utilisé dans le cas d'un milieu hétérogène, à coefficient de diffusion non constant, à condition d'introduire un teme correcteur dans les déplacements aléatoires des particules. Les essais numériques ont montré que cette technique améliore considérablement les résultats obtenus par un schéma MC classique. Simulation de la coagulation par la méthode stratifiée Sudoku

Introduction

La coagulation est le processus d'adhésion ou de fusion de particules entre elles, après qu'elles soient entrées en contact. Ce contact peut être dû à un mouvement aléatoire des particules ou causé par différentes forces internes moléculaires ou d'attraction (magnétique ou électrique) qui existent entre ces particules. En plus du processus de coagulation, le nombre et la taille des particules peuvent évoluer sous l'effet de plusieurs phénomènes tels que la nucléation, la condensation, la fragmentation et la déposition.

Les modèles de coagulation apparaissent dans les domaines des sciences et technologies comme l'étude des aérosols et des polymères réactifs et en ingénierie chimique, notamment dans la formation des sprays. Les sprays sont des fluides produits par l'éjection de gouttes de liquide dans un milieu gazeux, par des dispositifs appelés atomiseurs. Ils sont utilisés dans un grand nombre de processus industriels : bio-technologies, pharmacie, impression par jet d'encre, peinture, fabrication de piles à combustible. L'efficacité de l'atomisation dépend des caractéristiques de pulvérisation, parmi lesquelles la taille des gouttes éjectées constitue une des plus importantes caractéristiques du spray. Une méthode classique de modélisation de la distribution des tailles des gouttes est empirique : étant donnée une collection de distributions standards, on trouve une forme qui s'adapte aux données collectées pour une série d'atomiseurs. Cette approche est indépendante du temps, mais la distribution des tailles des gouttes peut changer suite à un choc ou une fusion entre des gouttes du spray.

La modélisation de la coagulation a commencé avec les travaux de von Smoluchowski [START_REF] Smoluchowski | Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen[END_REF][START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetic kolloider Lösungen[END_REF]. Il a proposé de la représenter par une infinité d'équations différentielles, qui décrivent l'évolution de la taille des particules (supposée discrète). L'article de Müller [START_REF] Müller | Zur allgemeinen Theorie der raschen Koagulation[END_REF] étend le travail précédent en considérant une distribution continue des tailles, ce qui conduit à une modélisation par une équation intégro-différentielle. S'il existe des méthodes déterministes de résolution de ces équations [START_REF] Filbet | Numerical simulation of the Smoluchowski coagulation equation[END_REF], la part des méthodes de Monte Carlo est importante [START_REF] Sabelfeld | Stochastic algorithms for solving Smolouchowsky coagulation equation and applications to aerosol growth simulation[END_REF][START_REF] Kolodko | Convergence of a Nanbu type method for the Smoluchowski equation[END_REF][START_REF] Kolodko | A stochastic method for solving Smoluchowski's coagulation equation[END_REF][START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF][START_REF] Eibeck | An efficient stochastic algorithm for studying coagulation dynamics and gelation phenomena[END_REF][START_REF] Debry | A stochastic approach for the numerical simulation of the general dynamics equation for aerosols[END_REF].

Dans ce chapitre un algorithme de simulation de l'équation de coagulation est présenté. Il suit les étapes principales de celui de Babovsky [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF] (lui-même inspiré du schéma de Nanbu [START_REF] Nanbu | Direct simulation scheme derived from Boltzmann equation, I. Monocomponent gases[END_REF] pour l'équation de Botzmann des gaz raréfiés) : la densité initiale de masse est approchée par une somme de mesures de Dirac, le temps est discrétisé et à chaque pas de temps, une simulation numérique de la coagulation est effectuée. Les tailles des particules évoluent selon la dynamique décrite par le système d'équations différentielles de Smoluchowski [START_REF] Smoluchowski | Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen[END_REF][START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetic kolloider Lösungen[END_REF] dans le cas discret ou par une équation intégro-différentielle [START_REF] Müller | Zur allgemeinen Theorie der raschen Koagulation[END_REF] dans le cas continu ; dans les deux cas on utilise un ensemble de points stratifié Sudoku. Une rénumérotation des particules est essentielle à chaque étape de la simulation. Nous appelons cette méthode SS. Des essais numériques montrent que les résultats obtenus par cet algorithme sont plus précis que ceux obtenus par une méthode de Monte Carlo classique, cela dans les deux cas discret et continu. L'avantage de la simulation SS est de pouvoir garder constant le nombre de particules numériques pendant la simulation (ce qui n'est pas le cas des méthodes de [START_REF] Sabelfeld | Stochastic algorithms for solving Smolouchowsky coagulation equation and applications to aerosol growth simulation[END_REF][START_REF] Kolodko | Convergence of a Nanbu type method for the Smoluchowski equation[END_REF][START_REF] Kolodko | A stochastic method for solving Smoluchowski's coagulation equation[END_REF]), avec un taux de convergence meilleur que celui de la méthode de [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF].

Nous présentons d'abord l'équation de coagulation discrète et son schéma de simulation numérique de type SS. Nous testons cette méthode par des essais numériques pour des données initiales particulières. Nous proposons ensuite un algorithme de simulation de l'équation de coagulation continue et il est également validé par des tests numériques. Une partie des résultats de ce chapitre a été publiée dans [START_REF] El Haddad | Simulation of coalescence with stratified sampling[END_REF].

Algorithme MC stratifié Sudoku pour l'équation de coagulation discrète

La modélisation des processus de coagulation a été abordée pour la première fois en 1916 par Smoluchowski dans ses ouvrages sur la coagulation des colloïdes [START_REF] Smoluchowski | Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen[END_REF][START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetic kolloider Lösungen[END_REF]. Il a considéré l'évolution d'un système de particules de tailles discrètes i = 1, 2, 3, . . . (la masse d'une particule étant supposée proportionnelle à sa taille) qui interagissent par coagulation. On désigne par c(i, t) la concentration des particules de taille i à l'instant t.

Le nombre total de particules à l'instant initial est noté N 0 et le nombre de particules de taille i à l'instant t est N 0 c(i, t). La rapidité du phénomène de coagulation est modélisée par un noyau de coagulation K(i, j) qui exprime le taux (ou la probabilité) de coagulation entre une particule de masse i et une autre de masse j, pour former une particule de masse i + j. La concentration c est solution du système d'équations différentielles :

∂c ∂t (i, t) = 1 2 i-1 j=1 K(i -j, j)c(i -j, t)c(j, t) - ∞ j=1 K(i, j)c(i, t)c(j, t), i ∈ N * , t > 0, (5.2.1)
qui décrit l'évolution, en fonction du temps, de la densité c(i, t) du nombre de particules de taille i à l'instant t.

Nous supposons que les données initiales satisfont :

c(i, 0) = c 0 (i) ≥ 0, i ≥ 1 et i≥1 c(i, 0) = 1. (5.2.2)
Nous supposons dans la suite que le noyau de coagulation K(i, j) est positif et symétrique

K(i, j) = K(j, i) ≥ 0, i, j ≥ 1.
À chaque instant, la concentration c(i, t) des particules de masse i peut :

• augmenter par coagulation de deux particules de tailles ij et j < i.

• diminuer par coagulation d'une particule de taille i avec d'autres particules.

En multipliant l'équation (5.2.1) par i et en sommant sur i ∈ N * , on obtient

d dt ∞ i=1 ic(i, t) = 1 2 ∞ i=1 i-1 j=1 (i -j + j)K(i -j, j)c(i -j, t)c(j, t) - ∞ i=1 ∞ j=1 iK(i, j)c(i, t)c(j, t), soit d dt ∞ i=1 ic(i, t) = ∞ i=1 i-1 j=1 (i -j)K(i -j, j)c(i -j, t)c(j, t) - ∞ i=1 ∞ j=1 iK(i, j)c(i, t)c(j, t) = ∞ h=1 ∞ j=1 hK(h, j)c(h, t)c(j, t) - ∞ i=1 ∞ j=1 iK(i, j)c(i, t)c(j, t) = 0.
Cela exprime la conservation de la masse ∞ i=1 ic(i, t) du système. Une approche intuitive de simulation MC consiste à considérer un grand nombre de particules numériques qui évoluent comme les particules physiques et peuvent donc coaguler selon la dynamique de l'équation. Un désavantage de cette méthode est qu'au cours de la simulation, le nombre de particules diminue à cause de la coagulation, ce qui détériore la convergence du schéma. L'injection de nouvelles particules de simulation s'avère alors nécessaire pour résoudre ce problème [START_REF] Sabelfeld | Stochastic algorithms for solving Smolouchowsky coagulation equation and applications to aerosol growth simulation[END_REF]. L'approche de [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF] permet de garder un nombre de particules constant. On note

C 1 := i≥1 ic(i, 0)
et on change la fonction inconnue

f (i, t) := ic(i, t) C 1 . Alors ∂f ∂t (i, t) = 1 2 i-1 j=1 (i -j + j) 1 C 1 K(i -j, j)c(i -j, t)c(j, t) - ∞ j=1 i 1 C 1 K(i, j)c(i, t)c(j, t) = i-1 j=1 C 1 K(i -j, j)f (i -j, t) f (j, t) j - ∞ j=1 C 1 K(i, j)f (i, t) f (j, t) j .
1. Renumérotation des particules. On effectue un tri des particules de l'instant t p selon leurs tailles croissantes :

j p 1 ≤ j p 2 ≤ • • • ≤ j p N .
Cette renumérotation est essentielle pour la convergence de la méthode. En effet, l'algorithme peut être décrit, à chaque étape, par une suite d'intégrations numériques et le tri vise à minimiser l'amplitude des sauts des fonctions à intégrer.

Ce type de tri a déjà été utilisé pour la simulation de la diffusion au chapitre 4. Dans le contexte de la coagulation, il a été employé pour la résolution QMC de l'équation de coagulation discète [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF].

2. Intégration numérique. On discrétise l'équation (5.2.4) à l'aide d'un schéma d'Euler explicite et on définit une mesure transitoire g p+1 par :

1 ∆t ∞ i=1 g p+1 (i)s E (i) - ∞ i=1 f p (i)s E (i) = ∞ i=1 ∞ j=1 K(i, j)f p (i)f p (j) (s E (i + j) -s E (i)) ,
pour tout E ⊂ N * . En remplaçant f p par son expression, on obtient :

∞ i=1 g p+1 (i)s E (i) = 1 N N k 1 =1 1 - ∆t N N k 2 =1 K(j p k 1 , j p k 2 ) s E (j p k 1 ) + ∆t N 2 N k 1 =1 N k 2 =1 K(j p k 1 , j p k 2 )s E (j p k 1 + j p k 2 ). ( 5
I 3 = [0, 1) 3 . Soit 1 k 1 ,k 2 la fonction indicatrice de k 1 -1 N , k 1 N × k 2 -1 N , k 2 N et I p k 1 ,k 2 la fonction indicatrice de 0, ∆t K(j p k 1 , j p k 2 ) . Pour tout E ⊂ N * , on note Γ p+1 E (u) := N k 1 =1 N k 2 =1 1 k 1 ,k 2 (u 1 , u 2 ) 1 -I p k 1 ,k 2 (u 3 ) s E (j p k 1 ) + I p k 1 ,k 2 (u 3 )s E (j p k 1 + j p k 2 ) , pour u = (u 1 , u 2 , u 3 ) ∈ I 3 . L'équation (5.2.6) s'écrit alors i≥1 g p+1 (i)s E (i) = I 3
Γ p+1 E (u)du.

(5.2.7)

On considère un ensemble de points stratifié Sudoku : (P1) il y a un point unique de l'ensemble dans chaque intervalle de la forme

W = {W ℓ : ℓ = (ℓ 1 , ℓ 2 , ℓ 3 ), 1 ≤ ℓ 1 , ℓ 2 , ℓ 3 ≤ n} défini par W ℓ = (W ℓ,1 , W ℓ,2 , W ℓ,3 ) := ℓ 1 -1 n + σ 1 (ℓ 2 , ℓ 3 ) -1 n 3 + ξ 1 ℓ n 3 , ℓ 2 -1 n + σ 2 (ℓ 1 , ℓ 3 ) -1 n 3 + ξ 2 ℓ n 3 , ℓ 3 -1 n + σ 3 (ℓ 1 , ℓ 2 ) -1 n 3 + ξ 3 ℓ n 3 , où σ 1 , σ
ℓ 1 -1 n , ℓ 1 n × ℓ 2 -1 n , ℓ 2 n × ℓ 3 -1 n , ℓ 3 n , pour 1 ≤ ℓ 1 ≤ n, 1 ≤ ℓ 2 ≤ n, 1 ≤ ℓ 3 ≤ n ;
(P2) il y a un point unique de l'ensemble dans chaque intervalle de la forme

I 2 × k -1 n 3 , k n 3 , I × k -1 n 3 , k n 3 × I, k -1 n 3 , k n 3 × I 2 . pour 1 ≤ k ≤ n 3 .
En effectuant une quadrature dans le membre de droite de (5.2.7) à l'aide de l'ensemble W, l'approximation f p+1 est déterminée par :

∀E ⊂ N * i≥1 f p+1 (i)s E (i) = 1 N ℓ Γ p+1 E (W ℓ ).
À chaque pas de temps, un nouvel ensemble de points stratifié Sudoku indépendant des précédents est utilisé.

À u ∈ I, on associe l'indice k(u) := ⌊N u⌋.

Les tailles des particules à l'étape p + 1 sont alors définies par :

j p+1 k(W ℓ,1 ) = j p k(W ℓ,1 ) + j p k(W ℓ,2 ) si W ℓ,3 < ∆t K(j p k(W ℓ,1 ) , j p k(W ℓ,2 ) ), j p k(W ℓ,1 )
sinon.

(5.2.8)

pour tout ℓ = (ℓ 1 , ℓ 2 , ℓ 3 ), avec 1 ≤ ℓ 1 ≤ n, 1 ≤ ℓ 2 ≤ n, 1 ≤ ℓ 3 ≤ n. Les composantes W ℓ,1 et W ℓ,
P c := ∆t K j p k(W ℓ,1 ) , j p k(W ℓ,2 ) .
La troisième composante W ℓ,3 est utilisée pour déterminer un événement :

• si 0 < W ℓ,3 ≤ P c , les deux particules fusionnent pour former une particule de masse j p k(W ℓ,1 ) + j p k(W ℓ,2 ) ; • si P c < W ℓ,3 < 1, il n'y a pas de coagulation.

Remarquons que les propriétés de répartition de l'ensemble stratifié Sudoku assurent que chaque particule sera considérée une fois comme particule primaire et une fois comme particule partenaire et cela à chaque pas de temps. L'approche Monte Carlo consiste à prendre un système de particules qui interagissent et forment des particules plus grandes en fonction de la dynamique décrite par l'équation de coagulation. À chaque pas de temps et pour chaque particule, un partenaire de coagulation est choisi de manière aléatoire dans l'ensemble des particules et une variable aléatoire de loi uniforme dans (0, 1) est utilisée pour déterminer si la coagulation a lieu ou non. Dans cette approche, aucun tri des particules n'est effectué. Plus précisément, la détermination des tailles des particules à partir des tailles précédentes se fait par :

j p+1 k = j p k + j p L k si U k < ∆t K(j p k , j p L k ), j p k sinon,
où L 1 , . . . , L N sont des variables aléatoires indépendantes uniformément distribuées dans {1, 2, . . . , N } et U 1 , . . . , U N sont des variables aléatoires indépendantes uniformément distribuées dans

[0, 1]. Les variables {L k , U k : 1 ≤ k ≤ N } sont supposées indépendantes.
Un algorithme QMC a été proposé dans [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF], où une borne d'erreur théorique est prouvée. Pour des entiers b ≥ 2 et m ≥ 1 on pose N := b m le nombre de particules. Soit X = {x 0 , x 1 , . . .} ⊂ I 3 une suite à faible discrépance ; on suppose que X est une suite-(t, 3) en base b pour t ≥ 0 (voir [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). On note

• X p = {x q : pN ≤ q < (p + 1)N } ; • π ′ la projection définie par π ′ (x 1 , x 2 , x 3 ) := (x 1 , x 2 ).
On suppose que l'ensemble des points π ′ (X p ) est un réseau-(0, m, 2) en base b.

Soit J 0 := {j 0 1 , j 0 2 , . . . , j 0 N } un ensemble de particules tel que la densité f 0 soit approchée par la somme de mesures de Dirac

f 0 := 1 N N k=1 δ j 0 k ≃ f 0 .
Soit J p l'ensemble de particules à l'instant t p , obtenu en faisant des quadratures quasi-Monte Carlo dans l'algorithme précédent : on remplace les ensembles de points stratifiés Sudoku W par les segments X p de la suite X. La notion de variation inférieure est l'analogue de la variation au sens de Hardy et Krause pour des suites à ensemble d'indices multi-dimensionnels. La discrépance à l'origine par rapport à une suite de termes positifs de somme 1 est une variante de la discrépance à l'origine (par rapport à la mesure de Lebesgue). Le résultat suivant est établi dans [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF],

Proposition 5.1 On suppose que

• pour tout i ≥ 1, la fonction t → f (i, t) est deux fois continûment differentiable sur (0, T ) et les séries i≥1 T 0 |f (i, t)|dt, i≥1 T 0 ∂f ∂t (i, t) dt, i≥1 T 0 ∂ 2 f ∂t 2 (i, t) dt convergent,
• K est à variation inférieure bornée V ⋆ ( K).

Alors la discrépance à l'origine de J p par rapport à f p vérifie :

D ⋆ N (J p ; f p ) ≤ e ctp D ⋆ N (J 0 ; f 0 ) + ∆t i≥1 tp 0 e c(tp-t) ∂ 2 f ∂t 2 (i, t) dt + 4V ⋆ ( K) + 3 K ∞ + 2 ∆t e ctp cb ⌊(m-t)/3⌋ , avec c := sup i≥1 V ( K(i, .)) + sup j≥1 V ( K(., j)) + 3 K ∞ .
Ce résultat indique une majoration (déterministe) de convergence en O(1/N 1/3 ), plus mauvaise que l'ordre (stochastique) de la méthode de Monte Carlo. Les résultats des expériences numériques repris dans [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF] montrent un ordre en O(1/N 0,67 ), ce qui pourrait suggérer un ordre théorique O(1/N 2/3 ). De plus la majoration de la proposition précédente croît linéairement avec le nombre de pas de temps. Ce n'est pas le comportement de la discrépance calculée dans les mêmes expériences numériques.

Essais numériques : cas de l'équation discrète

Dans ce paragraphe, nous testons la validité du schéma SS présenté précédemment en comparant les résultats obtenus par cette approche à ceux obtenus par une méthode MC classique.

Les noyaux de coagulation classiques de la littérature sont les noyaux constant, linéaire et multiplicatif donnés respectivement par :

K 0 (i, j) = 1, K 1 (i, j) = i + j et K 2 (i, j) = ij,
puisque dans ces cas, si la distribution intitale est une mesure de Dirac, les solutions de l'équation de Smoluchowski peuvent être explicitées : voir l'article de Aldous [START_REF] Aldous | Deterministic and stochastic models for coalescence (aggregation and coagulation) : a review of the mean-field theory for probabilists[END_REF].

Nous considérons ici les noyaux K 0 et K 1 avec une distribution initiale f 0 = δ 1 . Remarquons que le noyau modifié K 1 n'est pas borné. On a C 1 = 1. Le moment d'ordre µ ≥ 0 de la densité c est défini par

C µ (t) := i≥1 i µ c(i, t) = i≥1 i µ-1 f (i, t).
• Si K 0 (i, j) = 1 alors la solution exacte de (5.2.3) est donnée par :

f (i, t) = 4i (t + 2) 2 t t + 2 i-1 , i ≥ 1, t ≥ 0, ce qui entraîne C 0 (t) = 2 2 + t et C 2 (t) = 1 + t.
• Si K 1 (i, j) = i + j alors la solution exacte de (5.2.3) est donnée dans [START_REF] Golovin | The solution of the coagulating equation for cloud droplets in a rising air current[END_REF] par :

f (i, t) = i i i! (1 -e -t ) i-1 e -i(1-e -t )-t , i ≥ 1, t ≥ 0, et par suite C 0 (t) = e -t et C 2 (t) = e 2t . Comme C 0 (t) = i≥1 c(i, t) = i≥1 f (i, t) i et C 2 (t) := i≥1 i 2 c(i, t) = i≥1 if (i, t),
le schéma numérique permet de les approcher respectivement par

C 0 (t p ) ≈ 1 N N k=1 1 j p k et C 2 (t p ) ≈ 1 N N k=1 j p k .
On calcule les moments d'ordre 0 et 2 au temps T = 1, 0 avec N = (4m) 3 particules, 1 ≤ m ≤ 8 et un pas de temps ∆t = T /P où P = 4 000, par les schémas Monte Carlo usuel et Monte Carlo stratifié Sudoku. Afin d'estimer les variances de ces deux approches, nous avons répété les calculs des moments 1 000 fois puis nous avons calculé la variance empirique des valeurs ainsi obtenues. Les résultats sont représentés dans la Figure 5.1 pour le noyau constant et 5.2 pour le noyau linéaire.

Il est clair que les variances de la méthode SS sont inférieures à celles de la méthode MC, pour les deux noyaux considérés.

Si l'on suppose que la variance est de l'ordre de N -β , la valeur de β peut être estimée par une régression linéaire. Les résultats sont indiqués dans laTable 5.1.

Ces valeurs montrent que le nouvel algorithme est plus précis qu'une simulation stochastique classique.

Algorithme MC stratifié Sudoku pour l'équation de coagulation continue

Nous considérons à présent la version continue de l'équation de coagulation, introduite par Müller [START_REF] Müller | Zur allgemeinen Theorie der raschen Koagulation[END_REF] :

∂c ∂t (x, t) = 1 2 x 0 K(x -y, y)c(x -y, t)c(y, t)dy - +∞ 0 K(x, y)c(x, t)c(y, t)dy, x > 0, t > 0, (5.3.1)
avec une distribution initiale c(x, 0) = c 0 (x) ≥ 0, x > 0, et un noyau de coagulation K(x, y) qui, comme dans le cas discret, exprime la probabilité que deux particules de tailles x et y fusionnent en une particule de taille x + y. On suppose que ce noyau est positif et symétrique ∀x, y > 0, K(x, y) ≥ 0 et K(x, y) = K(y, x).

Cette modélisation de la coagulation est intéressante, quand on considère aussi des phénomènes de fragmentation : on ajoute alors des termes intégraux linéaires au membre de droite. La concentration des particules de taille dans [x, x + dx] à l'instant t est représentée par c(x, t)dx. La densité c(x, t) peut augmenter suite à la coagulation de particules de tailles xy et y < x ou diminuer par la coagulation de particules de tailles x avec d'autres. Le nombre total de particules On connaît peu de solutions exactes de l'équation ; elles se limitent à des conditions initiales spécifiques et des noyaux simples. Comme dans le cas discret, on transforme l'équation en changeant de fonction inconnue. Un algorithme quasi-Monte Carlo de simulation de l'équation continue a été proposé par Lécot et Tarhini [START_REF] Lécot | A quasi-stochastic simulation of the general dynamics equation for aerosols[END_REF][START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF]. Il est utilisé pour la modélisation des sprays dans [START_REF] Tembely | Prediction and evolution of drop-size distribution for a new ultrasonic atomizer[END_REF][START_REF] Lécot | Numerical simulation of the drop size distribution in a spray[END_REF].

Nous présentons dans ce qui suit un algorithme MC stratifié Sudoku pour la simulation de l'équation continue (5.3.1) semblable à celui de la section précédente. On introduit la densité de masse

f (x, t) := xc(x, t) C 1 , où C 1 := C 1 (0)
. En multipliant l'équation (5.3.1) par x et en effectuant le changement de fonction, on obtient

∂f ∂t (x, t) = x 0 K(x -y, y)f (x -y, t)f (y, t)dy - +∞ 0 K(x, y)f (x, t)f (y, t)dy, x > 0, t > 0, ( 5

Description du schéma

On suppose que le noyau modifié est borné et on note

K ∞ := sup x,y>0 K(x, y).
Comme dans le cas discret, le schéma peut fonctionner sans cette condition, en adaptant le pas de temps. Les étapes de l'algorithme sont les mêmes que dans le cas de l'équation discrète ; nous les rappelons.

Initialisation

Soit N un entier de la forme N := n 3 , avec n ≥ 1. On construit un ensemble de N points

X 0 = {x 0 1 , . . . , x 0 N } ⊂ R +
tel que la densité de masse initiale f 0 (x) soit approchée par :

f 0 (x) := 1 N N k=1 δ(x -x 0 k ).
Ceci peut être effectué en posant

x 0 k := F -1 0 2k -1 2N , 1 ≤ k ≤ N, où F 0 (x) := x 0 f 0 (y)dy
est la fonction de répartition associée à f 0 .

Évolution

Soit ∆t = T /P un pas de temps fixé tel que ∆t K ∞ < 1. On note t p = p∆t et f p (x) := f (x, t p ) la densité de masse à l'instant t p . Pour 0 ≤ p < P , on cherche une approximation de f p+1 de la forme

f p+1 (x) := 1 N N k=1 δ(x -x p+1 k ),
à partir des tailles des particules x p 1 , . . . , x p N de l'étape p. 1. Rénumérotation des particules. Comme dans le cas discret, les particules de l'étape p sont triées selon leur masse croissante :

x p 1 ≤ x p 2 ≤ • • • ≤ x p N .
Cette stratégie a été validée pour la simulation de la diffusion dans le chapitre 4 et par Lécot et Tarhini [START_REF] Lécot | A quasi-stochastic simulation of the general dynamics equation for aerosols[END_REF][START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF] et Lécot, Tembely, Soucemarianadin et Tarhini [START_REF] Lécot | Numerical simulation of the drop size distribution in a spray[END_REF] dans le contexte de l'approche QMC de l'équation de coagulation continue.

2. Intégration numérique. Une discrétisation d'Euler de l'équation (5.3.2) permet de définir une mesure g p+1 telle que pour toute fonction-test σ ∈ S(R)

1 ∆t +∞ 0 g p+1 (x)σ(x) - +∞ 0 f p (x)σ(x) = +∞ 0 +∞ 0 K(x, y)f p (x)f p (y)(σ(x + y) -σ(x)), c'est-à-dire, en exprimant f p : 1 ∆t +∞ 0 g p+1 (x)σ(x) = 1 N N k 1 =1 1 - ∆t N N k 2 =1 K(x p k 1 , x p k 2 ) σ(x p k 1 ) + ∆t N 2 N k 1 =1 N k 2 =1 K(x p k 1 , x p k 2 )σ(x p k 1 + x p k 2 ). Soient 1 k 1 ,k 2 et I p k 1 ,k 2 les fonctions indicatrices respectives de k 1 -1 N , k 1 N × k 2 -1 N , k 2 N et 0, ∆t K(x p k 1 , x p k 2 ) .
À σ ∈ S(R + ) on associe une fonction notée Γ p+1 σ définie par On considère un ensemble de points stratifié Sudoku :

Γ p+1 σ (u) := N k 1 =1 N k 2 =1 1 k 1 ,k 2 (u 1 , u 2 ) (1 -I p k 1 ,k 2 (u 3 ))σ(x p k 1 ) + I p k 1 ,k 2 (u 3 )σ(x p k 1 + x p k 2 ) , pour u = (u 1 , u 2 , u 3
W = {W ℓ : ℓ = (ℓ 1 , ℓ 2 , ℓ 3 ), 1 ≤ ℓ 1 , ℓ 2 , ℓ 3 ≤ n} défini par W ℓ = (W ℓ,1 , W ℓ,2 , W ℓ,3 ) := ℓ 1 -1 n + σ 1 (ℓ 2 , ℓ 3 ) -1 n 3 + ξ 1 ℓ n 3 , ℓ 2 -1 n + σ 2 (ℓ 1 , ℓ 3 ) -1 n 3 + ξ 2 ℓ n 3 , ℓ 3 -1 n + σ 3 (ℓ 1 , ℓ 2 ) -1 n 3 + ξ 3
ℓ n 3 , où σ 1 , σ 2 et σ 3 sont des bijections aléatoires de {1, . . . , n} 2 dans {1, . . . , n 2 } et ξ 1 ℓ , ξ 2 ℓ , ξ 3 ℓ sont des variables aléatoires uniformément distribuées sur I. Toutes ces variables aléatoires sont supposées indépendantes. Cet ensemble vérifie les deux propriétés d'uniformité (P1) et (P2) rappelées dans l'étude de l'équation de coagulation discrète. En effectuant une quadrature dans le membre de droite de (5.3.5) à l'aide de l'ensemble W, l'approximation f p+1 est déterminée par :

+∞ 0 f p+1 (x)σ(x) = 1 N ℓ Γ p+1 σ (W ℓ ),
pour toute fonction-test σ ∈ S(R + ). Si l'on note k(u) := ⌊N u⌋, pour u ∈ I, les tailles des particules à l'étape p + 1 sont déterminées par :

x p+1 k(W ℓ,1 ) = x p k(W ℓ,1 ) + x p k(W ℓ,2 )
si W ℓ,3 < ∆t K(x p k(W ℓ,1 ) , x p k(W ℓ,2 ) ), x p k(W ℓ,1 ) sinon, pour tout ℓ = (ℓ 1 , ℓ 2 , ℓ 3 ), avec 1 ≤ ℓ 1 ≤ n, 1 ≤ ℓ 2 ≤ n, 1 ≤ ℓ 3 ≤ n. Les composantes W ℓ,1 et W ℓ,2 des points permettent de sélectionner les particules à l'instant t p : la particule d'indice k(W ℓ,1 ) peut coaguler avec la particule d'indice k(W ℓ,2 ) avec une probabilité de coagulation P c := ∆t K x p k(W ℓ,1 ) , x p k(W ℓ,2 ) . La troisième composante W ℓ,3 est utilisée pour déterminer un événement :

• si 0 < W ℓ,3 ≤ P c , les deux particules fusionnent pour former une particule de taille x p k(W ℓ,1 ) + x p k(W ℓ,2 ) ; • si P c < W ℓ,3 < 1, il n'y a pas de coagulation.

Comme dans le cas de l'équation discrète, les propriétés de répartition de l'ensemble stratifié Sudoku assurent que chaque particule sera considérée une fois comme particule primaire et une fois comme particule partenaire et cela à chaque pas de temps. L'approche Monte Carlo consiste à prendre un système de particules qui interagissent et forment des particules plus grandes en fonction de la dynamique décrite par l'équation de coagulation. À chaque pas de temps et pour chaque particule, un partenaire de coagulation est choisi de manière aléatoire dans l'ensemble des particules et une variable aléatoire de loi uniforme dans (0, 1) est utilisée pour déterminer si la coagulation a lieu ou non. Dans cette approche, aucun tri des particules n'est effectué. La détermination des tailles des particules à partir des tailles précédentes se fait par : [START_REF] Niederreiter | Point sets and sequences with small discrepancy[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). On note Y p = {y q : pN ≤ q < (p + 1)N } et π ′ la projection définie par π ′ (y 1 , y 2 , y 3 ) := (y 1 , y 2 ).

x p+1 k = x p k + x p L k si U k < ∆t K(x p k , x p L k ), x p
On suppose que l'ensemble des points π ′ (Y p ) est un réseau-(0, m, 2) en base b.

Soit X 0 := {x 0 1 , x 0 2 , . . . , x 0 N } un ensemble de particules tel que la densité f 0 soit approchée par la somme de mesures de Dirac

f 0 := 1 N N k=1 δ(x -x 0 k ) ≃ f 0 .
Soit X p l'ensemble de particules à l'instant t p , obtenu en faisant des quadratures quasi-Monte Carlo dans l'algorithme précédent : on remplace les ensembles de points stratifiés Sudoku W par les segments Y p de la suite Y . La notion de variation de fonction sur R + est l'analogue de la variation au sens de Hardy et Krause. La discrépance à l'origine par rapport à une densité de probabilité est une variante de la discrépance à l'origine (par rapport à la mesure de Lebesgue). Le résultat suivant est établi dans [START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF] :

Proposition 5.2 On suppose :

• pour tout x > 0, la fonction t → f (x, t) est deux fois continûment differentiable sur (0, T ) et les fonctions f, ∂f ∂t , ∂ 2 f ∂t 2 sont intégrables sur R * + × (0, T ),

• K est à variation bornée au sens de Hardy et Krause. Comme dans le cas de l'équation discrète, ce résultat indique une majoration (déterministe) de convergence en O(1/N 1/3 ), plus mauvaise que l'ordre (stochastique) de la méthode de Monte Carlo. Les résultats des expériences numériques repris dans [START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF] montrent un ordre en O(1/N 0,64 ), ce qui pourrait suggérer un ordre théorique O(1/N 2/3 ). De plus la majoration de la proposition précédente croît linéairement avec le nombre de pas de temps. Ce n'est pas le comportement de la discrépance calculée dans les mêmes expériences numériques.

Essais numériques : cas de l'équation continue

Nous reprenons les essais numériques de [START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF] pour l'équation de coagulation continue. Nous comparons la performance du nouvel algorithme SS à celle du schéma MC classique en calculant les variances des estimations. On considère les deux noyaux • Noyau de coagulation constant. On considère le noyau de coagulation K 0 (x, y) = 1 et la distribution initiale c 0 (x) := e -x , x ≥ 0.

Alors C 1 = 1. Dans ce cas, la solution exacte de l'équation (5. , avec des points x P k qui peuvent être proches de 0. Cependant, l'écart entre les résultats des deux méthodes augmente avec le nombre de particules, ce qui suggère que le schéma stratifié est asymptotiquement meilleur que le schéma MC.

• Noyau de coagulation linéaire Nous considérons l'équation de coagulation avec le noyau K 1 (x, y) = x + y, et la distribution initiale c 0 (x) = e -x , x ≥ 0.

Ici encore C 1 = 1. Sous ces conditions, la solution exacte du problème s'écrit [START_REF] Ramabhadran | Dynamics of aerosol coagulation and condensation[END_REF] :

c(x, t) = 1 x √
1e -t e -x(2-e -t )-t I 1 2x √ 1e -t , x ≥ 0, t ≥ 0, (5.3.6) où I 1 est la fonction de Bessel modifiée de première espèce et du premier ordre [START_REF] Gradshteyn | Table of integrals, Series, and Products[END_REF].

Remarquons qu'ici, K 1 n'est pas borné et donc l'inégalité ∆t K 1 (x, y) < 1 n'est pas assurée. Pour que cette inégalité soit toujours vérifiée nous effectuons les calculs jusqu'à l'instant T = 0, 1.

Dans la Figure 5.4 sont représentées les variances en fonction du nombre de particules N pour les calculs des moments C 0 (T ) et C 2 (T ). Nous avons effectué M = 1 000 répétitions indépendantes de chacun des calculs, pour les deux méthodes avec un Ici encore, la variance de l'estimateur stratifié reste inférieure à celle de la méthode MC simple.

Pour estimer la vitesse de convergence des méthodes MC et SS, nous calculons l'ordre des variances en effectuant une régression linéaire. Les résultats sont regroupés dans la Table 5.2 pour 1 000 estimations. Pour le noyau constant, T = 1, 0 et pour le noyau linéaire, T = 0, 1.

Il est clair que la substitution des points aléatoires par des points stratifiés et le tri des particules de simulation améliorent la convergence.

Conclusion

Dans ce chapitre, des algorithmes de Monte Carlo stratifié Sudoku ont été développés pour la simulation des équations de coagulation discrète et continue.

Les schémas sont itératifs et consistent à approcher la distribution initiale par un ensemble de particules qui vont évoluer au cours du temps en pouvant s'assembler entre elles pour former des particules de tailles plus grandes, selon la dynamique de l'équation. L'évolution des particules se fait en utilisant des ensembles de points stratifiés Sudoku. Pour profiter des propriétés de régularité de répartition de ces ensembles stratifiés, les particules de simulation doivent être ordonnées suivant leur masse croissante à chaque pas de temps. Un avantage de ces algorithmes est que le nombre des particules reste fixe tout au long de la simulation.

Dans les essais numériques, nous avons remarqué que les méthodes proposées sont plus précises que la méthode MC classique puisque leurs variances sont inférieures (pour un même nombre de particules). Cela permet d'affirmer qu'il s'agit bien d'une technique de réduction de la variance de la méthode MC.

Conclusions et perspectives

L'objectif de ce travail était d'explorer les possibilités de méthodes de Monte Carlo stratifiées, utilisant des échantillons proches des réseaux et suites digitales. Il s'agissait de développer des schémas dont les vitesses de convergence soient meilleures que celles des méthodes MC classiques et qui permettent d'obtenir des intervalles de confiance plus accessibles que les bornes d'erreur des méthodes quasi-Monte Carlo.

Dans le premier chapitre nous avons mené une étude bibliographique sur les méthodes de Monte Carlo, les méthodes quasi-Monte Carlo et des méthodes d'hasardisation et nous les avons illustrées à l'aide d'essais numériques.

Ensuite, nous avons proposé des méthodes de Monte Carlo stratifiées. Le point de départ, au deuxième chapitre, était la méthode d'échantillonnage stratifiée simple qui consiste à diviser l'hypercube unité en des sous-cubes identiques et à choisir une variable aléatoire distribuée uniformément dans chaque sous-cube. Nous avons établi une borne de la variance de l'estimateur de la mesure d'un sous-ensemble de l'hypercube unité. L'ordre de convergence trouvé est meilleur que celui d'une méthode MC classique. Nous avons aussi analysé l'échantillonnage par hypercube latin, qui ne semble pas améliorer l'ordre de convergence de la variance.

Dans le troisième chapitre, nous avons proposé une méthode dite méthode de Monte Carlo stratifiée Sudoku, qui possède les propriétés de la méthode d'échantillonnage stratifiée simple et celles de la méthode d'échantillonnage par hypercube latin. Elle permet de construire un ensemble de points aléatoires vérifiant à la fois les propriétés de stratification multi-dimensionnelle et aussi uni-dimensionnelle. Elle retient certaines propriétés des réseaux digitaux. En considérant l'intégration de fonctions indicatrices d'ensembles, nous avons estimé la variance de cette nouvelle méthode. L'analyse théorique a montré une réduction effective de la variance de l'estimateur. La justesse de l'estimation a été validée par des expériences numériques. On a privilégié le cas des fonctions indicatrices, car on peut interpréter certaines simulations MC comme une suite de quadratures de fonctions indicatrices d'ensembles multi-dimensionnels à frontière compliquée.

Dans le quatrième et le cinquième chapitre, nous avons proposé des utilisations des techniques précédentes pour la résolution numérique de l'équation de diffusion puis de l'équation de Smoluchowski. Pour ces équations instationnaires, on discrétise le temps et on utilise des particules de simulation : à chaque pas de temps, les particules sont triées selon leurs positions ou tailles croissantes et on choisit les évènements subis par les particules à l'aide d'un ensemble stratifié Sudoku. On montre sur des exemples que cette méthode fournit des résultats plus précis qu'une méthode de marche aléatoire. Cette précision devient de plus en plus importante avec l'augmentation du nombre de particules utilisées dans la simulation.

Dans ces simulations, la propriété de stratification uni-dimensionnelle garantit la validité physique du schéma, en assurant que chacune des particules est considérée à chaque pas de temps ; la propriété de stratification multi-dimensionnelle garantit la précision numérique du schéma (de la quadrature qu'il contient).

Le travail effectué ici ouvre plusieurs perspectives de recherche. L'étude de la méthode de stratification développée mérite d'être complétée dans le but d'analyser sa convergence pour les problèmes d'équations différentielles ou intégrales abordées aux chapitres 4 et 5.

Elle a besoin d'être étendue aux cas d'équations en dimension plus élevée. Des algorithmes de tri multi-dimensionnels ont été mis au point dans l'approche quasi-Monte Carlo de la diffusion [START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF] ; une étude de différentes possibilités de tri pour la simulation QMC hasardisée a été faite dans [START_REF] L'ecuyer | On array-RQMC for Markov chains : mapping alternatives and convergence rates[END_REF].

Un développement de la méthode pour la simulation des sprays est aussi nécessaire. Il s'agira de suivre la démarche initiée dans le cadre de l'approche quasi-Montte Carlo, où l'on a montré que les fluctuations de la distribution des tailles de gouttes pouvait être réduite [START_REF] Tembely | Prediction and evolution of drop-size distribution for a new ultrasonic atomizer[END_REF][START_REF] Lécot | Numerical simulation of the drop size distribution in a spray[END_REF].

L'étude de l'adaptation de la méthode de stratification pour la simulation de la diffusion par résolution de l'équation de Feynman-Kac [START_REF] Tanré | Some new simulation schemes for the evaluation of Feynman-Kac representations[END_REF][START_REF] Tanré | Monte Carlo approximations of the Neumann problem[END_REF] peut aussi être abordée et semble prometteuse.

Un autre axe intéressant de recherche est l'adaptation de la méthode stratifiée Sudoku à la simulation des chaînes de Markov, suivant ce qui a été fait pour les méthodes QMC dans [START_REF] Lécot | Comparison of quasi-Monte Carlo-based methods for the simulation of Markov chains[END_REF][START_REF] Lécot | Quasi-Monte Carlo methods for estimating transient measures of discrete time Markov chains[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of discrete-time Markov chain on multidimensional state space[END_REF][START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] El Haddad | Quasi-Monte Carlo methods for Markov chains with continuous multi-dimensional state space[END_REF] ou bien pour les méthodes QMC hasardisées dans [START_REF] L'ecuyer | Randomized quasi-Monte Carlo simulation of Markov chains with an ordered state space[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF]. Il s'agira de proposer des algorithmes de simulation et d'étudier leur convergence.

Finalement, les travaux futurs incluront l'étude des techniques d'hasardisation des méthodes QMC, comme celles proposées par Owen [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF][START_REF] Owen | Variance with alternative scrambling of digital nets[END_REF]. Il s'agira d'estimer leur variance dans le cadre de l'intégration de fonctions indicatrices, puis de les appliquer à la résolution de problèmes d'équations différentielles ou intégrales et à la simulation de chaînes de Markov.

Introduction 1

 1 Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapitre 1 Méthodes

 1 de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique 1.1 Introduction Ce chapitre présente les méthodes de Monte Carlo (MC) pour le calcul approché d'une intégrale et introduit des techniques de réduction de la variance. Il propose également une approche élémentaire des méthodes quasi-Monte Carlo (QMC) et il décrit une technique d'hasardisation.

  der Corput dont nous détaillons la construction. Soit b ≥ 2 un entier : on note Z b := {0, 1, . . . , b -1}. À un entier n ≥ 0 on associe son développement en base b : n = ∞ j=0 a j (n)b j , où a j (n) ∈ Z b et a j (n) = 0 pour j assez grand. Définition 1.3 Si b ≥ 2 est un entier la fonction radicale inverse en base b est définie par

Définition 1 . 4

 14 Si b 1 ≥ 2, b 2 ≥ 2, . . . , b s ≥ 2 sont des entiers, la suite de Halton de bases b 1 , b 2 , . . . , b s est la suite (x n ) n≥0 où

Remarque 1 . 2

 12 des entiers. Un ensemble X de b m points dans I s est un réseau (t, m, s) en base b si A(J, X) = b t pour tout intervalle élémentaire J en base b tel que λ s (J) = b t-m . Soit t ≥ 0 un entier. Une suite infinie x 0 , x 1 , . . . de points de I s est une suite (t, s) en base b si, pour tous les entiers ℓ ≥ 0 et m ≥ t, l'ensemble {x n : ℓb m ≤ n < (ℓ + 1)b m } est un réseau (t, m, s) en base b. La Figure 1.1 représente les 16 points d'un réseau-(0, 4, 2) en base 2 ainsi que les différents intervalles élémentaires en base 2 de volume 2 -4 . On remarque que chaque intervalle contient exactement 1 point de réseau. Remarque 1.1 Si X est un réseau (t, m, s) en base b, alors pour N := b m , D N (J, X) = 0 pour tout intervalle élémentaire J en base b tel que λ s (J) = b t-m . Soit s ≥ 1 un entier.

Définition 1 . 6

 16 Soit b ≥ 2 un entier premier et s un entier tel que s ≤ b. Soient b 1 , b 2 , . . . , b s des éléments distincts de Z b . Une suite de Faure en base b en dimension s est la suite (x n ) n≥0 où

Proposition 1 . 5

 15 Soit b ≥ 2 un entier premier et s un entier tel que s ≤ b. Toute suite de Faure en base b en dimension s est une suite (0, s) en base b.

  une matrice triangulaire inférieure régulière dont les éléments sont dans F b . Ces suites sont des suites (0, s) en base b.

1 .Figure 1 . 2 -

 112 Figure 1.2 -Calcul de J 4 . Erreur des méthodes de MC, QMC et QMC hasardisées (échelles logarithmiques en base 5). Le nombre de points N varie entre 5 1 et 5 8 .

Figure 1 . 3 -

 13 Figure 1.3 -Calcul de J 8 . Erreur des méthodes de MC, QMC et QMC hasardisées (échelles logarithmiques en base 11). Le nombre de points N varie entre 11 1 et 11 5 .

Figure 1 . 4 -

 14 Figure 1.4 -Calcul de J 12 . Erreur des méthodes de MC, QMC et QMC hasardisées (échelles logarithmiques en base 13). Le nombre de points N varie entre 13 1 et .

Figure 1 . 5 -

 15 Figure 1.5 -Calcul de J 18 . Erreur des méthodes de MC, QMC et QMC hasardisées (échelles logarithmiques en base 19). Le nombre de points N varie entre 19 1 et .

Figure 1 . 6 -

 16 Figure 1.6 -Calcul de I. Erreurs par les méthodes de MC, QMC et QMC hasardisées (échelles logarithmiques en base 5). Le nombre de points N varie entre 5 1 et 5 8 .

  Figure 2.1. L'ensemble des N points ainsi construit possède la propriété d'indépendance et dans chacun des pavés de la subdivision, il y a un point unique. L'estimateur de Monte Carlo stratifié est défini par

Proposition 2 . 5

 25 On suppose que S f,y est fermée et que ∀i, j ∈ {1, 2, . . . , s} ∂f ∂x i et ∂ 2 f ∂x i ∂x j sont uniformément continues sur S f,y . On suppose qu'il existe λ > 0 et i ∈ {1, 2, . . . , s} tels que ∀x ∈ S f,y ∂f ∂x i (x) = max 1≤j≤s ∂f ∂x j (x) ≥ λ. Si la famille {C k : 1 ≤ k ≤ N } est formée d'hypercubes, alors

Lemme 2 . 1

 21 Soit f une fonction définie sur I s et à variation bornée au sens de Hardy et Krause. Soient n 1 , . . . , n s des entiers. Pour k = (k 1 , . . . , k s ) ∈ N s avec 1 ≤ k i ≤ n i , on note

Figure 2 . 4 -

 24 Figure 2.4 -Courbes des erreurs et ajustements linéaires par un calcul MC (tracés en pointillés), QMC (tracés continus) et MCS (tracés en tirets) en fonction du nombre de points N . Dimension s = 2 pour N = 10 2 , . . . , 400 2 (en haut) ; dimension s = 3 pour N = 10 3 , . . . , 200 3 (au milieu) ; dimension s = 4 pour N = 10 4 , . . . , 40 4 (en bas). Les échelles sont logarithmiques en base 2.

Figure 2 . 5 -

 25 Figure 2.5 -Variances des estimateurs de J en fonction du nombre de points N , pour M = 100 répétitions (tracés en pointillés) et M = 1 000 répétitions (tracés continus). Dimension s = 2 pour N = 10 2 , . . . , 400 2 (en haut) ; dimension s = 3 pour N = 10 3 , . . . , 200 3 (au milieu) ; dimension s = 4 pour N = 10 4 , . . . , 40 4 (en bas). Les variances des estimateurs MC sont les courbes supérieures, celles des estimateurs MCS sont les courbes inférieures ; les échelles sont logarithmiques en base 2.

s ≥ 2 .

 2 Soient s ≥ 2 un entier et N ∈ N * . Une méthode d'échantillonnage LHS est la suivante : On note N le nombre de points de l'échantillon et N = {1, 2, . . . , N } l'ensemble des indices qu'on utilisera. On note S N le groupe symétrique d'indice N : c'est l'ensemble des permutations de N .

Figure 2 . 6 -

 26 Figure 2.6 -Échantillon LHS de N = 6 2 points en dimension s = 2. Il y a un point unique dans chacun des rectangles de la forme I × [(ℓ -1)/N, ℓ/N ) (en haut) ou de la forme [(ℓ -1)/N, ℓ/N ) × I (en bas), pour 1 ≤ ℓ ≤ N .

Figure 2 . 7 -

 27 Figure 2.7 -Courbes des variances des estimateurs de J en fonction du nombre de points N pour M = 100 répétitions (tracés en pointillés) et M = 1 000 répétitions (tracés continus). Dimension s = 2 pour 20 2 ≤ N ≤ 400 2 (en haut) ; dimension s = 3 pour 10 3 ≤ N ≤ 100 3 (au milieu) ; dimension s = 4 pour 6 4 ≤ N ≤ 40 4 (en bas). Les variances des estimateurs MC sont les courbes supérieures, celles des estimateurs LHS sont les courbes inférieures ; les échelles sont logarithmiques en base 2.

1 .

 1 n ∈ N * un entier, N = n s , L = {ℓ = (ℓ 1 , . . . , ℓ s ) : ℓ 1 , . . . , ℓ s ∈ {1, . . . , n}} et M = {m = (m 1 , . . . , m s ) : m 1 , . . . , m s ∈ {1, . . . , n s-1 }};

Figure 3 . 2 -

 32 Figure 3.2 -Variance de l'estimateur SS de λ s (Q) en fonction du nombre de points N , pour M = 100 répétitions (tracés en pointillés) et M = 1 000 répétitions (tracés continus). Dimension s = 2 et 10 2 ≤ N ≤ 400 2 (en haut), dimension s = 3 et 10 3 ≤ N ≤ 200 3 (au milieu), dimension s = 4 et 10 4 ≤ N ≤ 40 4 (en bas). Échelles logarithmiques à base 2.

  (dimension s = 2), 3.4 (dimension s = 3) et 3.5 (dimension s = 4). L'estimation des variances se fait avec M = 1 000 répétitions. Nous avons indiqué dans les mêmes figures l'évolution des temps de calcul en fonction du nombre N de points de la simulation.

  .

2 .

 2 Pour terminer, nous reprenons le calcul de λ s (Q) par les méthodes MC, LHS, MCS, SS et une méthode quasi-Monte Carlo hasardisée par la technique d'Owen (MCO) introduite au chapitre 1. Dans ces essais numériques, le nombre N de points utilisés doit être une puissance de la base b, b étant le plus petit nombre premier supérieur ou égal à la dimension s. Le résultat du calcul de λ 2 (Q) par les méthodes MC, LHS, MCS, SS et MCO est tracé en fonction du nombre N de points sur la Figure 3.7 (en haut), en échelle logarithmique à base 10. Le nombre points de simulation est N = 2 m+6 , 1 ≤ m ≤ 13, pour la méthode MCO et N = (10m) 2 , 1 ≤ m ≤ 74, pour les autres méthodes. De nouveau, les quadratures MC et LHS donnent les résultats les moins bons, les méthodes MCS et MCO donnent des résultats proches et la méthode SS montre une convergence très rapide. Un inconvenient de l'approche MCO, est le temps de calcul. Sur la même Figure 3.7 (en bas), on compare les résultats de l'estimation de λ 3 (Q) par ces différentes méthodes en fonction du nombre de particules N = (10m) 3 , 1 ≤ m ≤ 9 pour les méthodes MC, LHS, MCS et SS et N = 3 (m+5) pour la simulation d'Owen avec 1 ≤ m ≤ 7.

Figure 3 . 3 -

 33 Figure 3.3 -Dimension s = 2. Variances (en haut) et temps de calcul (en bas) pour 1 000 estimations indépendantes de λ 2 (Q) en fonction du nombre de points N (N = , . . . , 400 2 ). Comparaison des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆).

Figure 3 . 4 -

 34 Figure 3.4 -Dimension s = 3. Variances (en haut) et temps de calcul (en bas) pour 1 000 estimations indépendantes de λ 3 (Q) en fonction du nombre de points N (N = 10 3 , . . . , 100 3 ). Comparaison des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆).

Figure 3 . 5 -

 35 Figure 3.5 -Dimension s = 4. Variances (en haut) et temps de calcul (en bas) pour 1 000 estimations indépendantes de λ 4 (Q) en fonction du nombre de points N (N = 6 4 , . . . , 40 4 ). Comparaison des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆).

Figure 3 . 6 -

 36 Figure 3.6 -Efficacités des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆) en fonction du nombre de points N en dimension s = 2 (en haut), s = 3 (au milieu) et s = 4 (en bas).

Figure 3 . 7 -

 37 Figure 3.7 -Valeurs des variances obtenues en répétant 1 000 fois de manière indépendante le calcul de λ 2 (Q) (en haut) et λ 3 (Q) (en bas) en fonction de N : comparaison des méthodes MC (+), LHS (△), MCS ( ), SS (⋆) et MCO (•).

.1. 1 )

 1 c(x, 0) = c 0 (x), x ∈ R, (4.1.2) où c est la concentration de matière en x à l'instant t et D = D(x, t) est le coefficient de diffusion. On suppose que la concentration initiale c 0 est une fonction intégrable et positive. Dans le cas où le domaine est borné, il faut ajouter des conditions aux limites. La matière est conservée au cours du temps et cela se traduit par l'égalité ∀t ≥ 0 R c(x, t)dx = R c 0 (x)dx.

Figure 4 . 1 -

 41 Figure 4.1 -Évolution de la concentration dans les intervalles [40, 60] (en haut) et [160, 180] (en bas). Comparaison de la solution exacte (ligne régulière) et de la simulation MC (ligne oscillante) ; coefficient de diffusion constant. .

Figure 4 . 2 -

 42 Figure 4.2 -Variances (en haut) et temps de calcul (en bas) pour 5 000 estimations indépendantes de l'approximation de J = a 0 c(x, T )dx en fonction de N , pour 10 2 ≤ N ≤ 200 2 . Comparaison des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆).

Figure 4 . 4 -

 44 Figure 4.4 -Valeurs des variances obtenues en répétant 1 000 fois de manière indépendante le calcul de J = a 0 c(x, T )dx en fonction de N . Comparaison des méthodes MC (+), LHS (△), MCS ( ), SS (⋆) et MCO (•).

  8u(x)), avec D 0 = 4 000 et u(x) := 1 2 (tanh(A(bx)) + 1), A = 0, 035 et b = 200.

Figure 4 . 5 -Figure 4 . 6 -

 4546 Figure 4.5 -Coefficient de diffusion constant. Évolution de la concentration dans les intervalles [40, 60] (en haut) et [100, 120] (en bas) par les méthodes MC (à gauche) et SS (à droite)

Figure 4 . 7 -Figure 4 . 8 -

 4748 Figure 4.7 -Coefficient de diffusion variable. Évolution de la concentration dans les intervalles [40, 60] (en haut) et [100, 120] (en bas) par les méthodes MC (à gauche) et SS (à droite)

Figure 5 . 1 -Figure 5 . 2 -

 5152 Figure 5.1 -Noyau de coagulation constant. Valeurs des variances et ajustements linéaires des estimations SS (⋆) et MC (•) des moments C 0 (T ) (en haut) et C 2 (T ) (en bas) en fonction de N (4 3 ≤ N ≤ 32 3 ), pour 1 000 répétitions indépendantes. Les échelles sont logarithmiques à base 2.

  décroît au cours du temps à cause du phénomène de coagulation, alors que la masse totale des particules C 1 (t) := +∞ 0 xc(x, t)dx reste constante. Cela se démontre comme pour l'équation discrète et exprime la conservation de la masse dans le système.

  .3.2) où K(x, y) := C 1 K(x, y) y , x, y > 0 est le noyau de coagulation modifié. On note f 0 (x) := xc 0 (x)/C 1 la densité de masse initiale. On a donc ∀t > 0, R + ) l'ensemble des fonctions définies sur R + , positives, mesurables et dont l'image est un sous-ensemble fini de R + . On considère la formulation faible de l'équation (5.3.2) suivante : pour toute fonction-test σ ∈ S(R + ) d dt +∞ 0f (x, t)σ(x)dx = , y)f (x, t)f (y, t)(σ(x + y)σ(x))dydx.(5.3.4)

) ∈ I 3 . 3 Γ

 33 Le problème d'approximation se transforme alors en un problème d'intégration sur I 3 puisque +∞ 0 g p+1 (x)σ(x) = I p+1 σ (u)du. (5.3.5)

  k sinon, où L 1 , . . . , L N sont des variables aléatoires indépendantes uniformément distribuées dans {1, 2, . . . , N } et U 1 , . . . , U N sont des variables aléatoires indépendantes uniformément distribuées dans [0, 1]. Les variables {L k , U k : 1 ≤ k ≤ N } sont supposées indépendantes. Un algorithme QMC a été proposé dans [123], où une borne d'erreur théorique est prouvée. Pour des entiers b ≥ 2 et m ≥ 1 on pose N := b m le nombre de particules. Soit Y = {y 0 , y 1 , . . .} ⊂ I 3 une suite à faible discrépance ; on suppose que Y est une suite-(t, 3) en base b pour t ≥ 0 (voir

  p ; f p ) ≤ e ctp D ⋆ N (X 0 ; f 0 ) + ∆t m-t)/2⌋ e ctp -1 c , où c K := 4V ( K) + 3 K ∞ et c := sup x>0 V ( K(x, .)) + sup y>0 V ( K(., y)) + 3 K ∞ .

K 0

 0 (x, y) = 1 et K 1 (x, y) = x + y, ∀x, y > 0 pour lesquels les solutions exactes sont connues (pour certaines données initiales).Le moment d'ordre µ ≥ 0 défini parC µ (t) := +∞ 0 x µ c(x, t)dx = +∞ 0 x µ-1 f (x, t)dxpeut être appoché numériquement, à l'instant t p par :

Figure 5 . 3 -

 53 Figure 5.3 -Noyau constant : variances de M = 1 000 estimations des moments C 0 (T ) (en haut) et C 2 (T ) (en bas), en fonction de N . Comparaison des approches MC (•) et SS (⋆) avec ajustements linéaires.

Figure 5 . 4 -Table 5 . 2 -

 5452 Figure 5.4 -Noyau linéaire : variances de M = 1 000 estimations des moments C 0 (T ) (en haut) et C 2 (T ) (en bas), en fonction de N . Comparaison des approches MC (•) et SS (⋆) avec ajustements linéaires.

  s , où N = m s est le nombre d'évaluations de la fonction f dans la formule. Pour diviser le majorant de l'erreur par deux, il faut multiplier le nombre de points d'intégration par 2

s/2 . Ce facteur croît exponentiellement avec la dimension s ; cela rend cette approche inutilisable en dimension élevée. Pour résoudre cette difficulté, on a développé les méthodes de Monte Carlo dont la vitesse de convergence est indépendante de la dimension du problème considéré

[START_REF] Hammersley | Monte Carlo Methods[END_REF][START_REF] Davis | Méthods of numerical integration[END_REF][START_REF] Fishman | Concepts, Algorithms, and Applications[END_REF][START_REF] Evans | Approximating Integrals via Monte Carlo and Deterministic Methods[END_REF][START_REF] Madras | Lectures on Monte Carlo Methods[END_REF]

. L'idée principale de la méthode de Monte Carlo est d'exprimer une intégrale comme l'espérance d'une variable aléatoire.

Soit π une loi de probabilité sur (R s , B(R s )). Soit f une fonction mesurable sur cet espace, telle que

  1. Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique Théorème 1.1 (Loi forte des grands nombre) Soit (Y n ) n≥1 une suite indépendante de variables aléatoires réelles équidistribuées. La suite (Y N ) N ≥1 converge presque sûrement si et seulement si les variables Y n sont intégrables. La limite est alors la valeur commune µ des espérances des Y n et la suite (Y N ) N ≥1 converge vers µ dans L 1 .

	Revenant à notre problème d'estimation, (f • X n ) n≥1 est une suite indépendante de variables aléatoires réelles équidistribuées intégrables, d'espérance I π (f ). D'après la loi forte des grands nombres, on a f • X N p.s.

Table 1 .

 1 1 -Niveau de confiance.

	a	1	1.96	2	2.5758	3
	2Φ(a)-1 0.6826 0.9500 0,9544 0.9900 0.9974
	On reprend notre problème d'intégration : π est une probabilité sur

  1. Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

	Comme on a	S N S N	p.s. → σ P

  Si X est un ensemble de N points de I s et si B ⊂ I s est mesurable, la discrépance locale de X pour B est définie par

	1. Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique
	1.3.1 Intégration numérique quasi-Monte Carlo
	On examine d'abord l'intégration de fonctions indicatrices de boréliens. L'erreur d'intégra-
	tion des fonctions indicatrices des rectangles permettra de majorer l'erreur d'intégration
	d'une classe de fonctions assez générale.
	Si X := {x 1 , x 2 , . . . , x N } est un ensemble de N points de I s et B est un borélien inclus dans I s , on définit
	A(B, X) := #{n : x n ∈ B}, où #E note le cardinal de l'ensemble E.
	Définition 1.1

Table 1 .

 1 2 -Variance des estimateurs de J s . Nombre de points N 5 8 = 390 625 11 5 = 161 051 13 5 = 371 293 19 4 = 130 321 Matousěk 1, 37 × 10 -22 1, 37 × 10 -18 1, 75 × 10 -18 2, 60 × 10 -19 I-Binomial 2, 45 × 10 -14 2, 81 × 10 -18 5, 77 × 10 -17 6, 25 × 10 -21 Owen 2, 07 × 10 -14 6, 13 × 10 -21 6, 37 × 10 -18 1, 80 × 10 -20

	Dimension s	4	8	12	18
	Méthode
	Monte Carlo	1, 56 × 10 -7	1, 97 × 10 -7	9, 31 × 10 -7	1, 23 × 10 -7

1. Méthodes de Monte Carlo et quasi-Monte Carlo pour l'intégration numérique

Table 1 .

 1 3 -Variance des estimateurs de I pour N = 5 8 = 390 625.

	Méthode	Variance
	Monte Carlo 1, 22 × 10 -7 Matousěk 8, 23 × 10 -16

I-Binomial 3, 70 × 10 -16 Owen 6, 15 × 10 -16

  {U 1 , . . . , U N } est une famille indépendante de variables aléatoires de loi uniforme sur I

s . L'expression (2.2.1) donne

Table 2 .

 2 1 -Ordres α des erreurs dans le calcul de J par les méthodes MC, MCS et QMC.

	Dimension s MC MCS QMC
	2	0,44 0,78	0,86
	3	0,57 0,57	0,80
	4	0,63 0,69	0,68

Table 2 .

 2 2 -Ordres β des variances des estimateurs MC et MCS de J .

	Dimension MC MCS Borne théorique
	2	0,99 1,49	1,50
	3	1,00 1,33	1,33
	4	1,00 1,25	1,25

Table 2 .

 2 3 -Ordres γ des variances des estimateurs MC et LHS de J .

	Dimension MC LHS
	2	0,99 1,00
	3	1,00 1,00
	4	1,01 1,00

Table 3 .

 3 1 -Ordre de convergence des variances des estimateurs MC et SS de J Dans cette section, nous testons numériquement la méthode de Monte Carlo stratifiée Sudoku (SS) étudiée dans ce chapitre. Nous reprenons le calcul du volume du secteur de la boule unité déjà considéré au chapitre 2. Soit

	Dimension MC SS Borne théorique
	2	0,99 1,50	1,50
	3	1,00 1,32	1,33
	4	1,00 1,25	1,25
	3.4 Essais numériques

Table 3 .

 3 2 -Ordre de convergence des variances des estimateurs MC, LHS, MCS et SS N = 6 4 , 8 4 , . . . , 40 4 = 2 560 000 points. Nous représentons les variances des estimateurs des différentes méthodes pour le calcul de λ s (Q) dans les Figures 3.3

	dimension MC LHS MCS SS
	2	0,99 1,00 1,48 1,50
	3	1,00 1,00 1,34 1,33
	4	1,01 1,00 1,26 1,24
	• pour s = 4,

Table 3 .

 3 3 -Ordre de convergence des variances des estimateurs MC, LHS, MCS, SS et MCO Dans ce chapitre, nous avons analysé une technique de stratification appelée Sudoku qui utilise des points aléatoires ayant une meilleure répartition dans l'hypercube unité que ceux d'une méthode de Monte Carlo classique. Elle possède les propriétés de stratification de la méthode stratifiée simple et aussi celles de la méthode LHS du chapitre précédent. Dans le cadre de l'évaluation des intégrales de fonctions indicatrices, l'analyse de la variance a confirmé qu'il s'agit bien d'une technique de réduction de la variance. Des essais numériques menés ont montré que cette approche donne des résultats très satisfaisants en ce qui concerne la précision et le temps de calcul.Dans les chapitres qui suivent, nous allons montrer que cette technique de stratification peut être utilisée dans des simulations Monte Carlo d'équations différentielles ou intégrales simples.

	dimension MC LHS MCS SS MCO
	2	0,99 1,00 1,49 1,50 1,43
	3	1,00 1,00 1,33 1,32 1,32
	3.5 Conclusion

  des bijections aléatoires de {1, . . . , n} dans {1, . . . , n} et ξ 1 ℓ , ξ 2 ℓ sont des variables aléatoires uniformément distribuées sur I. Toutes ces variables sont supposées indépendantes. Cet ensemble vérifie les deux propriétés d'uniformité suivantes :

Table 4 .

 4 1 -Ordre β de la variance du calcul de J par les méthodes MC, LHS, MCS et SS

	MC LHS MCS SS
	1,00 1,00 1,44 1,43
	et avec un coefficient de diffusion constant D = 1. On se propose de calculer l'intégrale

  Nous représentons les courbes d'évolution de la concentration dans différents intervalles. Nous considérons trois émissions successives de particules. La période des instants d'émission est égale à 10. Le temps final est T = 30 et le nombre de pas de temps est P = 187 500 (soit 62 500 pas de temps par période), donc ∆t := T /P = 1, 6 10 -4 . Le nombre de particules libérées durant chaque période d'émission est N = 212 . Afin de diminuer la fluctuations des courbes, nous faisons des moyennes sur des intervalles temporels d'amplitude 0, 16. Nous comparons les résultats des méthodes MC et Monte Carlo stratifié Sudoku (SS). Pour les simulations SS, les groupes de N particules évoluent séparément. L'évolution de la concentration des ions est calculée dans les intervalles :[START_REF] Haber | A modified Monte-Carlo quadrature[END_REF][START_REF] L'ecuyer | A randomized quasi-Monte Carlo simulation method for Markov chains[END_REF],[START_REF] Niederreiter | Constructions of (t, m, s)-nets and (t, s)-sequences[END_REF][START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], [160, 180] et [180, 200], dans un milieu où le coefficient de diffusion est D = 4 000. Les résultats des simulations sont représentés dans les Figures 4.5et 4.6.Nous constatons que la méthode SS atténue nettement les oscillations des résultats en comparaison avec les calculs MC. Les précisons atteintes sont presque inaccessibles à une méthode MC standard avec le même nombre de particules de simulation.

  2 et σ 3 sont des bijections aléatoires de {1, . . . , n} 2 dans {1, . . . , n 2 } et ξ 1 ℓ , ξ 2 ℓ , ξ 3 ℓ sont des variables aléatoires uniformément distribuées sur I. Toutes ces variables aléatoires sont supposées indépendantes. Cet ensemble vérifie les deux propriétés d'uniformité suivantes :

  2 des points permettent de sélectionner les particules à l'instant t p : la particule d'indice k(W ℓ,1 ) peut coaguler avec la particule d'indice k(W ℓ,2 ) avec une probabilité de coagulation

Table 5 .

 5 1 -Équation de coagulation discrète : estimation des ordres β des variances comme fonctions de N pour les calculs de C 0 (T ) et C 2 (T ).

	C 0 (T ) Noyau	MC	SS	C 2 (T ) Noyau	MC	SS
	constant 1, 01 1, 22	constant 0, 99 1, 28
	linéaire 1, 06 1, 15	linéaire 0, 78 1, 14

  3.1) a été calculée par Ramabhadran, Peterson et Seinfeld dans[START_REF] Ramabhadran | Dynamics of aerosol coagulation and condensation[END_REF] ; elle est donnée par : Les valeurs exactes des moments d'ordre 0 et 2 sont les suivantes :C 0 (t) = 2 2 + t et C 2 (t) = t + 2, t ≥ 0.Pour estimer les variances des deux approches, nous effectuons M estimations indépendantes des moments d'ordre 0 et 2 à l'instant T = 1, 0 : C 0 (T ) et C 2 (T ) par chacune des méthodes, puis nous calculons la variance empirique des échantillons.Les différences entre les estimations pour M = 100, . . . , 1 000 sont minimes ; nous nous limitons au cas M = 1 000. La Figure5.3 montre les variances en fonction du nombre de particules N , qui est ici de la forme N = (4m) 3 , avec 1 ≤ m ≤ 12. Le nombre de pas de temps est P = 4 000. La variance du schéma SS est inférieure à celle du schéma MC pour un même nombre de particules. À titre d'exemple, dans l'approximation de C 2 (T ), la variance obtenue avec 21 952 points stratifiés Sudoku est comparable à celle obtenue avec 46 656 points aléatoires. L'amélioration est moins visible dans l'estimation de C 0 (T ) ; il s'agit ici de calculer la somme

		C 0,P :=	1 N	N k=1	1 x P k
	c(x, t) =	4 (2 + t) 2 exp -	2x 2 + t	, x ≥ 0, t ≥ 0.

Méthode de Monte Carlo stratifiée Sudoku pour l'intégration numérique Figure 3.1 -Un exemple de 4 2 points stratifiés Sudoku en dimension s = 2.

Méthode de Monte Carlo stratifiée Sudoku pour l'intégration numérique

Figure 4.3 -Efficacités des méthodes MC (+), LHS (△), MCS ( ) et SS (⋆) en fonction du nombre de points N .

Cas général d'un coefficient variable

Dans ce cas on fait une approximation linéaire de D(x) en utilisant un développement de Taylor :

Pendant un intervalle de temps ∆t, une particule en x 0 fera un déplacement d'amplitude ∆x + ǫ, où ∆x = 2D(x 0 )∆tZ, Z de loi N (0, 1), (4.4.7)

où l'on a noté

Algorithme de simulation Monte Carlo stratifié en milieu hétérogène

Dans ce paragraphe, nous développons un algorithme pour la simulation de la diffusion dans un milieu hétérogène. Il reprend les étapes principales de celui présenté à la section 4.3 en introduisant le terme correcteur qu'il convient d'ajouter au pas de l'approche naïve. Une approche QMC a été étudiée dans [START_REF] El Haddad | Méthodes quasi-Monte Carlo de simulation des chaînes de Markov[END_REF][START_REF] El Haddad | Quasi-Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium[END_REF]. Soit n un entier naturel et N = n 2 ; soit ∆t un pas de temps. On approche la solution de l'équation de diffusion : On pose

Ces changements conduisent à approcher la densité de masse f (i, t) au lieu de la densité du nombre des particules c(i, t), c'est-à-dire à résoudre l'équation :

Pour tout sous-ensemble E de N * := {1, 2, 3, . . .}, on note s E la suite définie par :

En multipliant l'équation (5.2.3) par s E (i) et en sommant sur i ∈ N * , on obtient

pour tout E ⊂ N * . On note f 0 la densité de masse initiale. L'égalité

Un schéma de Monte Carlo de résolution de l'équation (5.2.3) a été proposé dans [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF]. La démonstration de sa convergence est également fournie dans [START_REF] Tarhini | Analyse numérique des méthodes quasi-Monte Carlo appliquées aux modèles d'agglomération[END_REF]. Plus tard, Lécot et Wagner en ont donné une version quasi-Monte Carlo dans [START_REF] Lécot | A quasi-Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF]. Nous rappelons plus loin un résultat de convergence pour cette approche.

Description du schéma

Nous proposons ici un schéma stratifié Sudoku pour la simulation numérique de (5.2.3) ; il comporte plusieurs étapes que nous abordons successivement.

Initialisation des particules

On considère un entier n ≥ 1 et on note N := n 3 . La première étape de l'algorithme consiste à choisir un ensemble de N particules J 0 := {j 0 1 , . . . , j 0 N } qui échantillonnent la distribution initiale f 0 , c'est-à-dire telles que

où δ j est la mesure de Dirac au point j définie par :

Ici, le nombre de particules de simulation de même valeur j 0 k représente non pas le nombre de particules de masse j 0 k , mais la masse totale des particules de masse j 0 k dans le système. Dans le cas d'une distribution initiale ponctuelle f 0 (1) = 1, f 0 (2) = f 0 (3) = . . . = 0, l'intitialisation se fait en prenant :

Évolution des tailles des particules

On suppose dans la suite que le noyau K est borné et on note

K(i, j).

L'algorithme proposé reste valable pour les noyaux de coagulation non bornés comme le montrent les essais numériques de la section 5.2.2 : on adapte le pas de temps. On choisit un pas de temps fixe ∆t vérifiant ∆t K ∞ < 1. Pour tout entier p, on pose t p := p∆t et f p (i) := f (i, t p ). Supposons qu'à l'instant t p soit déjà construit un ensemble J p = {j p 1 , . . . , j p N } tel que

(5.2.5) L'approximation de la solution à l'instant t p+1 se fait en deux étapes :