Méthodes de Simulation
Séances 5 et 6
28 avril 2025
Intégration de Monte-Carlo
▶ Dans la section précédante, on a vu les techniques de
simulation permettant de générer des variables aléatoires par
ordinateur.
▶ la présente section introduit les concepts majeurs des
méthodes de Monte-Carlo.
▶ Il s’agit de tirer avantage de la disponibilité des variables
aléatoires générées par ordinateur pour approcher des
intégrales unidimensionnelles et multidimensionnelles.
▶ Nous introduisons les notions de base des approximations de
Monte-Carlo comme un sous-produit de la Loi des Grands
Nombres
▶ Nous présentons l’universalité de cette approche en insistant
sur la versatilité de la représentation d’une intégrale comme
espérance.
▶ Par la suite on traitera la résolution des problèmes
d’optimisation par des techniques de simulation.
Inférence
▶ En inférence statistique, les 2 caractéristiques des problèmes
numériques sont l’optimisation et l’intégration.
▶ De nombreux exemples montrent qu’il n’est pas toujours
possible de calculer analytiquement les estimateurs associés à
un paradigme donné (que ce soit par maximum de
vraisemblance, par des techniques bayésiennes ou par la
méthode des moments).
▶ Par conséquent, quel que soit le type d’inférence statistique,
nous sommes souvent amenés à considérer des solutions
numériques.
Inférence
▶ On a déjà vu une variété de méthodes pour la génération par
ordinateur de variables aléatoires de distributions données
▶ Une des solutions les plus standard est bien sûr d’utiliser la
simulation, soit de la distribution originale, soit d’une
distribution de substitution, pour calculer la quantité qui nous
intéresse.
▶ Notons que la possibilité de produire une quantité presque
infinie de variables aléatoires distribuées selon une distribution
donnée ouvre l’accès à l’utilisation de résultats fréquentistes et
asymptotiques beaucoup plus facilement que dans le cas de
l’inférence usuelle, où la taille de l’échantillon est souvent
fixée.
▶ On peut donc appliquer des résultats probabilistes comme la
Loi des Grands Nombres ou le Théorème de la Limite
Centrale, puisqu’ils autorisent tous deux une évaluation de la
convergence des méthodes de simulation (équivalente aux
bornes déterministes utilisées par les approches numériques).
Intégration numérique
▶ Avant de nous engager dans la description des techniques de
Monte-Carlo, notons que l’alternative, a priori naturelle, à
l’utilisation des méthodes de simulation pour les intégrales de
la forme Z
h(x)f (x)dx,
χ
où f est une densité de probabilité, serait de se reposer sur des
méthodes numériques comme celle de Simpson et les règles du
trapèze.
▶ Par exemple, R propose deux fonctions apparentées qui
dirigent l’intégration unidimensionnelle, area dans la librairie
MASS et integrate.
▶ Cependant, area ne s’applique pas aux intégrales avec des
bornes infinies et nécessite donc des connaissances préalables
de la région d’intégration.
▶ L’autre fonction, integrate, accepte des bornes infinies
mais malheureusement, elle est très fragile et peut produire
des résultats faux.
Précision de la fonction integrate
▶ Nous comparons l’utilisation de integrate sur l’intégrale
Z ∞
xλ−1 exp(−x)dx = Γ(λ)
0
avec le calcul de Γ(λ) fourni par la fonction gamma.
▶ En implémantant cette comparaison avec le code suivant :
Précision de la fonction integrate
▶ Nous obtenons une suite qui montre qu’il n’y a aucune
divergence entre les deux évaluations même pour une très
faible valeur de λ.
▶ Une des difficultés des méthodes d’intégration numérique
comme integrate est qu’elles manquent souvent la région
la plus significative pour la fonction à intégrer.
▶ Au contraire, les méthodes de simulation ciblent
naturellement cette région en exploitant l’information
apportée par la densité de probabilité associée aux intégrales.
integrate versus area
▶ Considérons un échantillon de 10 variables aléatoires de
Cauchy xi (1 ≤ i ≤ 10) avec le paramètre de position θ = 350.
▶ La marginale de l’échantillon est alors
Z ∞Y 10
1 1
m(x) = dθ.
−∞ π 1 + (xi − θ)2
i=1
▶ Cependant, integrate retourne une valeur numérique
fausse
integrate versus area
▶ integrate ne réussit même pas à signaler la difficulté car
l’erreur d’évaluation est ridiculement petite.
▶ De plus, le résultat n’est pas comparable à area :
integrate versus area
▶ Utiliser area dans ce cas produit une évaluation plus sûre car
la commande area(lik,-a,a) reste relativement stable
quand a augmente, mais cette évidence requiert des
connaissances préalables à propos de la position du mode de
l’intégrande.
Figure: Comparaison de integrate et de area pour l’intégrale d’une
probabilité de Cauchy à l’échelle logarithmique (le résultat de area
correspond à la courbe en plein dans le graphe ci-dessus).
integrate versus area
▶ Les avancées continuelles de l’intégration numérique ne
peuvent pas faire face aux intégrales hautement (voire
modérément) multidimensionnelles qui sont la règle dans les
problèmes statistiques.
▶ Elaborer des astuces fines d’intégration spécifiques pour ces
problèmes serait difficile, d’autant que nous pouvons tirer
profit de la nature probabiliste de ces intégrales.
L’intégration classique de Monte-Carlo
▶ Avant de proposer des simulations comme réponse à
l’intégration, rappelons les propriétés qui justifient leur usage.
▶ Le problème générique est d’évaluer l’intégrale
Z
Ef [h(X)] = h(x)f (x)dx,
χ
où χ dénote l’espace où la variable aléatoire X prend ses
valeurs, qui est d’habitude égal au support de la densité f .
▶ Le principe de la méthode de Monte-Carlo pour approcher
cette intégrale est de générer (par ordinateur) un échantillon
(X1 , · · · , Xn ) suivant la densité f et de proposer comme
approximation la moyenne empirique
n
1X
hn = h(xj )
n
j=1
car hn converge presque sûrement, c’est-à-dire pour presque
toutes les suites générées, vers Ef [h(X)] par la Loi Forte des
Grands Nombres.
L’intégration classique de Monte-Carlo
▶ Quand h2 (X) a une variance finie sous f , la convergence de
√
hn se produit à la vitesse O( n) et la variance asymptotique
de l’approximation est
Z
V(hn ) = (h(x) − Ef [h(X)])2 f (x)dx
χ
▶ Cette variance peut aussi être estimée à partir de l’échantillon
(X1 , · · · , Xn ) par
n
1 X 2
vn = 2
h(xj ) − hn ) .
n
j=1
▶ Grâce au TCL, pour n suffisamment grand,
hn − Ef [h(X)]
√
vn
est asymptotiquement distribué comme une N (0, 1) ce qui
conduit à un test de convergence et aux bornes de confiance
sur l’approximation de Ef [h(X)]
Exemple de Convergence de Monte Carlo
▶ Considérons l’évaluation de l’intégrale sur [0, 1] de la fonction
artificiellement choisie
h(x) = [cos(50x) + sin(20x)]2
▶ Ce problème peut être vu comme une intégration uniforme.
▶ Par conséquent, nous générons des variables U1 , U2 , · · · , Un
R P
i.i.d. U[0,1] et approchons h(x)dx par h(Ui )/n.
Exemple de Convergence de Monte Carlo
▶ La Figure montre comment la convergence se produit, ainsi
que les bornes dérivées des erreurs standard estimées en
fonction du nombre n de simulations.
Figure: Approximation de l’intégrale: (haut) la fonction (bas) moyenne
± deux écarts types en fonction des itérations pour une seule série de
simulations.
Exemple de Convergence de Monte Carlo
▶ Notons que la bande de confiance produite dans cette figure
est bien une bande de confiance de 95% au sens classique
▶ C’est-à-dire qu’elle ne correspond pas à une bande de
confiance sur le graphe de convergence des valeurs estimées,
mais plutôt à l’intervalle de confiance que vous pouvez
produire pour tout nombre d’itérations sur lequel vous voulez
vous arrêter.
▶ Alors que le bonus apporté par l’évaluation simultanée de
l’erreur par l’estimation de Monte-Carlo est indiscutable, vous
devez rester attentif au fait que cette évaluation n’est valide
que lorsque vn est elle-même une estimation convergente de la
variance de hn .
▶ Dans des situations critiques où vn ne converge pas du tout
ou même ne converge pas assez vite pour appliquer un
Théorème de la Limite Centrale, cette estimation et la région
de confiance asssociée à cette dernière ne peuvent pas être
vérifiées.
Exemple de Convergence de Monte Carlo
▶ Quand vous évaluez la convergence des méthodes de
Monte-Carlo, problème qui sera traité en détail plus tard, la
commande R cumsum est assez pratique pour ce qui est de
calculer toutes les sommes partielles d’une suite à la fois et
pour autoriser une représentation immédiate de la suite
correspondante des estimateurs.
▶ La méthodologie de Monte-Carlo illustrée par l’exemple
ci-dessus peut être implémentée dans un grand nombre de cas
où les distributions impliquées dans le modèle peuvent être
simulées.
▶ Par exemple, nous pourrions utiliser les sommes de
Monte-Carlo pour calculer une fonction de répartition normale
(même si cette fonction est maintenant disponible dans tous
les logiciels de calcul et sur la plupart des calculatrices).
Approximation de la fonction de répartition normale
▶ Avec un échantillon normal N (0, 1) de taille n, (x1 , · · · , xn ),
l’approximation de
Z t
1 2
Φ(t) = √ e−y /2 dy
−∞ 2π
par la méthode de Monte-Carlo est
n
1X
Φ(t)
b = 1xi ≤t ,
n
i=1
avec une variance (exacte) de Φ(t)[1 − Φ(t)]/n
▶ car les fonctions 1xi ≤t sont des variables de Bernoulli
indépendantes avec probabilité de succès Φ(t).
Approximation de la fonction de répartition normale
▶ L’implémentation R est la suivante :
▶ Pour des valeurs de t autour de t = 0, la variance est alors
approximativement 1/4n et pour avoir pune précision de quatre
décimales, nous avons besoin de 2 × 1/4n ≤ 10−4 , soit
environ n = (104 )2 = 108 simulations.
▶ Notons qu’une plus grande précision (absolue) est atteinte
dans les queues et que des méthodes de simulation
(beaucoup) plus efficaces pourraient être utilisées.
Approximation de la fonction de répartition normale
▶ L’approximation de Monte-Carlo d’une fonction de répartition,
a des applications originales.
▶ Par exemple, elle peut être utilisée lors de la construction de
la distribution d’un test statistique comme un test du rapport
de vraisemblance sous une hypothèse nulle, aussi bien que
pour approcher sa puissance sous l’alternative.
▶ On pourrait penser à ce stade que la méthode de Monte-Carlo
introduite dans cette section est suffisante pour approcher les
intégrales de manière totalement contrôlée.
Approximation de la fonction de répartition normale
▶ En vérité, même si la méthode de Monte-Carlo directe fournit
de bonnes approximations dans les cas les plus communs, il
existe des alternatives plus efficaces qui n’évitent pas
seulement une simulation directe de f mais peuvent aussi être
utilisées plusieurs fois de suite pour plusieurs intégrales.
▶ L’utilisation peut de plus être recyclée soit pour une famille de
fonctions h ou pour une famille de densités f .
▶ De plus, les problèmes de queues de simulations illustrés dans
l’exemple précédent peuvent être résolus beaucoup plus
efficacement qu’en simulant suivant f car simuler des
événements de très petite probabilité nécessite un très grand
nombre de simulations sous f pour obtenir une précision
(relative) donnée.
Echantillonnage préférentiel
▶ La méthode que nous étudions à présent est appelée
échantillonnage préférentiel car elle relie une fonction
instrumentale, qui est une densité, à la distribution d’intérêt,
f , aussi appelée distribution cible.
▶ En fait, une évaluation de
Z
Ef [h(X)] = h(x)f (x)dx,
χ
fondée sur des simulations suivant f n’est presque jamais
optimale au sens où utiliser des distributions alternatives pour
la simulation peut diminuer la variance de l’estimateur
correspondant.
▶ La méthode d’échantillonnage préférentiel est fondée sur une
représentation alternative de Ef [h(X)].
Un changement arbitraire de mesure de référence
▶ Pour une densité arbitraire donnée g qui est strictement
positive quand h × f est différent de zéro, nous pouvons
réécrire
Z
f (x) h(X)f (X)
Ef [h(X)] = h(x) g(x)dx = Eg
χ g(x) g(X)
comme espérance sous la densité g.
▶ Notons que χ est encore l’espace où X ∼ f (x) prend sa valeur
et il peut donc être plus petit que le support de la densité g.
▶ Cette identité fondamentale de l’échantillonnage préférentiel
justifie l’usage de l’estimateur
n
1 X f (Xj )
h(Xj ) −→ Ef [h(X)]
n g(Xj )
j=1
fondé sur un échantillon X1 , · · · , Xn généré par g et non par
f.
Un changement arbitraire de mesure de référence
▶ Bien sûr, comme Ef [h(X)] peut aussi être écrit comme
espérance sous g, la dernière convergence est due pour la
même raison que l’estimateur régulier de Monte-Carlo hn
converge, quel que soit le choix de la distribution g (aussi
longtemps que supp(h × f ) ⊂ supp(g).
▶ Cette propriété universelle est reliée au fait que l’intégrale
Ef [h(X)] peut être représentée par une infinité de paires
(h, f ) et donc qu’une intégrale donnée n’est pas
intrinsèquement associée à une distribution donnée.
Un changement arbitraire de mesure de référence
▶ A l’opposé, il y a une liberté presque absolue dans sa
représentation en tant qu’espérance.
▶ La contrainte sur le support de g est fondamentale dans le
sens qu’utiliser un support plus petit tronque l’intégrale et
produit donc un résultat biaisé.
▶ Ceci signifie, en particulier, que quand on considère une
solution non paramétrique pour g, le support du noyau ne doit
pas être restreint.
Exemple d’approximation de la probabilité de queue
▶ Les probabilités de queue sont très mal approchées en utilisant
des sommes standard de Monte-Carlo une fois rendu assez
loin dans la queue.
▶ Par exemple, si Z ∼ N (0, 1) et si nous sommes intéressés par
la probabilité P (Z > 4.5), qui est très petite, simuler
Z (i) ∼ N (0, 1) produit une seule fois un succès sur environ 3
millions d’itérations !
▶ Bien sûr, le problème est que nous sommes maintenant
intéressés par la probabilité d’un événement très rare et donc
une simulation naı̈ve de f nécessitera un grand nombre de
simulations pour avoir une réponse stable.
▶ Mais, grâce à l’échantillonnage préférentiel, nous pouvons
grandement améliorer notre précision et donc diminuer le
nombre de simulations par de nombreux ordres de magnitude.
Exemple d’approximation de la probabilité de queue
▶ Par exemple, si nous considérons une distribution avec un
support restreint à [4.5, ∞[, l’additionnelle et inutile variation
de l’estimateur de Monte-Carlo due aux zéros simulés (quand
x < 4.5) disparaı̂t.
▶ Un des choix possibles est de prendre pour g la densité d’une
distribution exponentielle Exp(1) tronquée à 4.5,
Z ∞
−y
g(y) = e / e−x dx = e−(y−4.5) ,
4.5
▶ L’estimateur par échantillonnage préférentiel de la probabilité
de queue est
n n 2
1 X f (Y (i) ) 1 X e−Yi /2+Yi −4.5
= √ ,
n g(Y (i) )
i=1
n
i=1
2π
où les Yi sont des générations suivant g.
Exemple d’approximation de la probabilité de queue
▶ Le code correspondant est
▶ On obtient 3.312 × 10−6 versus la vraie valeur 3.398 × 10−6 .
Figure: Convergence de l’approximation de l’échantillonnage préférentiel
de la probabilité de queue normale P (Z ≥ 4.5). La ligne droite
correspond à la vraie valeur de l’intégrale.