RICM4 Probabilités et Simulation – Devoir 1
Antoine Saget
06/12/2019
Discussions avec : Claire Velut, Adrien Artaud, Tom Graugnard
Exercice 1 : Question préliminaire à propos d’estimation
Question 1.1
Intuition. On s’attend à ce que l’estimateur soit plus précis lorsque \(n\) est grand. On décide de tirer des réels, donc mon hypothèse est que \(a\) n’aura pas d’influence.
Variation de l’espérance en fonction de \(a\) et \(n\)
n = 10
i = 1
as = c(1, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000)
EM = c()
ecart = c()
for (a in as) { # Pour différentes valeurs de a
EM[i] <- mean(replicate(10000, max(runif(n, min = 0, max = a)))) # Faire la moyenne de 10000 répétitions du jeu avec ces paramètres a et n
ecart[i] <- EM[i] / a # Calculer le quotient entre la moyenne obtenue et la valeur attendue
i = i + 1
}
plot(as, ecart, ylim = c(0, 1), log = "x", main = "Quotient entre espérance trouvée empiriquement et valeur attendue \n pour différentes valeurs de a et n = 10", xlab = "a", ylab = "Quotient")
Il n’y a donc pas de variations de la qualité de l’estimateur pour différentes valeurs de \(a\). Le quotient entre l’espérance calculée et la valeur attendue est toujours le même. Cela fais sens puisqu’on tire des valeurs réels : tirer entre \(0\) et \(1\) est équivalent à tirer entre \(0\) et \(1000\).
On s’intéresse maintenant à l’espérance de \(M\) pour différentes valeurs de \(n\) :
a = 10
i = 1
ns = c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000)
EM = c()
ecart = c()
for (n in ns) { # Pour différentes valeurs de n
EM[i] <- mean(replicate(1000, max(runif(n, min = 0, max = a)))) # Faire la moyenne de 1000 répétitions du jeu avec ces paramètres a et n
i = i + 1
}
df <- data.frame(ns, EM, (ns / (ns + 1)) * a)
ggplot(df, aes(ns, EM)) +
geom_point(aes(y=EM), colour="red") +
geom_line(aes(y=(ns / (ns + 1)) * a), colour="green") +
scale_x_continuous(trans='log10') +
xlab("n") +
ylab("E[M]") +
ggtitle("E[M] pour différentes valeurs de n", "En vert les points obtenus avec la formule, en rouge les points obtenu empiriquement")
Ici on peut voir que plus n augmente plus le quotient entre l’esperance de l’estimateur et la valeur attendue ce rapproche de \(1\). L’estimateur est plus précis pour des grandes valeurs de \(n\).
On trouve empiriquement pour \(E[M] = {n \over n + 1} \times a\).
Correction de l’estimateur
Pour corriger M il suffit d’ajouter \({1 \over n + 1} \times a\) : \({1 \over n + 1} \times a + {n \over n + 1} \times a = a \times ({n + 1 \over n + 1}) = a\)
df <- data.frame(ns, EM + a * (1 / (ns + 1)), (ns / (ns + 1)) * a)
ggplot(df, aes(ns, EM)) +
geom_point(aes(y = EM + a * (1 / (ns + 1))), colour="blue") +
geom_point(aes(y = EM), colour="red") +
geom_line(aes(y = a), colour="green") +
scale_x_continuous(trans='log10') +
ylim(0, 11) +
xlab("n") +
ylab("E[M]") +
ggtitle("E[M] pour différentes valeurs de n", "En vert les valeurs attendues, en rouge les valeurs obtenues une fois corrigées, \n en bleu les valeurs obtenues sans correction")
Variance
On avais \(E[M] = {n \over n + 1} \times a\). Donc \(Var[M] = E[M^2] - E[M]^2 = ({n \over n + 2} - {n^2 \over (n + 1)^2})) \times a^2\)
a = 10
i = 1
ns = c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000)
VM = c()
for (n in ns) { # Pour différentes valeurs de n
VM[i] <- var(replicate(1000, max(runif(n, min = 0, max = a)))) # Faire la moyenne de 1000 répétitions du jeu avec ces paramètres a et n
i = i + 1
}
df <- data.frame(ns, VM, ((ns / (ns + 2)) - ((ns**2) / ((ns + 1)**2))) * (a**2))
ggplot(df, aes(ns, VM)) +
geom_point(aes(y=VM), colour="red") +
geom_line(aes(y=((ns / (ns + 2)) - ((ns**2) / ((ns + 1)**2))) * (a**2)), colour="green") +
scale_x_continuous(trans='log10') +
xlab("n") +
ylab("Var[M]") +
ggtitle("Var[M] pour différentes valeurs de n", "En vert les points obtenus avec la formule, en rouge les points obtenu empiriquement")
Question 1.2
Intuition. \({1 \over n }\sum_i X_i\) correspond à l’espérance de \(X\). Donc \(M'\) est égal à deux fois l’espérance. Puisque on tire des nombre uniformément entre \(0\) et \(a\), on s’attend à ce que \(E[X] = {a \over 2}\) et donc \(M' = a\). Empirquement on s’attend à observer des valeurs moins précises pour des petites valeurs de \(n\).
Espérance
a = 10
i = 1
ns = c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000)
EMM = c()
ecart = c()
for (n in ns) { # Pour différentes valeurs de n
EMM[i] <- mean(replicate(1000, (2 / n ) * sum(runif(n, min = 0, max = a)))) # Faire la moyenne de 1000 répétitions du jeu avec ces paramètres a et n
i = i + 1
}
df <- data.frame(ns, EMM)
ggplot(df, aes(ns, EMM)) +
geom_point(aes(y=EMM), colour="red") +
scale_x_continuous(trans='log10') +
xlab("n") +
ylab("E[M']") +
ggtitle("E[M'] pour différentes valeurs de n") +
ylim(0, 11)
Preuve de l’intuition :
\[ \begin{array}{lcl} E[MM] & = & E[{2 \over n} \sum_i X_i]\\ & = & {2 \over n} \sum_i E[X_i]\\ & = & {2 \over n} \sum_i {a \over 2}\\ & = & {2 \over n} \times {n \times a \over 2} \\ & = & a \end{array} \]
Variance
a = 10
i = 1
ns = c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000)
VMM = c()
ecart = c()
for (n in ns) { # Pour différentes valeurs de n
VMM[i] <- var(replicate(1000, (2 / n ) * sum(runif(n, min = 0, max = a)))) # Faire la moyenne de 1000 répétitions du jeu avec ces paramètres a et n
i = i + 1
}
df <- data.frame(ns, VMM)
ggplot(df, aes(ns, VMM)) +
geom_point(aes(y=VMM), colour="red") +
scale_x_continuous(trans='log10') +
xlab("n") +
ylab("Var[M']") +
ggtitle("Var[M'] pour différentes valeurs de n")
On remarque une variance très élevée pour \(n < 100\), l’estimateur n’est donc pas si bon. Dans le cadre du jeu entre Alice et Bob, il semble réaliste qu’ils choisissent un \(n\) relativement petit (\(< 100\)) donc il vaut peut-être mieux utiliser \(M\) que \(M'\).
Exercice 2 : Un deuxième jeu à base de max
Question 2.3
Intuition. On gagne lorsque l’on tire \(n\) nombre inférieur à \(A \over 1.1\). On peu aussi dire qu’on gagne lorsqu’on tire \(0\) nombre supérieur à \(A \over 1.1\). La probabilité qu’un nombre soit supérieur à \(A \over 1.1\) lors d’un tirage uniforme entre \(0\) et \(A\) est \(1 - {1 \over 1.1}\). Donc on va gagner \(0.1 \times M\) avec une probabilité \(P(X = 0)\) où \(X\) suit une loi binomiale de paramètre \(10\), \(1 - 1/1.1\). Pour \(n = 10\) :
dbinom(0, size=10, prob=1 - 1/1.1) # chances de gagner quelques chose## [1] 0.3855433On gagne peu souvent. La stratégie n’est peut-être pas la meilleure. Mais on ne peut pas perdre d’argent dans ce jeu donc toute stratégie est bonne.
N = 100000
W = sample(0, N, replace = TRUE)
n = 10
for(i in seq(1, N)) # Répéter n fois le jeu
{
a = runif(1)
X = runif(n, min = 0, max = a)
M = max(X)
R = M * 1.1
if (R < a) {
W[i] = R - M
}
}
mean(W)## [1] 0.01594129Pgain = sum(W != 0) / N # La proportion de parties que Bob a gagné
Pgain## [1] 0.38424Bob gagne donc en moyenne \(\approx 1.6\) centimes. On retrouve la probabilité de gain calculé avec la loi binomiale précédemment : \(Pgain = 0.38655 \approx 03855 = P(X = 0)\).
Question 2.4
Intuition. On s’attend à des \(P(A = a, M = m)\) qui diminue lorsque a augmente et m diminue : il est moins probable que la maximum des trois valeurs tirées soit petites lorsqu’on les tire dans un intervalle de plus en plus grand.
n = 3
N = 100000
res = matrix(sample(0, 11*11, replace = TRUE), nrow = 11, ncol = 11)
for(i in seq(1, N)) # On répète N fois l'expérience
{
a = sample(seq(0, 1, 0.1), 1) # A chaque fois on tire un nouveau a et m
X = sample(seq(0, a, 0.1), n, replace = TRUE)
m = max(X)
res[a*10 + 1, m*10 + 1] <- res[a*10 + 1, m*10 + 1] + 1 # On incrémente la case a,m du tableau
}
res = res / N # On divise par N : chaque case a,m correspond à P(A = a, M = m)
res## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.09101 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [2,] 0.01123 0.07730 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [3,] 0.00310 0.02330 0.06479 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [4,] 0.00152 0.00951 0.02823 0.05279 0.00000 0.00000 0.00000 0.00000 0.00000
## [5,] 0.00067 0.00486 0.01433 0.02664 0.04530 0.00000 0.00000 0.00000 0.00000
## [6,] 0.00043 0.00306 0.00768 0.01618 0.02520 0.03785 0.00000 0.00000 0.00000
## [7,] 0.00024 0.00173 0.00490 0.01026 0.01663 0.02417 0.03389 0.00000 0.00000
## [8,] 0.00017 0.00122 0.00330 0.00631 0.01043 0.01598 0.02273 0.02976 0.00000
## [9,] 0.00018 0.00097 0.00226 0.00454 0.00778 0.01175 0.01569 0.02113 0.02661
## [10,] 0.00006 0.00065 0.00168 0.00336 0.00580 0.00870 0.01180 0.01493 0.01997
## [11,] 0.00013 0.00049 0.00131 0.00272 0.00371 0.00603 0.00848 0.01112 0.01528
## [,10] [,11]
## [1,] 0.00000 0.00000
## [2,] 0.00000 0.00000
## [3,] 0.00000 0.00000
## [4,] 0.00000 0.00000
## [5,] 0.00000 0.00000
## [6,] 0.00000 0.00000
## [7,] 0.00000 0.00000
## [8,] 0.00000 0.00000
## [9,] 0.00000 0.00000
## [10,] 0.02423 0.00000
## [11,] 0.01859 0.02335La colonne 6 de la matrice correspond à \(M = 0.5\), on peut donc retrouver pour case de cette colonne une estimation de \(P[A = a | M = 0.5]\). La probabilité la p lus forte si \(M = 0.5\) est avec \(A = 0.5\). Donc on a tout intérêt à lui proposer cette valeur de A pour gagner.
Question 2.5
Sans simulation je n’arrive pas à avoir l’intuition de dire si la stratégie est bonne ou pas. Mais la stratégie est valable car \(M^\alpha > M\) et \(M^\alpha < 1\; \forall M,\alpha \in [0, 1]\).
N = 100000
n = 2
j = 1
Pgain = c()
for(alpha in c(0.3, 0.5, 0.7))
{
W = sample(0, N, replace = TRUE)
for(i in seq(1, N)) # Répéter n fois le jeu
{
a = runif(1)
X = runif(n, min = 0, max = a)
M = max(X)
R = M ** alpha
if (R < a) {
W[i] = R - M
}
}
Pgain[j] <- mean(W)
j = j + 1
}
Pgain## [1] 0.05602628 0.06619387 0.05618128La meilleure valeur semble être avec \(\alpha = 0.5\) et un gain moyen de \(\approx 6.6\) centimes.
N = 10000
n = 2
j = 1
Pgain = c()
for(alpha in seq(0.1, 0.9, 0.1))
{
W = sample(0, N, replace = TRUE)
for(i in seq(1, N)) # Répéter n fois le jeu
{
a = runif(1)
X = runif(n, min = 0, max = a)
M = max(X)
R = M ** alpha
if (R < a) {
W[i] = R - M
}
}
Pgain[j] <- mean(W)
j = j + 1
}
plot(seq(0.1, 0.9, 0.1), Pgain, xlab = "alpha", ylab = "Gain moyen")
Au vu de ce graphe qui test d’autres valeurs de \(\alpha\), le gain maximal semble être atteint en \(\alpha = 0.5\).
N = 10000
n = 10
j = 1
Pgain = c()
for(alpha in seq(0.1, 0.9, 0.02))
{
W = sample(0, N, replace = TRUE)
for(i in seq(1, N)) # Répéter n fois le jeu
{
a = runif(1)
X = runif(n, min = 0, max = a)
M = max(X)
R = M ** alpha
if (R < a) {
W[i] = R - M
}
}
Pgain[j] <- mean(W)
j = j + 1
}
plot(seq(0.1, 0.9, 0.02), Pgain, xlab = "alpha", ylab = "Gain moyen")
seq(0.1, 0.9, 0.02)[Pgain == max(Pgain)]## [1] 0.84Avec \(n = 10\), le maximum semble être au alentours de \(\alpha = 0.84\).
Exercice 3 : Question BONUS : un dernier jeu à base de max et de pile ou face
Question 3.6
Il a 50% de chance de gagner puisqu’il y a 50% de chance que le nombre choisie par la pièce soit le plus grand des deux.
Question 3.7
Soit \(M\) le nombre indiqué par Alice. On peut fixer un nombre \(a \in [0, 1]\), si \(M > a\) alors on dit oui, sinon non. Si \(A_1 < A_2 < a\) alors on a une chance sur deux de donner la bonne réponse, si \(a < A_1 < A_2\) alors on a une chance sur deux de donner la bonne réponse, si \(A_1 < a < A_2\) alors on a toujours la bonne réponse. On cherche à donc choisir \(a\) pour que \(P(A_1 < a < A_2)\) soit maximal. Si le tirage est uniforme alors la meilleure valeur pour \(a\) est \(0.5\).
Si Alice choisis \(A_1 = 0.4\) et \(A_2 = 0.6\) alors on a bien \(A_1 < 0.5 < A_2\) donc on gange dans 100% des cas.
Si Alice tire uniformément des nombres entres \(0\) et \(1\) :
n = 10000
B1 = runif(n, 0, 1) # Tirages des A1 et A2
B2 = runif(n, 0, 1)
A1 = c()
A2 = c()
for(i in seq(1, n)) {
if(B1[i] > B2[i]) {
A2[i] <- B1[i]
A1[i] <- B2[i]
} else {
A2[i] <- B2[i]
A1[i] <- B1[i]
}
}
choice = sample(c(0, 1), n, replace = TRUE) # Tirage des lancé de pièces
res = sample(0, n, replace = TRUE)
for(i in seq(1, n)) { # Vérification pour chaque partie si elle a été gagné ou perdue
if ((choice[i] == 1 && A2[i] > 0.5) || choice[i] == 0 && A1[i] <= 0.5) {
res[i] <- 1
}
}
mean(res) # Moyenne des partie gagnée## [1] 0.7504On gagne donc de \(\approx\) 75% des cas.