TD 3 - 20/03/20 - Simulation et Monte Carlo

• Exercice 4 : Loi géométrique
• Exercice 5 : Variables de contrôle, variables antithétiques, QMC

Exercice 4 : Loi géométrique

• Proposer un algorithme pour simuler une loi géométrique
• Le mettre en oeuvre
• Proposer une méthode pour vérifier les résultats

soit $X\sim \mathcal{G}(p)$ qui suit une loi géométrique de paramètre $p$.

$$ \forall k \in \mathbb{N}^*, \ \ P(X = k) = (p-1)^{k-1}p = q^{k-1}p $$

C'est la loi d'attente d'un premier succès dans une suite d'épreuves répétées indépendantes ayant chacun la probabilité de succès $p$

Par conséquent si

$$(U_i)_{i\geq 1} \sim^{iid} \mathcal{U}(0,1)$$

méthode 1 : la variable

$$ Y := \min\left\{k\in\mathbb{N}^*; U_k<p\right\} \sim \mathcal{G}(p) $$

méthode 2 :

code :

proposition de générateur de la loi géomatrique

# méthode 1
rand_geom <- function(n, p){
    rg = rep(0,n)
    for(i in 1:n){
        ru <- runif(1)
        k <- 1
        while(ru > p){
            k <- k+1
            ru <- runif(1)
        }
        rg[i] <- k
    }
    return(rg)
}

# méthode 2
rand_geom2 <- function(n, p){
  ru=runif(n)
  X=log(ru)/log(1-p)
  Y=floor(X)+1
  return(Y)
}

générons un échantillon de taille $n = 10000$

n = 10000
p = 0.3
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)

affichons son histogramme ainsi que la vraie distribution de la loi géométrique

proba_geom <- function(k, p){
    return((1-p)^(k-1) * p)
}

hist(sample_geom, col = "steelblue", proba = TRUE,
     breaks = seq(0,max(sample_geom),1),
     main = "histogramme et pdf de la loi géométrique",
     xlab = "k",
     ylab = "probabilité")
hist(sample_geom2, prob = TRUE, col = rgb(1,0,0,0.4), add=T, breaks = seq(0,max(sample_geom2),1))
k = seq(1,20)
lines(k, proba_geom(k, p), col = "orange", lwd = 3, type = "o")
legend("topright", c("Notre générateur", "Notre générateur", "Vraie proba"), col=c("steelblue", "red", "orange"), lwd=6)

Regardons les deux extrêmes : $p = 0.1$ et $p = 0.9$

• $p = 0.1$

n = 10000
p = 0.1
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)
hist(sample_geom, col = "steelblue", proba = TRUE,
     breaks = seq(0,max(sample_geom),1),
     main = "histogramme et pdf de la loi géométrique",
     xlab = "k",
     ylab = "probabilité")
hist(sample_geom2, prob = TRUE, col = rgb(1,0,0,0.4), add=T, breaks = seq(0,max(sample_geom2),1))
k = seq(1,60)
lines(k, proba_geom(k, p), col = "orange", lwd = 3, type = "o")
legend("topright", c("Notre générateur", "Notre générateur", "Vraie proba"), col=c("steelblue", "red", "orange"), lwd=6)

• $p = 0.9$

n = 80000
p = 0.9
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)
hist(sample_geom, col = "steelblue", proba = TRUE,
     breaks = seq(0,max(sample_geom),1),
     main = "histogramme et pdf de la loi géométrique",
     xlab = "k",
     ylab = "probabilité")
hist(sample_geom2, prob = TRUE, col = rgb(1,0,0,0.4), add=T, breaks = seq(0,max(sample_geom2),1))
k = seq(1,5)
lines(k, proba_geom(k, p), col = "orange", lwd = 3, type = "o")
legend("topright", c("Notre générateur", "Notre générateur", "Vraie proba"), col=c("steelblue", "red", "orange"), lwd=6)

Exercice 5 : Control variates, variables antithétiques, QMC

5.1.

code R :

Nous allons utiliser la librarie "mvtnorm" pour simuler un vecteur gaussien de dimension $d\geq 2$

library("mvtnorm")
library("scatterplot3d")
library("latex2exp")
d = 3 # d >= 2
n = 5000
mu = rep(0,d)
Sigma = diag(d)

rn <- rmvnorm(n, mean = mu, sigma = Sigma)

Jettons un coup d'oeil à l'échantillon généré

scatterplot3d(rn[,1], rn[,2], rn[,3],
              main="3D Scatter Plot",
              xlab = TeX("$x_1$"),
              ylab = TeX("$x_2$"),
              zlab = TeX("$x_3$"), angle = 45, color = "steelblue")