Probabilités

La méthode de Monte-Carlo

Calculer avec le hasard : tirer au sort, compter, faire la moyenne — et savoir de combien on se trompe.

L'idée en une phrase

Pour connaître une quantité difficile à calculer, on la déguise en moyenne d'une expérience aléatoire, on refait l'expérience un grand nombre de fois sur ordinateur, et l'on fait la moyenne des résultats : le hasard, répété, devient un instrument de mesure.

1. TirerOn produit des nombres au hasard X₁, X₂, …, X_N, tous de même loi, indépendants les uns des autres.

2. CalculerPour chacun, on calcule une valeur f(X_i) : « dedans ou dehors », une hauteur de courbe, un gain…

3. MoyennerLa moyenne des f(X_i) approche l'espérance E[f(X)], c'est-à-dire la valeur moyenne « à l'infini ».

E[f(X)] ≈ (1/N) · Σ f(X_i) avec une erreur typique σ / √N

Deux théorèmes font tout le travail. La loi des grands nombres garantit que la moyenne finit par s'approcher de l'espérance. Le théorème central limite dit de combien elle s'en écarte : l'erreur suit à peu près une courbe en cloche de largeur σ/√N, où σ, l'écart-type, mesure combien les valeurs f(X_i) sont dispersées. On peut donc annoncer un résultat avec sa marge d'erreur, par exemple ±1,96 σ/√N pour être juste 95 fois sur 100. Toute la page tourne autour de cette petite formule : ce qu'elle promet, ce qu'elle coûte, comment la battre.

Planche 1

Viser au hasard pour mesurer π

Comment des points jetés au hasard peuvent-ils mesurer une aire ?

On lance des points uniformément dans un carré de côté 1. Le quart de disque de rayon 1 occupe une part π/4 ≈ 78,5 % du carré : la proportion de points qui tombent dedans approche donc π/4, et quatre fois cette proportion approche π. On ne mesure rien à la règle : on compte.

π ≈ 4 × (points dans le quart de disque) / (points lancés)

dedans dehors — les 20 000 premiers points sont dessinés, les suivants sont seulement comptés.

Estimation courante (bleu) selon le nombre de points, en échelle logarithmique. En vert, π et l'entonnoir π ± 1,96 σ/√N prévu par le théorème central limite.

Vitesse

Points lancés N

—

Dans le quart de disque

—

Estimation 4·k/N

—

Erreur (estimation − π)

—

Marge à 95 % annoncée ± 1,96·σ̂/√N — contient π ?

—

σ·√N mesuré sur 4 000 séries de 1 000 points (théorie 4√(p(1−p)) = 1,6422)

calcul…

Pente log-log de l'erreur quadratique moyenne (N = 100 → 10 000)

calcul…

À retenir. L'erreur décroît comme 1/√N : la pente mesurée est proche de −0,5. Pour gagner une décimale, il faut cent fois plus de points. La méthode est lente, mais elle ne demande presque rien : savoir dire « dedans » ou « dehors ». Et elle livre sa propre marge d'erreur, calculée sur les mêmes points.

Planche 2

Une intégrale tirée au sort

Comment calculer une intégrale par le hasard, et quelle façon de tirer est la meilleure ?

Une intégrale, c'est une aire ; c'est aussi une hauteur moyenne multipliée par la largeur. Si U est tiré uniformément entre a et b, alors (b − a)·E[f(U)] est exactement l'intégrale. Deux recettes s'offrent à nous : la méthode du rejet (lancer des points dans un rectangle et compter ceux qui sont sous la courbe, comme en planche 1) et la méthode de la moyenne (tirer des abscisses et faire la moyenne des hauteurs).

∫_a^b f(x) dx ≈ (b − a) · (1/N) Σ f(U_i)

Fonction Méthode Tirages N

—

Estimation de l'intégrale

—

Valeur exacte

—

Erreur commise

—

Écart-type théorique — moyenne

—

Écart-type théorique — rejet

—

Variance rejet ÷ variance moyenne : théorie / mesure (4 000 essais de 100 tirages)

—

Contrôle : Simpson (2 000 intervalles) − valeur exacte

—

À retenir. Les deux recettes visent juste en moyenne, mais la moyenne des hauteurs est toujours plus précise : le rejet jette une information (la hauteur exacte) pour ne garder qu'un « oui / non ». Sur 4/(1+x²), sa variance est 6,5 fois plus grande : il faudrait 6,5 fois plus de tirages pour la même précision. Curiosité : l'écart-type du rejet y vaut 1,6422, exactement celui de la planche 1 — c'est la même expérience déguisée.

Planche 3

Tricher honnêtement : réduire la variance

Puisque l'erreur vaut σ/√N, peut-on diminuer σ plutôt qu'augmenter N ?

On cherche ∫₀¹ eˣ dx = e − 1 ≈ 1,71828 avec n évaluations de la fonction. Chaque méthode est lancée des milliers de fois ; l'histogramme montre où tombent ses estimations. Les quatre restent centrées sur la bonne valeur (aucune n'est biaisée), mais leur largeur change du tout au tout.

Simple	moyenne de e^U : la référence.
Antithétique	on tire U et on utilise aussi 1 − U : quand l'un donne trop, l'autre donne trop peu, et les erreurs se compensent.
Variable de contrôle	on corrige par une fonction voisine dont on connaît la moyenne : e^U − c·(U − ½), avec c = 1,6903 (valeur optimale).
Échantillonnage préférentiel	on tire plus souvent là où eˣ est grand (densité 2(1+x)/3) et l'on divise par cette densité pour rester juste.

Évaluations par estimation n

Une rangée par méthode, toutes sur la même échelle horizontale centrée sur e − 1 (trait vert). En pointillé vert, la cloche prévue par la variance théorique. Estimations faites : 0 par méthode.

Méthode	variance × n théorique	variance × n mesurée	gain théorique (à travail égal)	gain mesuré

À retenir. Diviser la variance par 61, c'est obtenir avec 100 évaluations la précision que la méthode simple atteindrait avec 6 100. Les variances théoriques viennent de calculs exacts (par exemple Var e^U = (e² − 1)/2 − (e − 1)² = 0,24204) et les mesures les retrouvent à quelques pour cent. La réduction de variance consiste à injecter dans le tirage ce que l'on sait déjà du problème.

Planche 4

Là où le hasard bat la grille : la dimension

Si le hasard est si lent, pourquoi s'en sert-on ?

Parce qu'il ne craint pas la dimension. On mesure le volume de la boule de rayon 1 en dimension d avec un même budget de B points. La grille place m points régulièrement espacés sur chaque axe, soit m^d points en tout : avec un million de points, cela fait 1 000 points par axe dans le plan, mais 3 seulement en dimension 12. Le hasard tire ses B points dans le cube [−1 ; 1]^d, et son erreur σ/√B décroît toujours au même rythme 1/√B, quel que soit le nombre d'axes : la dimension n'agit que sur σ.

Volume exact : V_d = π^d/2 / Γ(d/2 + 1) (V₂ = π, V₃ = 4π/3, …)

Budget B Dimension examinée d

Part du cube occupée par la boule, V_d/2^d (échelle logarithmique). Presque tout le cube finit dans les coins.

erreur relative de la grille ; erreur relative du hasard mesurée (points) et prévue σ/√B (pointillé).

Volume exact V_d

—

Grille : points par axe, total, estimation

—

Grille : erreur relative

—

Hasard : estimation

—

Hasard : erreur relative mesurée / prévue

—

Le hasard l'emporte à partir de la dimension

—

À retenir. En dimension 2, la grille gagne largement ; dès la dimension 3 à 5 selon le budget, le hasard passe devant, et l'écart se creuse ensuite. L'erreur d'une grille dépend de la finesse par axe, qui s'effondre comme B^1/d ; celle du hasard garde le rythme 1/√B. Ici σ grandit un peu avec d (la boule devient minuscule dans le cube, on « touche » rarement), mais sans commune mesure avec l'effondrement de la grille. C'est pourquoi la physique des particules, la finance ou l'imagerie, où l'on intègre sur des dizaines ou des milliers de variables, calculent par Monte-Carlo.

Planche 5

Le hasard qu'on fabrique

Un ordinateur ne sait pas jouer aux dés. D'où vient son hasard, et en faut-il vraiment ?

Les nombres « aléatoires » d'un ordinateur sont produits par un générateur pseudo-aléatoire : une formule déterministe qui, partie d'une graine, enchaîne des nombres qui ressemblent au hasard. Toute la méthode repose sur cette ressemblance. Deux expériences la mettent à l'épreuve : un générateur célèbre qui trompait son monde, et un « faux hasard » volontaire qui fait mieux que le vrai.

RANDU, distribué par IBM dans les années 1960, calcule x_k+1 = 65 539 · x_k modulo 2³¹, puis u = x/2³¹. Chaque nombre pris seul paraît uniforme. Mais comme 65 539 = 2¹⁶ + 3, on a toujours x_k+2 = 6 x_k+1 − 9 x_k (modulo 2³¹) : trois tirages consécutifs vérifient 9u₁ − 6u₂ + u₃ = un nombre entier. Les points (u₁, u₂, u₃) ne remplissent pas le cube : ils sont rangés sur quelques plans parallèles. Faites tourner le cube à la souris ou au doigt.

Générateur Triplets

Chaque point est un triplet de tirages consécutifs (u₁, u₂, u₃) placé dans le cube unité.

Valeurs distinctes de 9u₁ − 6u₂ + u₃ parmi les triplets

—

Écart moyen de 9u₁ − 6u₂ + u₃ à l'entier le plus proche (hasard idéal : 0,25)

—

Volume d'une tranche de cube entre deux plans : vrai / estimé par Monte-Carlo

—

Écart entre deux plans voisins, 1/√(9² + 6² + 1²)

0,0921

La « tranche » est l'ensemble des points où la partie décimale de 9u₁ − 6u₂ + u₃ est entre 0,45 et 0,55 : un dixième du cube. Avec RANDU, aucun point n'y tombe jamais ; un calcul de Monte-Carlo qui en dépendrait serait faux sans le moindre signal d'alerte.

À l'inverse, un hasard trop « honnête » forme des amas et laisse des trous : rien n'empêche deux points de tomber côte à côte. La suite de Halton renonce au hasard : elle écrit le rang i en base 2 et en base 3, renverse les chiffres derrière la virgule, et obtient des points qui se répartissent de façon délibérément régulière, chaque nouveau point comblant le plus grand vide. C'est le quasi-Monte-Carlo.

Intégrande

Pseudo-hasard : amas et trous.

Halton (bases 2 et 3) : remplissage régulier.

erreur du pseudo-hasard erreur de Halton, selon N (échelles logarithmiques). Pointillés gris : pentes −1/2 et −1.

Points N

—

Pseudo-hasard : estimation (erreur)

—

Halton : estimation (erreur)

—

Valeur exacte

—

Pentes mesurées, N = 100 → 102 400 : hasard / Halton

calcul…

À retenir. Le hasard de la machine n'est qu'une imitation, et une imitation ratée fausse les calculs en silence : la tranche de RANDU reste vide, quel que soit le nombre de points. À l'opposé, le quasi-hasard mesuré ici fait décroître l'erreur presque en 1/N sur la fonction lisse (au lieu de 1/√N), un peu moins vite sur le disque au bord net. Mais il perd ce que le vrai hasard offrait gratuitement : une marge d'erreur calculée sur les points eux-mêmes. Ce dont Monte-Carlo a besoin, au fond, ce n'est pas d'imprévisible, c'est d'uniforme.

Repères historiques

1777Buffon publie le problème de l'aiguille : jetée sur un parquet, elle coupe une rainure avec une probabilité qui fait intervenir π. Laplace (1812) remarque qu'on pourrait ainsi estimer π en lançant l'aiguille.
1901Mario Lazzarini annonce π ≈ 355/113 après 3 408 lancers d'aiguille : un résultat trop beau, que l'on soupçonne aujourd'hui arrangé : pour espérer une erreur de 3·10⁻⁷, il faudrait plus de cent mille milliards de lancers.
1908William Gosset (« Student ») tire au hasard des fiches de mensurations pour vérifier empiriquement la loi qui porte son nom.
1946À Los Alamos, Stanisław Ulam, en convalescence, se demande quelle est la chance de réussir une patience ; avec John von Neumann, il applique l'idée à la diffusion des neutrons sur l'ordinateur ENIAC.
1949Nicholas Metropolis et Ulam publient « The Monte Carlo Method » ; le nom, choisi par Metropolis, évoque le casino de Monaco.
1953Algorithme de Metropolis : tirer non plus des points indépendants, mais une marche aléatoire qui explore une loi compliquée (point de départ des méthodes MCMC).
1960John Halton propose ses suites à faible discrépancie ; Ilya Sobol en donne d'autres en 1967.
1968George Marsaglia montre que les générateurs congruentiels « tombent surtout dans les plans » ; RANDU en est l'exemple le plus caricatural.
1998Makoto Matsumoto et Takuji Nishimura publient le Mersenne Twister, de période 2¹⁹⁹³⁷ − 1, longtemps le générateur par défaut des logiciels scientifiques.

En résumé

Une moyenne qui mesure

Toute quantité qui s'écrit comme une espérance — aire, intégrale, probabilité, prix — peut se calculer en faisant la moyenne de tirages au sort. L'erreur décroît en σ/√N et s'estime sur les mêmes tirages.

Diminuer σ plutôt qu'augmenter N

Variables antithétiques, variables de contrôle, échantillonnage préférentiel : injecter ce que l'on sait du problème peut diviser la variance par des dizaines, ce qui vaut des dizaines de fois plus de calcul.

L'arme des grandes dimensions

Lente en dimension 1 ou 2, la méthode devient irremplaçable quand les variables se comptent par dizaines, car son erreur ignore la dimension. Elle exige un générateur réellement uniforme ; le quasi-hasard, lui, peut faire mieux encore.

Pour aller plus loin

Chaînes de Markov (MCMC). Quand on ne sait pas tirer directement selon la loi voulue (en physique statistique, en statistique bayésienne), l'algorithme de Metropolis-Hastings construit une marche aléatoire dont la loi limite est la bonne ; les tirages ne sont plus indépendants et l'erreur se calcule autrement.

Quasi-Monte-Carlo randomisé. En décalant ou en brouillant au hasard une suite de Halton ou de Sobol (Owen, 1995), on garde la régularité du quasi-hasard tout en retrouvant une marge d'erreur, par répétition des brouillages.

Monte-Carlo multiniveau. Proposé par Michael Giles (2008) pour les équations différentielles stochastiques : on combine beaucoup de simulations grossières et peu de simulations fines, ce qui réduit fortement le coût total.

Finance. Phelim Boyle (1977) a introduit la méthode pour évaluer les options : le prix est l'espérance actualisée du gain, sur des milliers de trajectoires du cours simulées.

Tirer selon une loi. Tout part de tirages uniformes : la méthode d'inversion (utilisée ici pour l'échantillonnage préférentiel, x = −1 + √(1 + 3u)), Box-Muller pour la loi normale, le rejet pour les lois plus compliquées.