RICM4/PS 2014: Images Panini
Arnaud Legrand
14/10/2014
L’objectif de ce devoir était d’écrire une petite simulation autour du célèbre problème du coupon collector
Approche expérimentale
Q1) Simulation
Commençons par écrire la fonction qui permet de simuler le remplissage de l’album et au bout de combien de temps la collection est terminée.
library(ggplot2)
library(plyr)
set.seed(42);
new_collection <- function(album_size=100) {
album=rep_len(0, album_size)
nb_images = 0;
time = 0;
while(nb_images<album_size) {
time = time + 1
idx = floor(runif(1,1, album_size+1))
if(album[idx]==0) {
album[idx] = 1
nb_images = nb_images + 1
}
}
time
}Q2) Étude
Génération d’un ensemble de populations
Réalisons maintenant un ensemble de simulations et rassemblons tout ceci dans une data-frame. Pour chaque taille d’album entre 10 et 200 par pas de 20, je génère 100 scénarios différents et je calcule combien de temps a été nécessaire pour compléter l’album.
df=data.frame(M=c(),T=c())
for(M in seq(10,200,20)) {
for(samples in 1:100) {
df = rbind(df,data.frame(M=c(M),T=new_collection(M)))
}
}Lors de nos investigations initiales, nous avions choisi un plus grand nombre de scénarios et un échantillonnage plus fin du nombre d’images. Les temps de calculs étaient alors assez longs sans que les tendances soient plus claires pour autant. Le choix précédent nous a donc semblé un compromis raisonnable.
Étude de E(T)
Regardons la distribution de nos échantillons pour l’album de taille 150.
ggplot(data=df[df$M==150,], aes(x=T)) + geom_histogram(binwidth=50) + theme_bw() + geom_vline(xintercept=mean(df[df$M==150,]$T))
La variabilité est assez importante. La distribution a l’air relativement sympa (un seul mode). Comme on pouvait s’y attendre, c’est non symmétrique et pas mal tirée vers la droite (skew). Pas clair qu’une centaine d’échantillons soit suffisant pour avoir une bonne estimation de l’espérence.
Faisons une première visualisation de notre échantillon:
ggplot(data=df, aes(x=M,y=T)) + geom_point(alpha=.2) + theme_bw()
Comme on pouvait s’y attendre, le temps moyen nécessaire à collecter augmente avec le nombre d’images. La variabilité augmente également énormément et il serait probablement plus judicieux de considérer plus de scénarios quand le nombre d’images augmente. Calculons donc la moyenne et l’écart type pour chacune des tailles possibles:
df_stat = ddply(df,c("M"),summarise,mean=mean(T),var=var(T),se=2*sqrt(var(T))/sqrt(length(T)))
df_stat## M mean var se
## 1 10 30.41 143 2.391
## 2 30 121.68 1378 7.425
## 3 50 228.70 4315 13.138
## 4 70 337.24 7233 17.009
## 5 90 444.92 10938 20.917
## 6 110 565.59 15527 24.922
## 7 130 700.11 22078 29.718
## 8 150 798.14 32080 35.822
## 9 170 947.80 50689 45.029
## 10 190 1104.44 40484 40.241Puis ajoutons cette moyenne (en rouge) sur la courbe précédente.
ggplot(data=df, aes(x=M,y=T)) + geom_point(alpha=.2) +
geom_line(data=df_stat,aes(y=mean),color="red") +
theme_bw()
Tiens, on dirait que c’est linéaire en fait… Qu’est-ce que ça donne si j’essaye de coller une droite ?
E_T <- function(x) {x*5.5}
ggplot(data=df, aes(x=M,y=T)) + geom_point(alpha=.2) +
geom_line(data=df_stat,aes(y=mean),color="red") +
stat_function(fun = E_T, color="purple") +
theme_bw()
Mmmh, c’est presque linéaire… Ça croit un peu plus vite.
Modélisation
Q3) Hypothèses
L’énoncé suggère que les \(X_n\) sont des variables indépendantes de loi uniforme.
Q4) Calculs des lois de probabilités
Étude de \(Y_i\)
Lorsque l’on a déjà \(i\) vignettes, il reste \(M-i\) vignettes à obtenir et donc la probabilité que la prochaine vignette soit correcte est \(p_i=\frac{M-i} M\). Le nombre de vignettes nécessaires à acheter pour obtenir la suivante suit donc une loi géométrique de paramètre \(p_i\), c’est à dire \(P(Y_i=k)=(1-p_i)^{k-1}p_i\). De ceci, on déduit \[ E(Y_i) = \sum_{k=1}^{\infty} k.P(Y_i=k) = \frac 1 {p_i}=\frac M {M-i}.\] La variance s’obtient après calcul \[ Var(Y_i) = E(Y_i^2) - E(Y_i)^2 = \frac {1-p_i}{p_i^2}.\] Les variables \(Y_i\) sont indépendantes car la loi des tirages est uniforme et les tirages sont indépendants.
Étude de \(T\)
On a \(T(M)=\sum_{i=0}^{M-1} Y_i\), en prenant l’espérance on obtient \[ Esp(T(M)) = M \left(\frac 1 M + \frac 1 {M-1}+ \cdots +\frac 1 2 +1\right)\simeq M\log M \] Ou bien un peu plus précisément, en utilisant une bonne approximation de la somme de la série harmonique; \[ Esp(T(M)) = M (\log M + \gamma + O(1/M)) \] Le log étant petit, c’est pour ça que ça a l’air linéaire mais que mes tentatives pour trouver “la bonne droite” étaient infructueuses…
Faisons le même travail sur \(Var(T)\) car ce n’était visiblement pas clair pour tout le monde. On a \[ Var(T(M)) = Var(\sum_{i=0}^{M-1} Y_i) = \sum_{i=0}^{M-1} Var(Y_i) \] car les \(Y_i\) sont indépendants deux à deux. On a donc: \[ Var(T(M)) = \sum_{i=0}^{M-1} \frac{1-p_i}{p_i^2} = \sum_{i=0}^{M-1} \frac{1}{p_i^2} - \sum_{i=0}^{M-1} \frac{1}{p_i} \] soit \[ Var(T(M)) = M^2\sum_{i=0}^{M-1} \frac{1}{(M-i)^2} - M\sum_{i=0}^{M-1} \frac{1}{M-i} \sim (M^2 \frac{\pi^2}{6} - M\log(M)) \sim \frac{\pi^2}{6}M^2. \]
La variance augmente effectivement bien plus vite que l’espérance. Il va falloir faire de plus en plus de simulations pour avoir une estimation correcte.
Synthèse
Q5) Vérification de la cohérence des calculs
E_T <- function(x) {x*(log(x)-digamma(1))}
ggplot(data=df_stat, aes(x=M)) + geom_point(data=df, aes(x=M,y=T),alpha=.1) +
geom_point(aes(y=mean), size=3, shape=21, fill="white") +
stat_function(fun = E_T, color="red") +
theme_bw()
Wow, ça marche quand même super bien! :)
Et si on rajoutais les intervalles de confiance, on aurait:
ggplot(data=df_stat, aes(x=M)) + geom_point(data=df, aes(x=M,y=T),alpha=.1) +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), colour="black", width=6) + geom_point(aes(y=mean), size=3, shape=21, fill="white") +
stat_function(fun = E_T, color="red") +
theme_bw()
Q6) Et mon petit cousin ?
Attention, si on veut être précis, ici, il faut pour chaque simulation regarder combien de paquet il faut acheter et en faire la moyenne. C’est différent de faire la moyenne du nombre de cartes nécessaires et de prendre une partie entière supérieure après coup… Comme je ne vais tester que pour un seul album, je peux me permettre d’augmenter le nombre d’échantillons pour avoir une bonne estimation.
df=data.frame(M=c(),T=c())
M = 300
for(samples in 1:1000) {
df = rbind(df,data.frame(M=c(M),T=new_collection(M)))
}
mean(2*ceiling(df$T/10))## [1] 379.5Finir l’album coûtera donc en moyenne 379.462 euros… Ça fait quand même cher l’album!
Extension au cas non uniforme
Q7) Simulation du marchand malhonnête
Fonction de simulation
new_collection <- function(album_size=100, proba="unif") {
switch(proba,
"unif" = {
p = rep.int(x = 1, album_size);
},
"harmonic" = {
p = 1:album_size;
p = 1/p;
},
"squared" = {
p = 1:album_size;
p = 1/p**2;
});
p = p/sum(p); # on normalise la proba
p = c(0,cumsum(p)); # on somme les valeurs pour retrouver les indices des cartes
album=rep_len(0, album_size)
nb_images = 0;
time = 0;
while(nb_images<album_size) {
time = time + 1;
idx = findInterval(runif(1),p); # pourrait aussi s'écrire sum(p<runif(1)) mais findInterval est quand même plus lisible
if(album[idx]==0) {
album[idx] = 1;
nb_images = nb_images + 1;
}
}
time
}Vérification de la cohérence
Commençons par vérifier que c’est cohérent avec le cas précédent:
df=data.frame(M=c(),T=c())
for(M in seq(10,200,20)) {
for(samples in 1:100) {
df = rbind(df,data.frame(M=c(M),T=new_collection(M)))
}
}
df_stat = ddply(df,c("M"),summarise,mean=mean(T),var=var(T),se=2*sqrt(var(T))/sqrt(length(T))) E_T <- function(x) {x*(log(x)-digamma(1))}
ggplot(data=df_stat, aes(x=M)) + geom_point(data=df, aes(x=M,y=T),alpha=.1) +
geom_point(aes(y=mean), size=3, shape=21, fill="white") +
stat_function(fun = E_T, color="red") +
theme_bw()
C’est cohérent avec ce qu’on avait avant. Parfait!
Étude du cas harmonique
df=data.frame(M=c(),T=c())
for(M in seq(10,200,20)) {
for(samples in 1:100) {
df = rbind(df,data.frame(M=c(M),T=new_collection(M,proba = "harmonic")))
}
}
df_stat = ddply(df,c("M"),summarise,mean=mean(T),var=var(T),se=2*sqrt(var(T))/sqrt(length(T))) E_T2 <- function(x) {3.5*x*(log(x)-digamma(1))}
ggplot(data=df_stat, aes(x=M)) + geom_point(data=df, aes(x=M,y=T),alpha=.1) +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), colour="black", width=6) +
geom_point(aes(y=mean), size=3, shape=21, fill="white") +
stat_function(fun = E_T, color="red") +
stat_function(fun = E_T2, color="purple") +
theme_bw()
Ah, visiblement, c’est bien plus long (la courbe rouge correspond à la valeur théorique pour le cas uniforme). La courbe violette qui correspond à la courbe rouge multipliée par 3.5 semble indiquer que ce n’est pas juste un facteur multiplicatif. La variance est également bien plus élevée qu’avant et pour avoir un bon estimateur, il faudrait donc faire beaucoup plus d’expériences… :(
Étude du cas quadratique
df=data.frame(M=c(),T=c())
for(M in seq(10,100,20)) {
for(samples in 1:100) {
df = rbind(df,data.frame(M=c(M),T=new_collection(M,proba = "squared")))
}
}
df_stat = ddply(df,c("M"),summarise,mean=mean(T),var=var(T),se=2*sqrt(var(T))/sqrt(length(T)))Aouch, déjà, on sent que c’est beaucoup plus long à simuler. C’est pour ça que je diminue le nombre de valeurs pour \(M\) et même comme ça il faut être patient.
ggplot(data=df_stat, aes(x=M)) + geom_point(data=df, aes(x=M,y=T),alpha=.1) +
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), colour="black", width=6) +
geom_point(aes(y=mean), size=3, shape=21, fill="white") +
stat_function(fun = E_T, color="red") +
stat_function(fun = E_T2, color="purple") +
theme_bw()
Wow et là, ça explose complètement! C’est assez logique en même temps, la probabilité d’avoir les dernières cartes sont extrêmement faibles…
Q8) Le cas d’une distribution exponentielle
La probabilité d’avoir la carte \(M\) est donc de l’ordre de \(1/2^M\), ce qui signifie qu’on l’obtient en moyenne au boût de \(2^M\) tirages! Pour \(M=200\), ça fait 1.6069 × 1060 donc notre simulation ne serait pas prête de terminer.