freeboson: reverie

A (Different) Game of Dice

Sujeet Akula -

This week's Riddler Classic really stumped my intuition. The premise is hard to summarize so I will quote verbatim:

From Ed Carl comes a surprising game of dice:

We’re playing a game where you have to pick four whole numbers. Then I will roll four fair dice. If any two of the dice add up to any one of the numbers you picked, then you win! Otherwise, you lose.

I generally enjoy counting puzzles, but this had some fun quirks and the answer was not clear to me at all.

Totally misinterpreted this

Actually, much of what was fascinating about this to me was due to my own misinterpretation. I thought that the game's description meant that if you had anything other than one match, you lose. I.e., if there was more than one pair of dice that summed to a choice, or to more than one choice, you'd lose the game, and only a single pair of dice matching only a single choice would be a win. That was incorrect. The remainder of this post is basically for the game that was in my head, which I thought was pretty interesting anyway. I posted a corrected version which adjusts the code (and greatly simplifies everything): https://freeboson.org/blog/dice-game/.

Initial thoughts

When you're looking at sums of pair of dice, any Craps player will tell you that 7 is the most common, while 2 and 12 are the rarest. This is because of the relative number of 2-way partitions. From two D6's (regular cubic / 6-faced dice), you can obtain 7 via: $(1, 6)$, $(2, 5)$, $(3, 4)$ and $(4, 3)$, $(5, 2)$, $(6, 1)$. (This puzzle is only about D6's, so we'll just call them die/dice.) Since there are $6^2 = 36$ possible outcomes for rolling two dice, the mode, 7, has a 16.7% chance. Compare that to 2 or 12 which can only be the result of $(1,1)$ (aka snake eyes) and $(6, 6)$ (aka box cars), respectively, meaning a $1/36$ or 2.78% chance.

But in this game, we have to watch out for collisions. Some further immediate observations:

  • If you pick any repeats, that number will always fail
  • It only says to pick four whole numbers, not necessarily in $[2, 12]$
  • You may want to intentionally pick one or more numbers outside that range to avoid collisions

I tried thinking about how I could solve this (semi-)analytically, but with the parameters of the problem, it's just so easy to simulate. Besides, I did not even get to use Julia for my last post.

Basic coding

First the dice rolls. A single die roll would be just rand(1:6), and rolling 4 dice would be rand(1:6, 4). But, I additionally want to see all distinct pairs of the 4 dice rolls. This function performs the dice roll and then provides a generator expression covering those pairs. I'll keep it general for now and allow the number of dice to vary as dice.

function diceroll(dice)
    v = rand(1:6, dice)
    ((v[i], v[j]) for i in 1:dice for j in (i+1):dice)
end

Now, a single trial. Given a selection of numbers choices, perform the dice roll. Returns 1 for success and 0 for failure (so that they can be tallied).

function trial(dice, choices)
    sums = map(sum, diceroll(dice))
    matched = false
    for choice in choices
        if choice in sums
            if matched
                return 0
            end
            matched = true
        end
    end
    Int(matched)
end

This is simple to wrap in count trials:

function trials(count, dice, choices)
    wins = 0
    for n in 1:count
        wins += trial(dice, choices)
    end
    wins
end

So I can try e.g. 100 trials of 4 dice with the choices [2, 5, 6, 9]:

julia> @btime trials(100, 4, [2, 5, 6, 9])
  23.245 μs (401 allocations: 29.78 KiB)
40

For these choices we won 40 out of 100 trials. Seems plausible.

Final setup

We just need to figure out how we make our choices. We probably could just consider all possible choices, but there's some obvious things to note.

  • As already mentioned, though the sum of two dice is in [2, 12], we need not choose only in that range.
  • But if we choose out of the range, it doesn't matter which number we pick, so choosing 1 as explicitly out of the range each time is sufficient.
  • When we do want to make a match, we should never repeat that number within a choice.
  • The ordering of the choices do not matter. E.g. [2, 5, 6, 9] is the same as [2, 5, 9, 6], so no point in repeating those.

It makes sense then to define the pool of elements for the choices as:

pool = [1; 1; 1; 2:12]

This is the list [1, 1, 1, 2, 3, 4, …, 12]. We explicitly include 1 three times, because it is meaningful to repeat 1, and not the other numbers. However, there's no point to have all 1's, since that can never succeed. Our list of choices is then

n = length(pool)
all_choices = ([pool[i], pool[j], pool[k], pool[l]]
            for i in 1:n
                for j in (i+1):n
                    for k in (j+1):n
                        for l in (k+1):n)

Now we can just put it all together, iterating over these choices, and storing the relevant results.

function run(count, dice)
    pool = [1; 1; 1; 2:12]
    n = length(pool)
    choices = ([pool[i], pool[j], pool[k], pool[l]]
               for i in 1:n
                   for j in (i+1):n
                       for k in (j+1):n
                           for l in (k+1):n)
    best = (0.0, 0.0, Vector{Int}())
    worst = (2.0, 0.0, Vector{Int}())
    for choice in choices
        wins = trials(count, dice, choice)
        # Beta(1,1) = uninformed prior
        d = Beta(wins + 1, count - wins + 1)
        p = mean(d)
        if p > best[1]
            best = (p, var(d), choice)
        end
        if p < worst[1]
            worst = (p, var(d), choice)
        end
    end
    (best, worst)
end

The function run does that, but additionally, I could see that the worst non-zero choice would be picking [1, 1, 1, 2] or [1, 1, 1, 12], since that reduces to the least common sum of the pairs of dice rolls, so I wanted to also track the worst choice as a sanity check.

Also, to get an idea of whether our choice of count, i.e. the number of trials is meaningful, it's helpful to understand the distribution in how we estimate the win probability, rather than simply noting the proportion.

Since we are performing repeated Bernoulli trials, we wish to estimate the parameter $\theta$ of the Binomial model. The conjugate prior to the Binomial distribution is the Beta distribution. If we initially consider all choices of the Binomial parameter $\theta \in [0, 1]$ equally likely, then we can express this as the distribution $\text{Beta}(1,1)$.

Once we perform a set of trials for a given choice, if we observe $x$ successes and $y$ failures and a priori considered $\theta \sim \text{Beta}(1,1)$, then our posterior distribution on the parameter $\theta$ is given by \begin{equation} p(\theta|(x,y)) = \text{Beta}(x+1,y+1) \ . \end{equation}
Thus in the function, we compare the mean of the resulting posterior distribution, and also store the variance.

Results

We find after 1 million runs:

julia> data = run(1000000, 4)
((0.659233681532637, 2.2464396073368056e-7, [1, 2, 7, 12]),
(0.13182973634052733, 1.1445031360597356e-7, [1, 1, 1, 12]))

This means that the optimal choice of numbers is $[1, 2, 7, 12]$, which does indeed sacrifice a slot! Notice that it selected the mode of dice pair sums, 7, and then the other two are maximally (and equally) spaced from 7. I wonder if someone could have easily guessed that as the optimal choice from just reading the problem! We estimate that the success probability of that choice is 65.9%. The variance of 2 in 10 million suggests that this is a pretty precise estimate.

We also note that the worst choice is $[1, 1, 1, 12]$ — this one at least I was able to see ahead of time!