La simulation de Monte-Carlo

Le principe de la simulation, au sens commun du terme, est d’utiliser un modèle, c’est-à-dire une représentation abstraite d’un système ou d’un problème, et d’étudier l’évolution de ce modèle sans faire fonctionner le système réel. C’est par exemple le cas lorsque vous allez à la banque pour un prêt : on vous dit combien vous payeriez chaque mois si vous acceptiez l’offre, sans pour autant avoir à le faire nécessairement. Modéliser consiste donc à poser un problème sur papier ou sur ordinateur afin de l’étudier. C’est aussi ce qui est fait pour les prévisions météorologiques, l’étude d’un phénomène physique, la construction et la mise au point de nouveaux avions, etc. Pour le cas d’un prêt bancaire, il s’agit de calculs relativement élémentaires, mais d’autres exemples nécessitent d’utiliser des méthodes particulières.

La simulation de Monte-Carlo est une méthode d’estimation d’une quantité numérique qui utilise des nombres aléatoires. Stanisław Ulam et John von Neumann l’appelèrent ainsi, en référence aux jeux de hasard dans les casinos, au cours du projet Manhattan qui produisit la première bombe atomique pendant la Seconde Guerre mondiale. La simulation d’un prêt bancaire par exemple n’est pas une simulation de Monte-Carlo car il s’agit de calculs exacts à partir du nombre de mensualités et du taux d’intérêt ; aucun phénomène aléatoire n’intervient dans les calculs.

La simulation de Monte-Carlo présente le double avantage d’être simple d’utilisation et de pouvoir être appliquée à un très large éventail de problèmes. Elle est utilisée en finance, pour déterminer quand lever une option sur un bien financier ; en assurance, pour évaluer le montant d’une prime ; en biologie, pour étudier les dynamiques intra et intercellulaires ; en physique nucléaire, pour connaître la probabilité qu’une particule traverse un écran ; en télécommunications, pour déterminer la qualité de service ; ou de façon générale pour déterminer la fiabilité d’un système, sa disponibilité ou son temps moyen d’atteinte de la défaillance.

Pour cela, il faut cependant poser le problème, le modéliser de sorte que la quantité à rechercher s’exprime comme l’espérance d’une variable aléatoire X, notée E(X). Une variable aléatoire est le résultat d’une expérience soumise au hasard ; son espérance est schématiquement ce que l’on s’attend à trouver en moyenne si l’on répète l’expérience un grand nombre de fois. Trouver cette modélisation est très probablement la tâche la plus complexe. Dans certains cas, la quantité à calculer est naturellement une espérance, par exemple la moyenne du nombre de clients dans une file d’attente quand les temps d’arrivée et de traitement sont aléatoires. Mais la méthode peut également être utilisée pour estimer une quantité purement déterministe, par exemple une surface ou une intégrale, en construisant artificiellement une variable aléatoire pour se ramener au calcul d’une moyenne.

Estimation du nombre π

Prenons l’exemple de l’estimation du nombre π, qui est la surface délimitée par un cercle de rayon 1, et vaut approximativement 3,14159. Comment la déterminer par simulation de Monte-Carlo ? Il faut donc représenter cette valeur comme la moyenne d’une variable aléatoire. Une méthode parmi d’autres peut être de dessiner un cercle de rayon 1 et de l’englober dans le carré de côté 2, donc d’aire 4, qui touche le cercle en 4 points, comme sur la figure. Dans le carré, choisissez un point aléatoirement et uniformément, ce qui signifie que chaque point a la même chance d’être choisi. Alors la probabilité qu’il se trouve à l’intérieur du cercle est la surface du cercle (π) divisée par la surface du carré (4), et vaut donc (π/4). Cette probabilité est l’espérance de la variable aléatoire X qui vaut 1 si le point choisi aléatoirement dans le carré est à l’intérieur du cercle et 0 sinon. Pour estimer l’espérance de X,  on considère un grand nombre de tirages indépendants uniformes de points dans le carré, qui vont se répartir uniformément dans le carré. Quand le nombre de points va tendre vers l’infini, la proportion de points dans le cercle va converger vers l’espérance de X, ici π/4.  Le nombre  π peut donc être estimé par 4 fois la proportion de points dans le cercle.

Il suffit de compter au fur et à mesure le nombre de points qui tombent à l’intérieur du cercle. Un algorithme basique pour estimer π est alors le suivant :

1. Demander le nombre n de points à générer ;
2. Initialiser un compteur à 0 ;
3. Pour i de 1 à n
   
   • générer un point (x,y) uniformément dans le carré (pour cela il suffit de générer x uniforme dans [-1,1] et y uniforme dans [-1,1])
   • Si le point est dans le cercle (c’est-à-dire si la distance par rapport à l’origine est plus petite que le rayon 1), ajouter 1 au compteur.
4. L’estimateur de π/4 est la valeur du compteur divisée par le nombre n.

L’animation suivante illustre l’estimation de π selon cette méthode.

Toute surface inconnue (ou volume si l’on veut passer en dimension 3) peut être calculée de la même manière en l’englobant dans un domaine plus large où il est facile de générer des points uniformément : on peut, par exemple, générer indépendamment les deux coordonnées d’un rectangle. La surface inconnue est alors la proportion de points multipliée par la surface du rectangle. L’idée utilisée ici est donc d’exprimer une surface comme l’espérance d’une variable aléatoire X qui vaut 1 si un point tiré aléatoirement dans le rectangle se trouve dans la surface d’intérêt et 0 sinon.

Notons que tout cet exposé sous-entend que nous disposons d’une source d’aléatoire selon des lois de probabilité désirées, un élément clé pour la simulation de Monte-Carlo. Ceci fera l’objet d’un prochain article.

La simulation de Monte-Carlo a connu une expansion rapide avec l’avènement des ordinateurs, sur lesquels il était facile de générer des variables et de multiplier les expériences. Mais la simulation de Monte-Carlo remonte à bien avant cela, même à l’Antiquité. Un autre exemple d’estimation de π  est via l’aiguille de Buffon en 1777. Il s’agissait de lancer des aiguilles sur les lattes de parquet et de compter la proportion touchant les séparations entre lattes. Selon la position du centre des aiguilles et l’angle par rapport aux lattes, on peut retrouver π via une formule  mathématique en fonction de la probabilité qu’une aiguille touche les lattes, estimée grâce à la proportion sur un échantillon.

Tout problème à résoudre par la méthode de Monte-Carlo doit être modélisé par le calcul de la moyenne d’une variable aléatoire X. En général, un algorithme de simulation de Monte-Carlo pour estimer la moyenne d’une variable aléatoire X,  est donc le suivant :

1. Demander le nombre n de points à générer ;
2. Initialiser un compteur sum à 0 ;
3. Pour i de 1 à n
   •     générer une copie X_i de X
   •     ajouter la valeur X_i à sum
4. L’estimateur est  sum/n, la moyenne des valeurs générées.

Estimation de l’erreur

Mais obtenir une estimation n’est pas tout. En effectuant deux expériences similaires indépendantes l’une de l’autre, nous devrions obtenir deux résultats différents — puisque des nombres aléatoires différents ont été utilisés pour la génération des variables. Se pose alors la question de l’erreur d’estimation commise. La simulation de Monte-Carlo permet en fait également d’estimer l’erreur de manière particulièrement simple, à l’aide d’un théorème statistique qui peut paraître assez magique, le théorème central limite.

Ce théorème dit schématiquement que quelle que soit la loi de probabilité d’un phénomène, si l’on répète un grand nombre de fois l’expérience d’observation de ce phénomène, la distribution des valeurs observées va s’approcher d’une loi dite gaussienne ou normale.

Cette convergence est très souvent illustrée par la planche de Galton : lancez des billes le long d’une planche où sont plantés des clous de sorte que la bille passe à chaque niveau soit à gauche, soit à droite du clou ; si vous lancez de nombreuses billes, la distribution observée en bas de la planche sera approximativement celle d’une loi gaussienne.

L’animation suivante simule une planche de Galton.

La loi gaussienne limite est caractérisée par deux paramètres uniquement, de manière remarquable, pour toute variable aléatoire X simulée. Ces paramètres sont l’espérance de X, qui est justement ce que l’on cherche à estimer, et  sa variance, qui détermine la dispersion par rapport à la moyenne. La variance est plus précisément  l’espérance d’une nouvelle variable aléatoire Y, qui vaut le carré de la différence entre X et l’espérance de X. La variance est souvent inconnue, mais peut cependant être facilement estimée en même temps que la moyenne dans l’algorithme que nous avons décrit.

Avec cette estimation, on peut obtenir une idée de l’erreur produite par la méthode de Monte-Carlo. Du fait de la nature statistique de l’estimation, on ne peut pas donner une réponse du type « mon estimation est x et je sais que mon erreur est assurément inférieure à e de sorte que je peux garantir une solution dans l’intervalle (x-e, x + e) ». Par contre, une réponse sera : « avec une probabilité souhaitée, la valeur cherchée est dans l’intervalle   (x-e, x + e) », alors appelé intervalle de confiance. On ne peut garantir d’avoir la solution dans l’intervalle, mais on peut l’avoir avec une probabilité de 90% , 99%, etc., cette probabilité de confiance étant spécifiée par l’utilisateur. Grâce au théorème central limite, on calcule le nombre e, qui est la demi-largeur de l’intervalle de confiance. La probabilité de confiance est associée à un nombre c que l’on sait calculer à partir de la loi gaussienne et que l’on peut trouver dans des tables : par exemple, si l’on veut un intervalle de confiance à 95%,  c=1,96. Pour un intervalle de confiance à 99%, c=2,58. On remarque que c augmente avec la probabilité de confiance. Le nombre e vaut \(c\sqrt{\dfrac{Var}{n}}\) où Var est la variance de la variable X, que l’on peut donc facilement estimer, et n le nombre d’expériences indépendantes. Ainsi, l’intervalle est d’autant plus grand que la demande de confiance est élevée.

L’animation suivante illustre l’estimation de π avec un intervalle de confiance.

On peut observer que plus on fait d’expériences, plus l’estimation est précise et la demi-largeur e de l’intervalle de confiance diminue. En raison de la racine carrée du nombre n au dénominateur, réduire la taille de l’intervalle par 2 demande 4 fois plus d’expériences, réduire par 10 demande 100 fois plus d’expériences , etc.

Ceci peut paraître laborieux, avec un temps de calcul proportionnel au nombre n d’expériences, qui doit être très grand. Toutefois, la simulation de Monte-Carlo est souvent la méthode la plus rapide pour résoudre certains problèmes, voire la seule méthode possible.

Simuler un modèle mathématique peut également être réalisé de manière plus ou moins subtile. La méthode standard est de simuler les lois du modèle telles quelles. Mais on peut également chercher à modifier le modèle de sorte que l’estimation soit plus précise pour un même temps de simulation. Un exemple  typique est lors de la simulation d’un événement rare, la défaillance d’un système par exemple. Comme cet événement est rare par définition, il faut un grand nombre d’expériences pour le voir intervenir. Il existe des techniques accélérant l’occurrence de l’événement sans perdre de vue la probabilité originale qu’il se produise.

Conclusion

Pour résumer, la simulation de Monte-Carlo est un outil statistique pour estimer la moyenne d’une variable aléatoire. Comme nous l’avons décrit pour le calcul d’une surface ou d’un volume, la quantité que l’on souhaite calculer n’a pas nécessairement de composantes aléatoires, mais peut être transformée sous cette forme. Cette première étape, la modélisation, est la plus importante. Ensuite, la simulation du modèle consiste à effectuer des expériences successives, à utiliser la moyenne obtenue comme estimation, et à évaluer la précision de l’estimateur. Mais ceci nécessite une source d’aléa dont les bonnes propriétés statistiques ont été prouvées.  De nombreux algorithmes ayant un comportement difficilement différenciable du hasard ont été ainsi mis au point.

Référence scientifique

B. Tuffin. La simulation de Monte Carlo. Edition Hermès, Février 2010.

Animations HTML5/JS réalisées par Ouest INSA.