Statistique et Probabilités [M3201]
TP 3 : Loi forte des grands nombres et théorème centrale limite
Introduction
Pour commencer, lisez la partie du cours relative à ces deux théorèmes fondamentaux.
Quelques distributions usuelles:
Loi | Notation | Esperance | Variance | Support | Densité/masse | sur \(\texttt{R}\) | |
---|---|---|---|---|---|---|---|
Binomiale | \(\mathcal B(n, p)\) | \(np\) | \(np(1-p)\) | \(k\in\{0, ..., n\}\) | \(C_n^k p^k(1-p)^{n-k}\) | \(\texttt{rbinom, dbinom, pbinom, qbinom}\) | |
Poisson | \(\mathcal P(\lambda)\) | \(\lambda\) | \(\lambda\) | \(k\in\mathbb N\) | \(\frac{e^{-\lambda}\lambda^k}{k!}\) | \(\texttt{rpois, dpois, ppois, qpois}\) | |
Normale | \(\mathcal N(\mu, \sigma^2)\) | \(\mu\) | \(\sigma^2\) | \(x\in\mathbb R\) | \(\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac 1 2\left(\frac{x- \mu}{\sigma}\right)^2\right)\) | \(\texttt{rnorm, dnorm, pnorm, qnorm}\) | |
Exponentielle | \(\mathcal E(\lambda)\) | \(\frac 1\lambda\) | \(\frac{1}{\lambda^2}\) | \(x\in[0, +\infty[\) | \(\lambda e^{-\lambda x}, x\geq 0\) | \(\texttt{rexp, dexp, pexp, qexp}\) |
Loi forte des grands nombres et Théorème Centrale Limite
Pour chacune des lois suivantes :
\(\mathcal N(0,1)\): Loi normale centrée réduite.
\(\mathcal N(4,9)\): Loi normale de moyenne 4 et de variance 9.
\(\mathcal B(5,1/4)\) : Loi binomiale de taille 5 et de paramètre 1/4 (\(\color{blue}{\texttt{?rbinom}}\))..
\(\mathcal B(1/3)\): Loi de Bernoulli de paramètre 1/3 (\(\color{blue}{\texttt{?rbinom, size=1}}\)).
\(\mathcal E(2)\): Loi exponentielle de paramètre 2 (\(\color{blue}{\texttt{?rexp}}\)).
\(\mathcal P(2)\): Loi de Poisson de paramètre 2 (\(\color{blue}{\texttt{?rpois}}\)).
Ecrire une fonction en R permettant de générer ces lois de taille \(N\) (échantillon de taille \(N\)).
Vérifiez pour chaque distribution la loi forte des grands nombres. \[\bar X_n = \dfrac{1}{n} \sum_{i=1}^n X_i \xrightarrow[n \to + \infty]{} \mathbb E[X_1] \]
On veut maintenant simuler \(M\) échantillons de taille \(N\). Ecrire une fonction permettant de le faire.
Vérifiez pour chaque distribution le théorème centrale limite. \[\bar X_n=\frac 1 n\sum_{i=1}^nX_i \xrightarrow[n\to +\infty]{} \mathcal N(\mu,\sigma^2/n)\]
Méthodes de Monte Carlo
Ces méthodes visent à calculer une intégrale par des simulations de variables aléatoires : \[ \rm I = \int g(x)f(x)\,dx \] où \(f\) est la fonction densité d’une variable aléaoire \(X\). Ainsi, on a \[ \rm I=\mathbb E[g(X)] \] D’après la loi forte des grands nombres, il suffit de simuler un échantillon i.i.d \(X_1,\ldots,X_n\), de loi à densité \(f\) et de calculer \[ \frac 1n \sum_{i=1}^n g(X_i), \] cette quantité tendant vers \(I=\mathbb E[g(X)]\) quand \(n\to +\infty\).
exemple :
On cherche à calculer \[ \rm I=\int_{0}^1 \frac{\ln(1+x)}{\sqrt{1-x^2}}\,dx \]
En R, cela donne
n<-1000000
ech<-runif(n)
mean(log(1+ech)/sqrt(1-ech^2))
## [1] 0.7424894
integrate(function(x) log(1+x)/sqrt(1-x^2),0,1)
## 0.7431381 with absolute error < 2.2e-06
On peut également utiliser une autre loi que la loi uniforme : par exemple la loi beta, notée \(\mathcal B(a,b)\). La densité de cette loi est \[ \forall x\in [0,1],\, f(x) = \frac{1}{B(a,b)} x^{a-1}(1-x)^{b-1} \] où \(B(a,b)=\dfrac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\). Pour rappel, la fonction gamma, notée \(\Gamma(x)\) vaut \[ \Gamma(x) = \int_0^{+\infty} t^{x-1}\exp(-t)\,dt \]
On peut donc écrire \(I\) sous la forme \[ \rm I = \int_0^1 \frac{g(x)}{f(x)}f(x)\,dx \] où \(f\) est la fonction densité de la loi Beta.
En R, cela donne
n<-1000000
a<-5 # 1er paramètre de la loi Beta
b<-1 # 2ème paramètre de la loi Beta
ech<-rbeta(n,a,b)
g<-function(x) log(1+x)/sqrt(1-x^2)
mean(g(ech)/dbeta(ech,a,b))
## [1] 0.7403863
integrate(g,0,1)
## 0.7431381 with absolute error < 2.2e-06
TP
Approximation de \(\pi\)
Soit \(f\) la fonction définie sur \([0,1]\) par \[ \forall x\in [0,1],\, f(x) = \sqrt{1-x^2} \] Faire une procédure permettant de calculer \[ I=\int_0^1 f(x)\,dx \]
En déduire une approximation de \(\pi\).
N<- 10000 # nombre de simulations
x<-runif(N) # échantillon iid de loi Uniforme sur [0,1]
res<-mean(sqrt(1-x^2))
4*res # Approximation de Pi
## [1] 3.130869
Une autre façon d’approcher \(\pi\) est de simuler des points dans le carré unité et de compter le nombre de points qui sont à une distance de l’origine plus petite que \(1\).
N<- 1000000 # nombre de simulations
x<-runif(N) # échantillon iid de loi Uniforme sur [0,1]
y<-runif(N)
res<-mean(sqrt(x^2+y^2) <= 1)
4 * res # Approximation de Pi
## [1] 3.144816
Calculer également par Monte Carlo l’intégrale suivante : \[ \int_0^1 \dfrac{1}{1+x^2}\,dx \]
Analytiquement avec \(f(x) =\arctan x, f'(x) = \frac{1}{1+x^2}\) \[ \int_0^1 \dfrac{1}{1+x^2}\,dx = \left[\arctan(x)\right]_0^1 = \pi/4 \]
Calcul d’intégrales
Exercice 1
Nous considérons une variable aléatoire \(X\) de loi \(\Gamma(a,b)\) de densité de probabilité : \[ f(x) = \frac{b^a}{\Gamma(a)}x^{a-1}\exp(-bx) \] Dans la suite, on prend \(a= 4\) et \(b= 1\).
Simuler un échantillon de 1000 réalisations de \(X\) (\(\texttt{?rgamma}\)).
Donner, par des simulations, des approximations de l’espérance et de la variance de \(X\).
Calculer l’intégrale suivante par Monte Carlo \[ \int_0^{+\infty} \cos(2x)x^{7/3}\exp(-3.2x)\,dx \]
Exercice 2
On cherche à calculer \[ \rm I = \int_0^{+\infty} \exp(-x^2/2)\,dx \]
Montrer grâce à la loi normale et par le calcul que \(\rm I=\sqrt{\dfrac{\pi}{2}}\).
En effectuant un changement de variables \(y=x^2/2\) dans l’intégrale, prouvez que \[ \rm I= \int_0^{+\infty} \frac{1}{\sqrt{2y}}\exp(-y)\,dy. \]
- Montrer que cette intégrale peut être vue comme l’espérance d’une fonction d’une variable aléatoire \(X\) qui suit une loi exponentielle. Donner le paramètre de cette loi exponentielle.
- Mettre en oeuvre un programme permettant de calculer cette intégrale par la Monte Carlo. Comparer avec la fonction integrate.
Exercice 3 : En dimension 2
On s’intéresse à l’intégrale \[ \rm I = \iint_{\mathcal D} \sqrt{r^2-x^2-y^2}\,dxdy \] où \(r>0\) et \(\mathcal D=\{(x,y)\in \mathbb R^2,\, x^2+y^2\leq r^2\}\).
Mettre en oeuvre un programme permettant de calculer cette intégrale double.