Throwing a die for a cookie; Acceptance-rejection policies

We have a class with \(N=19\) children and one biscuit. We would like to select one child to give the biscuit, but fate should decide which child. Can we use a fair six-sided die to select a winner, with multiple throws, such that all children have an equal probability of winning? And if so, what strategy would be best, in the sense that it minimizes the expected number of throws?

To answer these questions, we start with the analysis of simple cases to learn, hopefully, how to handle the case with \(19\) children. In this process we design a simple and intuitive strategy. As this policy is somehwat hard to analyze mathematically, we'll use simulation to assess its quality. It turns out that this policy is unfair: some children will have a larger probability of winning than others. Then we retract a bit and provide a better algorithm, which is fair and optimal in the sense that is requires the least number of throws on average to select the winner. {{{sidenote{(I discussed this problem with Henk Tijms and he pointed out an answer on Quora on how to handle a class with 20 pupils. Bases on this, I generalized this answer to an algorithm that can deal with an arbitrary number of children.)}}}

In this section we introduce some general concepts that are generally useful. The first is recursion, the second is to use acceptance-rejection policies to decide whether a game should continue or stop.

The first aim is to find a policy that selects in a fair way a winner. But before embarking on how to do this for general classes, we start with considering some simple problems. Perhaps we spot a pattern!

Obviously, when there is \(N=1\) child, there is no need to choose. The next simplest case is when we have \(N=6\) children: number the children from 1 to 6, and throw the die.

When \(N=4\), we can use the following policy. Throw the die for each child. Each child whose outcome is below \(5\) `survives' to the next round. When there are no survivors after the first round (all threw \(5\) or \(6\) ) let them all throw once again. When there are \(1\), \(2\), or \(3\) survivors, use the suitable strategy we developed above. Thus, the idea of this algorithm is like this: use the die to reduce the group such that a set of children remains for which we already have a policy. However, while this is a fool-proof strategy, it requires at least \(4\) throws. We'd better search for more efficient policies.

Note that the first throw splits the group into two subgroups \(\{1,2\}\) and \(\{3, 4\}\) of the same size. Then we have two smaller groups, and for these smaller group we already know how to find a fair solution. Clearly, the pattern of recursion is to try to express the solution in terms of solutions of problems we solved earlier.

This second strategy is known as an acceptance-rejection policy: accept an outcome of an experiment when the outcome lies in a certain subset of the state space, otherwise reject the outcome, and sample again.¹

We can compute the expected number of throws for the second policy with, so-called, first-step analysis. Write \(Y\) for the number of throws we need for the die to hit the set \(\{1,2,3,4\}\). The probability to hit this set is \(p=2/3\), while the probability to hit \(\{5,6\}\) is \(q=1-p\).

In probability terms, we say that \(Y\) has the first success distribution with failure probablity \(q\). This is rv is similar to a geometric rv, but slightly different. The first success rv counts the number of throws up and including the number of throws until a success; the geometric rv counts the number of failures until the success, so one less.² From this,³

So, for \(N=4\), \(\E Y = 3/2\).

Interestingly, the first strategy, based on recursion, always needs exactly \(2\) throws. Instead, the acceptance-rejection method requires less throws on average, but there is no formal bound on \(Y\). We remark in passing an important take away in the design of random algorithms. Sometimes it is required that an algorithm finishes in a guaranteed fixed amount time, rather than that it finishes fast. When time constraints are relevant, it's good to have a deterministic algorithm, even though it can be slower on average.

Observe that with the strategy discussed in 1 for \(N=19\), we always have to throw at least three times. However, the strategy of 1 sometimes completes in just two throws. Might there be a optimal policy that needs somewhere between \(2\) and \(3\) throws on average?

One promising idea to improve on the strategy of using complete rounds of two throws for \(N=19\) is as follows. If there is no winner after the first and second throw, just throw the die a third time. Now check whether the result of the second and third throw is in the list of winning outcomes. To illustrate, if \((X_1, X_2) = (4, 3)\), then, by 1, the outcome \((4,3)\) is not a winning sequence, so that we have to continue with a third throw. After seeing \(X_3=4\), there is winning sequence, since \((X_2, X_3) = (3,4)\), which corresponds to child \((3-1)\cdot 6 + 4 = 16\). If, however, \((X_1, X_2, X_3) = (4, 5, 1)\) there is not yet a winner, because neither \((4,5)\) nor \((5,1)\) correspond to a child. Supposing that \(X_{4}=3\), the result of the last two throws is \((X_3, X_4) = (1, 3)\) so child \(3\) wins.

Notwithstanding that this is an appealing strategy, I was not convinced that this strategy will result in all children having an equal probability to win. The reason is that it reminded me of the, so-called, `monkey-typing problem'. In a later chapter we will explore this problem, but since its analysis requires some new mathematical ideas, we will use here simulation to see how our this new `monkey' policy performs when \(N=19\).

Before we discuss the code below, you should know that in Python (like in C and many other programming languages) arrays start at index 0. This will turn out very practical later in this chapter when we use modulo computions. Hence, in all maths and the code below, children start couting at \(0\), and if we have, for example \(4\) children, the associated set will be \(\{0, 1, 2, 3\}\), and not \(\{1, 2, 3, 4\}\).

The main parts of the next python code work as follows. The tuple last_two corresponds to the outcomes of the last two throws. The rejection set reject contains the tuples that do not correspond to a winner, hence require at least one further throw; hence, reject contains the outcome pairs \((3,1), \ldots, (5,5)\). If last_two is an element of reject, we throw the die again and update last_two to include the last throw. We continue with throwing until last_two is accepted. Once we have a winner, we add one win for this child. Finally, the count object counts how often each outcome occurs.⁴

import collections
import itertools

import numpy as np

rng = np.random.default_rng(1)

reject = set().union((3, i) for i in range(1, 6))
for i, j in itertools.product([4, 5], range(6)):
    reject.add((i, j))

count = collections.Counter()
num_runs = 100_000
for _ in range(num_runs):
    last_two = (rng.integers(6), rng.integers(6))
    while last_two in reject:
        last_two = (last_two[1], rng.integers(6))
    count[last_two] += 1
most = count.most_common()[0][1]
least = count.most_common()[-1][1]
mean = num_runs // 19

print(f"{len(reject)=}, {len(count)=}")
print(f"{least=}, {mean=}, {most=}, {count.total()=}")

len(reject)=17, len(count)=19
least=4464, mean=5263, most=5677, count.total()=100000

We print three checks: the size of the rejection set is \(17\), as it should, the counter counted \(19 = 35 - 17\) different outcomes, and we see \(10^5\) samples in total. The difference between the children that win the most and the least is about 900, i.e., about \(1\%\) for a sample size of \(10^{5}\) throws. Is this large? To provide some context, we can compare this with a simulation of uniform numbers on \(\{0, 1, \ldots, 18\}\).

import numpy as np

rng = np.random.default_rng(1)
N = 100_000
X = rng.integers(19, size=N)
bins = np.bincount(X)
print(f"{min(bins)=}, {max(bins)=}")

min(bins)=5144, max(bins)=5363

This variation between the numbers that occur the most and least often is quite a bit smaller than the variation under the monkey strategy. We can do further statistical tests, but for the moment I simply conclude that, under the monkey-typing strategy, the children don't have uniformly distributed winning chances.

Still I find it interesting to analyze how many throws this monkey strategy needs on average tof inish, in particular how it compares to the acceptance-rejection policy designed in 1. If it slower, then the policy to use (multiple) rounds of three throws times is better in all respects, but if the monkey policy is faster, we might make a trade off between fairness and speed. Note that this is a common theme in optimization: nearly always there are multiple objectives that we like to minimize or maximize, but we cannot perform optimally on all dimenstions at the same time. If so, we need to make a choice which objective we prefer over another.⁵

Returning to the analysis of the average number of throws it takes the monkey to select a winner, consider the code below. The next exercise asks you to explain some parts of it.

T = 0
for _ in range(N):
    last_two = (rng.integers(6), rng.integers(6))
    T += 2
    while last_two in reject:
        last_two = (last_two[1], rng.integers(6))
        T += 1

print(f"{T/N=}")

T/N=2.94335

As T counts the total number of throws for \(N\) games, T/N is the average number of throws per game. Apparently, this lies just a bit below \(3\). In summary, although unfair, the monkey policy needs less throws on average than the best acceptance-rejection policy so far, as this already needs three throws for the first round.

In the above we learned that when the class contains \(N=19\) children, the monkey strategy, which prescribes to continue throwing the die until a winning pattern emerges, does not give a fair outcome. Hence, for the moment, by 1 the best fair algorithm for \(N=19\) seems to throw three times: if that gives a winner, stop; otherwise throw three times again, and continue with rounds until there is a winner.

However, there is a strategy that is fair and minimizes the expected number of throws at the same time. This strategy we discuss now. Although a bit harder to understand, it's an elegant interplay between probability and number theory.

In the sequel we write \(Z\sim\Unif{\{0, 1, \ldots, n\}}\) to say that a rv \(Z\) is uniformly distributed on the set \(\{0, 1, \ldots, n\}\). We also write \(Z|A\) to denote a rv \(Z\) given that the event \(A\) materialized.

Let us first estimate the number of throws needed for a class with \(N=20\) children; this is somewhat simpler than the case with \(N=19\). Consider the following code, where \(X\) is an iid uniform rv on \(\{0,\ldots,5\}\) corresponding to the throw of the die. The line t % 20 computes the remainder of t after dividing by \(20\).

import numpy as np

rng = np.random.default_rng(1)

t = 0
while t < 16:
    X = rng.integers(6)
    t = 6 * t + X
t = t % 20

The above two exercises explain the threshold \(16\) of the algorithm for \(N=20\). Now that we understand this, we can try to generalize the algorithm to other values for \(N\), such \(N=19\). Consider the following update rule for \(i=1, 2, \ldots\), \(t_0=0\), and \(X_1\sim\Unif{\{0, \ldots, 5\}}\),

Notice that for \(N=19\) the threshold is not the same for each throw, i.e., 17 for the second and 7 for the third throw.

Clearly, to deal with cases such as \(N=19\) we need to modify the above algorithm. Here is a generalized version.⁶

def select_uniformly_with_a_die(N):
    n = r = 1
    t = 0
    while t < r:
        X = rng.integers(6)
        t = 6 * t + X
        n = 6 * r
        r = n % N
    return t % N

Let us test this generalized algorithm numerically. Here is the code.

count = collections.Counter()
for _ in range(num_runs):
    s = select_uniformly_with_a_die(19)
    count[s] += 1

most = count.most_common()[0][1]
least = count.most_common()[-1][1]
mean = num_runs // 19

print(f"{len(count)=}, {count.total()=}")
print(f"{least=}, {mean=}, {most=}")

len(count)=19, count.total()=100000
least=5144, mean=5263, most=5434

If we compare this output to our earlier results, the difference between the least and most winning child is reasonable. In fact, comparing the occurrences, we see that our algorithm performs very well!

The last step is to estimate the expected number of throws \(\E Y\) required for the case \(N=19\). First, however, we consider \(N=20\) because in this case we can obtain a simple closed-form expression for \(\E Y\).

To compute \(\E Y\) for general \(N\) we can slightly modify the above algorithm.

N = 20
n = r = 1
alpha, eps = 1, 1e-3
EY = 0
while alpha >= eps:
    EY += alpha
    n = 6 * r
    r = n % N
    alpha *= r / n
    print(f"{n=:3d}, {r=:3d}, {r/n=:2.4f}, {alpha=:2.4f}, {EY=:2.4f}")

n=  6, r=  6, r/n=1.0000, alpha=1.0000, EY=1.0000
n= 36, r= 16, r/n=0.4444, alpha=0.4444, EY=2.0000
n= 96, r= 16, r/n=0.1667, alpha=0.0741, EY=2.4444
n= 96, r= 16, r/n=0.1667, alpha=0.0123, EY=2.5185
n= 96, r= 16, r/n=0.1667, alpha=0.0021, EY=2.5309
n= 96, r= 16, r/n=0.1667, alpha=0.0003, EY=2.5329

This seems to converges to \(2.533\). The last task is to run the code for \(N=19\). This gives \(\E Y=2.509\). Comparing this to the monkey strategy of 1 we see that this algorithm does better on all respects: it is fair and requires less throws on average.

And now two final questions.

By studying the above toy problem we covered many nice and useful ideas. We developed acceptance-rejection strategies to select a winner. Then we tried to speed up this algorithm to the monkey algorithm, but this algorithm turned out to be unfair. However, by improving the idea, we were able to find a fair and optimal algorithm. This final algorithm is not quite trivial, but effective and suprisingly short.

In the next section we'll look at the monkey-typing strategy. This will provide us with a number of useful tools, in particular first-step analysis.