Ce devoir fut discuté avec les élèves suivants : Alan GUIVARCH, Romain PASDELOUP et Raphaël AUDIN.
Durant ce TP mon intuition n’a pas toujours été un très bon guide, c’est pourquoi mes réponses sont souvent déduites de l’observation des graphes et des valeurs, que je n’ai pas toujours inclus, et c’est également pourquoi je ne mentionne pas toujours mon intuition.
set.seed(255)
# Dans le cadre du jeu de ALice et Bob a peut être modélisé par le code suivant toutefois dans le cadre de mes tests j'ai préféré fixé moi-même a
a = sample(x = 0:100, size = 1,replace=T)
# Fonction pour générer les X1, ..., Xn tirage indépendants et uniformes
generator = function(a,n){
if (a >= 0 && n >=2)
return(runif(n, 0, a))
else return(-1)
}
# On cherche à trouver l'espérance de M
nbexp =100
tabmax = c(1:nbexp)
for (i in tabmax) {
tabmax[i] = max(generator(100, 10))
}
cat("Test de base : \n")## Test de base :cat("E(M) =", mean(tabmax), "; a =", a, " et n = ",10,"\n\n")## E(M) = 91.91309 ; a = 52 et n = 10# On teste l'espérance pour différentes valeurs de a et n fixées
a1 = 2
a2 = 20
a3 = 200
n0 = 2
n1 = 10
n2 = 100
n3 = 1000
nbexp =1000
tabmax0 = c(1:nbexp)
tabmax1 = c(1:nbexp)
tabmax2 = c(1:nbexp)
tabmax3 = c(1:nbexp)
tabmax4 = c(1:nbexp)
tabmax5 = c(1:nbexp)
tabmax6 = c(1:nbexp)
for (i in tabmax1) {
tabmax0[i] = max(generator(a2, n0))
tabmax1[i] = max(generator(a1, n1))
tabmax2[i] = max(generator(a1, n2))
tabmax3[i] = max(generator(a1, n3))
tabmax4[i] = max(generator(a1, n3))
tabmax5[i] = max(generator(a2, n3))
tabmax6[i] = max(generator(a3, n3))
}
cat("TEST pour un n le plus petit possible : \n")## TEST pour un n le plus petit possible :cat("E(M) =", mean(tabmax0), "; a =", a2, " et n = ",n0,"\n\n")## E(M) = 13.43733 ; a = 20 et n = 2cat("TESTS concernant l'impact de la valeur de n: \n" )## TESTS concernant l'impact de la valeur de n:cat("E(M) =", mean(tabmax1), "; a =", a1, " et n = ",n1,"\n")## E(M) = 1.819809 ; a = 2 et n = 10cat("E(M) =", mean(tabmax2), "; a =", a1, " et n = ",n2,"\n")## E(M) = 1.980702 ; a = 2 et n = 100cat("E(M) =", mean(tabmax3), "; a =", a1, " et n = ",n3,"\n\n")## E(M) = 1.998072 ; a = 2 et n = 1000cat("TESTS concernant l'impact de la valeur de A : \n")## TESTS concernant l'impact de la valeur de A :cat("E(M) =", mean(tabmax4), "; a =", a1, " et n = ",n3,"\n")## E(M) = 1.998003 ; a = 2 et n = 1000cat("E(M) =", mean(tabmax5), "; a =", a2, " et n = ",n3,"\n")## E(M) = 19.98028 ; a = 20 et n = 1000cat("E(M) =", mean(tabmax6), "; a =", a3, " et n = ",n3,"\n")## E(M) = 199.7971 ; a = 200 et n = 1000Rappelons que l’espérance de M correspond à sa moyenne d’où l’utilisation de la fonction mean pour la calculer.
D’après les tests n possède une influence sur le résultat de la E(M), toutefois celle-ci ne se voit pas forcément bien avec la moyenne, le mieux est de faire un histogramme où on peut y voir que plus le n est grand plus les valeurs tirées ont de chance d’être proche de a. Sur E(M) cela se traduit par une moyenne plus proche de a.
Pour ma part j’avais observé graphiquement que : si on multiplie 10 fois n la longueur de l’intervalle des plus grandes barres de l’histogramme est divisée par 10 donc j’en avais déduit que n serait potentiellement au dénominateur dans E(M).
D’après les tests a possède une influence sur le résultat de la E(M). En effet si on multiplie a par 10, E(M) sera multiplié par 10. Il semblerait donc que a puisse avoir une place dans E(M) en tant qu’une partie du numérateur.
D’après les observations ci-dessus nous avions:
-a présent au numérateur
-n présent au dénominateur
-E(M) < a
Or si E(M) est inférieur à a nous savions que ce par quoi a est multiplié devait être inférieur à 1. Nous avons donc a * Z avec Z < 1. Nous savions que Z n’était pas inférieur de beaucoup à 1 et que n était au dénominateur. C’est avec le test de a = 20 et n = 2 qu’il a été le plus facile de définir que : \[E(M) = \frac{a \cdot n}{n+1} \underset{n\to +\infty}\to{a} \] Cette formule trouvée empiriquement a été vérifiée par calcul.
Afin que E(M) puisse être égale à a il faudrait multiplier E(M) par (n+1)/n.
Or d’après le cours : \[E(aX) = a\cdot E(X)\]
J’en déduis que :
\[\frac{n+1}{n} \cdot E(M) = E(M\cdot \frac{n+1}{n})\]
Donc la correction de M est: \[M = max_{i} X_{i} \times \frac{n+1}{n}\]
cat("Test de base : \n")## Test de base :cat("Var(M) =", var(tabmax), "; a =", a, " et n = ",10,"\n\n")## Var(M) = 58.63002 ; a = 52 et n = 10cat("TEST pour un n le plus petit possible : \n")## TEST pour un n le plus petit possible :cat("Var(M) =", var(tabmax0), "; a =", a2, " et n = ",n0,"\n\n")## Var(M) = 22.76948 ; a = 20 et n = 2cat("TESTS concernant l'impact de la valeur de n :\n" )## TESTS concernant l'impact de la valeur de n :cat("Var(M) =", var(tabmax1), "; a =", a1, " et n = ",n1,"\n")## Var(M) = 0.02796178 ; a = 2 et n = 10cat("Var(M) =", var(tabmax2), "; a =", a1, " et n = ",n2,"\n")## Var(M) = 0.0003490647 ; a = 2 et n = 100cat("Var(M) =", var(tabmax3), "; a =", a1, " et n = ",n3,"\n\n")## Var(M) = 3.706551e-06 ; a = 2 et n = 1000cat("TESTS concernant l'impact de la valeur de a : \n")## TESTS concernant l'impact de la valeur de a :cat("Var(M) =", var(tabmax4), "; a =", a1, " et n = ",n3,"\n")## Var(M) = 3.831428e-06 ; a = 2 et n = 1000cat("Var(M) =", var(tabmax5), "; a =", a2, " et n = ",n3,"\n")## Var(M) = 0.0003745123 ; a = 20 et n = 1000cat("Var(M) =", var(tabmax6), "; a =", a3, " et n = ",n3,"\n")## Var(M) = 0.04039645 ; a = 200 et n = 1000On sait que la formule de la variance est la suivante : \[Var(M) = E(M^{2}) - [E(M)]^{2}\] Donc si on applique la formule de E(M) trouvée précédemment on obtient : \[Var(M) =\frac{a^{2} \cdot n}{n+2} \cdot (\frac{a \cdot n}{n+1})^{2}\] Cette formule empirique a été vérifiée par le calcul.
set.seed(255)
a = sample(x = 0:100, size = 1,replace=T)
generator = function(a,n){
if (a >= 0 && n >=2)
return(runif(n, 0, a))
else return(-1)
}
# Définition de M'
Mprime = replicate(1000, (2/100)* sum(generator(a, 100)))
cat("E(M') =", mean(Mprime), "; a =", a, " et n = ",100,"\n\n")## E(M') = 51.8746 ; a = 52 et n = 100#TESTS
a1 = 2
a2 = 20
a3 = 200
n0 = 2
n1 = 10
n2 = 100
n3 = 1000
Mtest0 = replicate(1000, (2/n0)*sum(generator(a2, n0)))
Mtest1 = replicate(1000, (2/n1)*sum(generator(a1, n1)))
Mtest2 = replicate(1000, (2/n2)*sum(generator(a1, n2)))
Mtest3 = replicate(1000, (2/n3)*sum(generator(a1, n3)))
Mtest4 = replicate(1000, (2/n3)*sum(generator(a1, n3)))
Mtest5 = replicate(1000, (2/n3)*sum(generator(a2, n3)))
Mtest6 = replicate(1000, (2/n3)*sum(generator(a3, n3)))On remarque que M’ semble être un estimateur de a au même titre que M. On sait que :
\[ E(X) =\frac{\sum_{i} X_{i}}{n} = \frac {a}{2} \] Or d’après la linéarité de l’espérance on obtient : \[E(M') = E(\frac{2}{n} \sum_{i} X_{i}) = 2 \cdot E(\frac{\sum_{i} X_{i}}{n}) = 2\cdot E(X)\]
Conclusion : \[ E(M') = \frac {a}{2} \times 2 = a\]
cat("Var(M') =", var(Mprime), "; a =", a, " et n = ",100,"\n\n")## Var(M') = 9.463486 ; a = 52 et n = 100En testant sur différentes valeurs de a et de n on peut voir que a et n jouent toujours un rôle dans la variance comme dans la question précédente.
On avait Var(M) = 0.001725208 avec M comme estimateur de a et maintenant nous avons Var(M’) = 9.463486 pour M’ comme estimateur de a. Pour un même a on a une variance presque 10 000 fois plus grande (9,46/0.001 = 9460).
Ce qui s’explique par la formule de la variance : \[Var(M') = E(M'^2) - (E(M'))² \]
Qui lorqu’elle est écrite avec la linéarité de l’espérance et suivant l’espérance que nous avons déterminé précedemment donne :
\[Var(M') = E((2X)^2) - (E(2X))^2 = 4\cdot E(X^2) - 4\cdot(E(X))^2 = 4\cdot [E(X^2) - (E(X))^2] \]
Et même si l’espérance semble être plus proche de a, avec une variance aussi élevée M’ ne semble pas être un très bon estimateur comparé à M.
Cette stratégie ne semble pas très bonne. En effet l’espérance du gain sera assez faible puisque même si il gagne dans tous les cas il gagnerait :
\[Gain = r(M) - M = 1,1\cdot M - M = 0,1\cdot M \] Soit au mieux 1 centime d’euro si M = 1. Sachant que M est compris entre O et 1 cela est peu mais justement sachant que le gain se trouve entre 0 et 1 euros il faut aussi relativiser.
set.seed(255)
n = 10
# On génère A et les Xi
A = runif(1, 0,1)
X = runif(n, 0, A)
# On répète plusieurs fois l'expérience
nbexp =1000
M = c(1:nbexp)
rM = c(1:nbexp)
Gain = c(1:nbexp)
for (i in M) {
M[i] = max(X)
rM[i] = 1.1 * M[i]
if (rM[i] > A)
Gain[i] = 0
else
Gain[i] = rM[i] - M[i]
}
mean(M)## [1] 0.4428638mean(rM)## [1] 0.4871501mean(Gain)## [1] 0.04428638cat("Pour A = ", A,", Bob gagne en moyenne : ", mean(Gain), "euros")## Pour A = 0.5242901 , Bob gagne en moyenne : 0.04428638 eurosVia la simulation, on voit bien que son espérance de gain moyenne est faible puisqu’elle est de l’ordre de 0,04 euros.
set.seed(255)
n = 3
#Fonction qui va générer A, les Xi, M et retourner si A = a && M = m
est_egal = function(n, a, m){
A = sample(x=(0:10)/10, 1, replace=T)
X = sample(x=(0:(10*A))/10, n, replace=T)
M = max(X)
return(A == a && M == m)
}
#TEST de base
a = 0.8
m = 0.6
test = replicate(10000, est_egal(n,a,m))
#TESTS pour différentes valeurs de a et M
aVecteur = (0:10)/10
mVecteur = (0:10)/10
for (i in 1:10){
for (j in 1:10){
test = mean(replicate(1000, est_egal(n,aVecteur[i],mVecteur[j])))
#cat("Pour a =", aVecteur[i], " et m = ", mVecteur[j]," P[A = a, M = m] = ", test,"\n")
}
#cat("\n")
}set.seed(255)
# TESTS différentes valeurs de A pour M = 0.5
aVecteur = (0:10)/10
mVecteur = 0.5
for (i in 1:10){
test = mean(replicate(1000, est_egal(n,aVecteur[i],mVecteur)))
cat("Pour a =", aVecteur[i], " P[A = a, M = 0.5] = ", test,"\n")
}## Pour a = 0 P[A = a, M = 0.5] = 0
## Pour a = 0.1 P[A = a, M = 0.5] = 0
## Pour a = 0.2 P[A = a, M = 0.5] = 0
## Pour a = 0.3 P[A = a, M = 0.5] = 0
## Pour a = 0.4 P[A = a, M = 0.5] = 0
## Pour a = 0.5 P[A = a, M = 0.5] = 0.041
## Pour a = 0.6 P[A = a, M = 0.5] = 0.023
## Pour a = 0.7 P[A = a, M = 0.5] = 0.013
## Pour a = 0.8 P[A = a, M = 0.5] = 0.013
## Pour a = 0.9 P[A = a, M = 0.5] = 0.011Mes tests m’ont montré, comme me l’indiquait mon intuition, que si a < m la probabilité va être nulle. Donc si l’ont veut respecter l’esprit du jeu et avoir une probabilité positive, pour un m = 0.5 il faut prendre a >= 0.5.
Sachant que M <= 1 et alpha < 1, au maximum on peut espérer un gain de 1 euros ce qui semble mieux que le 1 centime de la question 2.3. De plus r(M) ne dépendra d’un pourcentage de M, qui est une petite valeur, mais d’une puissance de M.
set.seed(120)
n = 2
f = function (n){
Galpha = 0
# Calcul du gain pour différentes valeurs de alpha
for (alpha in c(0.7, 0.5, 0.3)){
# On génère A et Xi
A = runif(1, 0,1)
X = runif(n, 0, A)
nbexp =1000
M = c(1:nbexp)
rM = c(1:nbexp)
Gain = c(1:nbexp)
for (i in M) {
M[i] = max(X)
rM[i] = M[i]**alpha
if (rM[i] > A)
Gain[i] = 0
else
Gain[i] = rM[i] - M[i]
}
G = mean(Gain)
cat("pour alpha = ", alpha," G = ", G,"\n")
}
}
# On appelle la fonction pour n = 2
f(n)## pour alpha = 0.7 G = 0.1111848
## pour alpha = 0.5 G = 0.235369
## pour alpha = 0.3 G = 0.1113063La meilleure valeur de alpha semble être 0.5. Je ne pense pas que l’on puisse trouver de meilleur alpha. En effet si veut gagner plus on pourrait se dire qu’il faudrait augmenter alpha mais à 0.7 (qui est supérieur à 0.5) on obtient un gain moyen moins grand.
set.seed(120)
n = 10
# On appelle la fonction pour n = 10
f(n)## pour alpha = 0.7 G = 0
## pour alpha = 0.5 G = 0.04025185
## pour alpha = 0.3 G = 0Pour n = 10 on obtient la même conclusion. Toutefois je pense que la variance va être plus importante.
La probabilité que Bob a de gagner correspond à la probabilité qu’Alice face Pile.
Étant donné que la probabilité de gagner va dépendre uniquement du lancer de la pièce donc du résultat, soit Pile soit Face, si celle-ci n’est pas biaisée la probabilité devrait être de 0,5.
set.seed(255)
n=1000
# Lancer de la pièce
coin = sample(x=c("P","F"), n , replace=T)
mean(coin=="P")## [1] 0.481Au début j’avais pensé à une stratégie assez complexe en me disant que Bob répondrait forcément “non” si le chiffre est <= 0.1 et forcément “oui” si il est >= 0.9. Tandis que pour les nombres entre ]0.1,0.9[ il aurait répondu “oui” avec une probabilité de “1/2” et non avec une probabilité de 1/2.
Cependant ma stratégie était un peu complexe et dure à écrire.
Et c’est en relisant la question précédente que je me suis rendue compte qu’il suffisait de garder les mêmes probabilités que l’on avait trouvé mais que la réponse de Bob ne dépendrait pas de la valeur de la pièce mais plutôt de la valeur du nombre donné.
Bob répondra donc “oui” si le nombre donné par Alice est >= 0.5 et “non” si il est < 0.5. En effet un nombre plus grand a plus de chance d’être A2 que A1, et inversement.
Le gain de notre stratégie pour A1 = 0.4 et A2 = 0.6 est que Bob a 100% de chance de gagner.
set.seed(100)
n= 1000
# On génère A1 et A2
A1 = runif(n,0.0,1.0)
A2 = runif(n,A1,1.0)
nbex= 1000
# On lance la pièce
coin = sample(x=c("P","F"),size=nbexp,replace=T)
res = c(1:nbexp)
# Bob gagne si Alice donne A2 avec A2 >=0.5 ou bien si Alice donne A1 avec A1 <0.5
for(i in res) {
res[i] = (coin[i] == "P" && A2[i] >= 0.5 || coin[i] == "F" && A1[i] < 0.5)
}
cat("Bob a une espérance de gagner de : ", mean(res == TRUE) , "\n")## Bob a une espérance de gagner de : 0.67Avec notre stratégie Bob a 67% de chance de gagner ce qui est bien supérieur à la première stratégie de Bob.