# Rolling and Banking

Infront of you are three regular dice (3d6), and you start to play a game. The rules are simple: You initially roll all three dice. Any dice that lands on a six gets ‘banked’ and placed to one side. You continue to re-roll all the unbanked dice until all the dice are banked. What is the expected number of rolls it will take you to bank all three dice?
What is the expected number of rolls to bank all three dice?
We’ll build up to the answer.

## One Die

If we only had one die, what is the expected number of rolls it would take? It’s easy to use common sense and simply state the (correct) answer of six rolls, but let’s derive an expression as we will leverage this principle for more dice.
When we roll one die, one of two things could happen. Either we get a six (probability 1/6), and this takes us one roll, or we do not get a six (probability 5/6), and this takes us one roll, and we are back to where we start. If we define R1 as the expected number of rolls to bank one die, we can write a self-referential equation:
R1 = 5/6(R1+1) + 1/6(1)
Rearranging this, and collecting terms, gives the answer we expect:
6R1 = 5(R1+1) + 1

R1 = 6

## Two Dice

Now let’s consider what happens if we have two dice. There is a 1/36 chance that we roll two sixes (and this takes a single roll). There is a 25/36 chance we roll no sixes (takes a roll, and we are back to where we started). Finally, there is a 10/36 chance that we’ll roll one six (takes one roll), and then we’ll end up with just one die to roll and we already know what the expected number of rolls for one die is (we just calculated it, R1=6).
R2 = 1/36(1) + 25/36(R2+1) + 10/36(R1+1)

36R2 = 1 + 25R2 + 25 + 10(6+1)

R2 = 96/11
It will take approximately 8.7 expected rolls to bank two dice.

## Three Dice

Now we have this down. If we roll three dice, there is a 1/216 chance we roll all three sixes. There is a 15/216 chance we’ll roll two sixes (adding one roll, then converting the problem to a one roll problem). There is a 75/216 chance we’ll roll one six (taking one roll and reducing it to a two dice problem) and, finally, there is 125/216 chance we will roll no sixes (adding a roll, and returning us to the same problem).
Putting this is equation form:
R3 = 1/216(1) + 15/216(R1+1) + 75/216(R2+1) + 125/216(R3+1)

216R3 = 1 + 15(6+1) + 75(96/11+1) + 125(R3+1)

91R3 = 9606/11

R3 = 10556/1001 ≈ 10.555

So, the expected number of rolls to finish the three dice game is ≈10.555

## More Dice

How does the length of the game change if we add more dice? Interestingly, the curve flattens quite nicely. The more dice you have, the more likely you will be to roll a six. Approximately one sixth of all dice get banked each roll (meaning 5/6 survive). What we've just looked at is a simplified* view of radioactive decay. If we started with a bucket of a thousand dice, and each roll we discarded all the sixes, this is similar to how isotopes decay. Radioactive decay is a spontaneous process, can't be predicted individually, but can be described statistically. The number of dice we 'bank' each roll is analogous to the number of atoms that have decayed in that time period. Some isotopes decay rapdily, and some decay slowly. In our analogy, this is like the difference between rolling different sided dice. For a slower decay, you could roll twenty sided dice (or even a hundred sided). For faster decay, you could, for instance, flip a coin and discard all the heads.
Physicists like to talk about decay using a decay constant (typically given the greek symbol λ), which represents the fraction of nuclei that will decay per unit time. This exponential decay manifests as the disappearance of a constant fraction of those surviving per unit time.
Another way to quote the decay is by reference to the Half-Life. This is the time required for a sample to decay to 50% of its initial state. The relationship between the half-life and the decay constant is that the half-life is ln2/λ
Returning to our dice banking experiments, this equates to a half life of approximately four rolls. If we started with a large bucket of dice, and banked all the sixes on each roll, approximately every four rolls we'd be reduced to half of those we had four rolls ago.

## *Simplified View

(Feel free to skip this section)
I said this was a simplified view, because it is. It's also one of those confusing things that is often taught in schools incorrectly (that was disingenuous; I mean in a simplfied way). Radioacitve decay is a stochastic process at the level of a single atom, and happens on a continuous basis, not the quantized analogy modelled by the dice. I've not seen a school text book that describes the discrepancy between the answer usually achieved by the dice experiment (corresponding to the geometric model), and that expected from exponential decay. Usually getting a half-life of "around four" is considered correct within experimental error.
The more accurate description of the half-life of an atom is the time period over which it might have a 50% chance of decaying.
In the dice model, let's say we start with a thousand dice, what we're doing is rolling all thousand of the dice, and acting on these instantaneously for the next time step, then the next, then the next …
We're simplifying things to a geometric progression and assuming all decay happens (or not) at the end of the time step. In this simplified view, on the first time step we have 1000(1-1/6) of the initial thousand remaining. After two time steps, we have 1000(1-1/6)(1-1/6) remaning. After n steps, we have 1000(1-1/6)n remaining.
You can see that, in this model, we initially start with 1,000 dice, then on the next time step we start with 833. We are applying all this decay as though it occured at the end of the time step, whereas, in reality, the decay of these atoms might have occured at any point in this time sample.
Using this simple model, the half life can be calculated by determining the time it takes for half the sample to decay (from 1,000 to 500). From the formula (1-1/6)n = 500/1000 = 0.5, we need to determine n. Taking the log of both sides we get n × ln(1-1/6) = ln(1/2). This gives the answer as n ≈ 3.802 – If you perform this dice experiment in school (or model it in software), this is the answer you will get. This is different to the half-life calculated from using the classic formula ln2/λ ≈ 4.16 with the decay constant λ=1/6. Why the difference?
The more accurate definition of the half-life is that it is the time period over which an atom has a 50:50 chance of decaying. Let's imagine that you roll the thousand dice once an hour. Now imagine that, instead of rolling all 1,000 at the same time once an hour, you roll each dice individually, and continuously, ever 3.6 seconds. As soon as the dice lands six, you discard it. Taken to the limit on an immeasurably large number of atoms, you can see there is a difference.

## Numbers

Here, in tabular form, is a comparison of the two systems. The left column in the time step. We are assuming a starting condition of 1,000 'atoms'. The middle column shows the value calculated by exponential decay using the formula N=N0e-λt, with a decay constant λ=1/6
The right column shows the value calculated by the quantized geometric progression N=N0(1-λ)n. Each is rounded to the nearest integer.
Time
Step
Exponential
Decay
Geometric
Progression
010001000
1846833
2717694
3607579
4513482
5435402
6368335
7311279
8264233
9223194
10189162
The right column shows the expected number of dice remaining if the dice were all rolled at one instant, and the sixes removed. The middle column shows what the expected number of dice would be if each dice just decided to spontaneously evaporate based on when it seemed appropriate to itself based on the decay constant!
This is more visible in this graph. The tabular data above is plotted below. You can see the stair-step caused by the quantized Geometric Progression, compared to the smooth curve of continuous exponential decay. The exponential decay is always higher than the overaly aggresive geometric estimate at the end of each time step because the geometric model assumes all atoms can potentially decay anywhere in the time sample.

## Explanation

For a little bit more explanation, recall the binomial expansion for (1+x)n For the geometric progression of the dice rolls, we set x=-1/6, to give an expansion of (1-1/6)n which gives the following result for the number surving at stage n, from a starting population of N0: Compare this to the exponential decay. The expansion for the exponent function is below, and we know the formula for exponential decay and plugging in the decay constant shows is that for this expansion that x=-t/6 You can see the terms look very similar, but deviate because t2 != (n)(n-1), and so on …
The larger n becomes, the more it diverges.
We can mitigate for this by using a smaller decay constant. If, instead of using a d6, we used a d20, then we can reduced the percentage error down to 2.5% (cf. 8.6% with the d6) as this increases the half-life of the geometric solution to 13.513 cycles and the exponential solution to 13.863
Of course, a downside of this is that an expirement using d20 will take a lot longer to complete than one using a d6.
Here is a table comparing various different sized dice and the effect this has on the two models.
dNGeometric
Half-Life
Exponential
Half-Life
Delta
21.0001.386-27.865%
42.4092.773-13.099%
63.8024.159-8.586%
85.1915.545-6.389%
106.5796.931-5.088%
127.9668.318-4.227%
2013.51313.863-2.521%
10068.96869.315-0.501%
1000692.801693.147-0.050%