# Two dice problem

Two players are throwing each one die. The one that has higher value, receives a number of points equal to the difference of the values on both dice.

How do you estimate probability for a winner to collect more than 100 points in 100 throws?

• Is this a homework problem? Sure sounds like it, so, it should have the self-study tag Commented Jul 8, 2013 at 19:15
• The answer--if I'm correct this time ;-)--will be around 0.713.
– whuber
Commented Jul 8, 2013 at 20:37
• Ivana - the point with the CLT is that you avoid the $n^{th}$ order convolution. Can you work out the probability that a given person (without knowing they won or lost overall) collects > 100 points in 100 throws? More basically, but more to the point, can you work out the (approximate) distribution of their total score? Failing that, can you work out the mean and variance of the distribution of the number of points they gain from a single throw? Commented Jul 8, 2013 at 22:38
• It's not too difficult to also work out the covariance of the per-round score of the two participants, and hence of the total score. From there, it's not difficult to compute the relevant probability - but there is probably an easier way to do it than that. Commented Jul 8, 2013 at 22:45
• @COOLSerdash If you treat the pair of scores as bivariate normal, could you work out $P(\text{both players score} < 100.5)$? Commented Jul 9, 2013 at 8:50

Begin with an array of the possible (equiprobable) outcomes of the two dice, with my die indexing the rows and your die indexing the columns. In this table, the entries are the payoffs to each of us: the first is my payoff $X$ and the second is your payoff $Y$.

$$\begin{array}{c|cccccc} & 1 & 2 & 3 & 4 & 5 & 6 \\ \hline 1 & \left( \begin{array}{cc} 0 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 1 \end{array} \right) & \left( \begin{array}{cc} 0 & 2 \end{array} \right) & \left( \begin{array}{cc} 0 & 3 \end{array} \right) & \left( \begin{array}{cc} 0 & 4 \end{array} \right) & \left( \begin{array}{cc} 0 & 5 \end{array} \right) \\ 2 & \left( \begin{array}{cc} 1 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 1 \end{array} \right) & \left( \begin{array}{cc} 0 & 2 \end{array} \right) & \left( \begin{array}{cc} 0 & 3 \end{array} \right) & \left( \begin{array}{cc} 0 & 4 \end{array} \right) \\ 3 & \left( \begin{array}{cc} 2 & 0 \end{array} \right) & \left( \begin{array}{cc} 1 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 1 \end{array} \right) & \left( \begin{array}{cc} 0 & 2 \end{array} \right) & \left( \begin{array}{cc} 0 & 3 \end{array} \right) \\ 4 & \left( \begin{array}{cc} 3 & 0 \end{array} \right) & \left( \begin{array}{cc} 2 & 0 \end{array} \right) & \left( \begin{array}{cc} 1 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 1 \end{array} \right) & \left( \begin{array}{cc} 0 & 2 \end{array} \right) \\ 5 & \left( \begin{array}{cc} 4 & 0 \end{array} \right) & \left( \begin{array}{cc} 3 & 0 \end{array} \right) & \left( \begin{array}{cc} 2 & 0 \end{array} \right) & \left( \begin{array}{cc} 1 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 1 \end{array} \right) \\ 6 & \left( \begin{array}{cc} 5 & 0 \end{array} \right) & \left( \begin{array}{cc} 4 & 0 \end{array} \right) & \left( \begin{array}{cc} 3 & 0 \end{array} \right) & \left( \begin{array}{cc} 2 & 0 \end{array} \right) & \left( \begin{array}{cc} 1 & 0 \end{array} \right) & \left( \begin{array}{cc} 0 & 0 \end{array} \right) \end{array}$$

There are only eleven combinations of payoffs $(X,Y)$. By counting them up, we can make a table of them and their probabilities. To simplify the presentation, the next table expresses all probabilities as parts per $36$, so that for example $5$ represents a probability of $5/36$:

$$\begin{array}{l|llllll} & 0 & 1 & 2 & 3 & 4 & 5 \\ \hline 0 & 6 & 5 & 4 & 3 & 2 & 1 \\ 1 & 5 & 0 & 0 & 0 & 0 & 0 \\ 2 & 4 & 0 & 0 & 0 & 0 & 0 \\ 3 & 3 & 0 & 0 & 0 & 0 & 0 \\ 4 & 2 & 0 & 0 & 0 & 0 & 0 \\ 5 & 1 & 0 & 0 & 0 & 0 & 0 \end{array}$$

Applying the (elementary) definitions of mean, variance, and correlation to this tabulated distribution of $(X,Y)$, we find that the common mean of $X$ and $Y$ is $35/36$, their common variance is $2555/1296$ (almost $2$), and their correlation is $-35/73$. (The strong negative correlation makes this question a little trickier than one might expect.)

By virtue of the basic additive properties of means and variances of independent random variables, our payoffs after $100$ independent plays therefore have a joint distribution with means $3500/36$, variances $255500/1296$, and the same correlation $-35/73$. Let's approximate this by a bivariate normal distribution having those parameters. We can readily visualize it: it must be centered around $97$ (which is about $100$ times $35/36$) on both axes, have spreads around $\sqrt{255500/1296}\approx 14$, and have good negative correlation.

In this contour plot of the bivariate normal density, the region where $\text{round}(\max(X,Y)) \gt 100$ is shaded. Its integral approximates the answer to the problem. To the eye, the missing integral looks like about a quarter of the total, giving a rough preliminary estimate of $1 - 1/4 = 0.75$.

This much was completely elementary and can be done with mental arithmetic. Numeric integration over the shaded region (from 100.5 out to infinity in both variables) yields $0.723$ as a possibly more accurate approximation; a computer helps with this calculation ;-). The value must be a little high because the positive skewness of the distributions of $X$ and $Y$ will translate to positive skewness in the marginal distribution of both our scores after $100$ rounds, pulling a little of the probability out of the unshaded area of the figure into the shaded area.

For reference, here is an exact solution. The $100$ outcomes in this game can be obtained sequentially as follows:

1. Determine the number of times we tied. From the tables it is clear the chance of a tie is $6/36 = 1/6$. This number therefore has a Binomial$(100, 1/6)$ distribution.

2. Conditional on the number of ties $k$, the number of times I got a positive payoff, $j$, has a Binomial$(100-k, 1/2)$ distribution. The number of times you got a positive payoff equals $100-k-j$.

3. Conditional on $j$, the distribution of my payoffs is $1$ with $5/15$ chance, $2$ with $4/15$ chance, ..., and $5$ with $1/15$ chance. Its probability mass generation function (pmf) is $p(x) = (5x + 4x^2 + 3x^3 + 2x^4 + x^5)/15$.

The probability that my payoff is $100$ or less therefore equals the polynomial of degree $100$,

$$p(x)^j \mod x^{101},$$

evaluated at $x=1$. That's a reasonably fast calculation. Replacing $j$ by $100-k-j$ gives the probability that your payoff also is $100$ or less. The product of these two probabilities is the chance, conditional on $k$ and $j$, that we both had payoffs of $100$ or less. Therefore, computing the double sum over $k$ and then $j$, we can obtain the chance that the maximum does not exceed $100$. Subtracting that from $1$ gives the answer, equal approximately to

$$0.7131559594164915144059723.$$

We see the bivariate normal approximation was about $1.4$% too high.

### Edit

For those who like details, these might help:

• Define $t(j,m)$ to be the chance of a payoff of $m$ or less conditional on having received exactly $j$ positive payoffs. Its pmf $q(x;j,m)$ is the sum of all monomials of degrees up to and including $m$ in the expansion $p(x)^j$, the $j^\text{th}$ power of the fifth degree polynomial $p$. Therefore $t(j,m) = q(1;j,m)$. Using Mathematica I computed it as

(Expand[((5x + 4x^2 + 3x^3 + 2x^4 + x^5)/15)^j] + O[x]^(m+1) // Normal) /. {x -> 1}

• Conditional on $k$, the chance that both of us score $m$ points or less when there is some positive score in $100-k$ throws is obtained by summing over all the possibilities:

$$r(k,m) = \sum_{j=0}^{100-k} \binom{100-k}{j}2^{k-100}\ t(j,m)\ t(100-k-j,m).$$

• Summing these over all possibilities of $k$ gives the answer,

$$\sum_{k=0}^{100} \binom{100}{k}\left(\frac{1}{6}\right)^k\left(\frac{5}{6}\right)^{100-k} r(k,100).$$

### Reality check

It's easy to make mistakes with calculations like these--I often do. As a check, let's simulate the 100 throws a million times and track the winner's score. This R code will do it in less than a half minute:

set.seed(17)
s <- replicate(10^6, {
delta <- floor(runif(100, 0, 6)) - floor(runif(100, 0, 6));
max(sum(delta[delta > 0]), -sum(delta[delta < 0]));
})
sum(s > 100) / length(s)


The first line within the replicate function simulates the payoffs--positive for me and negative for you--and then the second line calculates which payoff was greater in size. The result is an array of a million winners' payoffs. The final line reports what fraction of those exceeded $100$.

The result, $0.713436$, is only $0.62$ standard errors greater than the exact result previously obtained, providing support for the correctness of this answer.

• Fantastic - a clear, detailed answer, covering all the bases - but also a very beautiful answer, wonderful tables and diagrams. A pleasure to read. Commented Jul 9, 2013 at 23:12
• @whuber Thanks again for this wonderful answer. Forgive my ignorance, but I'm having trouble to understand the exact part. How exactly can we sum over $k$ and $j$? And what do we sum up, the polynomial at $x=1$? Does $p(x)^j$ mean that we take the whole polynomial to the $j$th power? Commented Jul 10, 2013 at 7:26
• @COOL I added some equations for you. I hope there aren't too many typos in them :-).
– whuber
Commented Jul 10, 2013 at 13:28
• @whuber I'm really sorry to bug you further with this, but the expression for which you provided the Mathematica code seems to be 1 whenever $x=1$, irrespective of the values of $j$ and $m$. What am I missing here? Another thing I can't reproduce is the numerical integration: When I use the following Mathematica code, the result is 0.0924: NIntegrate[Integrate[PDF[MultinormalDistribution[{3500/36,3500/36},{{(255500/1296),(255500/1296)*(-35/73)},{(255500/1296)*(-35/73),(255500/1296)}}],{x, y}],{x, 100.5,Infinity}],{y,100.5,Infinity}]. Commented Jul 17, 2013 at 10:50
• @COOL You need to give m and j actual values, rather than symbolic ones, in order to have that Mathematica expression properly evaluated.
– whuber
Commented Jul 17, 2013 at 13:09

If you want to check the solutions given above using a simulation study, here is some useful code for it:

B = 100000
counter = 0

for(i in 1:B){

player1 = sample(1:6,100,replace=TRUE)
player2 = sample(1:6,100,replace=TRUE)

dice.diff = abs(player1 - player2)

player1.wins = ifelse(player1>player2,TRUE,FALSE)
player2.wins = ifelse(player2>player1,TRUE,FALSE)

player1.score = sum(dice.diff[player1.wins])
player2.score = sum(dice.diff[player2.wins])

if(player1.score > 100 | player2.score > 100){
counter = counter + 1
}

}

prob = counter/B


After running the code we obtain the following probability:

> prob
[1] 0.71306

• (1) For a 4X faster simulation use s <- replicate(10^5, { delta <- floor(runif(100, 0, 6)) - floor(runif(100, 0, 6)); max(sum(delta[delta > 0]), -sum(delta[delta < 0])) }); sum(s > 100)/ length(s). (2) I believe the answer is 19024559499426570037611458829573021458074826188755762634838975523692331755102046755611268880392383888911950488196749907935871387917415599926170449751285047 / 26676576488251712549810930719923358271436796584655589343467442433444558096851233688692008708797318581050800966667258114434654077532349702649000974494990336, approximately equal to $0.713155959.$
– whuber
Commented Jul 9, 2013 at 19:05