Méthode de Monte Carlo

[Sylviel]

Qu'est-ce que la méthode de Monte Carlo ?

La méthode de Monte-Carlo est une méthode pour estimer une espérance par simulation.

1) Qu'est-ce que la méthode de Monte Carlo précisément ?

Soit une fonction et une variable aléatoire telle que existe, i.e. telle que la variable aléatoire est intégrable.

La méthode de Monte Carlo consiste à simuler un grand nombre de fois, de manière indépendante, X, fournissant une liste de valeur . Alors l'espérance , inconnue, est approximativement égale à

2) Pourquoi est-ce que cela converge ?

Pour justifier proprement la méthode il faut voir chacun des N tirages comme une variable aléatoire de même loi que X. Ainsi j'appelle le premier tirage, le second...
Il s'agit bien ici de variables aléatoires et non de valeurs : peut prendre toutes les valeurs possible de , de même pour ...

La liste étant une réalisation possible de mon échantillon (plus de détails ici)

La loi des grands nombre garantit que la moyenne empirique
converge presque sûrement vers l'espérance inconnue .

"Converge presque sûrement" signifie que, toutes les simulations possibles donnerons une moyenne qui tends vers la valeur voulue sauf un ensemble négligeable (au sens mathématique du terme.

3) Quelle est l'erreur commise ?

La convergence de la méthode de Monte Carlo est garantie par la loi des grands nombre. Mais on peut aller plus loin et avoir une estimation de l'erreur commise. Comment ? Grâce au Théorème Central Limite !

Nous avons vu que les sont des variables aléatoires indépendantes et identiquement distribuées. Supposons qu'elles soit de variance finie. Ainsi le TCL nous garantit que "ressemble" à une loi normale.
Plus précisément converge en loi vers une loi Normale centrée de variance celle de que l'on note .

Cela signifie que, quand N tends vers l'infini, où est la fonction de répartition d'une loi normale centrée réduite.

En conséquence, pour N grand, on a


Cela est très classiquement utilisé de la manière suivante :

En clair, pour N grand, pour environ 95% des tirages possible, on aura la moyenne empirique qui à moins de de l'espérance voulue.

Reste une difficulté : est une valeur généralement non connue. On peut cependant l'estimer par l'écart-type des . Les résultats évoqués ci dessus reste vrais grâce au lemme de Slutsky.

4) Exemples d'utilisation

4.1) estimation de valeur d'une politique de gestion
Supposons que je gère un stock de produit face à une demande aléatoire, sur un an, au pas de temps hebdomadaire. Si j'ai de la demande et du stock je la sers, si je n'ai plus de stock la demande est supposée perdu.
Avoir des produits en stock a un coût (assurance, casse, cash immobilisé...), ne pas servir une demande a un coût (mécontentement client), passer une commande de renouvellement a un coût...

J'ai décidé d'une politique de gestion a priori (typiquement : passer commande d'une quantité fixe quand j'atteint un certain stock minimum) et je voudrais évaluer le coût de cette stratégie.

Supposons la demande indépendantes d'une semaine à l'autre, et pouvant prendre 10 valeurs par semaine.
Alors calculer l'espérance exacte de mon coût nécessite de faire la moyenne de 10^52 valeurs ! Ce nombre est titanesque (largement supérieur au nombre d'atome dans l'univers). Donc même dans un cas purement discret calculer l'espérance exacte est numériquement hors d'atteinte.

Heureusement on va pouvoir utiliser Monte Carlo pour estimer l'espérance du coût. Et ce de manière très simple :
- simuler une année et calculer le coût cumulé sur l'année, qu'on note
- recommencer fois.
- calculer la moyenne des N couts ainsi que leur écart-type
- Il y a environ 95% de chance que la vraie espérance, inconnue, soit dans


4.2) Calcul d'intégrale par Monte Carlo

L'une des grandes réussite de la méthode de Monte Carlo est la calcul d'intégrale en grande dimension.
En effet, si est une variable aléatoire de densité alors l'espérance
est donnée par .

Pour calculer ce genre d'intégrale, quand une formule analytique n'existe pas, on peut utiliser des méthodes d'intégration numériques qui discrétise l'espace et approxime "l'aire sous la courbe" par des rectangle, trapèze ou autres approche plus compliquées. Malheureusement toutes les méthodes d'intégration numériques sont sensible à la dimension : pour atteindre une erreur voulue le nombre de point de discrétisation augmente (dramatiquement) avec le nombre de dimension. Au contraire Monte Carlo est insensible à la dimension, et est donc la méthode numérique par défaut dès qu'on dépasse la dimension 4 ou 5...

Exemple pratique : calculer le volume d'une boule en dimension n.
a) initialiser un compteur à count 0
b) tirer n valeurs uniforme sur [0,1], lnoté ce qui nous donne un vecteur de dimension n tiré uniformément dans .
c) calculer la norme au carré du vecteur : , si elle est inférieure à 1 ajouter 1 au compteur
d) répéter b) et c) N fois (au moins 1000).
p=Count/N nous donne l'approximation du volume du quartier de boule ayant des coordonnées positive
et il y a 95% de chance que cette approximation soit à moins de du vrai volume
(on peut améliorer par ).
Le volume total estimé de la sphère étant naturellement

Un petit code qui simule cela https://repl.it/repls/WhichDeafeningPro ... ro#main.py
On peut y observer :
1) que l'estimation est proche de la valeur réelle
2) que l'estimation est aléatoire (à chaque fois que je relance une simulation j'ai une valeur différente)
3) que la vraie valeur peut être hors de l'intervalle d'incertitude (cela arrive en gros 5% du temps). Si on veut moins d'incertitude il suffit de prendre un intervalle plus grand, par exemple 3 écarts-type donnerons 1% environ (voire la section 2 plus haut).

P.S: si toute cette explication est de mon cru elle est néanmoins fort classique. Si quelqu'un voit une erreur qu'il n'hésite pas à la signaler.

[Sylviel]

Petit complément autour du Lemme de Slutsky.

La loi des grands nombre nous garantit que tends vers l'espérance recherchée .
Le TCL nous donne une information sur l'erreur (inconnue) entre l'estimation et la vraie valeur
ressemble, pour N grand, à une loi normale centrée d'écart-type .

Pour préciser mathématiquement cette ressemblance on dira que
converge en loi vers une Gaussienne centrée réduite G.

Malheureusement est très souvent une valeur inconnue (puisque l'espérance elle même est inconnue...).

Heureusement, tout comme on peut estimer par sa moyenne empirique, on peut également estimer le vrai à partir des données. Un estimateur classique est le suivant
.
On montre aisément que converge presque sûrement vers .

Le Lemme de Slutsky nous permet de montrer que le couple , dont la première coordonnée converge en loi vers G et la seconde converge presque sûrement vers converge en loi vers le couple .

La convergence en loi du couple, et la continuité de la multiplication, permet de garantir que, si on remplace la vraie valeur de par son estimation dans la définition de on conserve le résultat de convergence en loi.

Tout ceci pour dire tout simplement que le résultat, asymptotique (i.e. "quand N tends vers l'infini") est le même qu'on utilise la vraie valeur de l'écart-type ou un estimateur convergeant.

Ces détails techniques sont très souvent glissé sous le tapis dans les cours, et je dirais plutôt à raison : c'est un peu délicat à suivre et n'apporte pas grand chose.

Il est néanmoins bon de bien comprendre que le contrôle de l'erreur via le TCL exploite le vrai écart-type (inconnu) et qu'il faut travailler un peu plus pour utiliser l'estimateur de l'écart-type.

[ijkl]

Merci Sylviel

Je ne sais pas si je vais tout comprendre mais en tout cas un grand merci

...de toute façon j'ai pas le choix

[Sylviel]

Je ne sais pas si je vais tout comprendre

Tu peux poser des questions si tu veux.

...de toute façon j'ai pas le choix

[ijkl]

je suis punk donc ça explique que j'ai pas le choix

pour contrebalancer le délire du punk , j'ai besoin de comprendre les maths (qui sont le contraire d'un délire)

Merci à vous Sylviel

[ijkl]

NB : je ne vous poserais pas de questions Sylviel (mais vraiment merci)

Dans Kill Bill 2 (si mes souvenirs sont bons : elle doit creuser avec ses propres mains dans son cercueil mais elle n'aurait pas été capable de le faire si son Maître l'aurait encouragé à demander de l'aide)

et en ce qui me concerne je n'ai pas envie de devoir demander de l'aide si jamais je me retrouve dans un cercueil

en tout cas je ne peux que vous remercier (vous avez déjà fait le principal pour moi)

[Sylviel]

Utiliser Monte-Carlo pour calculer une intégrale sur [0,1]

Détaillons un peu comment utiliser la méthode de Monte-Carlo pour calculer l'intégrale d'une fonction continue de R sur un [0,1]: .
Bien entendu il s'agit ici d'un exemple pédagogique, la méthode n'est numériquement intéressante qu'en dimension supérieure (généralement a partir de la dimension 7-8).

Soit une fonction continue. Soit un segment.
Soit une variable uniforme sur , donc de densité
sur et 0 ailleurs.

Comme est continue elle est bornée sur et donc l'espérance de est bien définie.
Par définition de l'espérance (dans le cas d'une variable à densité) on a


Maintenant que nous avons exprimé notre intégral comme un calcul d'espérance on peut appliquer Monte Carlo. Pour mémoire :
- simuler un nombre N de , uniformément dans [0,1]
- pour chacun calculer l'image par f
- calculer la moyenne empirique sur le tirage en question

La loi des grands nombre nous garantit que converge vers l'espérance et donc vers l'intégrale à calculer.

Le TCL nous dis que si je répètes un grand nombre de fois l'estimation de l'intégrale, pour obtenir un grand nombre de , et que je fais l'histogramme des il aura une forme approximativement Gaussienne. Il ne dis bien entendu rien sur l'histogramme des d'un tirage.

Couplé au lemme de Slutsky le TCL permet de quantifier l'erreur entre l'estimation et la vraie valeur de l'intégrale :
soit l'estimateur de l'écart-type.
Alors on a que, quand N tends vers l'infini,
.

Cela se vérifie très bien numériquement :
- choisir une fonction f dont on peut calculer l'intégrale
- choisir N grand (au moins 1000)
- Calculer 10 000 fois (avec à chaque fois un tirage différent).
Compter le nombre de fois où ne sera pas dans : ce sera environ 500.

Voici le code pour vérifier cela : https://repl.it/repls/DefinitiveAwfulThing#main.py

P.S: le code est fait pour des bornes quelconques b>a

[Sylviel]

Bien entendu l'exemple ci-dessus est purement pédagogique, la méthode n'est numériquement intéressante qu'en dimension supérieure.

En petite dimension on utiliseras une méthode déterministe,
la plus mauvaise étant la méthode des rectangles (voir par exemple ici).

Que se passe-t-il en dimension supérieure ?

Et bien supposons que f soit une fonction continue de dans . Toutes les méthodes numériques exactes vont supposer une discrétisation de du cube, généralement en discrétisant chaque segment en n éléments, l'erreur étant contrôlée par 1/n.
Le problème est le suivant : si je veux discrétiser mon segment en 10 points, alors je discrétise le cube en 10^d points. On voit tout de suite la difficulté en dimension 20 ou 30...

A l'opposé la méthode de Monte Carlo a une précision insensible à la dimension, ce qui la rends indispensable en grande dimension. A noter qu'en dimension moyenne (disons de 5 à 20) il existe des méthodes de quasi-Monte Carlo (qui choisissent des xi "bien réparti" plutôt que tiré au hasard) qui sont plus efficace que Monte Carlo (et que les méthodes déterministes).

En dimension 1, pour des fonctions peu régulières, la méthode de Monte Carlo peut toutefois être plus intéressante qu'une méthode déterministe (voir par exemple :http://www.tangentex.com/MonteCarlo.htm).

[pascal16]

Comme matheux libre d'esprit et de trop grandes années de bourrage de crâne en math, après une recherche perso, sur la nullité de la méthode des rectangles :

1.1/ imaginons une fonction croissante, je peux faire dessus une méthode des rectangles aux n points "k/n", pour k=0 à n-1, j'ai une estimation par défaut, aux point k=1 à n, j'ai une valeur par exès. En faisant la moyenne, j'ai une valeur plus juste.
Si je met des coefficients à chacun des points de valuation f(k/n) pour k=0 à n, j'ai comme coefficient :
1/2-1-1-..... 1-1-1-1/2

faisons la méthode des trapèzes, on a comme coefficients : 1/2-1-1-..... 1-1-1-1/2

ie : améliorer la méthode des rectangle donne la méthode des trapèzes avec exactement la même erreur et le même nombre de calculs

NB : à l'heure des IA qui ne font que des gradients, je trouve normal de penser à améliorer une méthode par gradient (le haut du trapèze est une corde qui est l'approximation d'une dérivée) plutôt qu'avoir une convergence lente en augmentant n.

1.2/ si on va plus loin, à caque méthodes, on a des pages entières de calcul pour l'évaluation de l'erreur en fonction des dérivées pour les fonction dérivables. En approximent la dérivée, on a une méthode améliorée par gradient qui redonne une série de coefficients et on retombe systématiquement sur les méthodes dite "meilleurs". Même les amélioration de converge des suites redécoupage en 2 (n éval, 2n éval, 4n éval..) peuvent aboutir au même résultats.

1.3/ au final toute méthode utilisant des valuation f(k/n) pour k=0 à n ne sont toutes de la famille de la méthode des "rectangle'. La méthode des rectangle n'est pas nulle, c'est au contraire la base de toute une famille (c'est pour moi comme dire "les nombres entiers, c'est nul").

2/ MonteCarlo : Je n'ai pas vu ressortir le méthode où le choix de xi se fait de manière récursive en fonction de l'écart-type car comme dans le cas des rectangles, la répartition uniforme, c'est pas top. Dans le cadre du quart de cercle pour pi/4, ça peut permettre de comprendre ce qu'est un bon redécoupage

[Sylviel]

Bonjour Pascal,

j'ai un peu de mal à comprendre ton message.

Tout d'abord je ne dis pas que la méthode des rectangles est nulle mais que c'est la plus mauvaise des méthodes de calcul. Je pourrais préciser : c'est celle qui donne la plus mauvaise garantie sur l'erreur commise.

Je ne comprends pas vraiment ce que tu veux dire par "améliorer les rectangles par gradient". Mais si tu regardes le document en lien ci-dessus il montre bien qu'il y a tout une classe de méthodes de calcul d'intégrales (en approximant par des polynomes) et que les rectangles correspond à l'approximation par des polynomes constant.

Pour l'intérêt de Monte Carlo par rapport aux méthodes déterministe : le principal point c'est que les méthodes déterministes vont a priori discrétiser chaque dimension. Donc si d est grand (par exemple d=100) il faudra au moins 2^d points (ce qui est rédhibitoire). Alors que Monte Carlo pourras très bien faire le calcul.

Je n'ai

[pascal16]

La méthode des rectangle en pris en k/n pour k=0 à n-1 sous estime l'aire sous la courbe
La méthode des rectangle en pris en k/n pour k=1à n sur estime l'aire sous la courbe
En corrigeant l'erreur (la moyenne des deux) on retombe sur la formule des trapèzes

la méthode des trapèzes , sous-estime l'aire d'une portion convexe.
Si on compare l'erreur pour n points et 2n points, on peut faire une correction de cette erreur il me semble par (1/3)*(4*Aire(2n)-Aire(n)) qui est aussi un résultat sur l'accélération de convergence des suites.
Elle donnera la méthode de Simpson.
Donc en améliorant 2 fois la méthode des rectangles, on a celle de Simpson.

Donc, sur les méthodes utilisant la valuation en k/n pour k=0 à n, si on travaille intelligement en corrigeant les erreurs remarquées, on a à nombre de calculs identiques, une erreurs identiques pour toutes les méthodes une fois améliorées. C'est la façon par laquelle elles te sont présentées qui tu laisse penser qu'elle sont issues de méthodes totalement différentes.

Je ne comprend pas l'esprit mathématique très agrego-universito-capesien qui consiste à remplir une page entière de calculs pour trouver une erreur chiffrée que l'on a encadré 100 000 000 de fois dans l'histoire des mathématiques tout ça pour ne pas faire quelque chose de concret avec cette étude.

PS : si tu passes un concours, un oral, fait comme dans le bouquin, la page de calcul et rectangle < trapèze < Simpson < polynome de degré 3 < polynome de degré 4....


pour ce qui est des de "la valuation en k/n pour k=0 à n"
en fait quand tu fais le calcul sur ordinateur et que tu ne veux appeler qu'une seul fois f(k/n), tu dois rassembler les coefficients de chaque méthodes. Ca te donne une liste de coefficient pour chaque abscisses k/n.
rectangle [0 1 1 1 1 ... 1 1 ] ou [1 1 1 1.... 0]
trapeze [1 2 2 2 2 ...2 2 1 ] / 2
Simpson [ 1 4 2 4 2 4 2... 4 2 4 1] /3
du coup, elle paraissent moins éloignées les unes des autres