# Thought experiment for dice game

Two players in a dice game to see who can roll a total of 60 first, taking turns each rolling 2 dice.

For one player, one die is 4 sided and one is 6 sided. Therefore the average roll for this player will be a 6.

The other player has two 6 sided dice. So this players average roll will be a 7.

If they take turns, each rolling both dice, what % of the time will the first player (a 4-6 dice combo) get to 60 first and what % of the time will the second player (6-6 dice combo) get to 60 first. If they get to 60 on the same roll, the player who goes the furthest wins.

-
A draw is also possible in this game. Rough approximation gives: 83%,15%,2% (for draw). Hope I got no mistake. Will write more details when I'm back home tomorrow if you don't get an answer by then. – user3533 Dec 17 '10 at 17:51
cool, thanks! would be awesome if you could explain the method.. – jim_shook Dec 17 '10 at 18:00
Your question is slightly inconsistent - you say 'they take turns rolling dice', but then say 'get to 60 on the same roll'. The analysis below presumes that a 'turn' is both players rolling their respective dice and adding to their totals, but IMHO the question gets that much more interesting if the players alternate: player 1 first, then player 2, then player 1, etc, with no tiebreaker rule needed... – Steven Stadnicki Dec 17 '10 at 18:47
The program is not hard to understand, and I don't speak R. You just play the game as described. p1Roll and p2Roll return the sum of each player's dice. They are accumulated until one or both go to or over 60. Count one for the winner and repeat 10000 times. – Ross Millikan Dec 17 '10 at 19:01
@steven stadnicki yeah maybe wasnt entirely clear, meant that they simultaneously roll and move in tandem, so no first move advantage. but ya alternating rolls would add some intrigue/complexity – jim_shook Dec 17 '10 at 19:57

A rough answer, from Monte Carlo simulation, is:

"Player 1 wins  83.85 %, Player 2 wins  14.28 %, and  1.87 % are draws."


The following R code generates this output:

N <- 60

p1Wins  <- 0
p2Wins  <- 0
draws   <- 0

n.sim   <- 10000

p1Roll <- function() {
d1 <- sample(1:6,1)
d2 <- sample(1:6,1)
return(d1+d2)
}

p2Roll <- function() {
d1 <- sample(1:6,1)
d2 <- sample(1:4,1)
return(d1+d2)
}

for (i in 1:n.sim) {
p1Total <- 0
p2Total <- 0
while (p1Total < 60 & p2Total < 60) {
p1Total <- p1Total + p1Roll()
p2Total <- p2Total + p2Roll()
}
if (p1Total > p2Total) {
p1Wins <- p1Wins + 1
} else
if (p2Total > p1Total) {
p2Wins <- p2Wins + 1
} else {
draws <- draws + 1
}
}

print(paste("Player 1 wins ",p1Wins/n.sim*100,"%, Player 2 wins ",p2Wins/n.sim*100,"%, and ",draws/n.sim*100,"% are draws."))

-
looks great. any good resources where i can learn how to execute R code like this? – jim_shook Dec 17 '10 at 18:34
@jas225, download R from www.r-project.org, start it, and paste the code directly into the command prompt. – mpiktas Dec 17 '10 at 19:11
dl-ed it, thnks. – jim_shook Dec 17 '10 at 20:04

Here is the exact answer, which I used Mathematica to compute.

$P(player\ 1\ wins) = \frac{118598889714523902216022358617928917633636253614645787377890526452587094871641057161197}{778560366535929033488842048429259732340012411410736514701931793474234814991339289051136}$

$P(player\ 2\ wins) = \frac{858777828556435297558093986193070674512712644399583464104126713012616584348576395258605}{1038080488714572044651789397905679643120016548547648686269242391298979753321785718734848}$

$P(players\ tie) = \frac{63512421616314632416996800666111235287366697985612516983784929048741127433063741783941}{3114241466143716133955368193717038929360049645642946058807727173896939259965357156204544}$

Numerically, that's

.15233101351178373415...
.82727479987591082853...
.02039418661230543732...


Here are the Mathematica rules I used to compute the probability of player 1 winnng. The other two cases are similar.
(sorry for the image)

With these rules in place I then computed p[60,60]. I thought it was simpler to start the players at 60 and count down to 0, rather than the other way round.

The first rule just says if player 1 has reached 0 and he player 2 hasn't, then the probability of player 1 winning is 1. The second rule is similar for player 2.

The third rule handles the case where both players have crossed the finish line. Note that despite the game description, there aren't really turns involved—there are rounds. It doesn't matter which player goes first, so me might as well consider all four dice being rolled simultaneously.

The fourth and final rule is where all the action is. The recursive sum is over all possible outcomes of the four dice (two for each player). The constant $864$ is merely $4\cdot6\cdot 6\cdot 6$. The p[x,y] = thing in the center is a kind of memoization (caching). Without it, the simplistic recursion would take forever.

EDIT: Here is a plot of $P(n)$, the probability that the weaker player wins, when playing to a total of $n$ (instead of $60$), for $n = 1\dots80$. To my surprise, it's not quite monotonic! But a little thought will explain why.

-
Interesting! Could you add into your answer a brief description of the calculations you used to obtain this? – Mike Spivey Dec 19 '10 at 0:01
@Mike, I've added the Mathematica rules and some explanation. – I. J. Kennedy Dec 19 '10 at 1:45
+1: Nice work! And thanks for adding the details. :) – Mike Spivey Dec 19 '10 at 2:00

Here's how I approximated an answer:

We know player-2 usually wins. The expected value of player-2 after 9 games is 63, just over 60. So the expected number of turns the game will last is about 9. Let's consider another version of this game which is played for 9 turns exactly and the winner is the one who got the higher score.

How is the score of each player distributed after 9 games? We need to get the distribution of a sum of 18 random variables (because 9 turns correspond to 18 dice throws). The distribution of a sum of a random variable is exactly the convolution of the distributions. This can be efficiently computed by the convolution theorem which says: compute the Discrete Fourier Transform of each of the distributions, multiply the resulting distributions element by element, then compute the Inverse Discrete Fourier Transform of the resulting distribution.

I computed that in python. Now that we got the distribution of the scores for each player after 9 turns. To compute the probability player-2 wins, I computed the cumulative distribution of the score for each player (i.e: the probability of getting score >= i for each i). Now I just had to sum of each score i, the probability of player-1 getting score < i, multiplied by the probability of player-2 getting score >= i. The probability for a draw is computed similarly.

Repeating the calculation I get about: 81%,16%,3% (but it's just an approximation anyway).

To calculate the accurate numbers you should do this calculation for any number of turns, and then compute the average of these numbers, weighted by the probability of the game to end in any specific number of turns. But I expect the calculation for 9 turns to give a very good approximation, because the distribution of the number of turns is centered around 9 and I expect to get similar results of 8 or 10 turns.

Some notes in case you're not familiar with convolutions and Fourier Transform: Look at a random variable X over some finite probability space of n elements. Think of X as a vector of n real numbers (representing its distribution). Then the operator T, defined as T(Y)=X+Y for any Y (this is not element-wise addition. this is addition of random variables, also called convolution). Then T is a linear operator in $\mathbb{R}^n$. It happens to be diagonizable. The Discrete Fourier Transform is just a change of basis, after which T is a diagonal operator. It may be a little surprising, but this basis does not depend of X (although T does). The i,i element in the diagonal of the matrix representing T in this new basis happens to be Pr(X=i). Computing the compositions of diagonal operators in very easy - just multiply the diagonals element-wise. This is the justification for the convolution theorem. After you're done adding random variables in this new basis, just return to the original basis (i.e: do the inverse Fourier transform) and get the desired probability vector of the addition of all the random variables you added.

Here's how the probabilities are calculated in python (I'm supplying all the routines you need):

def dft(x):
N = len(x)
ret = [0]*N
for k in range(N):
for n in range(N):
ret[k] += x[n]*math.e**((-2*math.pi*(1j)/N)*k*n)
return ret
def idft(x):
N = len(x)
ret = [0]*N
for k in range(N):
for n in range(N):
ret[k] += x[n]*math.e**((2*math.pi*(1j)/N)*k*n)
ret[k] /= N
return ret
def mul_vec(x1,x2):
return [a1*a2 for (a1,a2) in zip(x1,x2)]
def pow_vec(x,n):
ret = [1]*len(x)
for i in range(n):
ret = mul_vec(ret, x)
return ret