Número esperado de vezes para rolar um dado até que cada lado apareça 3 vezes

15

Qual é o número esperado de vezes que você deve rolar um dado até que cada lado apareça 3 vezes?

Esta pergunta foi feita na escola primária da Nova Zelândia e foi resolvida usando simulações. Qual é a solução analítica para esse problema?

Edgar Santos
fonte
6
Como os resultados dos testes são aleatórios, não é possível saber antecipadamente quantos testes são necessários. Se a pergunta estiver procurando, por exemplo, o número esperado de jogadas antes de cada lado aparecer três vezes, isso deve ser declarado explicitamente. Nesse caso, stats.stackexchange.com/tags/self-study/info se aplica.
Juho Kokkala
3
Diga às crianças da Nova Zelândia que leiam Norman L. Johnson, Samuel Kotz, N. Balakrishnan "Distribuições multivariadas discretas" wiley.com/WileyCDA/WileyTitle/productCd-0471128449.html .
Mark L. Stone

Respostas:

28

Suponha que todos os lados d=6 tenham chances iguais. Vamos generalizar e encontrar o número esperado de rolagens necessárias até que o lado 1 apareça n1 vezes, o lado 2 apareça n2 vezes, ... e o lado d apareça nd vezes. Como as identidades dos lados não importam (todas têm chances iguais), a descrição desse objetivo pode ser condensada: suponhamos que i0 lados não precisam aparecer, i1 dos lados precisa aparecer apenas uma vez, ... e indos lados deve aparecer n=max(n1,n2,,nd) vezes. Seja

Eu=(Eu0 0,Eu1 1,...,Eun)
designe esta situação e escreva
e(Eu)
para o número esperado de rolagens. A pergunta pede e(0 0,0 0,0 0,6) : Eu3=6 indica que todos os seis lados precisam ser vistos três vezes cada.

Uma recorrência fácil está disponível. Na próxima rodada, o lado que aparece corresponde a um dos : isto é, ou não precisa vê-lo, ou que precisávamos para vê-lo uma vez, ..., ou que precisávamos para vê-lo n mais vezes. j é o número de vezes que precisávamos vê-lo.Eujnj

  • Quando , não precisamos vê-lo e nada muda. Isso acontece com probabilidade i 0 / d .j=0 0Eu0 0/d

  • Quando , precisávamos ver esse lado. Agora, há um lado a menos que precisa ser visto j vezes e outro lado que precisa ser visto j - 1 vezes. Assim, i j se torna i j - 1 e i j - 1 se torna i j + 1 . Deixe esta operação nos componentes de i ser designada ij , para quej>0 0jj1ijij1ij1ij+1iij

    ij=(i0,,ij2,ij1+1,ij1,ij+1,,in).

    Isso acontece com a probabilidade .ij/d

Nós apenas temos que contar esse teste e usar a recursão para nos dizer quantos mais testes são esperados. Pelas leis da expectativa e probabilidade total,

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

(Let's understand that whenever ij=0, the corresponding term in the sum is zero.)

If i0=d, we are done and e(i)=0. Otherwise we may solve for e(i), giving the desired recursive formula

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

Notice that

|i|=0(i0)+1(i1)++n(in)
is the total number of events we wish to see. The operation j reduces that quantity by one for any j>0 provided ij>0, which is always the case. Therefore this recursion terminates at a depth of precisely |i| (equal to 3(6)=18 in the question). Moreover (as is not difficult to check) the number of possibilities at each recursion depth in this question is small (never exceeding 8). Consequently, this is an efficient method, at least when the combinatorial possibilities are not too numerous and we memoize the intermediate results (so that no value of e is calculated more than once).

I compute that

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

That seemed awfully small to me, so I ran a simulation (using R). After over three million rolls of the dice, this game had been played to its completion over 100,000 times, with an average length of 32.669. The standard error of that estimate is 0.027: the difference between this average and the theoretical value is insignificant, confirming the accuracy of the theoretical value.

The distribution of lengths may be of interest. (Obviously it must begin at 18, the minimum number of rolls needed to collect all six sides three times each.)

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

Implementation

Although the recursive calculation of e is simple, it presents some challenges in some computing environments. Chief among these is storing the values of e(i) as they are computed. This is essential, for otherwise each value will be (redundantly) computed a very large number of times. However, the storage potentially needed for an array indexed by i could be enormous. Ideally, only values of i that are actually encountered during the computation should be stored. This calls for a kind of associative array.

To illustrate, here is working R code. The comments describe the creation of a simple "AA" (associative array) class for storing intermediate results. Vectors i are converted to strings and those are used to index into a list E that will hold all the values. The ij operation is implemented as %.%.

These preliminaries enable the recursive function e to be defined rather simply in a way that parallels the mathematical notation. In particular, the line

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

(1)1R1 rather than 0.

O tempo mostra que é preciso 0.01 seconds to compute e(c(0,0,0,6)); its value is

32.6771634160506

Accumulated floating point roundoff error has destroyed the last two digits (which should be 68 rather than 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))

Finalmente, aqui está a implementação original do Mathematica que produziu a resposta exata. A memorização é realizada através da e[i_] := e[i] = ...expressão idiomática , eliminando quase todas as Rpreliminares. Internamente, porém, os dois programas estão fazendo as mesmas coisas da mesma maneira.

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
fonte
5
+1 Eu imagino que algumas das notações seriam difíceis de seguir para os alunos que fizeram essa pergunta (não que eu tenha alguma alternativa concreta para sugerir agora). Por outro lado, imagino o que eles pretendiam fazer com essa pergunta.
Glen_b -Reinstala Monica
11
@Glen_b Eles poderiam aprender muito rolando os dados (e contabilizando os resultados). Parece uma boa maneira de manter a classe ocupada por meia hora enquanto o professor descansa :-).
whuber
12

A versão original desta pergunta começou a vida perguntando:

quantos rolos são necessários até que cada lado apareça 3 vezes

Obviamente, essa é uma pergunta que não tem resposta como o @JuhoKokkala comentou acima: a resposta é uma variável aleatória com uma distribuição que precisa ser encontrada. A pergunta foi modificada para perguntar: "Qual é o número esperado de rolagens". A resposta abaixo procura responder à pergunta original: como encontrar a distribuição do número de rolos , sem usar simulação, e apenas usando técnicas conceitualmente simples que qualquer estudante da Nova Zelândia com um computador poderia implementar Why not? The problem reduces to a 1-liner.

Distribution of the number of rolls required ... such that each side appears 3 times

We roll a die n times. Let Xi denote the number of times side i of the die appears, where i{1,,6}. Then, the joint pmf of (X1,X2,,X6) is Multinomial(n,16) i.e.:

P(X1=x1,,X6=x6)=n!x1!x6!16n subject to: i=16xi=n

Let: N=min{n:Xi3i}. Then the cdf of N is: P(Nn)=P(Xi3|n)

i.e. To find the cdf P(Nn), simply calculate for each value of n={18,19,20,}:

P(X13,,X63) where (X1,,X6)Multinomial(n,16)

Here, for example, is Mathematica code that does this, as n increases from 18 to say 60. It is basically a one-liner:

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

... which yields the exact cdf as n increases:

1814889875110199605761928290762544079842304203111983875176319369216211168408491253173748645888223283142988125507799783342082361483465418375609359740010496

Here is a plot of the cdf P(Nn), as a function of n:

enter image description here

To derive the pmf P(N=n), simply first difference the cdf:

enter image description here

Of course, the distribution has no upper bound, but we can readily solve here for as many values as practically required. The approach is general and should work just as well for any desired combination of sides required.

wolfies
fonte