TP5 - Méthode de Monte Carlo¶

In [1]:
import numpy as np
import numpy.random as rd
import matplotlib.pyplot as plt

Exercice 1¶

  1. Utiliser la méthode de Monte-Carlo pour déterminer une valeur approchée de l'intégrale $\int_{0}^{1} e^{-t^{2}}\mathrm{~d} t.$
  2. Faire de même pour l'intégrale $\int_{0}^{+\infty} \frac{e^{-t}}{1+t^{4}} \mathrm{~d} t$ et pour la somme $ \sum\limits_{n=1}^{\infty} \frac{\sqrt{n}}{2^{n}}$.

Solution¶

  1. Si $U$ suit la loi $\mathscr{U}([0,1])$, alors le théorème de transfert donne $$\mathbb{E}\left(e^{-U^2}\right) = \int_0^1 e^{-t^2} \; \mathrm{d}t.$$ Donc si $(U_n)_n$ est une suite de v.a.r. mutuellement indépendantes qui suivent la loi $\mathscr{U}([0,1])$, on sait que $\frac 1 n \sum\limits_{i=1}^n e^{-U_i^2}$ converge en probabilité vers $\mathbb{E}\left(e^{-U^2}\right)$ qui vaut l'intégrale souhaitée.
In [2]:
#1.
n=10000
U=rd.rand(n)
print("Approximation de l'intégrale :",np.mean(np.exp(-U**2)))

Approximation de l'intégrale : 0.7449561824878091
  1. Si $X\hookrightarrow \mathscr E(1)$, alors par le théorème de transfert,on a $$\mathbb E\left(\frac 1{1+X^4}\right)\;=\; \int_0^{+\infty} \frac{1}{1+t^4}\, \mathrm e^{-t} \, \mathrm dt.$$ Si $(X_n)_n$ est une suite de v.a.r. indépendantes suivant la loi $\mathscr E(1)$, on sait que $\left(\frac 1 n \sum\limits_{i=1}^n \frac{1}{1+{X_i}^4}\right)$ converge en probabilité vers l'intégrale souhaitée.
In [3]:
#2.
n=10000
X=rd.exponential(1,n)
print("Approximation de l'intégrale :",np.mean(1/(1+X**4)))

Approximation de l'intégrale : 0.6291438933971263

Si $X\hookrightarrow \mathscr G(1/2)$, alors toujours par le théorème de transfert, on a $$\mathbb E(\sqrt{X}) \;=\; \sum_{n=1}^{+\infty} \sqrt{n}\, \mathbb P(X=n) \;=\; \sum_{n=1}^{+\infty} \frac{\sqrt{n}}{2^n}.$$ Ainsi, si $(X_n)_n$ est une suite de v.a.r. indépendantes suivant la loi $\mathscr G(1/2)$, alors $\left(\frac 1 n \sum\limits_{i=1}^n\sqrt{X_i}\right)$ converge en probabilité vers la somme souhaitée.

In [4]:
n=10000
X=rd.geometric(1/2,n)
print("Approximation de la somme :",np.mean(np.sqrt(X)))

Approximation de la somme : 1.3442954879763536
In [5]:
#Vérification
S=0; u=1
for i in range(1,10000):
    u=u/2
    S+=np.sqrt(i)*u
print(S)

1.3472537527357502

Exercice 2¶

Soit $n \geqslant 2$ un entier. On considère une urne contenant $n$ boules numérotées de 1 à $n$. On effectue des tirages successifs d'une boule avec remise. On suppose que tous les tirages sont équiprobables. On s'arrête dès que l'on obtient une boule déjà obtenue lors d'un précédent tirage. On note ${T}_{n}$ la variable aléatoire égale au nombre de tirages effectués.

On souhaite comparer graphiquement à l'aide de Python la loi de la variable aléatoire ${Y}_{n}=\frac{\mathrm{T}_{n}}{\sqrt{n}}$ avec la loi qui a pour densité

$$ f: x \longmapsto\left\{\begin{array}{cc} 0 & \text { si } x \leqslant 0 \\ x {e}^{-x^{2} / 2} & \text { si } x>0 \end{array} .\right. $$
  1. Écrire une function $\verb|simuleY(n)|$ qui simule la variable aléatoire $Y_{n}$.
  2. Comment approcher la fonction de répartition de $Y_{n}$ par la méthode de Monte-Carlo ? La représenter graphiquement sur $[0,4]$ pour $n=50$.
  3. Représenter sur $[0,4]$ la fonction de répartition de $Y$. Illustrer la convergence en superposant les représentations précédentes pour différentes valeurs de $n$.

Solution¶

In [6]:
#1.
def simuleY(n):
    A=np.zeros(n)
    t=0
    while np.max(A)<2:
        k=rd.randint(0,n)
        A[k]+=1
        t+=1
    return(t/np.sqrt(n))
In [7]:
#2. Vérification avec la méthode de Monte Carlo pour approcher F_{Y_n}(t)
def SimuleF(n,t):   # Approche F(t) pour une valeur de n donnée
    N=1000          # Grand nombre de simulations
    A=np.array([simuleY(n) for i in range(N)])
    return(np.mean(A<t))

La loi de densité $f$ a pour fonction de répartition $F$ telle que pour tout $t\in \mathbb{R}$, $$F(t) = \int_0^\infty x e^{-\frac{x^2}{2}} \, \mathrm{d}x = 1-e^{-\frac{t^2}{2}}.$$

In [8]:
# Comparaison de l'approximation de F_{Y_n} et de F
n=100
N=20
t=np.linspace(0,4,N)
y=np.zeros(N)
for i in range(N):
    y[i]=SimuleF(n,t[i])
plt.plot(t,y)
plt.plot(t,1-np.exp(-t**2/2))
plt.show()
In [9]:
# Autre possibilité : vérification avec un histogramme (cf TP précédent)
n=100
A=np.array([simuleY(n) for i in range(1000)])
plt.hist(A,100,cumulative=True, density=True)
x=np.linspace(0,4,100)
plt.plot(x,1-np.exp(-x**2/2))
plt.show()

Exercice 3¶

Que retourne le code suivant ? Justifier.

n=10**4
X=rd.rand(n)
print(4*np.mean(np.sqrt(1-X**2)))

Solution¶

In [10]:
n=10**4
X=rd.rand(n)
print(4*np.mean(np.sqrt(1-X**2)))

3.1423147979458683

Le programme calcule $\frac 4 n \sum_{i=1}^n \sqrt{1-X_i^2}$, où les v.a.r. $X_i$ suivent la loi $\mathscr{U}([0,1])$. Or d'après la loi des grands nombres, $$\frac 1 n \sum_{i=1}^n \sqrt{1-X_i^2} \text{ converge en probabilité vers } \mathbb{E}\left(\sqrt{1-X_1^2}\right) = \int_0^1 \sqrt{1-x^2}\, \mathrm{d}x = \frac \pi 4.$$ Par conséquent, le script calcule une approximation de $\pi$.

Exercice 4¶

On considère de nouveau la méthode utilisée dans l'exercice 3 pour approcher $\pi$.

  1. Quelle valeur de $n$ choisir pour que $\overline Y_n$ soit une valeur approchée de $\pi$ à la précision $\varepsilon=10^{-2}$ avec un risque $\alpha=5 \%$ ?
  2. Comment vérifier expérimentalement le résultat précédent ? Mettre cette méthode en pratique.

Solution¶

Pour simplifier les écritures, on va chercher à approcher $\frac \pi 4$ avec une précision $\varepsilon = 10^{-2}$.

Avec Bienaymé-Tchebychev¶

On cherche la plus petite valeur de $n$ pour assurer que $\mathbb{P}(|Y_n-m|\geq \varepsilon)\leq 5.10^{-2}$. Si on utilise l'inégalité de Bienaymé-Tchebychev, on cherche donc le plus petite entier $n$ tel que $$\frac{\sigma^2}{n\varepsilon^2} \leq 5. 10^{-2},~~~ \text{c'est-à-dire}~~~ n\geq \frac{\sigma^2}{5} \, 10^6 = 2\sigma^2 10^5,$$ du fait que $\varepsilon=10^{-2}$.

Pour calculer une approximation de $\frac \pi 4$, on pose $Y_i=\sqrt{1-X_i^2}$, où $X_i \hookrightarrow \mathcal{U}([0,1])$. On a donc, comme vu ci-dessus, $\mathbb E(Y_i)=\frac \pi 4$, et $$\mathbb V(Y_i) = \mathbb E(Y_i^2)-\mathbb E(Y_i)^2 = \mathbb E(1-X_i^2) - \frac{\pi^2}{16} = \frac 2 3 - \frac{\pi^2}{16},$$ car $\mathbb E(1-X_i^2)= 1-\mathbb E(X_i^2) = 1- \int_0^1 x^2 dx = 1- \frac 1 3=\frac 2 3$. On est donc ammené à poser $$n = \left\lceil 2 \sigma^2 10^5 \right\rceil,~ \text{où}~ \sigma^2 = \frac 2 3 - \frac{\pi^2}{16}.$$

Avec le théorème central limite¶

Le théorème central limite exprime que pour $n$ assez grand, la loi de la variable aléatoire centrée réduite $$\frac{\overline Y_n-m}{\frac{\sigma}{\sqrt{n}}}$$ peut être approchée par la loi $\mathcal N(0,1)$.

Comme on a $$\mathbb P\left( |Y_n-m|\geqslant \varepsilon \right) ~=~ \mathbb P\left( \left| \frac{\overline Y_n-m}{\frac{\sigma}{\sqrt{n}}} \right| \geqslant \frac{\varepsilon}{\frac{\sigma}{\sqrt{n}}}\right) ~=~ \mathbb P\left( \left| \frac{\overline Y_n-m}{\frac{\sigma}{\sqrt{n}}} \right| \geqslant \frac{\varepsilon \sqrt{n}}{\sigma}\right),$$ En faisant l'approximation donnée par le théorème central limite, on cherche donc à avoir $$\mathbb P \left(|Z| \geq \frac{\varepsilon \sqrt{n}}{\sigma}\right) \;\leqslant\; 0{,}05,$$ où $Z$ suit la loi $\mathcal N(0,1)$.

Si on choisit $z=1{,}96$, on a $\mathbb P(|Z|\geq z)\simeq 5.10^{-2}$. Ceci suggère de choisir $$\frac{\varepsilon \sqrt{n}}{\sigma} \;\geqslant \; 1{,}96,~~ \mathrm{soit}~~ $$

Pour pouvoir utiliser le résultat de convergence dans notre cas, on souhaite avoir $$\mathbb P(|Y_n-m|\geq \varepsilon) \leq \mathbb P\left(|Y_n-m|\geq z \frac{\sigma}{\sqrt{n}}\right).$$ Par conséquent, il faut avoir $z\frac{\sigma}{\sqrt{n}} \leq \varepsilon$, soit $n\geq \frac{z^2 \sigma^2}{\varepsilon^2}$.

Si on chosit $z=1{,}96$, on a $\mathbb P(|Z|\geq z)\simeq 5.10^{-2}$, donc le résultat de convergence suggère de choisir $$n=\lceil 1{,}96^2 \sigma^2 10^4\rceil,~~~ \text{où}~ \sigma^2 = \frac 2 3 - \frac{\pi^2}{16}.$$

In [11]:
#Fonction approchant pi/4 :
def approx(n):
    X=rd.rand(n)
    return(np.mean(np.sqrt(1-X**2)))
In [12]:
#Bienaymé-Tchebychev
v=2/3-(np.pi/4)**2
n=int(np.ceil(2*v*10**5))
N=10000 # nombre de simulations
A=np.array([approx(n) for i in range(N)])
E=abs(A-np.pi/4)
print("Estimation de la probabilité que l'erreur soit supérieure à 0.01 :", np.mean(E>10**(-2)))

Estimation de la probabilité que l'erreur soit supérieure à 0.01 : 0.0
In [13]:
#TCL
n=int(np.ceil(1.96**2*v*10**4))
N=10000 # nombre de simulations
A=np.array([approx(n) for i in range(N)])
E=abs(A-np.pi/4)
print("Estimation de la probabilité que l'erreur soit supérieure à 0.01 :", np.mean(E>10**(-2)))

Estimation de la probabilité que l'erreur soit supérieure à 0.01 : 0.0495