Thought experiment for dice game

Question

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.

Answer 1

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.

Answer 2

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."))

Answer 3

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
def pad(x,n):
     return x+[0]*(n-len(x))
p1_probs=idft(pow_vec(dft(pad([0]+[1./6]*6,200)),18))
p2_probs=idft( mul_vec(pow_vec(dft(pad([0]+[1./6]*6,200)),9) ,  pow_vec(dft(pad([0]+[1./4]*4,200)),9)    )  )