tp3.utf8.md

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\)).

# N(0,1) Normale
sample1 <- function(n){rnorm(n, mean=0, sd=1)}

# N(4,9) Normale
sample2 <- function(n){rnorm(n, mean=4, sd=3)}

# B(1/3) Bernoulli:
sample3 <- function(n){rbinom(n, size=1, prob=1/3)}

# B(5, 1/4) Binomiale:
sample4 <- function(n){rbinom(n, size=5, prob=1/4)}

# E(2) Exponentielle:
sample5 <- function(n){rexp(n, rate=2)}

# P(2) Poisson:
sample6 <- function(n){rpois(n, lambda=2)}

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] \]

N = 10000
ech <- sample1(N)   #E[X1] = 0
(barX <- mean(ech))

## [1] -0.008348434

ech <- sample2(N)  #E[X1] = 4
(barX <- mean(ech))

## [1] 4.048126

ech <- sample3(N)  #E[X1] = 1/3
(barX <- mean(ech))

## [1] 0.3327

ech <- sample4(N)  #E[X1] = 5/4
(barX <- mean(ech))

## [1] 1.2455

ech <- sample5(N)  #E[X1] = 1/2
(barX <- mean(ech))

## [1] 0.5059285

ech <- sample6(N)  #E[X1] = 2
(barX <- mean(ech))

## [1] 1.9927

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.7397828

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.74263

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.145652

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.141592

Calculer également par Monte Carlo l’intégrale suivante : \[ \int_0^1 \dfrac{1}{1+x^2}\,dx \]

N<- 10000  # nombre de simulations
x<-runif(N) # échantillon iid de loi Uniforme sur [0,1]
res<-mean(1/(1+x^2))
res

## [1] 0.788333

f <-function(x) 1/(1+x^2)
integrate(f,0,1)

## 0.7853982 with absolute error < 8.7e-15

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 \]

x <- rgamma(1000, shape=4, rate=1)
# Loi forte des grands nombres:
(meanGamma <- mean(x))  # a/b

## [1] 3.948939

(var <- var(x))  # a/b^2

## [1] 3.9726

N <- 100000  # nombre de simulations
x <-rgamma(N, shape=4, rate=1) 
g <-function(x) cos(2*x)*x^(7/3)*exp(-3.2*x)
res <-mean(g(x) / dgamma(x, shape=4, rate=1))
res

## [1] -0.009360291

integrate(g, 0, Inf)

## -0.009534073 with absolute error < 5.8e-08

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}}\). \[ \begin{align*} 1 & = \int_{\infty}^{+\infty} \frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}2) dx\tag{Densité de la loi normale centrée réduite}\\ & = \int_{-\infty}^0 \frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}2) + \int_0^{+\infty}\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}2)\\ & = \int_{+\infty}^0\frac{1}{\sqrt{2\pi}}\exp(-\frac{(-x)^2}2)(-dx) + \int_0^{+\infty}\frac{1}{\sqrt{2\pi}}\exp(-\frac12 x^2) \tag{changement de variable à -x}\\ & = 2 \int_0^{+\infty}\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}2) \\ & = \frac{2}{\sqrt{2\pi}} \int_0^{+\infty}\exp(-\frac{x^2}2) \\ & = \sqrt{\frac2\pi} \rm I \end{align*} \] D’où \(\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. \] Avec le changement de variable \(y =x^2/2\) sur \(\mathbb R^+\) on a \(x = \sqrt{2y}\) et \(dy=x\, dx\).

\[ \rm I =\int_0^{+\infty} \exp(-x^2/2)\,dx= \int_0^{+\infty} \exp(-y)\,dy/x = \int_0^{+\infty} \frac{1}{\sqrt{2y}}\exp(-y)\,dy. \] 3. 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.

Soit \(X\) une v.a qui suit la loit exponentielle de parametre 1 i.e \(X\sim\mathcal E(1)\). On a \(\,f(x) =\begin{cases}\exp(-x),\,\text{ si } x\geq 0\\ 0,\,\text{ sinon}\end{cases}\)

\[ I= \int_0^{+\infty} \frac{1}{\sqrt{2x}}\exp(-x)\,dx = \int_{-\infty}^{+\infty} g(x)f(x)dx = E[g(X)],\,\,\text{avec } g(x)=\frac{1}{\sqrt{2x}}. \]

Mettre en oeuvre un programme permettant de calculer cette intégrale par la Monte Carlo. Comparer avec la fonction integrate.

# Monte Carlo:
N<- 10000  # nombre de simulations
x <- rexp(N, rate=1) 
g <- function(x) 1/sqrt(2*x)
res<- mean(g(x))
res

## [1] 1.230544

# Avec integrate
f <- function(x) 1/sqrt(2*x) * exp(-x)
integrate(f,0,Inf)

## 1.253314 with absolute error < 1.4e-05

# La valeur analytique:
sqrt(pi / 2)

## [1] 1.253314

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\}\).

Rappelons que l’équation d’une sphere de rayon \(r\) est:

\[\mathcal S=\{(x,y, z)\in \mathbb R^3,\, x^2+y^2+z^2 = r^2\}\] D’où:

\[ \rm I = \iint_{\mathcal D} \sqrt{r^2-(x^2 + y^2)} \,dxdy = \iint_{\mathcal D} z_+(x,y) \,dxdy = \int_{S, z\geq 0} z dz \] Ainsi, \(\rm I\) est le volume de la moitié d’une sphère (la partie à z positifs):

Soient \(X\) la v.a des abscisses et \(Y\) celle des ordonnées sur \(\mathcal D\).

On échantillonne \(X_\text{square}\) et \(Y_\text{square}\) de loi uniforme sur \([-r, -r]\) et on n’accepte que les points de \(\mathcal D\). Soit \(\alpha\) le taux d’acceptation.

On considere \(f\) la densite de la loi jointe \(f_{X,Y}(x,y) = \dfrac{1}{2r}\dfrac{1}{2r} = \dfrac{1}{4r^2}\)

La densité de \(X,Y \in \mathcal D\) peut etre approximée avec \(\frac 1\alpha.\, f.\, 1_D\) (\(1_D\) est l’indicatrice de \(\mathcal D\), elle vaut 1 si \((x,y)\in \mathcal D\) et 0 sinon).

\[\begin{align*} \iint_{\mathcal D} \sqrt{r^2-x^2-y^2}\,dxdy &= \mathbb E_{X, Y}\left[\frac{\sqrt{r^2-x^2-y^2}}{f(x)g(y) / \alpha}\right] \\ & = \mathbb E_{X, Y}\left[4\alpha r^2\sqrt{r^2-x^2-y^2}\right] \\ &= 4\alpha r^2\mathbb E_{X, Y}\left[\sqrt{r^2-x^2-y^2}\right] \end{align*}\]

Mettre en oeuvre un programme permettant de calculer cette intégrale double.

N <- 1e4  # nombre de simulations
radius <- 3
x <- runif(N, min=-radius, max=radius) 
y <- runif(N, min=-radius, max=radius)
xys = cbind(x, y)
disk = xys[x^2+y^2 <= radius^2, ]

par(pty="s")
plot(x, y, cex=.3, pch=19, main='Echantillonnement avec la methode de rejet')
points(disk, col='blue', pch=19, cex=.3)

Maha Elbayad

Statistique et Probabilités [M3201]

TP 3 : Loi forte des grands nombres et théorème centrale limite