Introduction aux méthodes de Monte Carlo

Introduction

Deux grandes familles de problèmes numériques apparaissent en analyse de données : l’optimisation et l’intégration. Il n’est pas toujours possible de calculer analytiquement les estimateurs associés à un modèle donné, et l’on doit souvent recourir à des solutions numériques. Une façon de contourner ce problème consiste à utiliser la simulation. L’estimation de Monte Carlo consiste à simuler des tirages hypothétiques d’une loi de probabilité afin de calculer des quantités caractéristiques de cette distribution.

L’idée de base de Monte Carlo consiste à écrire l’intégrale comme une espérance par rapport à une certaine distribution de probabilité, puis à l’approximer à l’aide de la méthode des moments ($E[g(X)] \\approx \\overline{g(X)} = \\dfrac{1}{n}\\sum g(X_{i})$).

Si l’on a une fonction continue $g(\\theta)$ et que l’on souhaite l’intégrer sur l’intervalle (a,b), on peut réécrire l’intégrale comme une espérance sous une loi uniforme $U \\sim U[a,b]$, à savoir :

En utilisant l’estimateur des moments, notre approximation d’intégrale est :

Où les

sont des valeurs simulées d’une loi uniforme.

Exemple 1 : approximation d’une intégrale exponentielle

1) Étant donnée la fonction $f(x)= e^{x}$, l’intégrale sur l’intervalle [3,5] est :

L’approximation Monte Carlo de l’intégrale est :

#       Définition de la fonction cible\nf = function(x){return(exp(x))}\n#       Définition de la fonction d'erreur absolue\nerror = function(x,y){return(abs(x-y))}\n#       Valeur exacte de l'intégrale\nans = exp(5)-exp(3)\n\nset.seed(6971)\n#       nombre d'itérations\nn = 10^2\n#       données uniformes simulées\nx= runif(n,3,5)\n#       Approximation Monte Carlo\nMCa= (5-3)*mean(f(x))\n#       Erreur d'approximation\ne = error(ans,MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\nset.seed(6971)\nfor(k in 3:6){\n  n = 10^k\n  x = runif(n,3,5)\n  mca = (5-3)*mean(f(x))\n  rest = rbind(rest,c(n,mca,error(ans,mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Résultats de l'approximation Monte Carlo de l'intégrale\",\n      col.names =c(\"Nombre de simulations\",\"Approximation Monte Carlo\",\"Erreur d'approximation\"))

Approximation de Monte Carlo généralisée

Dans le cas général, l’approximation d’intégrale pour une distribution f donnée est :

Un algorithme pour construire $\\widehat{I}$ peut se décrire comme suit :

1) Générer depuis une distribution f

2) Calculer :

3) Obtenir la moyenne empirique : $$\\overline{I} = \\dfrac{1}{n}\\sum_{k=1}^{n}\\dfrac{g(\\theta_k)}{f(\\theta_k)}$$

Dans le bloc suivant, nous présentons une fonction simple d’approximation Monte Carlo pour illustrer le fonctionnement de l’algorithme, où a et b sont les paramètres de la densité uniforme, n le nombre de simulations souhaitées, et f la fonction que l’on souhaite intégrer.

# Fonction Monte Carlo simple\nMCaf = function(n,a,b,f){\nx = runif(n,a,b)\nMCa = (b-a)*mean(f(x))\nreturn(MCa)\n}

Méthodes de Monte Carlo en analyse bayésienne des données

L’idée centrale de l’analyse bayésienne consiste à ajuster un modèle ( par exemple une régression ou une série temporelle ) via une approche d’inférence bayésienne. Nous supposons que nos paramètres d’intérêt suivent une distribution théorique ; cette distribution (a posteriori) est mise à jour à partir de la distribution des données observées (vraisemblance) et des informations préalables ou externes sur nos paramètres (distribution a priori) à l’aide du théorème de Bayes.

$$p(\\theta / X) \\text{ } \\alpha\\text{ } p(X/\\theta)p(\\theta)$$ Où :

$p(\\theta/X)$ est la distribution a posteriori du paramètre

$p(X/\\theta)$ est la distribution d’échantillonnage des données observées ( vraisemblance ).

$p(\\theta)$ est la distribution a priori du paramètre.

Le principal défi de l’approche bayésienne est l’estimation de la distribution a posteriori. Les méthodes de Monte Carlo par chaînes de Markov ( MCMC ) génèrent un échantillon de la distribution a posteriori et approchent espérances, probabilités ou quantiles via des méthodes de Monte Carlo.

Dans les deux sections suivantes, nous proposons deux exemples pour approximer des probabilités et des quantiles d’une distribution théorique. La procédure décrite ci-dessus correspond aux méthodologies usuelles de l’approche bayésienne.

Exemple 2 : approximation de probabilité pour une loi gamma

Supposons que nous voulions calculer la probabilité qu’une variable aléatoire $\\theta$ soit comprise entre 0 et 5, $P(0 < \\theta < 5)$, où $\\theta$ suit une loi gamma de paramètres a = 2 et b = 1/3 ($\\theta \\sim Gamma(a = 2,b = 1/3)$). Cette probabilité s’écrit :

Où $I_{[0,5]}(\\theta) = 1$ si $\\theta$ appartient à l’intervalle [0,5]. L’idée de l’approximation de Monte Carlo est de compter le nombre d’observations appartenant à l’intervalle [0,5], puis de le diviser par le total des données simulées.

set.seed(6972)\n#       nombre d'itérations\nn = 10^2\n#       données simulées (gamma)\nx= rgamma(n,shape = 2,1/3)\n#       Approximation Monte Carlo\nMCa= mean(x <= 5)\n#       Erreur d'approximation\ne = error(pgamma(5,2,1/3),MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\n\nfor(k in 3:6){\n  n = 10^k\n  x= rgamma(n,shape = 2,1/3)\n  mca= mean(x <= 5)\n rest = rbind(rest,c(n,mca,error(pgamma(5,2,1/3),mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Résultats de l'approximation Monte Carlo de la probabilité\",\n      col.names =c(\"Nombre de simulations\",\"Approximation Monte Carlo\",\"Erreur d'approximation\"))

Exemple 3 : approximation d’un quantile d’une loi normale

Supposons que nous voulions calculer le quantile 0,95 d’une variable aléatoire $\\theta$ suivant une loi normale de paramètres $\\mu = 20$ et $\\sigma = 3$ ($\\theta \\sim normal(\\mu = 20,\\sigma^{2} = 9)$). Le quantile 0,95 est alors :

L’idée principale est de trouver la plus grande valeur d’échantillon donnant une probabilité inférieure ou égale à 0,95. L’approximation Monte Carlo du quantile s’estime via la fonction quantile() appliquée aux données simulées.

set.seed(6973)\n#       nombre d'itérations\nn = 10^2\n#       données simulées (normales)\n  x= rnorm(n,20,3)\n#       Approximation Monte Carlo\nMCa= quantile(x,0.95)\n#       Erreur d'approximation\ne = error(qnorm(0.95,20,3),MCa)\n\nrest =  data.frame(n = n,MCapprox = MCa,error = e)\n\nfor(k in 3:6){\n  n = 10^k\n  x= rnorm(n,20,3)\n  mca= quantile(x,0.95)\n rest = rbind(rest,c(n,mca,error(qnorm(0.95,20,3),mca) ) )\n}\n\nkable(rest,digits = 5,align = 'c',caption = \"Résultats de l'approximation Monte Carlo du quantile\",\n      col.names =c(\"Nombre de simulations\",\"Approximation Monte Carlo\",\"Erreur d'approximation\"),row.names = FALSE)

Discussions et conclusions

Les méthodes d’approximation de Monte Carlo offrent une alternative pour l’approximation d’intégrales et constituent un outil essentiel de l’inférence bayésienne, en particulier lorsque l’on travaille avec des modèles sophistiqués et complexes. Comme on l’observe dans nos trois exemples, les méthodes de Monte Carlo fournissent d’excellentes approximations, mais exigent un très grand nombre de simulations pour obtenir une erreur d’approximation proche de zéro.

Références

1) Introducing Monte Carlo methods with R, Springer 2004, Christian P. Robert and George Casella.

2) Handbook of Markov Chain Monte Carlo, Chapman and Hall, Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng.

3) Introduction to mathematical Statistics, Pearson, Robert V. Hogg, Joseph W. Mckean, and Allen T. Craig.

4) Statistical Inference An Integrated Approach, Chapman and Hall, Helio S. Migon, Dani Gamerman, Francisco Louzada.