Dice rolls, simulation vs. theory

Question

I just ran into something interesting, and was hoping someone could shed some light on it. I'm trying to calculate rolling at least N on multiple dice. I wrote a simulation script in R that uses the sample(8,1,TRUE) to simulate a d8 (8 sided die) roll and iterates and sums. I also calculated the probability in excel using a normal distribution with variance (for kdN) = k*(N^2-1)/12 and mean = k*(N+1)/2. I used excel function of: 1-Norm.dist(number.to.roll, mean, sd, TRUE). I consistently get a few percentage points higher on the simulation than on the calculation. For example, to roll 8 on a 2d8, simulation succeeds 67.6% of the time (+/- 0.3%). Excel predicts 62.1%. I think it has something to do with dice being integers and excel calculating based on a continuous distribution. I have tried using +1, or -1, or -.5 to get excel to match. It gets closer, but I can't find a constant that work across different dice sizes. Is there a discrete distribution that will model the rolls better? I tried the R "triangle" package, but it gives values very close to the Normal approximation in Excel. My question is less how to get the formula right, than can someone provide me insight as to why there is the difference? This is for my own edification only. Am fiddling with numbers while I wait for parts to show up at work. :) Thanks.

Answer 1

The distribution of 2d8 is discrete triangular.

     x       2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
64.P(X=x)    1    2    3    4    5    6    7    8    7    6    5    4    3    2    1    
       %   1.56 3.13 4.69 6.25 7.81 9.38 10.9 12.5 10.9 9.38 7.81 6.25 4.69 3.13 1.56 

If you need an algebraic expression for it, for 2d8 it's:

$p(x) = P(X=x) = \frac{1}{64} \min(x-1,17-x)$

As you add more dice, the cdf becomes closer and closer to a normal distribution, but if you want to use normal distributions to approximate probabilities for it, I'd suggest using a continuity correction. If you don't get too far into the tail, it should work pretty well for more than 3 dice. However, it's not all that hard to do the convolution - or even complete enumeration - by hand for small numbers of dice to get exact answers. e.g. here's me doing 3d8 in R:

 o2d8=c(outer(1:8,1:8,"+"))
 o3d8=c(outer(o2d8,1:8,"+"))
 table(o3d8)
o3d8
 3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
 1  3  6 10 15 21 28 36 42 46 48 48 46 42 36 28 21 15 10  6  3  1 

The table at the end shows the number of ways (out of $8^3=256$) of getting each result on 3d8. You convert them to probabilities by dividing by $8^3$:

> print(table(o3d8)/8^3,d=2)
o3d8
     3      4      5      6      7      8      9     10     11     12     13     
0.0020 0.0059 0.0117 0.0195 0.0293 0.0410 0.0547 0.0703 0.0820 0.0898 0.0938 

    14     15     16     17     18     19     20     21     22     23     24 
0.0938 0.0898 0.0820 0.0703 0.0547 0.0410 0.0293 0.0195 0.0117 0.0059 0.0020 

The algebraic formula becomes relatively complicated past 3 dice, and I wouldn't bother with it, but the probabilities are easy to work out in either R or Excel (or any other tool with the relevant ability to do the kind of calculations you need)

Let's say I want to compute probability of rolling at least 9 on 3d8 from a normal approximation (I suggested more than 3 dice, but let's try it anyway). The exact answer is easily computed by summing probabilities in the above table (it's 0.890625).

> pnorm(8.5,3*(8+1)/2,sqrt(3*(8+1)*(8-1)/12),lower.tail=FALSE)
[1] 0.896144

(The use of 8.5 rather than 9 is because of the continuity correction).

That's not too bad, a relative error of a little over half a percent. But the exact answer takes only a few seconds longer to generate.

---

Incidentally, a simulation script for rolling 2d8 in R is as simple as:

r2d8=replicate(10000,sum(sample(8,2,TRUE)))

And to display a table of the results as proportions:

table(r2d8)/length(r2d8)

The results of the simulation can be seen as red circles, compared with the exact values (black dots):

--

There's detailed instructions on how to use Excel to do similar calculations to the ones I did in R to compute the exact probabilities here