A Game of Dice

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.

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.

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))
    for choice in choices
        if choice in sums
            return 1
        end
    end
    0
end

This is simple to wrap in count trials:

trials(count, dice, choices) = sum(trial(dice, choices) for _ in 1:count)

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])
  30.062 μs (501 allocations: 39.17 KiB)
91

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

Final setup

We just need to figure out how we make our choices. There's not much to consider other than noting that 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. Additionally, we do not want to waste a slot by repeating a number within a choice. The pool of elements is simply collect(2:12), just the list [2, 3, 4, …, 12]. Our list of choices is then

pool = collect(2:12)
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 = collect(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

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> run(1_000_000, 4)
((0.9753700492599015, 2.402324419690625e-8, (4, 6, 8, 10)),
(0.5990128019743961, 2.401957444579457e-7, (2, 3, 11, 12)))

This means that the optimal choice of numbers is $(4, 6, 8, 10)$. We estimate a win probability of 97.5% for this, and with a tremendously low variance (as expected after a million Remarkably, it does not include 7! This must be because of how the integer partitions are more complementarily spread out for these.

We also note that the worst choice is $(2, 3, 11, 12)$ — this is of course because we forced no repeats. If we hadn't, $(2, 2, 2, 2)$ or $(12, 12, 12, 12)$ would each be the worst.