Illustration de la méthode de Monte-Carlo et de l'échantillonnage préférentiel¶
I) Méthode de Monte-Carlo¶
1) Loi des grands nombres¶
La méthode de Monte-Carlo repose sur la loi forte des grands nombres. L'idée est la suivante. On se donne une expérience aléatoire $X$, pour laquelle $\mathbb{E}[|X|]<+\infty$, qu'on va répéter de manière identique et indépendamment des autres. Ces répétitions nous donnent une suite d'observations $(X_n)_n$ qui sont des v.a i.i.d ayant pour loi celle de $X$ (un échantillon de $X$). La méthode de Monte-Carlo étudie le comportement moyen de ces observations.
Pour tout $N\in \mathbb{N}^\ast$, on note
$$ \bar X_N = \frac{1}{N}\sum_{j=1}^N X_j\qquad\text{et}\qquad m = \mathbb{E}[X]. $$
Alors, on a
$$ \mathbb{P}\left( \lim_{N \to \infty} \bar X_N = m \right) = 1. $$
Illustrons ce résultat numériquement, en illustrant l'évolution de $\bar X_N$ en fonction de $N$. Pour cela on va faire appelle à la fonction cumsum, de la librairie numpy, qui calcule les sommes cumulées des éléments d'un vecteur. Par exemple, si x=[a,b,c], la fonction numpy.cumsum(x) nous donnera un vecteur composé de [a, a+b, a+b+c]. On va aussi utiliser numpy.arange(u, v, d) qui donne un tableau numpy représentant une discrétisation de l'intervalle $[$ u,v $]$ avec un pas de d. Par défaut d$=1$.
# importons les libraires dont nous aurons besoin
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
On définit une loi pour $X$ (on pourra tester différentes lois.)
X = st.norm(0,1) # loi normale centrée réduite
#X = st.expon() # loi exponentielle de paramètre 1
On génère des nombres aléatoires indépendants de loi $X$, et on calcul $\bar X_N$ pour différent $N$.
N = 1000 # taille de l'échantillon
Xn = X.rvs(N) # génère un échantillon de X de taille N
barXN = np.cumsum(Xn) / np.arange(1,N+1) # calcul l'évolution de \bar X_N
On obtient le graphique suivant.
plt.plot(np.arange(1,N+1), barXN, label="$\overline{X}_N$") # dessine l'évolution de \bar X_N
plt.hlines(X.mean(), 0, N+1, color='orange', label="$\mathbb{E}[X]$") # représente la vraie valeur de E[X]
plt.legend() # affiche la légende sur le graphique
plt.xlabel(r'$N$', fontsize =15) # affiche une légende pour l'axe des abscisses
plt.show()
Ce graphique illustre la convergence de $\bar X_N$ vers $\mathbb{E}[X]$. On peut répéter cette expérience un certain nombre de fois pour illustrer le caractère "prèsque sûr" de la convergence.
cmap = plt.get_cmap("tab10") # définie des couleurs
K = 100 # nombre de répétition du calcul de \bar X_N
N = 1000 # taille de l'échantillon
Xn = X.rvs((K,N)) # matrice dont chaque ligne représente un échantillon diférent
barXN = np.cumsum(Xn, axis = 1) / np.arange(1,N+1) # calcul de l'évolution de \bar X_N pour chaque échantillon
for j in range(K-1): # répétitions du calcul de \bar X_N
plt.plot(np.arange(1,N+1), barXN[j], color=cmap(0), alpha = 0.3) # illustration de \bar X_N pour le jème échantillon
plt.plot(np.arange(1,N+1), barXN[K-1], color=cmap(0), label="$\overline{X}_N$") # illustration pour le dernier et faire la légende
plt.hlines(X.mean(), 0, N+1, color='orange', lw=3, label="$\mathbb{E}[X]$") # représente la vraie valeur de E[X]
plt.legend()
plt.xlabel(r'$N$', fontsize =15)
plt.show()
On peut remarquer que si $X$ suit une loi de Cauchy, on a $\mathbb{E}[|X|]=+\infty$, $X$ n'admet pas d'espérance, et la convergence n'a pas lieu.
X = st.cauchy(1)
K = 5
Xn = X.rvs((K,N)) # matrice dont chaque ligne représente un échantillon différent
barXN = np.cumsum(Xn, axis = 1) / np.arange(1,N+1) # calcul de l'évolution de \bar X_N pour chaque échantillon
for j in range(K): # N répétitions du calcul de \bar X_N
plt.plot(np.arange(1,N+1), barXN[j], color=cmap(0), alpha = 0.3) # illustration de \bar X_N pour le jème échantillon
plt.xlabel(r'$N$', fontsize =15)
plt.show()
2) Théorème centrale limite¶
Sous condition que $\mathbb{E}[X^2]<+\infty$, le théorème central limite nous guaranti que, pour $N$ assez grand, la v.a
$$ Z_N = \frac{\sqrt{N}}{\sigma}(\bar X_N - m), $$
où $m=\mathbb{E}[X]$ et $\sigma^2=Var(X)$, suit à peu près une loi $\mathcal{N}(0,1)$. Illustrons numériquement ce résultat avec des réalisons de $Z_N$ obtenues à partir de celles de $\bar X_N$. Une fois que nous aurons les réalisons de $Z_N$, nous comparerons l'histogramme normalisé de ces observations avec la densité de la loi $\mathcal{N}(0,1)$.
X = st.uniform(0,1) # loi normale centrée réduite
#X = st.expon() # loi exponentielle de paramètre 1
m, sigma = X.mean(), X.std()
K = 10000 # nombre de répétition du calcul de \bar X_N
N = 1000 # taille de l'échantillon
Xn = X.rvs((K,N)) # matrice dont chaque ligne représente un échantillon différent
barXN = np.mean(Xn, axis = 1) # calcul de \bar X_N pour chaque échantillon
ZN = np.sqrt(N)*(barXN - m)/sigma # calcul de ZN
x = np.linspace(-5, 5, 100)
plt.hist(ZN, bins=x, density=True, label="histogramme")
plt.plot(x, st.norm(0,1).pdf(x), lw=3, label="densité de la loi N(0,1)")
plt.legend()
plt.show()
On a donc une trés bonne adéquation de la loi de $Z_N$ avec celle de la loi $\mathcal{N}(0,1)$.
3) Approximation de la fonction de répartition et de la densité.¶
Pour la fonction de répartition, on souhaite évaluer numériquement la probabilité $\mathbb{P}(X\leq x)$, où $x$ parcours un intervalle de $\mathbb{R}$.
X = st.expon() # on définit une loi pour $X$
x = np.linspace(-1, 5, 100) # on définit une discrétisation de l'intervalle [-1,5]
N = 1000 # taille de l'échantillon
Xn = X.rvs(N) # génère l'échantillon
FN = np.array([np.mean( (Xn <= u)) for u in x]) # calcul une approximation de la fonction de réparition aux points de x
plt.plot(x, FN, label="aproximation")
plt.plot(x, X.cdf(x), label="fct de répartition de X")
plt.legend()
plt.xlabel(r'$x$', fontsize =15)
plt.show()
On observe une bonne adéquation entre la fonction de répartition théorique et celle évaluée numériquement à partir de l'échantillon. Pour la densité, on procède de la même manière.
h = 0.05 # valeurs de h
N = 10000 # taille de l'échantillon
Xn = X.rvs(N)
fNh = np.array([np.mean( (u-h <= Xn)*(Xn <= u+h)) for u in x])/(2*h)
# calcul une approximation de la densité aux points de x
plt.plot(x, fNh, label="aproximation")
plt.plot(x, X.pdf(x), label="fct de répartition de X")
plt.legend()
plt.xlabel(r'$x$', fontsize =15)
plt.show()
On observe encore une bonne adéquation entre la densité théorique et celle évaluée numériquement à partir de l'échantillon. Par contre, pour avoir une bonne approximation, la taille de l'échantillon a besoin d'être plus grand que pour la fonction de répartition.
II) Réduction de variance échantillonnage préférentiel¶
Illustrons numériquement l'exemple de la section 5.1.2 du cours. On rappelle qu'on cherchait à évaluer numériquement $I=\mathbb{E}[e^{X}]$ avec $X\sim \mathcal{N}(0,1)$ à l'aide d'un échantillon de $X$ et de
$$ \bar X_N = \frac{1}{N}\sum_{j=1}^N e^{X_j}. $$
On a vu en section 5.1.2 que la convergence vers $I$ devrait être plus rapide à l'aide de
$$ \bar Y_N = \frac{1}{N}\sum_{j=1}^N e^{Y_j(1-\mu) + \mu^2/2}, $$
où $(Y_n)_n$ est un échantillon de la loi $\mathcal{N}(\mu, 1)$ avec $\mu$ proche de $1$.
# définition des lois
X = st.norm(0,1)
mu = 0.8
Y = st.norm(mu,1)
N = 1000 # taille des échantillons
Xn = X.rvs(N)
Yn = Y.rvs(N)
# calcul des moyenne empirique
barXN = np.cumsum(np.exp(Xn)) / np.arange(1,N+1)
barYN = np.cumsum(np.exp(Yn*(1-mu) + mu**2/2)) / np.arange(1,N+1)
# illustration de l'évolution des moyennes empiriques
plt.plot(barXN, lw=2, label = 'N(0,1)')
plt.plot(barYN, lw=2, label = 'N(%s,1)' % mu)
plt.hlines(np.e**0.5, 1 , N, color = 'k',linestyle = 'dashed', label = 'valeur de I')
plt.xlabel(r'$N$', fontsize =15)
plt.title("Illustration de l'accélaration de convergence par IS")
plt.legend(loc = 'best')
plt.show()
On voit sur ce graphique qu'en modifiant la loi sous-jacente dans la méthode de Monte-Carlo, on peut obtenir une convergence plus rapide.