set.seed(10293482) # Rajouté à la toute fin pour figer les valeurs et une fois que je m'étais assuré que j'avais suffisemment d'échantillons pour avoir des estimations stables.

Exercice 1: Question préliminaire à propos d’estimation

Alice choisit \(a>0\) puis tire \(X_1, ..., X_n\) uniformément dans \([0,a]\).

Q1.1

On s’intéresse à la loi de \(M=\max_i X_i\). Réalisons une fonction permettant de simuler \(M\).

rM = function(size=1, a=1.5, n=10) {
  M = rep(x = 0, size = size);
  for(i in 1:n) {
    X=runif(n = size, min=0, max=a);
    M = pmax(M,X);
  }
  return(M);
}
rM(size = 20, a=1.7, n=10)
##  [1] 1.479919 1.434738 1.559035 1.633652 1.545177 1.646887 1.576863
##  [8] 1.644450 1.577479 1.550442 1.621595 1.567799 1.676014 1.508064
## [15] 1.547348 1.584508 1.693927 1.412147 1.535669 1.691097

Q1.1.a

Les échantillons de \(M\) sont effectivement proches de \(a\) mais ils sont toujours inférieurs. \(M\) est donc un estimateur biaisé (son espérance est trop faible) qu’il va falloir corriger.

Q1.1.b

A posteriori, il y avait plus direct que ce que je réponds ici mais je vous livre volontairement mon raisonnement dans l’ordre où j’y ai pensé en répondant à cette question.

Pour estimer l’espérance de \(M\) en fonction de \(a\) et de \(n\), il suffit d’en tirer un “grand” nombre et d’en faire la moyenne. Je vais tout mettre dans une data frame pour pouvoir facilement faire des graphiques

df = data.frame();
for(a in c(1.5, 2.5, 4.5)) {
  for(n in 1:15) {
    df = rbind(df, data.frame(a=a, n=n, mean = mean(rM(size=100000,a=a,n=n))));
  }
}
df
##      a  n      mean
## 1  1.5  1 0.7506455
## 2  1.5  2 1.0010705
## 3  1.5  3 1.1255092
## 4  1.5  4 1.1994222
## 5  1.5  5 1.2492377
## 6  1.5  6 1.2854289
## 7  1.5  7 1.3120378
## 8  1.5  8 1.3338932
## 9  1.5  9 1.3498784
## 10 1.5 10 1.3636676
## 11 1.5 11 1.3750940
## 12 1.5 12 1.3842913
## 13 1.5 13 1.3926145
## 14 1.5 14 1.4006110
## 15 1.5 15 1.4066918
## 16 2.5  1 1.2537473
## 17 2.5  2 1.6679730
## 18 2.5  3 1.8781674
## 19 2.5  4 1.9995346
## 20 2.5  5 2.0827056
## 21 2.5  6 2.1443313
## 22 2.5  7 2.1865779
## 23 2.5  8 2.2237706
## 24 2.5  9 2.2496375
## 25 2.5 10 2.2723896
## 26 2.5 11 2.2906595
## 27 2.5 12 2.3072543
## 28 2.5 13 2.3211342
## 29 2.5 14 2.3331044
## 30 2.5 15 2.3426556
## 31 4.5  1 2.2443923
## 32 4.5  2 2.9986367
## 33 4.5  3 3.3775081
## 34 4.5  4 3.6018919
## 35 4.5  5 3.7492933
## 36 4.5  6 3.8572074
## 37 4.5  7 3.9363317
## 38 4.5  8 3.9990782
## 39 4.5  9 4.0498844
## 40 4.5 10 4.0916775
## 41 4.5 11 4.1233086
## 42 4.5 12 4.1533587
## 43 4.5 13 4.1789715
## 44 4.5 14 4.1999118
## 45 4.5 15 4.2179005
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
ggplot(data=df, aes(x=n, y=mean, color=as.factor(a))) + geom_point() + theme_bw()

On retrouve comme on pouvait s’y attendre que \(E[M]\) tend vers \(a\) quand \(n\) devient grand. Cette forme avec une asymptote rappelle quelque chose en \(1/n\). Si je multiplie par \(n\), j’aurai peut-être quelque chose de plus clair.

library(ggplot2)
ggplot(data=df, aes(x=n, y=mean*n+1, color=as.factor(a))) + geom_point() + theme_bw()

Parfait, du coup, on dirait bien que \(n*E[M]\) est linéaire et la pente vaut forcément \(a\), il faudrait que j’enlève \(a.n\) pour voir où ça coupe:

library(ggplot2)
ggplot(data=df, aes(x=n, y=mean*n-a*n, color=as.factor(a))) + geom_point() + theme_bw()

Ah, encore une forme en \(1/n\), c’est que je n’ai pas du multiplier par la bonne valeur. Mmmh, reprenons, ça tends vers \(a\) avec une forme en \(1/n\). Si c’est bien le cas, en enlevant \(a\) et en regardant l’inverse, je verrais si c’est bien \(n\) qu’il faut mettre au dénominateur ou autre chose. Au passage, vu que tout est proportionnel à \(a\), je peux normaliser en re-multipliant par \(a\):

library(ggplot2)
ggplot(data=df, aes(x=n, y=a/(a-mean), color=as.factor(a))) + geom_point() + theme_bw()

Ah, ben là, c’est clair, pour \(n=4\), j’ai \(5\), pour \(n=8\), c’est \(9\), etc. C’est donc un \(n+1\) qu’il y a dans ma formule! Vérifions ça:

library(ggplot2)
ggplot(data=df, aes(x=n, y=mean*(n+1)-a*(n+1), color=as.factor(a))) + geom_point() + theme_bw()

Autrement dit, on dirait bien que \((n+1)E[M] = a(n+1) -a = an\) soit \(E[M]=\frac{n}{n+1}a\). Regardons si ça marche “bien” ou pas:

library(ggplot2)
ggplot(data=df, aes(x=n, y=mean, color=as.factor(a))) + geom_point() + 
    geom_line(aes(y=a*n/(n+1))) + theme_bw()

C’est presque trop beau. :) Et numériquement:

df$err = df$mean - df$a*df$n/(df$n+1)
df$err
##  [1]  6.455425e-04  1.070543e-03  5.091695e-04 -5.778265e-04 -7.622793e-04
##  [6] -2.853429e-04 -4.621560e-04  5.598442e-04 -1.216065e-04  3.119151e-05
## [11]  9.398247e-05 -3.240475e-04 -2.426103e-04  6.109831e-04  4.417986e-04
## [16]  3.747253e-03  1.306358e-03  3.167356e-03 -4.653594e-04 -6.277283e-04
## [21]  1.474202e-03 -9.220890e-04  1.548357e-03 -3.624811e-04 -3.377199e-04
## [26] -1.007123e-03 -4.380010e-04 -2.943937e-04 -2.289106e-04 -1.094403e-03
## [31] -5.607714e-03 -1.363256e-03  2.508084e-03  1.891865e-03 -7.067164e-04
## [36]  6.450653e-05 -1.168258e-03 -9.217653e-04 -1.156059e-04  7.683966e-04
## [41] -1.691395e-03 -4.874592e-04  4.000923e-04 -8.822593e-05 -8.495332e-04

C’est pas mal du tout. J’ai tiré \(100000=10^6\) échantillons et, si ma formule est correcte, l’erreur est de l’ordre de \(10^{-3}\), une convergence en racine du nombre d’échantillons, c’est logique.

Du coup, pour débiaiser \(M\), il suffit de le multiplier par \(\frac{n+1}{n}\). Ainsi, l’espérance de ce nouvel estimateur sera exactement \(a\). Certains m’ont proposé de rajouter \(\frac{a}{n+1}\) à \(M\) pour le débiaiser mais ce n’est pas possible si on cherche justement à estimer \(a\) que l’on ne connaît pas…

Q1.1.c

Estimons la variance maintenant:

df = data.frame();
for(a in c(1.5, 2.5, 4.5)) {
  for(n in 1:20) {
    M = rM(size=10000,a=a,n=n)
    df = rbind(df, data.frame(a=a, n=n, mean = mean(M), var = var(M)));
  }
}
head(df)
##     a n      mean        var
## 1 1.5 1 0.7503672 0.18921178
## 2 1.5 2 0.9998454 0.12449968
## 3 1.5 3 1.1234256 0.08546866
## 4 1.5 4 1.2043782 0.05755376
## 5 1.5 5 1.2499419 0.04415615
## 6 1.5 6 1.2850600 0.03496370

L’espérance étant proportionnelle à \(a\), la variance sera proportionnelle à \(a^2\) et devrait décroître avec \(n\) mais à quelle vitesse ?

library(ggplot2)
ggplot(data=df, aes(x=n, y=(var/a^2), color=as.factor(a))) + geom_point() + theme_bw()

Pas facile à estimer. Supposons que \(Var(M_{a,n)})\approx \alpha.a^2/n^\beta\) (ça parait raisonnable vu l’allure de la courbe précédente). En regardant ce qu’il se passe pour des puissances de 2, je pourrai estimer \(\beta\) puis \(\alpha\).

df2 = data.frame();
for(a in c(1)) {
  for(n in c(1,2,4,8,16,32,64,128,256,512)) {
    M = rM(size=100000,a=a,n=n)
    df2 = rbind(df2, data.frame(a=a, n=n, mean = mean(M), var = var(M)));
  }
}
X = df2$var
X[-1]/X[-length(X)]
## [1] 0.6652875 0.4743428 0.3697313 0.3185609 0.2759765 0.2635159 0.2611709
## [8] 0.2532809 0.2564493

Le ratio entre deux puissances de 2 successive tend vers 1/4, autrement dit quand on double \(n\), la variance est divisée par \(4\) et on a donc probablement quelque chose en \(1/n^2\) mais c’est une approximation de mauvaise qualité au début. Estimons \(\alpha\), maintenant:

df2$var*df2$n^2
##  [1] 0.08350082 0.22220820 0.42161145 0.62353179 0.79453138 0.87708784
##  [7] 0.92450630 0.96581671 0.97849178 1.00373430

Ça tend vers \(1\) donc \(Var(M) \approx a^2/n^2\) semble une approximation correcte (pour \(n\) grand).

Regardons ce que ça donne pour de plus petites valeurs de \(n\) (les points représentent les estimations et les lignes l’approximation basée sur la formule précédente).

library(ggplot2)
ggplot(data=df, aes(x=n, y=var, color=as.factor(a))) + 
  geom_point() + geom_line(aes(y=a^2/n^2)) + ylim(0,1) + theme_bw()
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 7 rows containing missing values (geom_path).

Bon, ça n’est pas terrible mais \(a^2/n^2\) a quand même l’air d’être une borne supérieure. Ce que je viens de faire, là, c’est une étude empirique et je l’ai faite sans connaître la véritable formule. C’est assez incroyable le nombre de personnes qui m’ont dit “oui, on trouve facilement grâce à la formule de la variance que \(Var[M] = E[M^2]-E[M]^2 =\frac{na^2}{n+2} - \left(\frac{na}{n+1}\right)^2\)”. Mais d’où est-ce que c’est sorti ?!? Ça n’a rien d’évident.

Q1.2

On s’intéresse maintenant à l’estimateur basé sur la moyenne. L’espérance de chaque \(X_i\) est \(a/2\) et donc l’espérance de \(M'\) est exactement \(a\) pour tout \(n\) grace à la la linéarité de l’espérance.

rM2 = function(size=1, a=1.5, n=10) {
  M = rep(x = 0, size = size);
  for(i in 1:n) {
    X=runif(n = size, min=0, max=a);
    M = M+X;
  }
  return(2*M/n);
}

df = data.frame();
for(a in c(1.5, 2.5, 4.5)) {
  for(n in 1:20) {
    M = rM(size=10000,a=a,n=n)
    M2 = rM2(size=10000,a=a,n=n)
    df = rbind(df, data.frame(a=a, n=n, mean = mean(M), var = var(M), mean2 = mean(M2), var2 = var(M2)));
  }
}
head(df)
##     a n      mean        var    mean2      var2
## 1 1.5 1 0.7483912 0.19073088 1.506643 0.7397509
## 2 1.5 2 0.9980844 0.12321513 1.506892 0.3754922
## 3 1.5 3 1.1241895 0.08451639 1.503398 0.2468918
## 4 1.5 4 1.2012198 0.05914434 1.503552 0.1867644
## 5 1.5 5 1.2493703 0.04448315 1.502696 0.1491913
## 6 1.5 6 1.2862442 0.03427571 1.495475 0.1262828

Vérifions que \(\frac{n+1}{n}M\) (les points) et \(M'\) (les lignes) sont de bons estimateur de \(a\).

library(ggplot2)
ggplot(data=df, aes(x=n, color=as.factor(a))) + geom_point(aes(y=mean*(n+1)/n)) + geom_line(aes(y=mean2))+ theme_bw()

Regardons leur variance:

library(ggplot2)
ggplot(data=df, aes(x=n, color=as.factor(a))) + geom_point(aes(y=var*((n+1)/n)^2)) + geom_line(aes(y=var2)) + theme_bw()

La variance de \(M'\) est systématiquement (enfin, sauf pour \(n=1\)…) bien plus grande que celle de \(M\). En fait, on peut aisément la calculer: la variance de \(X_i\) est \(\frac{a^2}{12}\) et donc la variance de \(M'\) est \(\frac{a^2}{3n}\) (on multiplie la somme par \(2\) donc la variance par \(4\)), ce qui est bien plus grand que le \(a^2/n^2\) que j’ai estimé précédemment. \(\frac{n+1}{n}M\) est donc un bien meilleur estimateur de \(a\) que \(M'\) puisque sa variance est bien plus faible.

Exercice 2: Un deuxième jeu à base de max

Commençons par écrire une fonction simulant ce que fait Alice:

Alice = function(size=1000, n=10) {
  a = runif(n = size);
  M = rep(0, size=size)
  for(i in 1:n) {
    U = runif(n = size);
    X = a*U; # pour avoir un nombre uniforme entre 0 et a
    M = pmax(M,X)
  }
  return(data.frame(a=a, M=M))
}
Alice(size=10, n=10)
##             a          M
## 1  0.98008174 0.87084437
## 2  0.36054854 0.34838306
## 3  0.92937986 0.85670980
## 4  0.80496263 0.80305321
## 5  0.42201378 0.37516639
## 6  0.08774429 0.08497445
## 7  0.32229151 0.29256453
## 8  0.06604859 0.06374569
## 9  0.70236379 0.59338853
## 10 0.48791369 0.41387302

Q2.3: Une stratégie naïve

Bob sachant comment débiaiser le \(M\) pour avoir un bon estimateur se propose de répondre \(r(M)=\frac{n+1}{n}M\). Cette réponse sera par contre assez souvent bien supérieure à \(a\) pour être vraiment bonne. En particulier, c’est le cas quand \(M\) est supérieur à \(1/1.1=0.9090..\) car on propose alors une valeur supérieure à \(1\)… Dans le graphique ci dessous, on n’obtient un gain que lorsque l’on répond entre les deux lignes en pointillé et la ligne noire représente la réponse proposée dans cette question. Pour les grandes valeurs de \(M\), on ne gagnera rien.

ggplot(data=data.frame()) + xlim(0,1) + ylim(0,1.1) + theme_classic() +
    geom_abline(slope=1,intercept = 0, linetype=2) + 
    geom_abline(slope=1.1,intercept = 0, linetype=1)+
    geom_hline(yintercept = 1, linetype=2) + xlab("M") + ylab("r(M)")

Bob1 = function(M,n=10) {
  return((n+1)/n*M)
}

Gain = function(size=1000,n=10, Bob) {
  game = Alice(size, n);
  R = Bob(game$M)
  ifelse(R>game$a,0,R-game$M)
}

Gain(size=100, n=10, Bob1)
##   [1] 6.191501e-02 0.000000e+00 1.110856e-02 2.262995e-02 0.000000e+00
##   [6] 0.000000e+00 1.463039e-02 0.000000e+00 6.627482e-02 0.000000e+00
##  [11] 8.584771e-02 0.000000e+00 0.000000e+00 0.000000e+00 5.004571e-02
##  [16] 2.560146e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 6.625065e-02
##  [26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [31] 0.000000e+00 0.000000e+00 2.439036e-02 0.000000e+00 0.000000e+00
##  [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [41] 0.000000e+00 3.127988e-02 0.000000e+00 0.000000e+00 6.797084e-02
##  [46] 0.000000e+00 0.000000e+00 1.637388e-02 5.992927e-02 0.000000e+00
##  [51] 0.000000e+00 6.988380e-03 5.158463e-02 2.721865e-02 0.000000e+00
##  [56] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 7.756021e-02
##  [61] 0.000000e+00 0.000000e+00 0.000000e+00 2.864102e-02 0.000000e+00
##  [66] 5.882687e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [71] 0.000000e+00 0.000000e+00 0.000000e+00 6.565391e-02 6.979059e-02
##  [76] 0.000000e+00 3.737816e-02 0.000000e+00 6.243301e-02 0.000000e+00
##  [81] 1.541149e-02 0.000000e+00 3.623192e-02 1.801094e-02 1.390667e-02
##  [86] 4.089539e-02 8.000225e-02 6.967415e-02 1.290431e-02 0.000000e+00
##  [91] 3.900927e-03 0.000000e+00 0.000000e+00 7.103273e-05 0.000000e+00
##  [96] 0.000000e+00 8.256991e-05 0.000000e+00 0.000000e+00 1.227643e-02

Le gain avec cette stratégie est d’environ:

mean(Gain(size=100000, n=10, Bob1))
## [1] 0.01581778

Q2.4 Simplifions le jeu

AliceDiscrete = function(size=1000, n=3) {
  A = sample(size = size, x=.1*(0:10), replace=T);
  M = rep(0, size=size)
  for(i in 1:n) {
    U = runif(n = size);
    X = .1*floor((10*A+1)*U);  
    M = pmax(M,X)
  }
  return(data.frame(A=A, M=M))
}
AliceDiscrete(size=10, n=3)
##      A   M
## 1  0.1 0.1
## 2  0.4 0.3
## 3  0.2 0.2
## 4  0.6 0.5
## 5  0.1 0.1
## 6  0.5 0.4
## 7  0.9 0.4
## 8  0.5 0.1
## 9  0.3 0.2
## 10 0.4 0.3

Estimons \(P[A=a,M=m]\)

Estimons \(P[A=a,M=m]\) en simulant un grand nombre de fois et en comptant chaque évènement.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
N=100000
AliceDiscrete(size=N, n=3) %>% group_by(A,M) %>% summarise(freq=n()/N)
## # A tibble: 66 x 3
## # Groups:   A [?]
##        A     M    freq
##    <dbl> <dbl>   <dbl>
##  1   0     0   0.0918 
##  2   0.1   0   0.0108 
##  3   0.1   0.1 0.0780 
##  4   0.2   0   0.00347
##  5   0.2   0.1 0.0235 
##  6   0.2   0.2 0.0642 
##  7   0.3   0   0.00127
##  8   0.3   0.1 0.00949
##  9   0.3   0.2 0.026  
## 10   0.3   0.3 0.0527 
## # … with 56 more rows

Voici les \(P[A|M=m]\) pour chaque valeur possible de \(m\). La valeur de \(A\) la plus probable est à chaque fois \(m\) mais le gain associé est nul donc ce n’est pas ce qu’il faut répondre…

N=100000
AD = AliceDiscrete(size=N, n=3);
AD %>% group_by(A,M) %>% summarize(count=n()) %>% ungroup() %>% mutate(freq=count/sum(count)) -> AD
AD %>% group_by(M) %>% mutate(freq=freq/sum(freq)) -> PcondM;

ggplot(data=PcondM, aes(x=A, y=freq, color=as.factor(M))) + geom_point() + geom_line() + theme_classic()

Intéressons nous au cas où Alice nous a indiqué que \(M\) valait \(0.5\). Je commence à conditionner par les situations où \(M=0.5\) pour renormaliser et obtenir \(P[A=a|M=0.5]\):

library(dplyr)
N=100000
AliceDiscrete(size=N, n=3) %>% filter(M==0.5) %>% group_by(A,M) %>% summarise(count=n()) -> BobDiscrete
N2 = sum(BobDiscrete$count)
BobDiscrete %>% group_by(A,M) %>% summarise(freq=count/N2) -> BobDiscrete
BobDiscrete
## # A tibble: 6 x 3
## # Groups:   A [?]
##       A     M   freq
##   <dbl> <dbl>  <dbl>
## 1   0.5   0.5 0.366 
## 2   0.6   0.5 0.236 
## 3   0.7   0.5 0.149 
## 4   0.8   0.5 0.110 
## 5   0.9   0.5 0.0769
## 6   1     0.5 0.0616

Si Alice m’indique \(M=0.5\) et que je répond \(0.5\), mon gain sera évidemment nul. Si je répond \(0.6\), je ne gagnerai quelque chose (\(0.1\)) que si \(0.6\leq A\), ce qui arrive avec probabilité \(\approx 1-0.37=0.62\). Mon espérance de gain en répondant \(0.6\) est donc de \(0.062\). Quel gain puis-je espérer selon la valeur répondue:

G0.5 = data.frame();
for(r in BobDiscrete$A) {
  G0.5 = rbind(G0.5, data.frame(r=r, gain = sum(ifelse(r>BobDiscrete$A,0,(r - BobDiscrete$M)* BobDiscrete$freq))));
}
G0.5
##     r       gain
## 1 0.5 0.00000000
## 2 0.6 0.06338803
## 3 0.7 0.07947876
## 4 0.8 0.07459459
## 5 0.9 0.05540541
## 6 1.0 0.03079151

Mon gain sera maximal en répondant \(0.7\). Regardons ce que ça donne pour toute valeur de \(M\):

N=100000
AD = AliceDiscrete(size=N, n=3);
G = data.frame();
for(m in unique(AD$M)) {
  AD %>% filter(M==m) %>% group_by(A,M) %>% summarise(count=n()) %>% ungroup() %>% mutate(freq=count/sum(count)) -> BD;
  GBD = data.frame() # Je stoque les gains possibles quand M=m 
  for(r in unique(BD$A)) {
     GBD = rbind(GBD, data.frame(M=m, R=r, Gain = sum(ifelse(r>BD$A,0,(r - BD$M)* BD$freq))));
  }
  G = rbind(G, GBD %>% filter(Gain==max(Gain))) # Je ne conserve que la meilleure stratégie
}
G %>% arrange(M)
##      M   R       Gain
## 1  0.0 0.1 0.01689789
## 2  0.1 0.2 0.03755180
## 3  0.2 0.4 0.05622737
## 4  0.3 0.5 0.07025468
## 5  0.4 0.6 0.07814685
## 6  0.5 0.7 0.08030676
## 7  0.6 0.8 0.07783961
## 8  0.7 0.9 0.06961913
## 9  0.8 0.9 0.05528508
## 10 0.9 1.0 0.04300138
## 11 1.0 1.0 0.00000000

Regardons à quoi ressemble la stratégie optimale en fonction de M:

ggplot(data=G, aes(x=M, y=R)) + geom_point() + ylim(0,1) + geom_line(aes(y=M)) + theme_classic()

Bon difficile de dire. Ça doit forcément rester au dessus de M et se terminer en 1 mais la discrétisation rend assez difficile de repérer une forme simple. Peut-être quelque chose du genre \(\sqrt{M}\)…

Q2.5 Retour au cas général

En répondant \(M^\alpha\) pour \(\alpha>1\) on a une stratégie qui renvoie toujours une valeur relativement proche de (mais plus grande que) \(M\) et toujours valide car plus petite que \(1\). Petite illustration:

ggplot(data=data.frame(x = c(0, 1)),aes(x)) + xlim(0,1) + ylim(0,1.1) + theme_classic() +
    geom_abline(slope=1,intercept = 0, linetype=2) + 
    stat_function(fun=function(M) {M**0.5}, geom="line") +
    geom_hline(yintercept = 1, linetype=2) + xlab("M") + ylab("r(M)")

C’est a priori un peu plus intelligent que renvoyer \(\frac{n+1}{n}M\), en particulier pour les grandes valeurs de \(M\). Le fait que ça soit assez optimiste avec des valeurs assez supérieures à \(M\) pour les valeurs entre 0.1 et 0.8 semble une bonne idée mais c’est peut-être un peu trop optimiste. Reste à trouver le bon \(\alpha\) (en fonction de \(n\)) et à voir si on obtient un gain substantiel…

Bob2 = function(M, alpha=.5) {
  return(M**alpha)
}

Gain = function(size=1000,n=10, Bob) {
  game = Alice(size, n);
  R = Bob(game$M)
  ifelse(R>game$a,0,R-game$M)
}

Galpha = data.frame()
for(alpha in (0:20)/20) {
  Galpha = rbind(Galpha, data.frame(alpha=alpha, Gain=mean(Gain(size=10000, n=2, function(M) {Bob2(M,alpha=alpha)}))))
}
ggplot(Galpha, aes(x=alpha, y=Gain)) + geom_point() + theme_bw()

Et bien, pour n=2, il semble que renvoyer \(\sqrt{M}\) soit une assez bonne stratégie.

Regardons ce que ça donne pour d’autres valeurs de \(n\):

Galpha = data.frame()
for(n in 1:10) {
  for(alpha in (1:20)/20) {
    Galpha = rbind(Galpha, data.frame(n=n, alpha=alpha, Gain=mean(Gain(size=100000, n=n, function(M) {Bob2(M,alpha=alpha)}))))
  }
}
ggplot(Galpha, aes(x=alpha, y=Gain, color=as.factor(n))) + geom_point() + geom_line() + theme_bw()

Plus \(n\) est grand plus il faut un \(\alpha\) proche de \(1\), ce qui est assez logique. Regardons est le meilleur \(\alpha\) (par incrément de \(0.1\)) en fonction de \(n\):

Galpha %>% group_by(n) %>% filter(Gain==max(Gain))
## # A tibble: 10 x 3
## # Groups:   n [10]
##        n alpha   Gain
##    <int> <dbl>  <dbl>
##  1     1  0.35 0.105 
##  2     2  0.5  0.0666
##  3     3  0.6  0.0479
##  4     4  0.65 0.0365
##  5     5  0.7  0.0293
##  6     6  0.7  0.0246
##  7     7  0.75 0.0208
##  8     8  0.8  0.0181
##  9     9  0.8  0.0160
## 10    10  0.8  0.0145

Pour \(n=10\), c’est \(\alpha\approx0.8\) qui donne le meilleur gain. Comparons avec la stratégie naïve du début:

mean(Gain(size=1000000, n=3, Bob1))
## [1] 0.02562411
mean(Gain(size=1000000, n=3, function(M) {Bob2(M,alpha=0.6)}))
## [1] 0.04739925
mean(Gain(size=1000000, n=10, Bob1))
## [1] 0.01593749
mean(Gain(size=1000000, n=10, function(M) {Bob2(M,alpha=0.8)}))
## [1] 0.01432032

Bon, pour \(n\) petit, on gagne beaucoup, mais pour \(n=10\), ça semble moins bon. C’est probablement dû au fait que pour les petites valeurs de \(M\), la réponse \(M^{0.8}\) n’est pas assez élevée pour amener à un gain substantiel. Ça pourrait aussi être dû au fait que \(0.8\) n’est pas si bon (après tout, on n’a testé que des multiples de \(0.1\)) mais mes tentatives pour des valeurs plus élevées ou plus faibles n’ont pas été concluante. C’est donc probablement un problème de forme.

mean(Gain(size=1000000, n=10, function(M) {Bob2(M,alpha=0.78)}))
## [1] 0.01426231
mean(Gain(size=1000000, n=10, function(M) {Bob2(M,alpha=0.82)}))
## [1] 0.01438536

Exercice 3: Question Bonus

Q3.6

En répondant systématiquement oui, Bob gagne si Alice lui renvoie la plus grande valeur, i.e., si la pièce est tombée sur Pile. Il a donc une chance sur deux de gagner et il est difficile de casser cette symétrie.

Q3.7

Intuitivement, si Alice renvoie une valeur élevée, il y a de bonnes chances pour que ça soit le maximum et si la valeur est basse, il y a de bonnes chances pour que ça ne soit pas la bonne. Bob pourrait donc choisir de répondre Vrai si la valeur est supérieure à \(0.5\) et Faux sinon. Évidemment, si Alice est au courant de cette stratégie, elle choisira volontairement des valeurs toutes deux élevées (auquel cas, Bob se trompera avec une chance sur deux) ou toutes deux basses (auquel cas, Bob se trompera avec une chance sur deux). En revanche, si Alice choisit 0.4 et 0.6, Bob aura toujours raison…

Mettons qu’Alice tire deux nombres au hasard uniformément dans \([0,1]\). S’ils sont tous deux plus petits que 0.5 ou plus grands que 0.5 (cas 1), Bob gagne avec une chance sur deux mais si l’un est plus petit que 0.5 et l’autre plus grand que 0.5 (cas 2) alors Bob gagne à tous les coups. La probabilité d’être dans chacun des 2 cas est d’une chance sur deux (faites un dessin avec \(x\in[0,1]\) et \(y\in[0,1]\) et regardez dans quelle situation l’un des deux est plus petit que \(0.5\) et l’autre plus grand que \(0.5\)). La probabilité de gain de Bob dans ces conditions est donc de \(3/4>1/2\).

Comment faire quand Alice choisit ses deux nombres? Pour que Bob gagne, il faut que le seuil qu’il utilise soit compris entre ces deux nombres. Il suffit que Bob choisisse son seuil au hasard, uniformément dans \([0,1]\) par exemple. La probabilité pour que Bob gagne à coup sûr est alors \(a_2-a_1\) et sa probabilité de gain est donc \(\frac{1}{2}+\frac{1}{2}(a_2-a_1)\). On retrouve bien le fait que si \(a_2\) est arbitrairement proche de \(a_1\), la probabilité de gain est d’une chance sur deux alors que si \(a_1=0\) et \(a_2=1\), Bob gagne à tous les coups.

Encore plus surprenant avec cette dernière solution, on n’utilise nulle part le fait que les valeurs doivent être comprises entre \(0\) et \(1\). Alice pourrait choisir deux nombres arbitrairement dans \(\mathbb{R}\). Si Bob tire un seuil selon une loi dont la densité est strictement positive en tout point (par exemple une loi normale), il augmentera strictement (bon, pas de grand chose, hein ?) sa probabilité de gain.