Número esperado de veces para lanzar un dado hasta que cada lado haya aparecido 3 veces

Question

¿Cuál es el número esperado de veces que hay que lanzar un dado hasta que cada lado haya aparecido 3 veces?

Esta pregunta se planteó en la escuela primaria de Nueva Zelanda y se resolvió mediante simulaciones. ¿Cuál es la solución analítica para este problema?

Answer 1

Supongamos que todos los $d=6$ los lados tienen las mismas oportunidades. Generalicemos y encontremos el número esperado de tiradas necesarias hasta que el lado $1$ ha aparecido $n_1$ veces, lado $2$ ha aparecido $n_2$ tiempos, ..., y lado $d$ ha aparecido $n_d$ veces. Como las identidades de los lados no importan (todos tienen las mismas posibilidades), la descripción de este objetivo se puede condensar: supongamos que $i_0$ los lados no tienen que aparecer en absoluto, $i_1$ de los lados tiene que aparecer sólo una vez, ..., y $i_n$ de los lados tienen que aparecer $n=\max(n_1,n_2,\ldots,n_d)$ tiempos. Dejemos que $$\mathbf{i}=(i_0,i_1,\ldots,i_n)$$ designe esta situación y escriba $$e(\mathbf{i})$$ para el número esperado de tiradas. En la pregunta se pide $e(0,0,0,6)$ : $i_3 = 6$ indica que los seis lados deben ser vistos tres veces cada uno.

Hay una recurrencia fácil. En la siguiente tirada, el lado que aparece corresponde a uno de los $i_j$ Es decir, o no necesitábamos verlo, o necesitábamos verlo una vez, ..., o necesitábamos verlo $n$ más veces. $j$ es el número de veces que necesitamos verlo.

• Cuando $j=0$ No necesitamos verlo y nada cambia. Esto sucede con probabilidad $i_0/d$ .

• Cuando $j \gt 0$ entonces sí necesitábamos ver este lado. Ahora hay un lado menos que necesita ser visto $j$ veces y un lado más que hay que ver $j-1$ tiempos. Así, $i_j$ se convierte en $i_j-1$ y $i_{j-1}$ se convierte en $i_j+1$ . Que esta operación sobre los componentes de $\mathbf{i}$ ser designado $\mathbf{i}\cdot j$ para que

  $$\mathbf{i}\cdot j = (\color{gray}{i_0, \ldots, i_{j-2}}, i_{j-1}+1, i_j-1, \color{gray}{i_{j+1},\ldots, i_n}).$$

  Esto ocurre con probabilidad $i_j/d$ .

Sólo tenemos que contar esta tirada y utilizar la recursividad para saber cuántas tiradas más se esperan. Por las leyes de la expectativa y la probabilidad total,

$$e(\mathbf{i}) = 1 + \frac{i_0}{d}e(\mathbf{i}) + \sum_{j=1}^n \frac{i_j}{d}e(\mathbf{i}\cdot j)$$

(Entendamos que siempre que $i_j=0$ el término correspondiente en la suma es cero).

Si $i_0=d$ , hemos terminado y $e(\mathbf{i}) =0$ . De lo contrario, podemos resolver para $e(\mathbf{i})$ , dando lugar a la fórmula recursiva deseada

$$e(\mathbf{i}) = \frac{d + i_1 e(\mathbf{i}\cdot 1) + \cdots + i_n e(\mathbf{i}\cdot n)}{d - i_0}.\tag{1}$$

Observe que $$|\mathbf{i}| = 0(i_0) + 1(i_1) + \cdots + n(i_n)$$ es el número total de eventos que deseamos ver. La operación $\cdot j$ reduce esa cantidad en uno para cualquier $j\gt 0$ proporcionó $i_j \gt 0$ que es siempre el caso. Por lo tanto, esta recursión termina en una profundidad de precisamente $|\mathbf{i}|$ (igual a $3(6) = 18$ en la pregunta). Además (como no es difícil de comprobar) el número de posibilidades en cada profundidad de recursión en esta pregunta es pequeño (nunca superior a $8$ ). En consecuencia, se trata de un método eficiente, al menos cuando las posibilidades combinatorias no son demasiado numerosas y memorizamos los resultados intermedios (para que ningún valor de $e$ se calcula más de una vez).

Calculo que $$e(0,0,0,6) = \frac{2\,286\,878\,604\,508\,883}{69\,984\,000\,000\,000}\approx 32.677.$$

Eso me pareció terriblemente pequeño, así que hice una simulación (usando R ). Después de más de tres millones de tiradas de dados, este juego se había jugado hasta su finalización más de 100.000 veces, con una duración media de $32.669$ . El error estándar de esa estimación es $0.027$ La diferencia entre esta media y el valor teórico es insignificante, lo que confirma la exactitud del valor teórico.

La distribución de las longitudes puede ser interesante. (Obviamente debe comenzar en $18$ el número mínimo de tiradas necesarias para recoger las seis caras tres veces cada una).

# 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))

---

Aplicación

Aunque el cálculo recursivo de $e$ es sencillo, presenta algunos retos en algunos entornos informáticos. El principal es el almacenamiento de los valores de $e(\mathbf{i})$ a medida que se calculan. Esto es esencial, ya que, de lo contrario, cada valor se calculará (redundantemente) un gran número de veces. Sin embargo, el almacenamiento potencialmente necesario para una matriz indexada por $\mathbf{i}$ podría ser enorme. Idealmente, sólo los valores de $\mathbf{i}$ que realmente se encuentran durante el cómputo deben ser almacenados. Esto requiere una especie de matriz asociativa.

Para ilustrar, aquí está el trabajo R código. Los comentarios describen la creación de una simple clase "AA" (associative array) para almacenar los resultados intermedios. Vectores $\mathbf{i}$ se convierten en cadenas y éstas se utilizan para indexar en una lista E que contendrá todos los valores. La página web $\mathbf{i}\cdot j$ se implementa como %.% .

Estos preliminares permiten la función recursiva $e$ para ser definida de forma bastante sencilla de manera paralela a la notación matemática. En particular, la línea

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

es directamente comparable a la fórmula $(1)$ arriba. Obsérvese que todos los índices se han incrementado en $1$ porque R comienza a indexar sus matrices en $1$ en lugar de $0$ .

El tiempo muestra que se necesita $0.01$ segundos para calcular e(c(0,0,0,6)) su valor es

32.6771634160506

El error de redondeo en coma flotante acumulado ha destruido los dos últimos dígitos (que deberían ser 68 en lugar de 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))

Por último, aquí está el original Mathematica implementación que produjo la respuesta exacta. La memoización se realiza a través de el idiomático e[i_] := e[i] = ... expresión, eliminando casi todo el R preliminares. Sin embargo, internamente, los dos programas hacen las mismas cosas de la misma manera.

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}]

$\frac{2286878604508883}{69984000000000}$

Answer 2

La versión original de esta pregunta comenzó preguntando:

cuántas tiradas son necesarias hasta que cada lado haya aparecido 3 veces

Por supuesto, esa es una pregunta que no tiene respuesta, como comentó @JuhoKokkala más arriba: la respuesta es una variable aleatoria con una distribución que hay que encontrar. La pregunta se modificó entonces para preguntar: "Cuál es el número esperado de tiradas". La respuesta que se ofrece a continuación trata de responder a la pregunta original planteada: cómo encontrar el distribución del número de rollos Sin recurrir a la simulación, y utilizando únicamente técnicas conceptualmente sencillas que cualquier estudiante neozelandés con un ordenador podría aplicar $\rightarrow$ ¿Por qué no? El problema se reduce a una sola línea.

Distribución del número de rollos necesarios ... de forma que cada lado aparezca 3 veces

Tiramos un dado $n$ tiempos. Dejemos que $X_i$ denotan el número de veces que el lado $i$ de la matriz aparece, donde $i \in \{1, \dots, 6\}$ . Entonces, el pmf conjunto de $(X_1, X_2,\dots, X_6)$ es $\text{Multinomial}(n,\frac16)$ es decir :

$$P\left(X_1=x_1,\ldots ,X_6=x_6\right) \; = \; \frac{n! }{ x_1! \cdots x_6!} \; \frac{1}{6^n} \quad \text{ subject to: } \quad \sum _{i=1}^6 x_i=n$$

Déjalo: $\quad N = \min\big\{n: \; {X_i \geq 3 \; \forall_i } \big\}. \;$ Entonces la fdc de $N$ es: $\quad P(N \leq n) \; = \; P\big(X_{\forall_i} \geq 3 \; \big| \; n\big) $

Es decir, para encontrar la fdc $P(N \leq n)$ simplemente calcula para cada valor de $n = \{18, 19, 20,\dots\}$ :

$$P(X_1 \geq3, \dots , X_6 \geq 3) \quad \text{ where } \quad (X_1, \dots, X_6) \sim \text{Multinomial}(n,\frac16)$$

Aquí, por ejemplo, está Mathematica código que hace esto, como $n$ aumenta de 18 a, digamos, 60. Se trata básicamente de una sola línea:

 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}]

... lo que da como resultado la fdc exacta como $n$ aumenta:

$$\begin{array}{cc} 18 & \frac{14889875}{11019960576} \\ 19 & \frac{282907625}{44079842304} \\ 20 & \frac{3111983875}{176319369216} \\ 21 & \frac{116840849125}{3173748645888} \\ 22 & \frac{3283142988125}{50779978334208} \\ 23 & \frac{61483465418375}{609359740010496} \\ \vdots & \vdots\\ \\ \end{array}$$

Este es un gráfico de la fdc $P(N\leq n)$ en función de $n$ :

$ $

Para derivar el pmf $P(N=n)$ Simplemente, primero hay que diferenciar la fdc:

Por supuesto, la distribución no tiene un límite superior, pero podemos resolverlo fácilmente para tantos valores como sea necesario en la práctica. El enfoque es general y debería funcionar igual de bien para cualquier combinación deseada de lados requeridos.