Approche expérimentale :

1. Ecrire un simulateur en R générant T en fonction de M.

T étant le nombre total de vignettes achetées pour avoir toute la collection. M étant les 300 places distinctes pour les vignettes autocollantes.

#choisi le num de la vignette entre 1 et 300, de façon uniforme, en considérant qu'il n'existe pas de vignette ayant pour numéro 0
achat_vignette <- function (M = 300){
  v = c();
  v = sample(1:M,1);
  v;
}
achat_vignette();
## [1] 236
# regarde si la vignette est dans l'album et l'ajoute si elle ne l'est pas.
complete_collection <- function (a = c(), v = 1, M = 10){
  if (a[v]== 0){
    a[v] = 1;
    a[M+1] = a[M+1] + 1; 
  }
  a;
}
#génère l'album complet et compte le nombre de vignette acheté
album_complet <- function (M = 300){
  a = c();
  T = c(0);
  v = c();
  i=1;
  for (i in 1:(M+1)){
    a[i] = 0;
  }
  while(a[M+1] < M){
    T[1] = T[1] + 1;
    v = achat_vignette(M)
    a = complete_collection(a, v[1], M);
  }  
  T;
}
album_complet();
## [1] 1577

2. Etudiez expérimentalement la valeur moyenne de T en fonction de M. On prendra des échantillons de l’ordre de 100 ou 1000 et une dizaine de valeurs M entre 10 et 200. Une représentation graphique sera appréciée. . . Que remarquez vous ?

#fonction permettant d'executer E fois la fonction album_complet
album_echantillion <- function (M = 10, E = 100){
  e = c();
  for(i in 1:E){
    e[i] = album_complet(M);
  }
  e;
}
album_echantillion();
##   [1] 34 21 43 20 25 37 26 27 29 50 19 20 19 15 27 67 27 44 21 20 27 22 53
##  [24] 15 14 47 15 36 24 23 30 36 26 34 29 46 31 41 18 21 29 44 31 31 36 38
##  [47] 27 12 21 42 25 20 25 23 29 13 45 23 16 22 38 27 28 18 45 31 57 26 39
##  [70] 21 26 32 43 51 18 43 29 28 25 28 25 27 26 37 31 32 33 58 43 35 25 36
##  [93] 57 50 37 28 39 60 26 24
# calcule la moyenne d'un echantillion e 
moyenne_T <- function (e = c()){
  mean(e);
}
n=album_echantillion(300,100);
moyenne_T(n);
## [1] 1921.66
#1839.02 pour 300 et 100 lancers
#fonction d'affichage de l'histograme et du plot de x 
affiche <- function(x)
{
  par(mfrow=c(1,2))
  hist(x)
  plot(x[1:length(x)-1],x[2:length(x)])
}
n_10=album_echantillion();
moyenne_T(n_10);
## [1] 27.89
affiche(n_10);

Echantillonnage :

Moyenne <- function (M = 10){
  m = c();
  
  m[1] = moyenne_T(album_echantillion(M,100));
  m[2] = moyenne_T(album_echantillion(M,250));
  m[3] = moyenne_T(album_echantillion(M,500));
  m[4] = moyenne_T(album_echantillion(M,750));
  m[5] = moyenne_T(album_echantillion(M,1000));
   
  m;
}
m = Moyenne(10);
affiche(m);

Les valeurs tendent à s’uniformiser. Nous avons pour une même valeur de M la même moyenne (à quelques 10ème près) quelque soit le nombre de tirage.

Modélisation

3. Quelles sont les hypothèses statistiques sur la suite des variables Xn que l’énoncé suggère ?

On peut supposer que, comme les vignettes sont imprimées dans les mêmes proportions, les variables aléatoires aurons la même probabilité d’apparaitre, qui est de 1/M. Le tirage se comporte donc comme un tirage aléatoire avec remise.

4. Lorsque la collection comporte déjà i-1 vignettes, on note Yi le nombre de vignettes à acheter pour avoir une vignette qui n’est pas dans la collection. Donner, en la justifiant, la loi de Yi (vous pourrez calculer la probabilité P(Yi = k) et vous en déduirez l’espérance et la variance de Yi). Les variables Yi sont-elles indépendantes ? En déduire E(T) et Var(T) ainsi que des équivalents quant M est grand.

P(Yi=k) = ((M-(i-1)/(M^k))*((i-1)^(k-1))

On veut montrer que P(Yi = k) = ((M-(i-1))/M^k) * (i-1)^(k-1) avec ((M-(i-1))/M^k) la probabilité qu’une combinaison, dont la k-ème vignette soit une vignette inédite, soit tirée et (i-1)k-1 le nombre total de ces combinaisons

Pour nos démonstrations, on aura M = 5 avec les vignettes comprises dans l’ensemble {A, B, C, D, E}

Exemple: P(Y2 = 2) => on a donc 1 vignette dans notre collection (on choisit A). On va donc lister toutes les combinaisons où A est en première position et une autre lettre en deuxième. On obtient alors: P(Y2 = 2) = P( (A,B) , (A,C) , (A,D) , (A,E) ) Sachant que chaque vignette a 1/5 chance d’être tirée, on peut en conclure que chaque combinaison de deux vignettes a une probabilité de 1/5^2 = 1/25 d’être tirée. Comme le fait de tirer une vignette est un évènement indépendant des autres, on peut noter: P(Y2 = 2) = 1/25 * 4 = 4/25 Or ((M-(i-1))/M^k) * (i-1)^(k-1) = ( (5-1)/5^2 ) * 1^1 = 4/25 * 1 = 4/25 => notre calcul est vérifié

Autre exemple: P(Y3 = 3) => on a 2 vignettes dans notre collection, A et B. En suivant la même logique que précédemment, on a donc la liste suivante: P(Y3 = 3) = P( (A, A, C) , (A, A, D) , (A, A, E) , (A, B, C) , (A, B, D) , (A, B, E) , (B, A, C) , (B, A, D) , (B, A, E) , (B, B, C) , (B, B, D) , (B, B, E) ) Cette fois, la probabilité d’apparition est de 1/5^3 = 1/125. De plus, on voit qu’on a 34 = 12 combinaisons différentes, on peut donc conclure que: P(Y3 = 3) = 1/125 12 = 12/125 Et ((M-(i-1))/M^k) * (i-1)^(k-1) = ((5-2)/5^3) * 2^2 = 3/125 * 4 = 12/125 => le calcul est vérifié

On remarque que P(Yi = k) = ((M-(i-1)/(M^k))*((i-1)^(k-1)) suit une loi géométrique de paramètre p = (M-(i-1)/(M^k)).

On peut donc en déduire que : E(Yi) = 1/p = 1/(M-(i-1)/(M^k)) V(Yi) = q/(p)² = ((i-1)(k-1))/(((M-(i-1)/(Mk))²)

On voit que Yi ne s’exprime pas en fonction des vignettes déjà achetées. Donc Les variables sont indépendantes.

5. Vérifiez graphiquement que les résultats obtenus par simulation sont cohérents avec vos calculs.

# Compte le nombre de vignette à acheter avant de trouver une nouvelle vignette 
Yi <- function(M = 10,a = c(1, 1, 1, 1)){
  T = c(0);
  v = c();
  i=1;
  a[M+1] = 0;
  for (i in 1:M){
    if(is.na(a[i])){
      a[i] = 0;
    }else{
      a[M+1] = a[M+1] + 1;
    }
  }
  temp = a[M+1];
  while(temp == a[M+1]){
    T[1] = T[1] + 1;
    v = achat_vignette(M)
    a = complete_collection(a, v[1], M);
  }  
  T;
}
Yi();
## [1] 1
# génère un vecteur contenant différente valeur de Y sur N tirages.
generateurYi <- function (N = 100){
  Y = c();
  for(i in 1:N){
    Y[i] = Yi();
  }
  Y;
}
generateurYi();
##   [1] 1 5 2 2 2 1 2 1 2 1 1 1 2 1 1 2 4 1 4 3 1 1 2 2 1 1 1 2 1 3 1 2 3 1 2
##  [36] 5 2 1 1 2 1 2 1 1 1 1 1 3 1 2 2 1 2 1 1 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1
##  [71] 1 1 1 1 2 1 3 1 1 3 3 1 1 1 2 2 1 2 5 1 2 1 3 1 2 1 3 1 4 3
# Calcule la probabilité de P(Yi = k)
proba_Yi <- function(y = c()){
  P = c();
  k = 1;
  autre = 1;# sert à savoir s'il y a une autre valeur de k que l'on n'a pas déjà traitée
  while(autre == 1){
    autre = 0;
    i = 1;
    P[k] = 0;
    while (i < length(y)){
      i = i + 1;
      if(y[i] == k){
        P[k] = P[k] + 1;
      }else{
         if(y[i] > k && autre != 1){
           autre = 1;
         }
      }
    }
    P[k] = P[k]/length(y);
    k = k + 1;
  }
  P;
}
proba_Yi(generateurYi(5000));
##  [1] 0.6046 0.2440 0.0886 0.0412 0.0116 0.0034 0.0038 0.0024 0.0000 0.0002
Y = generateurYi(5000);
affiche(Y);

#plot de la probabilité de Yi pour toutes les valeurs de k
P = proba_Yi(generateurYi(5000));
plot(P);

6. Au fait, combien votre petit cousin devra-t-il dépenser en moyenne sachant que chaque paquet de 2 euros contient 10 vignettes et que l’album comporte 300 places ?

On va calculer la dépense de notre petit cousin en se basant sur la quantité T moyenne qu’il devra acheter de vignettse. Nous savons que la loi est uniformes (les vignettes ont la même probabilité d’être tirées) et les tirages sont indépendants. On peut donc en déduire que le coût moyen peut être calculé en divisant le nombre de vignettes acheté par 10 (on aura ainsi le nombre de paquets), le multiplier par 2 et ajouter 4. Ce qui nous donne (T moyen/10) * 2 + 4. Ou plus précisément (E(T)/10) * 2 + 4.

La fonction simule la dépense en calculant par étape. On aurait pu simplement faire le calcul (T moyen/10 * 2 + 4).

# Simulation de la dépense de notre petit cousin par rapport au nombre T moyen
depenser_en_moyenne <- function(T = c()){
  Depense = c(0);
  i = 0;
  while( i < T[1]){
    Depense[1] = Depense[1] + 2;
    i = i + 10;
  }
  Depense[1] = Depense[1] + 4; 
  Depense;
}
depenser_en_moyenne(T = c(1839.02));# valeur que l'on avait trouvée précédemment lors du calcul de T moyen 
## [1] 372
# ce qui nous donne 372

Extension au cas non uniforme (ou au marchand malhonnête)

Q7) On s’intéresse maintenant au cas où la loi de probabilité d’obtenir chaque vignette n’est plus uniforme.

Pour adapter notre simulateur aux nouvelles probabilités on va utiliser l’aliasing

# calcul de alpha
alpha_cas1 <- function(M = 300){
  alpha = c(0);
  for(k in 1:M){
    alpha[1] = alpha[1] + 1/k
  }
  alpha;
}
alpha_cas1();
## [1] 6.282664
# création du tableau de probabilité
proba_cas1 <- function(M = 300, alpha = c(1)){
  proba = c();
  for(i in 1:M){
    proba[i] = 1/(alpha*i);
  }
  proba;
}
P=proba_cas1(300,alpha_cas1(300));
plot(P);

Aliasing

Table_Alias <- function(M = 300, P = c(1))
{
    Alias = data.frame();
    
    S = c();
    A = c();
    
    L = c(0);
    l = 0;
    U = c(0);
    u = 0;   
    
    for(k in 1:M){
      if(P[k] < 1/M){
          l = l + 1;
          L[l] = k;
      }else {
        if(P[k] > 1/M){
          u = u + 1;
          U[u] = k;
        } 
      }     
    }
        
    while(l > 0 && u > 0){
      i = L[l];
      l = l - 1;
      k = U[u];
      u = u - 1;
      S[i] = P[i];
      A[i] = k;
      P[k] = P[k] - (1/M - P[i]);
      if(P[k] < 1/M) {
        l = l + 1;
        L[l] = k; 
      } else {
          if(P[k] > 1/M){
            u = u + 1;
            U[u] = k;
          }
      }         
    }
    S[1] = P[1];
    A[1] = 0;    
    Alias = rbind(Alias,data.frame(S, A));
    Alias;
}
Table_Alias(5,proba_cas1(5,alpha_cas1(5)));
##            S A
## 1 0.20000000 0
## 2 0.10656934 1
## 3 0.14598540 1
## 4 0.10948905 1
## 5 0.08759124 2
Genere <- function(S = c(), A = c(), M = 300){
  k = sample(1:M,1);
  rand = (runif(1)*(1/M));
  if( rand < S[k]){
    k;
  } 
  else{
    A[k];
  }
}
alias = Table_Alias(300,proba_cas1(300,alpha_cas1(300)));
Genere(alias$S, alias$A, 300);
## [1] 1

Simulation de l’achat de toutes les vignettes de l’album dans le cas 1 :

album_complet_cas1 <- function (M = 300){
  a = c();
  v = c();
  T = c(0);
  i=1;
  for (i in 1:(M+1)){
    a[i] = 0;
  }
  proba = proba_cas1(M,alpha_cas1(M));
  alias = Table_Alias(M,proba);
  while(a[M+1] < M){
    T[1] = T[1] + 1;
    v = Genere(alias$S, alias$A, M);
    a = complete_collection(a, v[1], M);
  }  
  T;
}
album_complet_cas1();
## [1] 9849

Echantillion dans le cas1 :

album_echantillion_cas1 <- function (M = 10, E = 100){
  e = c();
  for(i in 1:E){
    e[i] = album_complet_cas1(M);
  }
  e;
}
e = album_echantillion_cas1(150,100);
moyenne_T(e);
## [1] 3308.34
# 8539.17 pour M = 300 et E = 100
affiche(e);

Deuxième cas :

alpha_cas2 <- function(M = 300){
  alpha = c(0);
  for(k in 1:M){
    alpha[1] = alpha[1] + 1/(k*k)
  }
  alpha;
}
alpha_cas2();
## [1] 1.641606
proba_cas2 <- function(M = 300, alpha = c(1)){
  proba = c();
  for(i in 1:M){
    proba[i] = 1/(alpha*i*i);
  }
  proba;
}
P = proba_cas2(300,alpha_cas2(300));
plot(P);

Simulation de l’achat de toutes les vignettes de l’album pour le cas2

album_complet_cas2 <- function (M = 300){
  a = c();
  v = c();
  T = c(0);
  i=1;
  for (i in 1:(M+1)){
    a[i] = 0;
  }
  proba = proba_cas2(M,alpha_cas2(M));
  alias = Table_Alias(M,proba);
  while(a[M+1] < M){
    T[1] = T[1] + 1;
    v = Genere(alias$S, alias$A, M);
    a = complete_collection(a, v[1], M);
  }  
  T;
}
album_complet_cas2();
## [1] 404943

Echantillion cas2 :

album_echantillion_cas2 <- function (M = 10, E = 100){
  e = c();
  for(i in 1:E){
    e[i] = album_complet_cas2(M);
  }
  e;
}
e2 = album_echantillion_cas2(30,100);
moyenne_T(e2);
## [1] 3234.45
# 589439.3 pour M = 300 et E = 100
affiche(e2);

Cas 1 : Pour M = 10

e1_M10_100 = album_echantillion_cas1(10);
moyenne_T(e1_M10_100);# 59.51
## [1] 59.66
affiche(e1_M10_100);

Pour M = 25

e1_M25_100 = album_echantillion_cas1(25);
moyenne_T(e1_M25_100);# 237.09
## [1] 239.04
affiche(e1_M25_100);

Pour M = 50

e1_M50_100 = album_echantillion_cas1(50);
moyenne_T(e1_M50_100);# 655.42
## [1] 689.35
affiche(e1_M50_100);

Pour M = 75

e1_M75_100 = album_echantillion_cas1(75);
moyenne_T(e1_M75_100);# 1211.13
## [1] 1295
affiche(e1_M75_100);

Pour M = 100

e1_M100_100 = album_echantillion_cas1(100);
moyenne_T(e1_M100_100);# 1785.75
## [1] 1913.84
affiche(e1_M100_100);

Pour M = 125

e1_M125_100 = album_echantillion_cas1(125);
moyenne_T(e1_M125_100);# 2491.87
## [1] 2604.71
affiche(e1_M125_100);

Pour M = 150, on trouve des valeurs approchant 3157.95.

Pour M = 175, on trouve des valeurs approchant 3924.86.

Pour M = 200, on trouve des valeurs approchant 4905.24.

Cas 2 : Pour M = 10

e2_M10_100 = album_echantillion_cas2();
moyenne_T(e2_M10_100);# 262.02
## [1] 253.07
affiche(e2_M10_100);

Pour M = 25

e2_M25_100 = album_echantillion_cas2(25);
moyenne_T(e2_M25_100);# 2243.36
## [1] 2120.64
affiche(e2_M25_100);

Pour M = 75, on trouve des valeurs approchant 26660.18.

Pour M = 150, on trouve des valeurs approchant 118358.

Dans les deux cas, la moyenne de T augmente avec M et augmente d’une façon plutôt régulière.

Pour le cas 1 , on a remarqué que lorsqu’on multiplie M par 2, T est multiplié par environ 3.

8 .Selon vous, que se passera-t-il pour une simulation avec pi = 1/(alpha)*2^i et alpha = sum(1/(2^i) ?

Les probabilités vont décroitre de façon exponentielle plus rapidement que dans le cas2 et T va augmenter beaucoup plus rapidement.