Nombre de fois prévu pour lancer un dé jusqu'à ce que chaque camp apparaisse 3 fois

15

Quel est le nombre de fois que vous devez lancer un dé jusqu'à ce que chaque camp apparaisse 3 fois?

Cette question a été posée à l'école primaire en Nouvelle-Zélande et elle a été résolue à l'aide de simulations. Quelle est la solution analytique à ce problème?

Edgar Santos
la source
6
Étant donné que les résultats des rouleaux sont aléatoires, il n'est pas possible de savoir à l'avance combien de rouleaux sont nécessaires. Si la question est à la recherche, par exemple, du nombre prévu de rouleaux avant que chaque camp ne soit apparu 3 fois, cela doit être indiqué explicitement. Dans ce cas, stats.stackexchange.com/tags/self-study/info s'applique.
Juho Kokkala
3
Dites à ces enfants néo-zélandais de lire Norman L. Johnson, Samuel Kotz, N. Balakrishnan «Discrete Multivariate Distributions» wiley.com/WileyCDA/WileyTitle/productCd-0471128449.html .
Mark L. Stone

Réponses:

28

Supposons que tous les côtés ont des chances égales. Généralisons et trouvons le nombre attendu de rouleaux nécessaires jusqu'à ce que le côté 1 apparaisse n 1 fois, le côté 2 apparaisse n 2 fois, ..., et le côté d apparaisse n d fois. Parce que les identités des côtés n'ont pas d'importance (elles ont toutes des chances égales), la description de cet objectif peut être condensée: supposons que i 0 côtés ne doivent pas du tout apparaître, i 1 des côtés doivent apparaître juste une fois, ... et i nd=61n12n2dndi0i1jendes côtés doivent apparaître fois. Soit i = ( i 0 , i 1 , , e ( 0 , 0 , 0 , 6 ) : i 3 =n=max(n1,n2,,n) désigner cette situation et écrire e ( i ) pour le nombre attendu de rouleaux. La question demande

je=(je0,je1,,jen)
e(je)
e(0,0,0,6) indique que les six côtés doivent être vus trois fois chacun.je3=6

Une récurrence facile est disponible. Au rouleau suivant, le côté qui apparaît correspond à l'un des : c'est-à-dire que nous n'avions pas besoin de le voir, ou nous devions le voir une fois, ..., ou nous devions le voir n plus fois. j est le nombre de fois que nous avons eu besoin de le voir.jejnj

  • Quand , nous n'avions pas besoin de le voir et rien ne change. Cela se produit avec la probabilité i 0 / d .j=0je0/

  • Lorsque nous devions voir ce côté. Maintenant, il y a un côté de moins qui doit être vu j fois et un autre côté qui doit être vu j - 1 fois. Ainsi, i j devient i j - 1 et i j - 1 devient i j +j>0jj-1jejjej-1jej-1 . Soit cette opération sur les composantes de i désignée ij , de sorte quejej+1jejej

    jej=(je0,,jej-2,jej-1+1,jej-1,jej+1,,jen).

    Cela se produit avec la probabilité .jej/

Nous devons simplement compter ce jet de dé et utiliser la récursivité pour nous dire combien de rouleaux de plus sont attendus. Par les lois de l'espérance et de la probabilité totale,

e(i)=1+i0de(i)+j=1nijde(ij)

(Comprenons que chaque fois que , le terme correspondant dans la somme est zéro.)ij=0

Si , nous avons terminé et e ( i ) = 0i0=de(i)=0 . Sinon, nous pouvons résoudre pour , en donnant la formule récursive souhaitéee(i)

(1)e(i)=d+i1e(i1)++ine(in)di0.

Notez que est le nombre total d'événements que nous souhaitons voir. L'opération j réduit cette quantité par une pour tout j > 0 fourni i j > 0 , ce qui est toujours le cas. Par conséquent, cette récursivité se termine à une profondeur de précisément | je | (égal à 3 ( 6 ) =

|i|=0(i0)+1(i1)++n(in)
jj>0ij>0|i| dans la question). De plus (comme cela n'est pas difficile à vérifier) ​​le nombre de possibilités à chaque profondeur de récursivité dans cette question est faible (ne dépassant jamais 8 ). Par conséquent, c'est une méthode efficace, du moins lorsque les possibilités combinatoires ne sont pas trop nombreuses et que nous mémorisons les résultats intermédiaires (de sorte qu'aucune valeur de e n'est calculée plus d'une fois).3(6)=188e

Je calcule que

e(0,0,0,6)=22868786045088836998400000000032.677.

Cela me semblait terriblement petit, alors j'ai exécuté une simulation (en utilisant R). Après plus de trois millions de lancers de dés, ce jeu avait été joué jusqu'à son terme plus de 100 000 fois, avec une longueur moyenne de . L'erreur type de cette estimation est32.669 : la différence entre cette moyenne et la valeur théorique est insignifiante, confirmant l'exactitude de la valeur théorique.0.027

La distribution des longueurs peut être intéressante. (Évidemment, il doit commencer à , le nombre minimum de rouleaux nécessaires pour collecter les six côtés trois fois chacun.)18

Figure

# Specify the problem
d <- 6   # Number of faces
k <- 3   # Number of times to see each
N <- 3.26772e6 # Number of rolls

# Simulate many rolls
set.seed(17)
x <- sample(1:d, N, replace=TRUE)

# Use these rolls to play the game repeatedly.
totals <- sapply(1:d, function(i) cumsum(x==i))
n <- 0
base <- rep(0, d)
i.last <- 0
n.list <- list()
for (i in 1:N) {
  if (min(totals[i, ] - base) >= k) {
    base <- totals[i, ]
    n <- n+1
    n.list[[n]] <- i - i.last
    i.last <- i
  }
}

# Summarize the results
sim <- unlist(n.list)
mean(sim)
sd(sim) / sqrt(length(sim))
length(sim)
hist(sim, main="Simulation results", xlab="Number of rolls", freq=FALSE, breaks=0:max(sim))

la mise en oeuvre

Bien que le calcul récursif de soit simple, il présente certains défis dans certains environnements informatiques. Le plus important d'entre eux est le stockage des valeurs de e ( i ) lors de leur calcul. Ceci est essentiel, sinon chaque valeur sera (redondante) calculée un très grand nombre de fois. Cependant, le stockage potentiellement nécessaire pour un tableau indexé par i pourrait être énorme. Idéalement, seules les valeurs de i réellement rencontrées lors du calcul devraient être stockées. Cela nécessite une sorte de tableau associatif.ee(i)ije

Pour illustrer, voici le Rcode de travail . Les commentaires décrivent la création d'une simple classe "AA" (tableau associatif) pour stocker les résultats intermédiaires. Les vecteurs sont convertis en chaînes et ceux-ci sont utilisés pour indexer dans une liste qui contiendra toutes les valeurs. L' ijiEij opération est implémentée comme %.%.

Ces préliminaires permettent la fonction récursive e de définir assez simplement d'une manière parallèle à la notation mathématique. En particulier, la ligne

x <- (d + sum(sapply(1:n, function(i) j[i+1]*e.(j %.% i))))/(d - j[1])

est directement comparable à la formule ci-dessus. Notez que tous les index ont été augmentés de 1 car commence à indexer ses tableaux à 1 plutôt qu'à 0(1)1R10 .

Le timing montre qu'il faut seconde pour calculer0.01e(c(0,0,0,6)) ; sa valeur est

32,6771634160506

L'erreur d'arrondi à virgule flottante accumulée a détruit les deux derniers chiffres (qui devraient être 68plutôt que 06).

e <- function(i) {
  #
  # Create a data structure to "memoize" the values.
  #
  `[[<-.AA` <- function(x, i, value) {
    class(x) <- NULL
    x[[paste(i, collapse=",")]] <- value
    class(x) <- "AA"
    x
  }
  `[[.AA` <- function(x, i) {
    class(x) <- NULL
    x[[paste(i, collapse=",")]]
  }
  E <- list()
  class(E) <- "AA"
  #
  # Define the "." operation.
  #
  `%.%` <- function(i, j) {
    i[j+1] <- i[j+1]-1
    i[j] <- i[j] + 1
    return(i)
  }
  #
  # Define a recursive version of this function.
  #
  e. <- function(j) {
    #
    # Detect initial conditions and return initial values.
    #
    if (min(j) < 0 || sum(j[-1])==0) return(0)
    #
    # Look up the value (if it has already been computed).
    #
    x <- E[[j]]
    if (!is.null(x)) return(x)
    #
    # Compute the value (for the first and only time).
    #
    d <- sum(j)
    n <- length(j) - 1
    x <- (d + sum(sapply(1:n, function(i) j[i+1]*e.(j %.% i))))/(d - j[1])
    #
    # Store the value for later re-use.
    #
    E[[j]] <<- x
    return(x)
  }
  #
  # Do the calculation.
  #
  e.(i)
}
e(c(0,0,0,6))

Enfin, voici l' implémentation originale de Mathematica qui a produit la réponse exacte. La mémorisation est réalisée via l' e[i_] := e[i] = ...expression idiomatique , éliminant presque tous les Rpréliminaires. En interne, cependant, les deux programmes font les mêmes choses de la même manière.

shift[j_, x_List] /; Length[x] >= j >= 2 := Module[{i = x},
   i[[j - 1]] = i[[j - 1]] + 1;
   i[[j]] = i[[j]] - 1;
   i];
e[i_] := e[i] = With[{i0 = First@i, d = Plus @@ i},
    (d + Sum[If[i[[k]] > 0, i[[k]]  e[shift[k, i]], 0], {k, 2, Length[i]}])/(d - i0)];
e[{x_, y__}] /; Plus[y] == 0  := e[{x, y}] = 0

e[{0, 0, 0, 6}]

228687860450888369984000000000

Whuber
la source
5
+1 J'imagine qu'une partie de la notation serait difficile à suivre pour les étudiants à qui on a posé cette question (pas que j'ai une alternative concrète à proposer pour le moment). D'un autre côté, je me demande ce qu'ils avaient l'intention de faire avec une telle question.
Glen_b -Reinstate Monica
1
@Glen_b Ils pourraient apprendre beaucoup en lançant les dés (et en comptant les résultats). Cela semble être un bon moyen de garder une classe occupée pendant une demi-heure pendant que l'enseignant se repose :-).
whuber
12

La version originale de cette question a commencé sa vie en demandant:

combien de rouleaux sont nécessaires jusqu'à ce que chaque côté apparaisse 3 fois

Pourquoi pas? Le problème se réduit à une doublure.

Répartition du nombre de rouleaux requis ... de sorte que chaque face apparaisse 3 fois

nXjejeje{1,,6}(X1,X2,,X6)Multinomial(n,16)

P(X1=X1,,X6=X6)=n!X1!X6!16n sujet à: je=16Xje=n

Laisser: N=min{n:Xje3je}.N est: P(Nn)=P(Xje3|n)

ie Pour trouver le cdf P(Nn), calculez simplement pour chaque valeur de n={18,19,20,}:

P(X13,,X63) où (X1,,X6)Multinomial(n,16)

Voici, par exemple, le code Mathematica qui fait cela, commen passe de 18 à 60. Il s'agit essentiellement d'un vol simple:

 cdf = ParallelTable[ 
   Probability[x1 >= 3 && x2 >= 3 && x3 >= 3 && x4 >= 3 && x5 >= 3 &&  x6 >= 3, 
       {x1, x2, x3, x4, x5, x6} \[Distributed] MultinomialDistribution[n, Table[1/6, 6]]],
    {n, 18, 60}]

... qui donne le cdf exact comme n augmente:

1814889875110199605761928290762544079842304203111983875176319369216211168408491253173748645888223283142988125507799783342082361483465418375609359740010496

Voici un tracé du cdf P(Nn), en tant que fonction de n:

enter image description here

Pour dériver le pmf P(N=n), commencez simplement par différencier le cdf:

enter image description here

Bien sûr, la distribution n'a pas de limite supérieure, mais nous pouvons facilement résoudre ici autant de valeurs que cela est pratiquement nécessaire. L'approche est générale et devrait fonctionner aussi bien pour toute combinaison désirée de côtés requise.

Wolfies
la source