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 :

4.jpg

code :

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

In [31]:
# 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$

In [32]:
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

In [33]:
proba_geom <- function(k, p){
    return((1-p)^(k-1) * p)
}
In [37]:
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$
In [8]:
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$
In [12]:
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.

5_1.jpg

code R :

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

In [73]:
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é

In [74]:
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")

Programmons cette méthode Monte Carlo :

In [43]:
phi <- function(x, c){
    return((abs(prod(x)) < c)*1)
}
In [44]:
monte_carlo_estimate <- function(n, d, c){
    # n : taille de l'échantillon
    # d : dimension
    # c le paramètre pour la région A
    
    rn <- rmvnorm(n, mean = rep(0,d), sigma = diag(d))
    
    I_hat = 0
    
    for(i in 1:n){
        I_hat = I_hat + phi(rn[i,], c)
    }
    I_hat = I_hat / n
    
    return(I_hat)
}
In [45]:
monte_carlo_estimate(n = 100, d = 3, c = 1)
0.9
In [46]:
monte_carlo_estimate(n = 100, d = 3, c = 0.5)
0.77
In [52]:
monte_carlo_estimate(n = 100, d = 3, c = 0.0001)
0
In [53]:
monte_carlo_estimate(n = 10000,  d = 3, c = 0.0001)
0.0033
  • Regardons l'évolution de l'estimateur en fonction de $n$ pour $c = 0.5$
In [56]:
x = seq(10, 3000, 10)
y = c()
for(n in x){
    y = c(y, monte_carlo_estimate(n = n, d = 3, c = 0.5))
}
plot(x, y, xlab = "n, la taille de l'échantillon", ylab = "I_hat")

Regardons l'évolution de l'estimateur en fonction de $n$ pour $c = 0.0001$

In [58]:
x = seq(10, 10000, 10)
y = c()
for(n in x){
    y = c(y, monte_carlo_estimate(n = n, d = 3, c = 0.0001))
}
plot(x, y, xlab = "n, la taille de l'échantillon", ylab = "I_hat")

5.2.

5_2.jpg

code R :

Jettons un coup d'oeil sur la forme de la fonction $\varphi$ :

In [28]:
phi_2d <- function(x, y, c){
    return((abs(x*y) < c)*1)
}
phi_cpetit <- function(x, y){
    return(phi_2d(x, y, 0.01))
}
phi_cgrand <- function(x, y){
    return(phi_2d(x, y, 0.9))
}
In [29]:
x <- y <- seq(-2, 2, 0.1)
z <- outer(x, y, phi_cpetit)
persp(x, y, z,
      main="Représentation de phi pour c = 0.01",
      xlab = "x_1",
      ylab = "x_2",
      zlab = "X in A",
      theta = 50, phi = 30,
      col = "steelblue", shade = 0.2)
In [30]:
x <- y <- seq(-2, 2, 0.1)
z <- outer(x, y, phi_cgrand)
persp(x, y, z,
      main="Représentation de phi pour c = 0.9",
      xlab = "x_1",
      ylab = "x_2",
      zlab = "X in A",
      theta = 50, phi = 30,
      col = "steelblue", shade = 0.2)

Idée : Prendre la statistique $h(X) = \sum\limits_{k=1}^{d} X^k $ comme variable de contrôle, pour deux raisons :

  • $h(X)$ est corrélé à $\{X\in A\}$
  • il est facile de calculer l'esperance de $h(X)$ : $\mu = \mathbb{E}[h(X)] = d\cdot\mathbb{E}[X^1] = d\cdot\sqrt{\frac{2}{\pi}}$

code R :

In [82]:
library("mvtnorm")

# rappel
phi <- function(x, c){
    return((abs(prod(x)) < c)*1)
}

monte_carlo_estimate <- function(n, d, c){
    # n : taille de l'échantillon
    # d : dimension
    # c le paramètre pour la région A
    
    rn <- rmvnorm(n, mean = rep(0,d), sigma = diag(d))
    
    I_hat = 0
    
    for(i in 1:n){
        I_hat = I_hat + phi(rn[i,], c)
    }
    I_hat = I_hat / n
    
    return(I_hat)
}
In [83]:
h <- function(x){
    return(sum(abs(x)))
}

control_variates_estimate <- function(n, d, c){
    # n : taille de l'échantillon
    # d : dimension
    # c le paramètre pour la région A
    
    rn <- rmvnorm(n, mean = rep(0,d), sigma = diag(d))

    fX = rep(0, n)
    for(i in 1:n){ fX[i] = phi(rn[i,], c) }
    fXbar = mean(fX)
    
    hX = rep(0, n)
    for(i in 1:n){ hX[i] = h(rn[i,]) }
    mean_hX = d*sqrt(2/pi)
    var_hX = var(hX)
    cov_hX_fX = cov(fX, hX)
    beta = - cov_hX_fX / var_hX
    
    Z =  fXbar + beta * mean(hX - mean_hX)
    
    return(Z)
}
In [91]:
n = 500
d = 3
c = 0.01
Imc = c()
Icv = c()
for(i in 1:n){
    Imc = c(Imc, monte_carlo_estimate(n, d, c))
    Icv = c(Icv, control_variates_estimate(n, d, c))
}
In [92]:
boxplot(Imc, Icv)

5.3.

5_3.jpg

Idée : Nous pouvons générer

  • $X_i = G(U_i)$ par Box-Muller avec les variables uniformes $U_1, U_2$ et
  • $-X_i = G(1-U_i)$ avec les variables uniformes $1 - U_1, 1 - U_2$