Tout d’abord, rappelons le problème posé : Nous avons un ensemble de vignettes à retrouver, en tout \(M=300\). Ces vignettes sont vendues par paquet de \(N=10\) et le paquet coûte \(coutN=2\) euros. L’album pour stocker les \(M\) vignettes coûte \(album=4\) euros.

Dans un premier temps, nous supposerons qu’il achète les vignettes une par une. Il achète donc chaque vignette \(cout=0.2\) euros.

Approche expérimentale

Question 1

Dans cette partie, nous réalisons un simulateur pour pouvoir avoir une meilleure intuition du problème.

Nous commençons par définir quelques paramètres du problème ainsi que la seed pour avoir des expériences reproductibles.

coutN <- 2
album <- 4
N <- 10
cout <- coutN/N
set.seed(42)

Ensuite, nous définissons le générateur qui nous donne pour un M donné, le nombre de cartes a acheter pour avoir l’ensemble de la collection (soit M cartes différentes).

Pour le générateur, nous gardons une table qui contient les cartes que nous avons (carte). Nous gardons également une variable reste qui se décrémente à chaque fois que nous obtenons une carte que nous n’avions pas (cela permet d’éviter de parcourir table à chaque iteration). Nous tirons une nouvelle carte uniformément et indépendament jusqu’à avoir reste égal à 0.

simulateur_carte = function (M=300) {
  carte <- logical(M)
  nbtir <- 0
  reste <- M
  while (reste > 0) {
    nbtir <- nbtir+1
    tirage <- ceiling(runif(1,0,M))
    if (!carte[tirage]) {
      carte[tirage] <- TRUE
      reste <- reste-1
    }
  }
  return(nbtir)
}

Dans ce bout de code, nous pouvons vérifier que l’instruction ceiling(runif(1,0,M)) nous génère bien uniformément un entier compris entre 1 et M par la commande suivante en choisissant M=10 et un nombre de générations grand (100000):

barplot(table(ceiling(runif(100000,0,10))), xlab="nombres générés", main="Répartition des nombres générés")

Ici, runif nous génère 100000 nombres réels compris entre 0 et 5 exclus et nous calculons un seuil de ces valeurs pour avoir les entiers allant de 1 à 5. L’histogramme ci-dessus nous montre que notre générateur à l’air d’être uniforme et génère bien des entiers de 1 à 5.

Question 2

Nous pouvons maintenant tester le nombre généré sur quelques valeurs : Pour cela, nous commençons par implémenter une fonction simu_carte_moy qui prend deux arguments : nbtest et M et qui fait la moyenne sur nbtest échantillons de simulateur_carte(M).

simu_carte_moy <- function(nbtest=100, M=300) {
  sum<-0
  for ( i in 1:nbtest)
    sum <- sum + simulateur_carte(M)
  return(sum/nbtest)
}

Nous pouvons donc écrire une fonction qui à partir d’un ensemble de valeurs de M applique simu_carte_moy et renvoie les résultats assossiés.

evolution_achats <- function(nbtest=100,valeursM=seq(10,200,20)){
  res<-valeursM
  for ( i in 1:length(valeursM) )
    res[i] <- simu_carte_moy(nbtest,valeursM[i])
  return(res)
}

Nous pouvons enfin tracer la courbe correspondante à l’évolution de \(T(M)\) :

x <- seq(10,200,20)
y <- evolution_achats(1000,x)
plot(x,y, xlab="M",ylab="T(M)", main="Évolution de T(M) en fonction de M")
lines(x,y)

D’après l’évolution de la coube ci-dessus, on a l’impression que \(T(M)\) est linéaire. Réalisons un test simple pour s’en convaincre : Nous essayons de diviser \(T(M)\) par \(M\). En effet, si \(T(M)\) est linéaire, elle est de la forme \(aM+b\). Or, nous savons que \(T(0)=0\) car pour avoir 0 cartes différentes, nous n’avons pas besoin d’en acheter. donc \(b=0\).

\(T(M)/M\) serait de la forme \((aM)/M = a\) (un terme constant).

plot(x,y/x,xlab="M",ylab="T(M)/M", main="Évolution de T(M)/M en fonction de M")
lines(x,y/x)

D’après cette courbe \(T(M)/M\) n’est pas une fonction constante. Lorsque \(M\) augmente, sa valeur croît de manière logarithmique. Cela reste à vérifier, mais nous pourrons nous attendre par la suite à ce que \(T(M)\) ne soit pas tout à fait linéaire.

Modélisation

Question 3

D’après l’énoncé, les vignettes sont imprimées dans les mêmes proportions. Cela se traduit mathématiquement par : \(\forall a \in \{1,..,M\}, \mathbb{P}(X=a)=1/M\)\(X\) représentent la vignette nouvellement achetée. On dit aussi que la loi de \(X\) est uniforme.

De plus, l’énoncé suggère que les \(X_i\) sont indépendants les uns avec les autres car le fait de tirer une carte n’influe pas le fait d’en tirer une autre.

Question 4

\(Y_i\) correspond au nombre de tirages qu’il faut effectuer pour tomber sur une carte que l’on n’a pas. On a donc la formule : \(P(Y_i=k) = (1-p)^{k-1}p \)

On remarque alors que \(Y_i\) suit une loi géométrique de paramètre \(p\) qui correspond à la chance de succès (d’obtenir une nouvelle carte). \(p\) correspond à la probabilité de tirer une nouvelle carte qui est de \(\frac{M-(i-1)}{M} = \frac{M-i+1}{M}\).

Vérifions sur les cas extrêmes :

  • \(i=1\) : Nous n’avons pas encore tiré de carte. Si la formule est correcte, nous devrions trouver \(1\) comme probabilité d’en tirer une nouvelle. \(\frac{M-i+1}{M}=\frac{M-1+1}{M} = \frac{M}{M} = 1\) Ok.
  • \(i=M\) : Nous avons tiré toutes les cartes sauf une. Si la formule est correcte, nous devrions trouver \(\frac{1}{M}\) comme probabilité de tirer la carte manquante. \(\frac{M-i+1}{M}=\frac{M-M+1}{M}=\frac{1}{M}\) Ok.

Nous avons donc :

  • \(\mathbb{E}(Y_i) = \frac{1}{\frac{M-i+1}{M}} = \frac{M}{M-i+1}\) (l’espérance d’une loi géométrique de paramètre \(p\) est \(\frac{1}{p}\))

  • \(Var(Y_i) = \frac{1-\frac{M-i+1}{M}}{(\frac{M-i+1}{M})^2} = \) \(\frac{\frac{i-1}{M}}{\frac{(M-i+1)^2}{M^2}}=\) \(\frac{M(i-1)}{(M-i+1)^2}\) (la variance d’une loi géométrique de paramètre \(p\) est \(\frac{1-p}{p^2}\))

Les variables \(Y_i\) sont indépendantes car elles ne varient qu’en fonction de i, c’est à dire le nombre de cartes qu’il manque à notre collection à un moment donné, et qu’elles ne sont pas influencées par les autres résultats.

  • \(\mathbb{E}(T) = \sum_{i=1}^{M} \mathbb{E}(Y_i) = \sum_{i=1}^{M} \frac{M}{M-i+1}=\) \(M \sum_{i=1}^{M} \frac{1}{M-i+1} = \) \(M \sum_{i=1}^{M} \frac{1}{i} \approx M \log M \) Rappelons que dans la question 2, la division du nombre de tirages par la longueur de la collection nous donnait une fonction à l’allure logarithmique, ce qui vient donc corroborer ce résultat calculatoire.

  • \(Var(T) = \sum_{i=1}^{M} Var(Y_i) = \sum_{i=1}^{M} \frac{M(i-1)}{(M-i+1)^2} \) Car les \(Y_i\) sont indépendants

\(Var(T)\) étant difficile a simplifier, implémentons un code R calculant sa valeur pour nous :

varT <- function(M){
  res <- 0
  for ( i in 1:M)
    res <- res + M*(i-1)/((M-i+1)*(M-i+1)) 
  return(res)
}

Traçons son évolution en fonction de \(M\) :

x2 <- 1:300
y2 <- x2
for(i in 1:length(x2))
  y2[i] <- varT(x2[i])

plot(x2,y2,xlab="M",ylab="Var(T)", main="Évolution de Var(T) en fonction de M", type="l")

Synthèse

Question 5

Dans la partie Approche expérimentale, nous avons réalisé une simulation du modèle qui nous intéresse. L’avantage est qu’elle est facile à réaliser, l’inconvéniant est qu’elle est imprécise.

Dans la partie Modélisation, nous avons réalisé un ensemble de calculs liés aux évolutions de notre modèle suivant le nombre de cartes à collectionner.

Nous allons dans cette partie nous intéresser à comparer nos deux résultats et vérifier qu’ils sont proches l’un de l’autre.

Note : Pour être plus précis, la somme des \(1/i\) correspond à la série harmonique \(H_n\). Pour s’en approcher d’avantage, nous ajoutons la constante d’euler (\(\approx 0.577\) calculée par -digamma(1)) au log car \(H_n = ln(n)+\gamma+o(1)\).

esperanceT <- function(M){
  return(M*(log(M)-digamma(1)))
}

Dans le code R suivant, nous traçons un graphe montrant les résultats de la simulation et de la théorie effectués dans les questions précédentes :

plot(x,y,xlab="M",ylab="T et E(T)", main="Évolution de T simulé et E(T) en fonction de M")
lines(x,y)
lines(x,esperanceT(x),col="red") # traçé de E(T)
legend("topleft", legend=c("simulation","théorie"), lwd=1, col=c("black","red"))

On remarque que notre modèle théorique semble s’approcher de notre expérimentation. De plus, On avait observé dans la question 2 une courbe presque linéaire et dans la question 4, nous avons trouvé une espérance de \(T\) de \(M \log M\), ce qui colle également de ce point de vue.

Question 6

Maintenant que nous avons notre modèle, nous pouvons donner une estimation du prix qu’il devra payer pour avoir sa collection complète :

Le nombre de cartes qu’il devra acheter est de :

esperanceT(300)
## [1] 1884.299

Il devra donc payer au total:

album+coutN*ceiling(esperanceT(300)/N)
## [1] 382

Question 6 Bis

Comme nous avons calculé la moyenne et la variance de notre \(T\), il nous semble intéressant de poser un problème différent (ce qui nous permettra de vérifier au passage que nos calculs de variance sont corrects) : supposons que notre cousin achète ses paquets d’un seul coup. Dans ce cas, il a un risque de ne pas avoir sa collection complète. Combien de paquets doit-il acheter pour avoir un risque faible de ne pas avoir sa collection complète ?

Nous allons donc déterminer une fonction qui aura en paramètre le risque maximal qu’il est prêt à prendre et donne en sortie le prix minimal qu’il doit payer pour ne pas avoir un risque au delà de ce qu’il souhaite.

Intuitivement, cette fonction donnera un prix bas si il est prêt a prendre du risque et des prix élevés si il souhaite vraiment avoir sa collection complète.

Pour l’instant, intéressons nous au nombre de cartes qu’il doit acheter en fonction du risque à prendre.

Notre variable aléatoire \(T\) suit une loi normale : \(T \rightsquigarrow \mathcal{N}(m,\sigma^2)\) où :

  • \(m = \mathbb{E}(T) = M \log M\)
  • \(\sigma^2 = Var(T) = \sum_{i=1}^{M} Var(Y_i) = \sum_{i=1}^{M} \frac{M(i-1)}{(M-i+1)^2}\)

Pour \(M=300\), nous avons :

m <- esperanceT(300)
sigma <- sqrt(varT(300))

Nous cherchons donc à déterminer \(t\) la valeur limite qui respecte le seuil \(s \in [0,1]\) :

\(\mathbb{P}(T \leq t) = 1-s\)

Réalisons maintenant la fonction R qui nous donne le nombre de carte a acheter pour garantir le seuil \(s\) :

nbCartes <- function(s){
  return(ceiling(max(300,qnorm(1-s,mean=esperanceT(300),sqrt(varT(300))))))
}

On note qu’on utilise un max entre 300 et la fonction quantile de la loi normale. Cela vient du fait que si notre cousin est prêt à prendre de gros risques avec des seuils très élevés (0.9999…), la modélisation via la loi normale pourrait nous donner des résultats en dessous de 300 qui est la valeur limite pour laquelle on a encore une chance d’avoir la collection complète.

Quelques chiffres :

nbCartes(0.05) # risque faible
## [1] 2513
nbCartes(0.5)
## [1] 1885
nbCartes(0.95) # risque élevé
## [1] 1257
nbCartes(0.9999999) # nombre de cartes minimal (M=300)
## [1] 300

Vérification de ces résultats avec la simulation :

On génère un ensemble de tirages de T avec notre générateur et on compare au résulat théorique. On vérifie que le nombre de rejets correspond à la valeur qu’on a entrée dans notre fonction nbCartes :

seuilCartes <- nbCartes(0.25)
rejets <- 0
nbTests <- 500
for ( i in 1:nbTests){
  if(simulateur_carte(300)>seuilCartes)
    rejets <- rejets+1
}
rejets/nbTests
## [1] 0.234

On remarque que notre résultat est proche de ce qu’on doit obtenir, ce qui est bon signe. Cela signifie que nos calculs de variance ont de fortes chances d’être corrects.

Nous pouvons donner un aperçu de la fonction nbCartes :

x3 <- seq(0,1,0.01)
y3 <- x3
for( i in 1:length(x3) )
  y3[i] <- nbCartes(x3[i])
plot(x3,y3,type="l",xlab="seuil",ylab="nombre de cartes a acheter",main="nombre de cartes a acheter en fonction du seuil")

Nous pouvons maintenant calculer la fonction de coût à partir de nos travaux précédents :

coutSeuil <- function(s=0.95){
  return(album+ceiling(nbCartes(s)/N)*coutN)
}

Quelques chiffres :

coutSeuil(0.05) # risque faible
## [1] 508
coutSeuil(0.5)
## [1] 382
coutSeuil(0.95) # risque élevé
## [1] 256
coutSeuil(0.9999999) # (prix minimal : $4+30 \times 2$)
## [1] 64

On peut aussi tracer cette fonction :

x4 <- seq(0,1,0.01)
y4 <- x4
for( i in 1:length(x4) )
  y4[i] <- coutSeuil(x4[i])
plot(x4,y4,type="l",xlab="seuil",ylab="prix à payer",main="prix à payer en fonction du seuil")

Et ainsi retrouver le résultat de la question 6 en prenant un seuil de 0.5 :

coutSeuil(0.5)
## [1] 382

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

Question 7

Intéressons nous maintenant au cas où le marchand ne serait pas honnête et ne distriburait pas ses cartes de façon uniforme. Intuitivement, on imagine que si on a des cartes communes et d’autres rares, à chaque fois qu’on tirera une nouvelle carte, les chances d’en avoir une nouvelle seront plus faibles que dans le cas uniforme.

Considérons deux lois de probabilité pour le tirage des cartes

alpha <- function(M,loi){
  a <- 0
  for (k in 1:M)
    a <- a+loi(k,1)
  return (a)
}

loi1 <- function(i,a){
  return(1/(a*i))
}

loi2 <- function(i,a){
  return(1/(a*i*i))
}

Nous pouvons les tracer pour observer les différences de probabilité :

x <- seq(0,300,1)
a1 <- alpha(300,loi1)
y <- loi1(x,a1)
plot(x,y,main="distribution des probabilités de loi1 et loi2", xlab="ième carte",ylab="probabilité",type="l",col="blue")
a2 <- alpha(300,loi2)
y2 <- loi2(x,a2)
lines(x,y2,col="red")
legend("topright", legend=c("loi1","loi2"), lwd=1, col=c("blue","red"))

On remarque que pour la loi 2, les probabilités sont d’avantage disproportionnées que dans la loi 1.

Maintenant que nous avons nos lois de probabilité, nous pouvons réaliser une simulation en reprennant notre générateur fait à la question 1 pour tester le nombre de cartes en moyenne qu’il faudra acheter par la simulation :

Pour tirer un nombre suivant une loi donnée, nous utilisons le code R suivant :

tirage2 <- function(loi,a){
  alea <- runif(1) # tirage uniforme d'un nombre entre 0 et 1
  somme <- 0
  i<-1
  while(alea>somme+loi(i,a)){
    somme <- somme+loi(i,a)
    i <- i+1
  }
  return(i)
}

Testons notre générateur avec la loi 1 :

y <- 1:10000
a1 <- alpha(10,loi1) 
for ( i in 1:length(y) )
  y[i] <-tirage2(loi1,a1)
barplot(table(ceiling(y)), xlab="nombres générés", main="Répartition des nombres générés")

Avec loi 2 :

y <- 1:10000
a2 <- alpha(10,loi2) 
for ( i in 1:length(y) )
  y[i] <-tirage2(loi2,a2)
barplot(table(ceiling(y)), xlab="nombres générés", main="Répartition des nombres générés")

Le code de la question 1 modifié. On ajoute la loi en paramètre du générateur et on modifie la façon d’obtenir tirage

simulateur_carte2 <- function (loi,M=300, a=alpha(M,loi)) { # On rajoute le paramètre loi
  carte <- logical(M)
  nbtir <- 0
  reste <- M
  while (reste > 0) {
    nbtir <- nbtir+1
    tirage <- tirage2(loi,a) # On modifie la façon d'obtenir tirage
    if (!carte[tirage]) {
      carte[tirage] <- TRUE
      reste <- reste-1
    }
  }
  return(nbtir)
}

De la même façon que la question 1, on réalise une moyenne sur un ensemble d’appels à simulateur_carte2 et pour un certain nombre de \(M\) :

simu_carte_moy2 <- function(nbtest=100, M=300,loi) {
  sum<-0
  for ( i in 1:nbtest) {
    a = alpha(M,loi)
    sum <- sum + simulateur_carte2(loi,M,a)
  }

  return(sum/nbtest)
}

evolution_achats2 <- function(nbtest=100,valeursM=seq(10,200,20),loi){
  res<-valeursM
  for ( i in 1:length(valeursM) )
    res[i] <- simu_carte_moy2(nbtest,valeursM[i],loi)
  return(res)
}

Forts de ces fonctions, on peut maintenant tracer \(T\) en fonction de \(M\).

x <- seq(5,55,10)
y <- evolution_achats2(10,x,loi1)
plot(x,y, xlab="M",ylab="T(M)", main="Évolution de T(M) en fonction de M avec la loi 1")
lines(x,y)

x <- seq(5,25,3)
y <- evolution_achats2(20,x,loi2)
plot(x,y, xlab="M",ylab="T(M)", main="Évolution de T(M) en fonction de M avec la loi 2")
lines(x,y)

On se rend compte que le fait d’avoir des probabilités qui ne sont pas uniformes augmente considérablement le nombre d’achats de cartes requis pour avoir la collection.

De plus, on remarque que les lois 1 et 2 font grandir plus vite \(T\) en fonction de \(M\).

Question 8

Si on effectue le même essai avec \(\alpha = \sum_{k=1}^{M} \frac{1}{2^k}\), nous aurons encore plus de mal qu’avec les lois 1 et 2 :

Admettons que l’on prenne \(M=300\) :

  • pour la loi 1, nous avons un rapport de 300 entre les probabilités d’avoir la carte la plus fréquente et la moins fréquente.
  • pour la loi 2, nous avons un rapport de \(300^2=90.000\) entre les probabilités d’avoir la carte la plus fréquente et la moins fréquente. (c’est déjà beaucoup !)
  • pour la loi 3, nous avons un rapport de \(2^{300}\) (un nombre avec a peu près 100 zéros) entre les probabilités d’avoir la carte la plus fréquente et la moins fréquente, ce qui n’est pas calculable avec les machines actuelles.