Simulation :

De la loi uniforme aux équations différentielles stochastiques

Fabien Campillo, Frédéric Cérou, David Miglior

SYSDYS^{∗}

Table des matières

1 Avertissement 2

2 Introduction 2 2.1 Génération de nombres pseudo–aléatoires . . . . . . . . . . . . . . . . . . . . 2 2.2 Méthodes numériques probabilistes . . . . . . . . . . . . . . . . . . . . . . . 2

3 Variables aléatoires 2 3.1 Loi uniforme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2 Loi normale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.3 Algorithmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

4 Processus stochastiques 5 4.1 Marche aléatoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4.2 Mouvement brownien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

4.2.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 4.2.2 Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4.3 Processus d’Ornstein-Uhlenbeck . . . . . . . . . . . . . . . . . . . . . . . . . 6

5 Quelques liens 7

6 Crédits 7

∗http://www.inria.fr/mefisto/sysdys-fra.html

1 Avertissement

Les noms précédés du signe “▶” sont des pointeurs vers la base de données sur l’histoire des mathématiques, MacTutor^{1}, de l’ Université de St Andrews^{2}.

2 Introduction

Le but de ce “petit cours” est de tenter de sensibiliser les communautés scientifiques et les étudiants à certaines méthodes numériques dites “probabilistes”.

Nous avons tenté, dans la mesure du possible, d’illustrer notre propos à l’aide de simulations interactives fondées sur Java (Sun Microsystem).

Il faut distinguer deux problèmes : celui des algorithmes de génération de nombres pseudo aléatoires et celui des méthodes numériques, que nous qualifierons de “probabilistes”, faisant appel à ces algorithmes.

2.1 Génération de nombres pseudo–aléatoires

Nous n’aborderons pas ici le problème (difficile) de la “nature” aléatoire d’une suite de nombres. L’article de vulgarisation Chaitin [2] propose un exposé intéressant et relativement abordable du sujet.

2.2 Méthodes numériques probabilistes

Nous citerons quelques exemples où des méthodes de “simulation” de variables aléatoires où de processus stochastiques sont présentées/utilisées :

– méthodes numériques pour les processus stochastiques : Bouleau–Lépingle [1]

– méthodes de Monte–Carlo : Lapeyre–Pardoux–Sentis [5],

– simulation d’équations différentielles stochastiques : Kloeden–Platen [4],

– simulation géostatistique : Deutsch–Journal [3].

3 Variables aléatoires

On considère un espace de probabilité (Ω, F, P) et X une variable aléatoire (à valeurs dans ℝ, ℝ^{2}...). La simulation de pseudo–réalisations de cette variable est un problème complexe qui trouve de nombreuses applications.

En effet, beaucoup de méthodes numériques utilises de tels outils. On pensera en premier aux méthodes de Monte–Carlo. Mais d’autres algorithmes sont apparus : recuit simulé et algorithmes génétiques en optimisation combinatoire ; algorithmes stochastiques comme le gradient stochastique.

1. http://www-groups.dcs.st-and.ac.uk:80/˜ history/index.html 2. http://www.st-and.ac.uk/

3.1 Loi uniforme

Soit X une variable aléatoire de loi uniforme U[0, 1].

loi(X) = 1_{[0,1]}(x) dx .

Applet : Densité empirique d’un échantillon d’une loi uniforme U[0, 1]

Les projectiles sont lancés de façon à retomber sur le sol en suivant un loi uniforme sur un intervalle. Les histogrammes représentent la densité empirique du point d’impact : sur chaque sous-intervalle la hauteur est égale au nombre d’impacts sur celui-ci divisé par le nombre total de jets. On vérifie que sur un grand nombre de jets, l’histogramme converge vers une constante. Le générateur utilisé, comme dans toute la suite, est le générateur fourni par Java. La barre de contrôle règle le nombre de sous-intervalles utilisés pour les histogrammes.

Applet : Pile ou face

Cette applet simule le jet d’une pièce. Les deux barres à gauche et à droite de la pièce indiquent respectivement la proportion de “piles” et de “faces”. On vérifie que sur un grand nombre de tirages les deux s’équilibrent à 1/2.

Applet : loi uniforme sur le carré unité

L’applet ci-dessus montre des tirages pseudo-aléatoire d’une variable uniforme dans le carré unité, et la façon dont elle “remplit” uniformément le carré. La barre de contrôl selectionne le nombre de tirages entre deux images.

3.2 Loi normale

Soit X une variable aléatoire X de loi normale N[0, 1] également appelée loi gaussienne (cf. ▶ Johann Carl Friedrich Gauss^{3}).

Applet : La planche de Galton La “planche de ▶ Francis Galton^{4}” illustre la propriété de convergence en loi de la loi binomiale B(1/2, n) vers la loi normale N(0, 1). Il s’agit d’un résultat du théorème central limite (désigné ainsi par le mathématicien hongrois, ▶ George Pólya^{5}).

Planche de Galton : Illustre la converge de la binomiale B(1/2, n) (i.e. n tirages indépendants d’une pièce de monnaie non biaisée) vers la loi normale centrée réduite. Le paramètre n, ici limité à 16, se règle à l’aide du curseur.

3. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Gauss.html 4. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Galton.html 5. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Polya.html

3.3 Algorithmes

3.3.1 Deux méthodes pour calculer π

Les deux applets ci-dessous montrent comment calculer une quantité, ici le nombre π, en utilisant deux méthodes de Monte-Carlo différentes.

Applet : L’Aiguille de Buffon Il s’agit d’un problème de probabilités géométriques posé (et résolu) par ▶ Buffon^{6}. Sur un parquet composé de planches parallèles, toutes de même largeur, on jette des aiguilles au hasard, et on cherche la probabilité pour qu’une aiguille tombe à cheval sur deux planches. Soit c la largeur des planches, et d la longueur de l’aiguille. Pour simplifier, on suppose que d ≤ c. On note r la distance entre le centre de l’aiguille et la jointure de deux planches la plus proche, et θ l’angle de l’aiguille avec la jointure. Il est alors facile de vérifier que l’aiguille intersecte une jointure si et seulement si

\frac{d}{2} sin θ > r.

Si on suppose, ce qui paraît raisonnable, que les jets sont faits de sorte que les lois de r et θ sont uniformes sur respectivement [0, c /2] et [0, π/2], la probabilité de tomber sur une rainure du parquet est donc de

\frac{∫0^{\frac{π}{2}} \frac{d}{2} sin θ dθ}{\frac{π}{2} × \frac{c}{2}} = \frac{2d}{πc}.

Ainsi, sur un grand nombre d’expériences, en estimant de façon empirique cette probabilité, on en déduit une estimation ▶ du nombre^{7} π.

Applet : Calcul de π par méthode de Monte–Carlo Voici une autre façon relativement simple, bien que certainement pas la plus efficace, de calculer ▶ le nombre^{8} π.

Si on tire selon une loi uniforme un point sur le carré, il est facile de voir que la probabilité de tomber dans le cercle est de π /4. Donc sur un grand nombre de réalisations, en multipliant par 4 la fréquence empirique des points intérieurs au cercle, on doit retrouver approximativement π. Notons que cette méthode peut s’appliquer au calcul d’une aire quelconque, pourvu que l’on sache décider facilement si un point donné est à l’intérieur ou à l’extérieur de la surface considérée.

3.3.2 Recuit simulé

Il est clair que pour minimiser une fonction F non convexe un algorithme de type gradient peut donner de très mauvais résultats : si la fonction comporte beaucoup de minima locaux, on est à peu près certain de rester bloqué sur l’un d’eux. La stratégie des algorithmes de recuit consiste, en effectuant une exploration aléatoire de l’espace d’état, à favoriser les descentes, mais sans interdire tout à fait les remontées. Plus précisément, on se donne une chaine de Markov sur l’espace d’état, et on accepte ou on refuse une transition avec une probabilité 1 si F décroît, et une probabilité égale à

6. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Buffon.html

7. http://www-groups.dcs.st-and.ac.uk/˜ history/HistTopics/Pi_through_the_ages.html

8. http://www.cecm.sfu.ca/pi/pi.html

exp [−\frac{ΔF}{T} ]

si F croît. Ici T est un paramètre, appelé température par analogie avec la mécanique statistique. Il est clair que plus la température est grande, plus seront facilitées les transitions ascendantes. En revanche, à la limite T = 0, on obtient un algorithme de descente. Tout au long de l’algorithme, on va donc faire décroître T, ni trop vite pour ne pas rester bloqué autour d’un minimum local, ni trop lentement si on veut avoir un résultat en un temps raisonnable.

Applet : Recuit simulé Crédit : Cette Applet est due à Matthew Chalwell^{9} de l’University Of Western Australia^{10}

Cette Applet résout un problème d’optimisation (non combinatoire) mais fait apparaître tous les paramètres de stratégie de recuit.

Applet : Problème du voyageur de commerce Etant donné un ensemble de villes (des points dans le plan), il s’agit de trouver le plus court chemin fermé passant une fois et une seule par chacun des points. La difficulté du problème vient de l’explosion combinatoire des chemins à explorer lorque l’on accroît le nombre de villes à visiter.

4 Processus stochastiques

4.1 Marche aléatoire

Le terme de “marche aléatoire” désigne le plus souvent un processus stochastique en temps discret à valeur dans ℤ^{d}. Un cas simple est celui du processus X_{n} à valeur dans ℤ^{2} qui, lorsqu’il se trouve à l’instant n au point (i, j), décide de passer à l’un de ses plus proche voisin :

(i + 1, j) , (i − 1, j) , (i, j + 1) , (i + 1, j − 1) ,

avec une probabilité 1/4 (indépendemment des positions passées). Naturellement il s’agit d’un processus de Markov [▶ Andrei Andreyevich Markov^{11}] il est de plus à accroissement indépendant.

9. http://docker.library.uwa.edu.au/matt/ 10. http://docker.library.uwa.edu.au 11. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Markov.html

On peut imaginer des marches aléatoires sur des espaces et/ou combiné à un système de voisinage plus complexes.

La simulation d’un tel processus est triviale.

Applet : marche au hasard

Le stylet trace une marche aléatoire simple. Après avoir lancé la première marche par “start”, on peut en lancer d’autres en cliquant sur la fenètre graphique. La barre de contrôle règle le nombre de pas à effectuer.

4.2 Mouvement brownien

4.2.1 Définition

Le mouvement brownien ou procesus de Wiener décrit un phénomène observé en 1827 par le botaniste écossais Robert Brown. Il observa que des grains de pollens en suspension dans l’eau suivaient un mouvement rapide et désordonné.

En 1905, Albert Einstien [▶ Albert Einstein^{12}] developpa cette théorie par une approche de type mécanique statistique.

À partir de 1920, Nobert Wiener [▶ Norbert Wiener^{13}] proposa une définition mathématique de ce phénomène.

– B_{0} = 0

– pour tout 0 ≤ s ≤ t, B_{t} − B_{s} est une variable gaussienne N(0, t − s) indépendant de la tribu engendrée par {B_{u}; u ≤ s}.

4.2.2 Approximation

La simulation la plus simple du mouvement brownien consiste à se donner un pas de discrétrisation en temps δ > 0 fixe et à poser : t_{n} = n δ. Le mouvement brownien standard unidimensionnel B_{t} est donné aux instants t_{n} par :

B_{t_{n}} = ∑_{k<n} ξ_{k}

avec B̄_{0} = 0 et où

ξ_{k} = B_{t_{k+1}} − B_{t_{k}} .

Par définition du mouvement brownien, ξ_{k} est une suite i.i.d. de variables gaussienne N(0, δ).

12. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Einstein.html 13. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Wiener_Norbert.html

Applet : Mouvement brownien dans le carré unité

Clicker sur “start” ou sur le carré unité pour lancer une simulation d’une trajectoire d’un mouvement brownien standard. On peut ensuite clicker plusieurs fois sur le carré et lancer ainsi des trajectoires différentes. n est le nombre de pas de l’approximation du mouvement brownien.

Applet : Mouvement brownien, aspect fractal Crédit : cette Applet est due à Charle (Charlie) Geyer^{14} de l’Université du Minnesota^{15} (cf. page originale^{16}).

Cette Applet présente un mouvement brownien 1d. Le mouvement brownien est auto–similaire en loi, donc toutes les images zoomées ont le même “aspect”. Plus précisément, si

{B_{t}, t ≥ 0}

est un mouvement brownien, alors pour tout c > 0

\{\frac{1}{c}B_{c^{2}t}, t ≥ 0\}

est aussi un mouvement brownien.

Clicker sur la fenêtre pour lancer le zoom. Clickez à nouveau pour arrêter.

4.3 Processus d’Ornstein-Uhlenbeck

▶ George Eugene Uhlenbeck^{17}

Applet : processus d’Ornstein-Uhlenbeck 2D

5 Quelques liens

– L’ Institut de Statistiques et des Sciences de la Décision^{18} de l’Université de Duke propose un site d’applets Java^{19} dans le domaine des probabilités et des statistiques.

– Un des sites les plus importants proposant des applications Java est celui de Gamelan^{20}.

14. http://www.stat.umn.edu/˜ charlie/ 15. http://www.stat.umn.edu/ 16. http://www.stat.umn.edu/˜ charlie/Stoch/brown.html 17. http://www-groups.dcs.st-and.ac.uk:80/˜ history/Mathematicians/Uhlenbeck.html 18. http://www.stat.duke.edu/ 19. http://www.stat.duke.edu/sites/java.html 20. http://www.gamelan.com

6 Crédits

– David Miglior (Java)

– Frédéric Cérou (Conception, Java, Mathématiques)

– Fabien Campillo (Conception, Mathématiques)

Références

[1] N. BOULEAU and D. LÉPINGLE. Numerical Methods for Stochastic Processes. John Wiley & Sons, 1994.

[2] G. CHAITIN. Les suites aléatoires. In Le Hasard, Pour la Science, Dossier, pages 68–73, avril 1996.

[3] C.V. DEUTSCH and A.G. JOURNEL. GSLib : Geaostatistical Software Library and User’s Guide. Oxford University Press, 1992.

[4] P.E. KLOEDEN and E. PLATEN. Numerical Solution of Stochastic Differential Equations, volume 23 of Applications of Mathematics. Springer Verlag, New York, 1992.

[5] B. LAPEYRE, É. PARDOUX, and R. SENTIS. Méthodes de Monte–Carlo. To appear.