[RICM4-PS] DM
Samuel COURTHIAL
05/12/2019
Exercice 1 : Question préliminaire à propos d’estimation
Q1.1
A première vue, M est bien un estimateur de a. Je suppose que en tirant suffisament de \(X_{i}\) , la valeur maximale de ces derniers se rapprochera de a.
set.seed(0)
a = 10
n = 10
X = runif(n,0,a)
print(X)## [1] 8.966972 2.655087 3.721239 5.728534 9.082078 2.016819 8.983897
## [8] 9.446753 6.607978 6.291140M = max(X)
print(M)## [1] 9.446753Le code précédent permet d’obtenir une valeur de M. Si nous souhaitons approximer l’espérance de M, la loi des grands nombres nous dit que pour une valeur n suffisament grande, la moyenne empirique de n valeurs de M nous donne une bonne approximation de son espérance.
a = 10
n = 10
m = 100000
M = c(0)
i = 0
for (i in 1:m) {
X = runif(n,0,a)
M[i] = max(X)
}
print(mean(M))## [1] 9.089103Nous pouvons tenter de faire varier a pour déterminer son impact sur l’espérance de M.
n = 10
M_a = c(0)
for (a in 1:100) {
m = 1000
M = c(0)
i = 0
for (i in 1:m) {
X = runif(n,0,a)
M[i] = max(X)
}
M_a[a] = mean(M)
}
plot(x=1:100,y=M_a)
print((M_a[100]-M_a[10])/90)## [1] 0.90373On observe une influence linéaire de a avec E(X) = 0.9*a
a = 10
M_n = c(0)
for (n in 1:100) {
m = 1000
M = c(0)
i = 0
for (i in 1:m) {
X = runif(n,0,a)
M[i] = max(X)
}
M_n[n] = mean(M)
}
plot(x=1:100,y=M_n)
Il me semble logique que le 0,9 trouvé précédemment dépend en fait de n, comme cela avait était intuité au début de l’exercice. On observe qu’il s’agit d’une courbe croissante tendant vers a. Cette courbe pourrait correspondre à une fonction de la forme \(a-\frac{1}{x}\). En sachant que pour n = 0, on a E(M) = 0, on peut supposer que \[E(M) = (1-\frac{1}{n+1}) \times a\] On remarque que cette espérance tend bien vers a quand n tend vers l’infini. Pour avoir E(M) = a, il faut que l’on enlève \(\frac{a}{n+1}\) à l’espérance actuelle. On peut donc prendre \[M = max(X_i) + \frac{a}{n+1}\]
Pour estimer la variance, on peut utiliser le code suivant
a = 10
n = 10
V = c(0)
i = 0
for (i in 1:m) {
X = runif(n,0,a)
V[i] = (max(X) + (a/n+1) - a)**2
}
print(mean(V))## [1] 1.900644La variance étant une espérance d’une variable bien spécifique (\((M - E(M))²\)), la méthode utilisée plus haut semble appropriée.
Q1.2
Commençons par faire varier n
a = 10
m = 100000
M = c(0)
i = 0
for (n in 1:10) {
for (i in 1:m) {
X = runif(n,0,a)
M[i] = (2/n)*sum(X)
}
print(mean(M))
}## [1] 10.01499
## [1] 9.997676
## [1] 10.00287
## [1] 9.999698
## [1] 9.99341
## [1] 9.994702
## [1] 10.00573
## [1] 9.989586
## [1] 9.996735
## [1] 9.988793On n’observe pas de réelle correlation entre E(M’) et n. Faisons varier a
n = 10
m = 100000
M = c(0)
i = 0
for (a in 1:10) {
for (i in 1:m) {
X = runif(n,0,a)
M[i] = (2/n)*sum(X)
}
print(mean(M))
}## [1] 0.9998151
## [1] 1.997724
## [1] 2.999171
## [1] 4.00001
## [1] 4.99395
## [1] 5.999819
## [1] 7.001686
## [1] 8.00019
## [1] 8.996291
## [1] 10.00651Cette fois-ci, on remarque une correlation entre E(M’) et a. M’ semble être un autre estimateur de a. En effet, en sachant que l’espérance d’une loi uniforme continue entre a et b est \(\frac{a+b}{2}\), \[ E(M') = \frac{2}{n} E(\sum_{i} X_i) \\ = \frac{2}{n} \sum_{i} E(X_i) (linéarité) \\ = \frac{2}{n} \sum_{i} \frac{a - 0}{2} \\ = \frac{2na}{2n} \\ = a\] On retrouve bien que M’ est un estimateur de a.
On peut estimer la variance de M’ de la même manière que dans l’exercice 1
a = 10
n = 10
V = c(0)
i = 0
for (i in 1:m) {
X = runif(n,0,a)
V[i] = ((2/n)*sum(X) - a)**2
}
print(mean(V))## [1] 3.333333On observe une variance empirique en moyenne 2 fois plus grande. Cela semble indiquer que M’ est un moins bon estimateur que M corrigé.
Exercice 2 : Un deuxième jeu à base de max
Q2.3
L’idée de choisir une valeur 1.1 fois plus grande que M me semble une bonne idée. Avec 10 tirages, le maximum de ces tirages devrait être généralement assez proche de A. Cependant, ce choix pour r(M) donne une valeur supérieur à 1 si M est trop grand (ce qui risque d’arriver assez souvent, compte tenu de l’espérance de M), ce qui n’est pas forcément pertinent.
Pour déterminer si cette stratégie est bonne, il faut la comparer à d’autre startégie. Dans mon cas, je vais comparer cette stratégie avec les autre stratégies linéaires selon M. Cela me donnera seulement la pertinence du coefficient devant M mais permettra d’effectuer un premier classement. Je vais donc la simuler, calculer le gain et faire de même avec plusieurs startégie de même type. On sait que prendre \(r(M) \leq M\) assure la défaite. Il n’y a ensuite pas de limite supérieur pour r(M).
M = c(0)
G = c(0)
EG = c(0)
n = 10
m = 1000
i = 0
for (j in 1:100) {
coef = 1 + j / 100
for (i in 1:m) {
A = runif(1,0,1)
X = runif(n,0,A)
M = max(X)
r = coef * M
if (r <= A) {
G[i] = r - M
} else {
G[i] = 0
}
}
EG[j] = mean(G)
}
plot(x=1:100, y=EG)
On remarque que l’estimation de l’espérance de gain est maximale pour une valeur de j proche de 10, impliquant que \(r(M) = 1.1 \times M\) est une stratégie très efficace. Pour cette stratégie, l’éspérance de gain estimé est
print(EG[10])## [1] 0.01524569Q2.4
Les probabilités sont présentés dans une matrice
E = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
M = c(0)
G = c(0)
P = matrix(data=0,nrow=11,ncol=11)
rownames(P) = c("P(A=0)", "P(A=1)", "P(A=2)", "P(A=3)", "P(A=4)", "P(A=5)", "P(A=6)", "P(A=7)", "P(A=8)", "P(A=9)", "P(A=10)")
colnames(P) = c("P(M=0)", "P(M=1)", "P(M=2)", "P(M=3)", "P(M=4)", "P(M=5)", "P(M=6)", "P(M=7)", "P(M=8)", "P(M=9)", "P(M=10)")
n = 3
nbTimes = 100000
i = 0
for (i in 1:nbTimes) {
A = sample(E,1)
X = sample(E[E <= A],n,replace=TRUE)
M = max(X)
P[10*A+1,10*M+1] = P[10*A+1,10*M+1] + 1
}
P = P/nbTimes
sum(P)## [1] 1print(P)## P(M=0) P(M=1) P(M=2) P(M=3) P(M=4) P(M=5) P(M=6) P(M=7)
## P(A=0) 0.09162 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## P(A=1) 0.01115 0.07878 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## P(A=2) 0.00346 0.02332 0.06414 0.00000 0.00000 0.00000 0.00000 0.00000
## P(A=3) 0.00144 0.00981 0.02662 0.05263 0.00000 0.00000 0.00000 0.00000
## P(A=4) 0.00069 0.00548 0.01351 0.02712 0.04442 0.00000 0.00000 0.00000
## P(A=5) 0.00040 0.00257 0.00786 0.01593 0.02559 0.03835 0.00000 0.00000
## P(A=6) 0.00025 0.00178 0.00519 0.00948 0.01582 0.02389 0.03357 0.00000
## P(A=7) 0.00021 0.00120 0.00359 0.00631 0.01063 0.01624 0.02256 0.02956
## P(A=8) 0.00012 0.00095 0.00236 0.00482 0.00785 0.01120 0.01557 0.02077
## P(A=9) 0.00010 0.00058 0.00172 0.00365 0.00563 0.00814 0.01104 0.01502
## P(A=10) 0.00006 0.00040 0.00139 0.00240 0.00447 0.00613 0.00927 0.01165
## P(M=8) P(M=9) P(M=10)
## P(A=0) 0.00000 0.00000 0.00000
## P(A=1) 0.00000 0.00000 0.00000
## P(A=2) 0.00000 0.00000 0.00000
## P(A=3) 0.00000 0.00000 0.00000
## P(A=4) 0.00000 0.00000 0.00000
## P(A=5) 0.00000 0.00000 0.00000
## P(A=6) 0.00000 0.00000 0.00000
## P(A=7) 0.00000 0.00000 0.00000
## P(A=8) 0.02768 0.00000 0.00000
## P(A=9) 0.01970 0.02467 0.00000
## P(A=10) 0.01607 0.01848 0.02294On sait que \(P(A = a, M = 0.5) = P(M) \times P(A = a | M = 0.5)\) d’où \(P(A = a | M = 0.5) = \frac{P(A = a, M = 0.5)}{P(M = 0.5)}\) A partir des estimations précédentes, on obtient pour chaque valeur de A (dans la ligne 1) une probabilité (dans la ligne 2)
P2 = matrix(data=0,nrow=2,ncol=11,byrow=TRUE)
i = 0
for (i in 1:11) {
P2[1,i] = i/10
P2[2,i] = P[i,6] * 10
}
print(P2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 0.1 0.2 0.3 0.4 0.5 0.6000 0.7000 0.8000 0.900 1.0000 1.1000
## [2,] 0.0 0.0 0.0 0.0 0.0 0.3835 0.2389 0.1624 0.112 0.0814 0.0613On remarque que la probabilité est maximale pour A = 0.5, ce qui représente donc le meilleur choix. Après avoir vérifié avec plusieurs autres valeur de M, il semble que la valeur la plus judicieuse est A = M.
Q2.5
Encore une fois, la stratégie semble pertinente. Pour un \(\alpha\) < 1, r(M) est un peu plus grand que M et ne dépasse jamais 1, contrairement à la stratégie de la question 2.3. En revanche, si A est petit, M le sera aussi, et alors r(M) sera assez grand relativement à M. Si A est assez grand (proche de 1), M peut aussi l’être et alors r(M) varie peu par rappot à M. Autrement dit, si A est petit, on prend un r plus éloigné de M que si A est grand, contrairement à la stratégie de la question 2.3. Ce choix me parait étrange car si l’intervale de tirage augmente en taille, il me semble plus cohérent de prendre un r plus éloigné du M trouvé, à nombre de tirages égal.
Simulons cette situation pour n = 2 et pour \(\alpha \in \{0.7, 0.5, 0.3\}\)
M = c(0)
G = c(0)
n = 2
m = 10000
i = 0
j = 0
alphas = c(0.7,0.5,0.3)
for (j in 1:3) {
alpha = alphas[j]
for (i in 1:m) {
A = runif(1,0,1)
X = runif(n,0,A)
M = max(X)
r = M**alpha
if (r <= A) {
G[i] = r - M
} else {
G[i] = 0
}
}
print(mean(G))
}## [1] 0.05606518
## [1] 0.06604396
## [1] 0.05493837On remarque que le gain est maximum pour \(\alpha = 0.5\)
Généralisons
M = c(0)
G = c(0)
EG = c(0)
n = 2
m = 1000
i = 0
j = 0
for (j in 1:100) {
alpha = j / 101
for (i in 1:m) {
A = runif(1,0,1)
X = runif(n,0,A)
M = max(X)
r = M**alpha
if (r <= A) {
G[i] = r - M
} else {
G[i] = 0
}
}
EG[j] = mean(G)
}
plot(x=1:100, y=EG)
Sur cette courbe, on peut voir l’espérance de gain empirique en fonction de \(\alpha\) sur [0,1]. On remarque que cette courbe présente un maximum de gain pour \(\alpha = 0.5\) environ.
Mettons maintenant n = 10
M = c(0)
G = c(0)
EG = c(0)
n = 10
m = 1000
i = 0
j = 0
for (j in 1:100) {
alpha = j / 101
for (i in 1:m) {
A = runif(1,0,1)
X = runif(n,0,A)
M = max(X)
r = M**alpha
if (r <= A) {
G[i] = r - M
} else {
G[i] = 0
}
}
EG[j] = mean(G)
}
plot(x=1:100, y=EG)
On observe assez clairement que le maximum s’est décalé vers \(\alpha = 0.8\) environ.