Enchantillonage préférentiel¶
Omar El Mellouki & Bertrand PATUREL
Professeur : Pierre GLOAGUEN
Cours ENSTA : https://papayoun.github.io/courses/2020_monte_carlo/
Dans ce problème nous étudions la probabilité qu'un joueur gagne le jackpot lors d'un jeu de hasard à espérance $\mu$ négative. Ce joueur s'arrete si il est ruiné (noté $ r<0 $) ou bien si la somme de ses gains dépasse le jackpot (noté $j > 0$). Si l'on note $X_n$ la variable aléatoire de ce que le joueur remporte à chaque fois qu'il joue, alors $S_n=\sum_{k=1}^nX_k$ est la somme de ses gains. Nous notons $\tau$ le temps d'arret (i.e. le dernier coup du joueur) qui vérifie $S_{\tau}\leq r ou S_{\tau}\geq j$.
Introduction :¶
Nous allons simuler une marche aléatoire à n pas avec condition d'arrêt dans le but d'étudier le problème de la ruine du joueur dans un jeu à espérance négative qui décide de s'arrêter quand il est ruiné ou quand il gagne le jackpot. Pour cela l'on va simuler la réalisations de variables aléatoires i.i.d gaussienne à moyenne négative et de variance 1, puis l'on va calculer la somme de ces valeurs propres tant que le joueur n'est pas ruiné ou qu'il n'a pas gagné le jackpot.
Calcul de l'estimateur de Monte Carlo de la probabilité que le joueur gagne.¶
On note $p∗=\mathbb{P}(S_{\tau}≥j)$, avec $\tau$ représentant le temps d'arrêt et j le jackpot. $p∗$ est donc la probabilité de gagner le jackpot.
Afin de calculer la probabilité par la méthode de Monte-Carlo, on l'exprime sous la forme de l'espérance d'une indicatrice.
On a alors $p∗=\mathbb{E}(\mathbb{1}_{S_\tau \geq j})$
Un estimateur de Monte-Carlo serait $\hat{p}=\frac{1}{M}\sum_{k=1}^{M}\mathbb{1}_{S_\tau(k) \geq j}$
Implémentation de l'estimateur de Monte Carlo en R.¶
mu = 1
r = -50
j = 5
M = 100000
#Simulation d'une trajectoire
S_t = 0
while((S_t<j) & (S_t>r)){
X = rnorm(1,-mu,1)
S_t = S_t + X
}
# La fonction indicatrice
indicatrice = function(vec,A)#On créée la fonction indicatrice ou A est le palier à dépasser
{
n=length(vec)
value = rep(x=0,n)
for(i in 1:n)
{
if(vec[i]>=A){value[i] = 1}
}
return(value)
}
#Le vecteur des phi
phi = c(1:M)
for(i in 1:M){
S_t = 0
while((S_t<j) & (S_t>r)){
X = rnorm(1,-mu,1)
S_t = S_t + X
}
phi[i] = indicatrice(S_t,j)
}
#L'estimateur de Monte Carlo
p_star = mean(phi)
p_star
I_hat <- mean(phi) # Estimateur de Monte Carlo de la proba
gamma2_hat <- var(phi) # Variance de Monte Carlo de la proba
z_975 <- qnorm(0.975, 0, 1) # Quantile à 0.975 de la loi normale centrée réduite
IC = I_hat + z_975 * sqrt(gamma2_hat / M) * c(-1, 1) #Intervalle de confiance
### On affiche finalement l'intervalle de confiance à 95 %
library(ggplot2)
library(dplyr)
my_samples = tibble(index = 1:M, sample = phi)
my_samples %>%
mutate(Estimate = cumsum(sample) / index,
IC_95_inf = Estimate - z_975 * sqrt(gamma2_hat / index),
IC_95_sup = Estimate + z_975 * sqrt(gamma2_hat / index)) %>%
ggplot(aes(x = index)) +
geom_ribbon(aes(ymin = IC_95_inf, ymax = IC_95_sup),
fill = "lightblue") +
geom_line(aes(y = Estimate)) +
labs(x = "M", title = "Estimation Monte Carlo de I",
y = expression(hat(I)[M])) +
theme_bw()
p_star
#On calcule la variance de l'estimateur :
gamma2_hat
Calcul de la densité du vecteur $S_{1:\tau}$.¶
Le vecteur $S_{1:\tau}$ correspond au vecteur $(S_i)=(\sum_{k=1}^iX_k)=AX$ où $X$ est le vecteur $(X_i)$ et $A$ est la matrice triangulaire inférieure $(\mathbb{1}_{i \leq j})$. Le vecteur $X$ suit la loi normale $N(-\mu,\Sigma)$ où $\Sigma$ est l'identité puisque $\Sigma = (Cov[X_i,X_j])$ avec $X$ vecteur aléatoire construit à partir des variables $(X_i)_{ 1\leq i \leq \tau}$ i.i.d. selon $N(-\mu,1)$, donc $S_{1:\tau}=AX$ suit la loi $N(-A \mu ,A \Sigma A^{\top})$.
La densité du vecteur $S_{1:\tau}$ est donc de la forme : $f(S_{1:\tau})=\frac{1}{(2\pi)^{\frac{\tau}{2}}\sqrt{det(AA^{\top})}}exp(-\frac{(S+A\mu)^{\top}(AA^{\top})^{-1}(S+A\mu)}{2})=\frac{1}{(2\pi)^{\frac{\tau}{2}}\sqrt{det(AA^{\top})}}exp(-\frac{(X+\mu)^{\top}A^{\top}(AA^{\top})^{-1}A(X+\mu)}{2})=\frac{1}{(2\pi)^{\frac{\tau}{2}}\sqrt{det(AA^{\top})}}exp(-\frac{(X+\mu)^{\top}(X+\mu)}{2})$.
- $det(AA^{\top})=det(A)^2=1$
$f(S_{1:\tau})=\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(x_k+\mu)^2}{2})$
Puisque $x_i=s_i-s_{i-1}$ où $s_0=0$, on a donc : $f(S_{1:\tau})=\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(s_k-s_{k-1}+\mu)^2}{2})=\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{s_1^2+(s_2-s_1)^2+(s_3-s_2)^2+..+(s_{k}-s_{k-1})^2+(s_{k+1}-s_{k})^2+..+(s_{\tau}-s_{\tau-1})^2+\mu(2s_{\tau}+\tau \mu)}{2})$
Reduction de la variance de notre estimation de la probabilité de remporter le jackpot.¶
On a en réutilisant ce que l'on a fait ci dessus (on change $-\mu$ en $\mu$) :
$g(S_{1:\tau})=\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(x_k-\mu)^2}{2})=\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(s_k-s_{k-1}-\mu)^2}{2})$
$\frac{f(S_{1:\tau})}{g(S_{1:\tau})}=exp(\frac{-\sum_{k=1}^{\tau}(s_k-s_{k-1}+\mu)^2 + \sum_{l=1}^{\tau}(s_l-s_{l-1}-\mu)^2}{2})=exp(-2\mu s_{\tau})$
$p*=\mathbb{P}(S_{\tau}\geq j)=\mathbb{P}((S_{1:\tau})_{\tau} \geq j)=\displaystyle \int_{-\infty}^{+\infty}\displaystyle \int_{j}^{+\infty}f(S_{1:\tau}) \, \mathrm{d}s_{\tau}\, \mathrm{d}S_{1:\tau-1}=\displaystyle \int_{-\infty}^{+\infty}\displaystyle \int_{j}^{+\infty}\frac{f(S_{1:\tau})}{g(S_{1:\tau})}g(S_{1:\tau}) \, \mathrm{d}s_{\tau}\, \mathrm{d}S_{1:\tau-1}=\displaystyle \int_{-\infty}^{+\infty}\displaystyle \int_{j}^{+\infty}exp(-2\mu s_{\tau})\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(s_k-s_{k-1}-\mu)^2}{2}) \, \mathrm{d}s_{\tau}\, \mathrm{d}S_{1:\tau-1}=\displaystyle \int_{-\infty}^{+\infty}\mathbb{1}_{S_{\tau}\geq j}exp(-2\mu s_{\tau})\frac{1}{(2\pi)^{\frac{\tau}{2}}}exp(-\frac{\sum_{k=1}^{\tau}(s_k-s_{k-1}-\mu)^2}{2}) \, \mathrm{d}S_{1:\tau}$
Ainsi un estimateur de Monte Carlo d'effort M de $p*$ est :
$\hat{p}_M^{IS}=\frac{1}{M}\sum_{k=1}^{M}\mathbb{1}_{(S_{\tau(k)}\geq j)}exp(-2\mu s_{\tau(k)})\leq\frac{1}{M}\sum_{k=1}^{M}exp(-2\mu j)=exp(-2\mu j)$ où $S_{1:\tau}$ suit la loi $N(A \mu ,A \Sigma A^{\top})$, $S_{\tau}$ la loi $N(\tau \mu , \tau)$ et $X_{1:\tau}$ suit la loi $N(\mu ,1)$.
Or $\hat{p}_M^{IS}$ tend vers $p*$. Donc $p*\leq exp(-2\mu j)$.
Implémentation en R.¶
set.seed(123) #Onfixe le germe pour avoir le même résultat
#Dans la suit M est l'effort de Monte Carlo, mu la moyenne mu, j et r sont définit comme dans le TD
S_tau=function(M,mu,j,r)#On créée la fonction S_tau
{
S_tau=1:M
for (i in 1:M)
{
S=0
while (S>r & S<j)
{
X=rnorm(n=1,mean=mu,sd=1)
S=S+X
}
S_tau[i]=S
}
return(S_tau)
}
IS=function(M,mu,j,r)#Voici notre fonction estimateur
{
S=S_tau(M,mu,j,r)
IS=mean(S*indicatrice(S,j)*exp(-2*mu*S))
return(IS)
}
IS(10000,1,5,-50)#7.7e-05
### On affiche finalement l'intervalle de confiance à 95 %
library(ggplot2)
library(dplyr)
#Calcul de la variance
S=S_tau(10000,1,5,-50)
gamma2_hat <- var(S*indicatrice(S,5)*exp(-2*1*S)) # Variance de Monte Carlo de la proba
z_975 <- qnorm(0.975, 0, 1) # Quantile à 0.975 de la loi normale centrée réduite
my_samples = tibble(index = 1:5000, sample = IS(5000,1,5,-50))
my_samples %>%
mutate(Estimate = cumsum(sample) / index,
IC_95_inf = Estimate - z_975 * sqrt(gamma2_hat / index),
IC_95_sup = Estimate + z_975 * sqrt(gamma2_hat / index)) %>%
ggplot(aes(x = index)) +
geom_ribbon(aes(ymin = IC_95_inf, ymax = IC_95_sup),
fill = "lightblue") +
geom_line(aes(y = Estimate)) +
labs(x = "M", title = "Estimation Monte Carlo de I",
y = expression(hat(I)[M])) +
theme_bw()
#Regardons si il y a bien une diminution de la variance.
gamma2_hat
Nous avons donc réussi à réduire la variance de l'estimateur de Monte - Carlo.
Remerciements¶
Ce problème nous a été donné par Mr Pierre Gloaguen à l'ENSTA Paris : https://papayoun.github.io/courses/2020_monte_carlo/ .