DM_PS
OTBA_ZERAMDINI
04/12/2019
R Markdown
set.seed(40)
#Ce DM a été réalisé en collaboration avec 2 de mes camarades : Raphaël Audin et Aleck Bilounga.
#Exercice 1 : Question préliminaire à propos d'estimation
#Q1.1
a=5# a correspond au nombre a choisi secrètement par Alice.
n=2# n correspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
#Fonction qui retourne l'espérence de M, et qui prend a en paramètre afin de pouvoir
#représenter E(M) en fonction de a.
EMa<-function(a){
TM=replicate(nb,max(runif(n)*a)) #TM contient les nb maximuns des n nombres que tire Alice, nb fois dans [0,a].
return (sum(TM)/nb) #on somme ces maximums obtenus par répition nb fois de l'expérience : Alice tire n nombres dans l'intervalle[0,a], et on
#divise cette somme par nb dans le but d'avoir une moyenne des maximums obtenus ce qui correspond à l'espérance de M.
}
#On trace ensuite la fonction Ema en fonction de a dans le but de faire une hypothèse sur l'espérence de M, E(M).
#Nous fixons alors n, et nous faisons varier a.(n est içi fixé à 2). Nous ferrons varier n et nous fixerons à plus tard.
x=0:20#Nous traçons l'espérance en fonction de a sur 0 20, pour voir comment celle-ci évolue en fonntion de a.
y=lapply(x,EMa)
plot(x,y,xlab="a",ylab="E(M)")
#Le résultat est trés parlant nous trouvons une fonction linéaire en a. Ce qui veut dire que l'esperence est lineaire en a. Ainsi l'espérence s'écrit donc de la manière suivante : E(M)= K*a avec K dans R.
#Reste maintenant à déterminer ce K. Pour cela nous changeons les rôles de a et de n. Nous fixons n, et nous faisons varier a.
#On crée alors une fonction qui retourne l'espérence de M, et qui prend n en paramètre afin de pouvoir représenterreprésenter E(M) en fonction de n.
EMn<-function(n){
TM=replicate(nb, max(runif(n)*a))#TM contient les nb maximums des n nombres que tire Alice, nb fois dans [0,a].
return(sum(TM)/nb)#on somme ces maximums obtenu par répétition nb fois de l'expérience Alice : tire n nombre dans l'intervalle[0,a], et on
#divise cette somme par nb dans le but d'avoir une moyenne des maximums obtenus ce qui correspond à l'espérance de M.
}
#On trace ensuite la fonction EMn en fonction de n dans le but de faire une hypothèse sur l'espérance de M, E(M). Nous fixons alors a, et nous faisons varier n.(a est içi fixé à 5)
x=2:20
y=lapply(x,EMn)
plot(x,y,xlab="n",ylab="E(M)")
#Nous voyons que E(M) tend vers a=5, on observe une convergence asymptotique de la courbe obtenue vers a=5. Ainsi d'après ce qui précède nous avons vu que E(M)=K*a. Maintenant nous avons fixé a et nous avons trouvé que E(M) tendait vers a, ainsi K tend vers 1.
#De plus K dépend de n. Après plsuieurs essais et après une étude théorique que nous exposerons pas la suite nous pensons que K=(n/n+1). Une représentation graphique de la fonction n->a*(n/(n+1)) permet d'en avoir l'esprit clair.
#Représentons cette focntion qui à n->a*(n/n+1) :
a=5
fn<-function(n){
return(a*(n/(n+1)))
}
x=2:20
y=lapply(x,fn)
plot(x,y,xlab="n",ylab="a*(n/n+1)")
#La ressemblance entre EMn et fn est tout à fait frappante. Nous faisons alors l'hypothèse que K=n/n+1.
#Ainsi nous faisons l'hypothèse que E(M)=a*K=a*(1/n+1).
#Démontrons-le théoriquement : 
Espérence de M
set.seed(44)
#Nous savons que V(M)=E(M²)-E(M)², nous allons utlisé cette propriété pour déterminier la variance :
#Je rappelle que nous avons pris les valeurs suivantes, continuons avec ces mêmes valeurs :
a=5# a correspond au nombre a choisi secrètement par Alice.
n=2# n corrspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
#J'implémente alors une fonction EMa2 qui prend a en paramètre et qui renvoie l'espérence M², c'est-à-dire E(M²).
EMa2<-function(a){
TM2=replicate(nb,max(runif(n)*a)^2)#TM2 contient les nb maximuns au carée, des n nombres que tire Alice, nb fois dans [0,a].
return(sum(TM2)/nb)#on somme ces maximums obtenu par répition nb fois de l'expérience Alice : tire n nombre dans l'intervalle [0,a], et on
#divise cette somme par nb dans le but d'avoir une moyenne des carrés des maximums obtenus ce qui correspond à l'espérance.
}
#Puisque nous avons calculer l'espérence précédemment, il suffit de rappeler cette espérence puisque nous en avons besoins.
EMa<-function(a){
TM=replicate(nb,max(runif(n)*a))
return (sum(TM)/nb)
}
#Nous créons alors une fonction qui va prendre en paramètre a et nous rentre la variance de M. C'est-à-dire la fonction qui réalise la propriété suivante V(M)=E(M²)-E(M) :
VMa<-function(a){
return (EMa2(a)-(EMa(a))^2)
}
#Représentons, alors graphiquement cette variance en focntion de a et en fixant n pour se faire un idée de la manière avec laquelle varie en fonction de a.
x=1:20
y=lapply(x,VMa)
plot(x,y,xlab="a",ylab="V(M)")
#On remarque alors de la courbe ressemble grandement à une parabole. Ainsi nous faisons l'hypthèse que V(M)=K*a². Reste maintenant à déterminer K qui appartient à l'ensemble des réels et dépend de n.
#on remarque que la variance augmente de façons polynomiale en a² selon a.
#Maintenant faisons l'inverse : on fixe a et onf ait varier n.
#Ainsi nous créons les même fonctions que précédemment sauf que ces nouvelles fonctions prennent en paramètre n et non a, pour que l'on puisse représenter V(M) en fonction de n.
#J'implémente alors une fonction EMn2 qui prend n en paramètre et qui renvoie l'espérence M², c'est-à-dire E(M²).
EMn2<-function(n){
TM2=replicate(nb,max(runif(n)*a)^2)#TM2 contient les nb maximums au carée, des n nombres que tire Alice, nb fois dans [0,a].
return(sum(TM2)/nb)#on somme ces maximums obtenus par répition nb fois de l'expérience Alice : tire n nombre dans l'intervalle[0,a], et on
#divise cette somme par nb dans le but d'avoir une moyenne des carrés des maximums obtenus ce qui correspond à l'espérance de M.
}
#Puisque nous avons calculé l'espérence précédemment, il suffit de rappeler cette espérence puisque nous en avons besoin.
EMn<-function(n){
TM=replicate(nb,max(runif(n)*a))
return (sum(TM)/nb)
}
#Nous créons alors une fonction qui va prendre en paramètre n et nous rentre la variance de M. C'est-à-dire la fonction qui réalise la propriété suivante V(M)=E(M²)-E(M) :
VMn<-function(n){
return (EMn2(n)-(EMn(n))^2)
}
#Représentons, alors graphiquement cette variance en fonction de a et en fixant n pour se faire un idée de la manière avec laquelle varie en fonction de a.
x=2:20
y=lapply(x,VMn)
plot(x,y,xlab="n",ylab="V(M)")
#La courbe obtenue de V(M) en fonction de M ressemble énormément à une hyperbole ainsi une étude théorique nous permet de connaitre le K(n) recherhé en fait grâce à la démonstration théorique suivante on en déduit que v(M)=a²*(n/((n+2)*(n+1)²))
#Démontrons le : 
Espérence de M
set.seed(45)
#On peut alors tracer la variance théorique pour la comparer à la variance empirique obtenue plus haut :
#on rappelle les variables fixées pour l'éxemple :
a=5# a correspond au nombre a choisi secrètement par Alice.
n=2# n correspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
VTHn<-function(n){
return((n*a^2)/((n+2)*(n+1^2)))
}
x=2:20
y=lapply(x,VTHn)
plot(x,y,xlab="n",ylab="V(M)")
#Le résultat est alors très similaire en la théorie et la pratique.#Q1.2
a=5# a correspond au nombre a choisi secrètement par Alice.
n=2# n correspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
#Fonction qui retourne l'espére de M', et qui prend a en paramètre afin de pouvoir
#représenter E(M') en fonction de a.
EMPa<-function(a){
TS=replicate(nb,(2/n)*sum(runif(n)*a)) #TS contient les nb : (2/n) fois la somme des n nombres que tire Alice, nb fois dans [0,a].
return (sum(TS)/nb) #on somme alors ces nombres et on divise cette somme par nb dans le but d'avoir une moyenne des maximums obtenus précedemment, ce qui correspond à l'espérance de M'.
}
#On trace alors cette espérence E(M') en fonction de a et en fixant n pour voir comment varie E(M') en fonction de a :
x=1:100
y=lapply(x,EMPa)
plot(x,y,xlab="a",ylab="E(M')")
#On obtient clairement d'après la courbe et les résultats obtenus, une fonction linéaire de coefficient directeur 1.
#Il est donc claire empiriquement, que E(M')=1*a=a.
#Pour en être certain faisons le contraire fixons a et faisons varier n :
#Comme précedemment on crée une fonction qui retourne l'espére de M', et qui prend n en paramètre afin de pouvoir représenter E(M') en fonction de n.
EMPn<-function(n){
TS=replicate(nb,(2/n)*sum(runif(n)*a)) #TS contient les nb : (2/n) fois la somme des n nombres que tire Alice, nb fois dans [0,a].
return (sum(TS)/nb) #on somme alors ces nombres et on divise cette somme par nb dans le but d'avoir une moyenne des maximums obtenus précédemment ce qui correspond à l'espérance de M'.
}
#On trace alors cette espérence E(M') en fonction de a et en fixant n pour voir comment varie E(M') en fonction de a :
x=2:100
y=lapply(x,EMPn)
plot(x,y,xlab="n",ylab="E(M')") 
#On remarque les valeurs de E(M') tourne autour de a=5, et cela pour tout valeurs de n, ce qui veut dire que E(M') ne dépend pas de n.
#Ainsi on peut conclure qu'empiriquement E(M')=a.
#Prouvons le théoriquement :
Espérence de M’
set.seed(46)
#Nous savons que V(M)=E(M'²)-E(M')², nous allons utlisé cette propriété pour déterminier la variance :
#Je rappelle que nous avons pris les valeurs suivantes continuons avec ces mêmes valeurs :
a=5# a correspond au nombre a choisi secrètement par Alice.
n=2# n corrspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
#J'implémente alors une fonction EMPa2 qui prend a en paramètre et qui renvoie l'espérence M'², c'est-à-dire E(M'²).
EMPa2<-function(a){
TM2=replicate(nb,(2/n*(runif(n)*a))^2)#TM2 contient les nb maximuns au carée, des n nombres que tire Alice, nb fois dans [0,a].
return(sum(TM2)/nb)#on somme ces maximum obtenu par répition nb fois de l'expérience Alice : tire n nombre dans l'intervalle[0,a], et on divise cette somme par nb dans le but d'avoir une moyenne des carrés des maximums obtenus ce qui correspond à l'espérance.
}
#Puisque nous avons calculer l'espérence précédemment, il suffit de rappeler cette espérence, puisque nous en avons besoins.
EMPa<-function(a){
TM=replicate(nb,max(runif(n)*a))
return (sum(TM)/nb)
}
#Nous créons alors une fonction qui va prendre en paramètre a et nous retourner la variance de M'. C'est-à-dire la fonction qui réalise la propriété suivante V(M')=E(M'²)-E(M') :
VMPa<-function(a){
return (EMPa2(a)-(EMPa(a))^2)
}
#Représentons, alors graphiquement cette variance en fonction de a et en fixant n pour se faire un idée de la manière avec laquelle varie V(M') en fonction de a.
x=1:20
y=lapply(x,VMPa)
plot(x,y,xlab="a",ylab="V(M')")
#On remarque alors de la courbe ressemble grandement à une parabole. Ainsi nous faisons l'hypothèse que V(M')=K*a². Reste maintenant à déterminer K qui appartient à l'ensemble des réels et dépend de n. Ainsi plus a est grand et plus la variance est grande. Et la variance croit de manière polynomiale selon a.
#A présent fixons a et faisons varier n.Nous implémentons alors les mêmes fonctions que précedemment. Mais cette fois çi elles auront n en paramètre n pour puvoir représenter V(M') en fonction de n et ainsi tirer des conclusions finales quant à la valeur de V(M') :
EMPn2<-function(n){
TM2=replicate(nb,(2/n*(runif(n)*a))^2)#TM2 contient les nb maximuns au carée, des n nombres que tire Alice, nb fois dans [0,a].
return(sum(TM2)/nb)#on somme ces maximum obtenu par répition nb fois de l'expérience Alice : tire n nombres dans l'intervalle[0,a], et on
#divise cette somme par nb dans le but d'avoir une moyenne des carrés des maximums obtenus ce qui correspond à l'espérance.
}
#Puisque nous avons calculé l'espérence précédemment, il suffit de rappeler cette espérence puisque nous en avons besoin.
EMPn<-function(n){
TM=replicate(nb,max(runif(n)*a))
return (sum(TM)/nb)
}
#Nous créons alors une fonction qui va prendre en paramètre a et nous rendre la variance de M'. C'est-à-dire la fonction qui réalise la propriété suivante V(M')=E(M'²)-E(M') :
VMPn<-function(n){
return (EMPn2(n)-(EMPn(n))^2)
}
#Représentons, alors graphiquement cette variance en fonction de n et en fixant a pour se faire un idée de la manière avec laquelle varie V(M') en fonction de n.
x=1:20
y=lapply(x,VMPn)
plot(x,y,xlab="n",ylab="V(M')")
#On remarque que ça à l'air d'évolue en 1/n. Mais il est très compliqué de déterminer ce coefficient empiriquement. Une étude mathématiques toute simple nous permet alors de conclure que V(M')=a²/3n.
#Prouvons-le :
Variance de M’
set.seed(47)
#Nous allons utilisé le risque quadratique, pour savoir lequel des estimateurs est meilleur que l'autre :
#C'est-à-dire que nous allons calculé le risque quadratique des deux estimateurs de M et M' puis les comparer. L'estimateur qui aura le risque quadratique les plus bas est meilleur que l'autre.
#Nous allons réalisé cela de manière totalement théorique :
#D'abord déterminons le risque quadratique de M comme estimateur de a :
Risque quadratique de M comme estimateur de a
#Nous trouvons un risque quadratique pour M comme estimateur de a de E((M-a)²)=2a²/((n+2)*n+1))
#Cherchons maintenant le risque quadratique de M' comme estimateur de a et comparons par la suite ces deux riques quadratique afin de savoir quel estimateur de a est meilleur entre M et M' :
Risque quadratique de M’ comme estimateur de a et comparaison des deux estimateurs
#Nous trouvons un risque quadratique pour M' comme estimateur de a de E((M'-a)²)=a²/3n.
#Comme on a pu le voir dans la preuve mathématique le risque quadratique pour M comme estimateur de a est plus petit que le risque quadratique pour M' comme estimateur de a. Aisni l'estimateur M est meilleur que l'estimateur M'.set.seed(49)
#Exercice 2
#Pour savoir si Bob a eu raison de chosir r(M)=1.1, et si la stratégie choisie par Bob est une bonne stratégie, nous allons représenter l'espérence du Gain, donc de G en fonction de q où R(M)=q*M, et voir si c'est bien en q=1.1 que l'espérence du gain est maxiamle.
a=runif(1)# a correspond au nombre a choisi secrètement et uniformément par Alice dans l'intervalle [0,1].
n=10# n corrspond au nombre, de nombre que tire Alice dans l'intervalle [0,a]
nb=100000#nb correspond au nombre de fois où Alice tire n nombre dans [0,a]
#Nous créons une fonction EG qui prend q en paramètre et retourne l'espérence du gain G.
EG<-function(q){
TM=replicate(nb,max(runif(n)*a))#TM contient les nb maximuns des n nombres que tire Alice, nb fois dans [0,a].
TR=q*TM#TR va contenir les nb, R(M), calculés pour les nb maximums, c'est-à-dire les q*M, pour les différents maximums.
gain=0#Nous initialisons une variable gain à 0, qui sera incrémenter lors du parcours du tableau TR.
for (i in 1:nb) {
if(TR[i]<a){#Si R(M)<a, Bob gagne
gain=gain+(TR[i]-TM[i])#Bob gagne R(M)-M
}
#Dans tous les autres cas Bob ne gagne pas.
}
return((gain)/nb)#Nous divison ensuite ce gain par le nombre de partie, afin d'obtenir l'espérence du gain.
}
x=seq(from=0, to=4, by=0.1)
y=lapply(x,EG)
plot(x,y,xlab="q",ylab="E(G)")
#Je trace une asymptote verticale à q=1.1, pour que la figure soit plus parlante.
abline(v = 1.1, col="red", lwd=3, lty=2)
#On remarque que pour n=10, l'épérence du gain est à la fois maximale et positive. Ainsi nous pensons que la stratégie de Bob est une bonne stratégie.#Q2.4
set.seed(55)
n=3
nb=10000
#a=sample(size=1,x=seq(from=0,to=1, by=0.1),replace=T)
#J'ai écrit plusieurs fonctions qui retourne la probabilité P(A=a,M=m), ces fonctions je les ai écrite pour des M différents. Ainsi il suffit seulement de changer la valeur de M==0.? dans le if. Pour trouver la valeur de la probabilite P(A=a|M=m)
#Cette prmière fonction nous permet de conclure quant à l'estimation de la probabilité de P(A=a|M=0.5). Puisqu'il s'agit de la fonction qui calcul P(A=a,M=0.5), puisque P(A=a|M=0.5)=P(A=a,M=0.5)/P(M=0.5). Ainsi le but de Bob c'est d'augmenter ça probabilité de réussite. Donc d'augmenter P(A=a,M=0.5), puisque P(M=0.5) est une constante.
PAM5<-function(a){
compteur=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.5){
compteur=compteur+1
}
}
return (compteur/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM5)
plot(x,y,xlab="a",ylab="P(M=0.5)")
#Ainsi on peut voir que c'est pour A=0.5 que la probabilité P(A=a|M=0.5) est maximale, donc je proposerai à Alice la valeur A=0.5.
PAM8<-function(a){
compteur=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.8){
compteur=compteur+1
}
}
return (compteur/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM8)
plot(x,y,xlab="a",ylab="P(M=0.8)")
PAM7<-function(a){
compteur=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.7){
compteur=compteur+1
}
}
return (compteur/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM7)
plot(x,y,xlab="a",ylab="P(M=0.7)")
PAM9<-function(a){
compteur=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.9){
compteur=compteur+1
}
}
return (compteur/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM9)
plot(x,y,xlab="a",ylab="P(M=0.9)")
PAM3<-function(a){
compteur=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.3){
compteur=compteur+1
}
}
return (compteur/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM3)
plot(x,y,xlab="a",ylab="P(M=0.3)")
PAM4<-function(a){
compt=0
for(i in 1:nb){
M=max(sample(size=n,x=seq(from=0,to=a, by=0.1),replace=T))
if(M==0.4){
compt=compt+1
}
}
return (compt/nb)
}
x=seq(from=0, to=1, by=0.1)
y=lapply(x,PAM4)
plot(x,y,xlab="a",ylab="P(M=0.4)")
#Q2.5
set.seed(60)
n=2
nb=10000
a=runif(1)
alpha=1;
#Comme précedemment, on crée une fonction EG qui retourne l'espére de G, et qui prend alpha en paramètre afin de pouvoir EG en fonction de alpha et savoir si cet alpha est optimal.
EG<-function(alpha){
M=replicate(nb,max(runif(n)*a))
R=M^alpha
gain=0
for(i in 1:nb){
if(R[i]<a){
gain=gain+(R[i]-M[i])
}
}
return (gain/nb)
}
x=seq(from=0.1,to=5,by=0.1) # alpha qui va de 0.1 jusqu'à 5 avec un pas de 0.1
y=lapply(x,EG)
plot(x,y,xlab="alpha",ylab = "E(G)")
#On peut conclure que c'est une bonne stratégie en vue de la courbe puisque si il prend un alpha plus grand que 1, on peut voir que l'espérence du gain devient négative pour un alpha plus grand que 1, donc s'il prend alpha plus grand que 1, il perd de l'argent. Par contre pour alpha plus petit que 1, l'espérence est positive, donc c'est une stratégie gagnante donc bonne.
#la meilleur valeur de alpha parmi {0.3,0.5,0.7} est 0.5 puisque parmi les 3 c'est pour cette valeur de alpha que l'on obtient la plus grande valeur de l'esperance du gain.
#De plus d'après la courbe on ne trouve pas d'alpha meilleur que 0.5. En fait le meilleure alpha se trouve entre 0.5 et 0.6.#On refait les mêmes questions pour n=10:
set.seed(15)
n=10
#De même que précedemment on peut conclure que c'est une bonne stratégie en vue de la courbe puisque si il prend un alpha plus grand que 1, on peut voir que l'espérence du gain devient négative pour un alpha plus grand que 1, donc s'il prend alpha plus grand que 1 il perd de l'argent. Par contre pour alpha plus petit que 1, l'espérence est positive, donc c'est une stratégie gagnante donc bonne.
EG2<-function(alpha){#E(G)
M=replicate(nb,max(runif(n)*a))
R=M^alpha
gain=0
for(i in 1:nb){
if(R[i]<a){
gain=gain+(R[i]-M[i])
}
}
return (gain/nb)
}
x1=seq(from=0.1,to=5,by=0.1)
y1=lapply(x1,EG2)
plot(x1,y1,xlab="alpha",ylab = "E(G) avec n=10")
#la meilleur valeur de alpha parmi {0.3,0.5,0.7} est 0.7 puisque parmi ces 3 valeurs c'est pour alpha=0.7 qu'on obtient la plus grande valeur de l'esperance.
#Après analyse de la courbe nous trouvons que le meilleur alpha qui est proche de 0.9 c'est-à-dire entre 0.9. Il se trouve en réalité entre 0.9 et 1.
#Finalement empiriquement nous pouvons conclure que plus n devient grand, plus le alpha donnant la meilleure stratégie, se rapproche de 1. #Exercie 3 Bonus.
#Q3.6
#Ici la réponse est assez simple si Bob répond toujours "oui". Alors si Alice fait Pile Bob à gagner, et si Alice fait Face Bob à perdu. Ainsi la probabilité de Bob de gagner est tout simplement la probabilité que Alice face pile au lancer de Pile ou face. C'est à dire que la probabilité de gagner pour Bob est de P(G)=1/2.