Trying to work out the probabilities of a dice game I used to play

Question

At college, my friends and I would sometimes waste time playing a game with dice. We would roll 25 dice, pick out all the dice that landed on a 6, then roll the rest. This would carry on until all the dice were removed. The winner was the person who managed to remove all the dice in the least amount of turns. (A turn counts as rolling the dice - i.e. the first turn is rolling all 25 dice, the second turn is rolling all the dice which didn't roll six on the first turn, etc.)

I want to find the expected amount of turns it would take to finish the game, but I am having difficulty modelling it. I know that if there was one dice then you could model it using the geometric distribution; and the number of sixes on each turn can be modeled using the binomial distribution, but I do not know how to combine the two (or if this is even useful).

Please help, this problem has been bugging me at the back of my head for several years.

Answer 1

More generally, let's say there are N dice, and X is the number of rolls required until every die has rolled a 6. The probability that $n$ or more rolls is required to get a 6 on a single die is $(5/6)^{n-1}$, so $\Pr(X \ge n) = 1 - [1 - (5/6)^{n-1}]^N$. Then

$$\begin{align} E(X) &= \sum_{n=1}^{\infty} \Pr(X \ge n) \\ &= \sum_{n=1}^{\infty} ( 1 - [1 - (5/6)^{n-1}]^N )\\ &= \sum_{n=1}^{\infty} \left( 1 - \left[ 1 - \sum_{i=0}^N (-1)^i \binom{N}{i} (5/6)^{i(n-1)} \right] \right) \\ &= \sum_{i=1}^N \left( (-1)^{i+1} \binom{N}{i} \sum_{n=1}^{\infty}(5/6)^{i(n-1)} \right) \\ &=\sum_{i=1}^N \left( (-1)^{i+1} \binom{N}{i} \frac{1}{1-(5/6)^i} \right) \end{align}$$

Answer 2

Sometimes insight can be gained by sampling the distribution. This doesn't prove anything but it helps validate and supplement other answers. I attached a short python script to sample $10^5$ of your games. In addition to plotting the mean, it is simple to not only estimate the higher order moments, but the actual shape of the distribution:

The results of a Monte-Carlo sampling indicate that @awkward's answer is correct:

Monte-Carlo, mean, std: 21.42546, 6.92985885775

@awkward's answer : 21.4298245245

Code below:

import scipy.misc
import pylab as plt
import numpy as np

N = 25
sample_rolls = 10**5

def sim(n, count=0):
    while n>0:
        count +=1
        n  = (np.random.randint(1,6+1,size=n)!=6).sum()
    return count

ROLLS = np.array([sim(N) for _ in xrange(sample_rolls)])

def awkward_answer(n):
    k = np.arange(1,n+1)
    return (((1.0/(1-(5/6.)**k)))*(-1)**(k+1)*scipy.misc.comb(n,k)).sum()

print "Monte-Carlo, mean, std: ", ROLLS.mean(), ROLLS.std()
print "@awkward's answer     : ", awkward_answer(N)

plt.hist(ROLLS,normed=True,bins=range(0,ROLLS.max()+1),
         histtype='stepfilled')
plt.xlabel("rolls")
plt.show()

Answer 3

For one die, the expected number of rolls is $6$. Let's calculate it to show the technique. Let $E(1)$ be the expected number of rolls for one die to get a $6$. We have $E(1)=\frac 16\cdot 1 + \frac 56(1+E(1))$ because we get one roll with probability $\frac 16$ and with probability $\frac 56$ we return to where we started with one roll under our belts.
Now $E(2)=\frac 1{36}\cdot 1+\frac {10}{36}(1+E(1))+\frac {25}{36}(1+E(2))$ because we have $\frac 1{36}$ chance to get double sixes and be done, $\frac {10}{36}$ chance to get exactly one six and go to the one die case having used one roll, and $\frac {25}{36}$ chance to get no sixes and be one roll along. This becomes $11E(2)=1+10+60+25=96, E(2)=\frac {96}{11}$
It gets messier as we go along. The next equation would be $E(3)=\frac {3 \choose 3}{216}+\frac {{3 \choose 2}5}{216}(1+E(1))+\frac {{3 \choose 1}5^2}{216}(1+E(2))+\frac{{3 \choose 0}5^3}{216}(1+E(3))$
I think you can organize it in a spreadsheet-Excel has the COMBIN function to do the combinations, and I wrote them that way to show how to use it.

Added: you can do what you suggested in your comment. It might be easier. The chance to end in one turn is $(\frac 16)^{25}$ The chance that one die has come up six at least once in two tries is $1-(\frac 56)^2$, so the chance they all have come up six is $\left(1-(\frac 56)^2\right)^{25}$. Generally, the chance you have won by turn $n$ is $\left(1-(\frac 56)^n\right)^{25}$, so the expectation is $$\sum_{n=1}^\infty n\left(\left(1-(\frac 56)^n\right)^{25}-\left(1-(\frac 56)^{n-1}\right)^{25}\right)$$