# Distribution of the sum of $N$ loaded dice rolls

I would like to calculate the probability distribution of the sum of all the faces of $N$ dice rolls. The face probabilities ${p_i}$ are know, but are not $1 \over 6$. I have found answers for the case of a fair dice (i.e. $p_i={1 \over 6}$) here and here

For large $N$ I could apply the central limit theorem and use a normal distribution, but I don't know how to proceed for small $N$. (In particular, $N=2,4, 20$)

• Are you sure you know what you are asked to calculate? Please, tell us what have you already done and what are you trying to do. – V-X Feb 8 '16 at 15:17
• I guess I simply want to know how the distribution in math.uah.edu/stat/apps/DiceExperiment.html is calculated. It's probably simple but I don't know how to do it. – Ramon Crehuet Feb 8 '16 at 16:14

You can use generating functions.

Let $P=p_1x+p_2x^2+p_3x^3+p_4 x^4+p_5 x^5 +p_6 x^6$ where $p_i$ is the probability of $i$ occurring when rolling the die once.

Then the coefficient of $x^k$ in $P^N$ gives the probability of rolling a sum of $k$ when rolling the die $N$ times and summing.

For example, suppose $P=\frac{1}{7}x+\frac{1}{7}x^2+\frac{1}{7}x^3+\frac{1}{7}x^4+\frac{1}{7}x^5 + \frac{2}{7}x^6$.

Then, using a computer algebra system (I like PARI/GP), we find $P^3 = \frac{8}{343} x^{18} + \frac{12}{343} x^{17} + \frac{18}{343} x^{16} + \frac{25}{343} x^{15} + \frac{33}{343} x^{14} + \frac{6}{49} x^{13} + \frac{40}{343} x^{12} + \frac{39}{343} x^{11} + \frac{36}{343} x^{10} + \frac{31}{343} x^9 + \frac{24}{343} x^8 + \frac{15}{343} x^7 + \frac{10}{343} x^6 + \frac{6}{343} x^5 + \frac{3}{343} x^4 + \frac{1}{343} x^3$.

From this, we can conclude, for instance, that the probability of a sum of $10$ when rolling $3$ times (or rolling once with three identical copies of this die) is $\frac{36}{343}.$

(Using Bruce's example, we get $P^3=\frac{1}{64} x^{18} + \frac{3}{64} x^{17} + \frac{3}{32} x^{16} + \frac{1}{8} x^{15} + \frac{9}{64} x^{14} + \frac{9}{64} x^{13} + \frac{25}{192} x^{12} + \frac{7}{64} x^{11} + \frac{5}{64} x^{10} + \frac{91}{1728} x^9 + \frac{19}{576} x^8 + \frac{11}{576} x^7 + \frac{1}{108} x^6 + \frac{1}{288} x^5 + \frac{1}{576} x^4 + \frac{1}{1728} x^3$, and so the probability of $10$ is $\frac{5}{64}=0.078125$ (exactly).)

• Nice approach. But are you sure of your answer for a sum of ten? After my initial simulation mean(t == 10) returns 0.078097, whereas $36/343 \approx 0.1049563.$ Too far away. – BruceET Feb 8 '16 at 18:21
• OK, we're using differently biased dice. With your bias, I get 0.104877 for a total of 10 and 0.02338 for a total of 18. Got rid of one of by Comments to simplify for readers. – BruceET Feb 8 '16 at 18:54
• As a curiosity, do you have any idea of how could this be implemented in this applet math.uah.edu/stat/apps/DiceExperiment.html I am pretty sure the algebra system is not part of it... – Ramon Crehuet Feb 9 '16 at 9:09
• @RamonCrehuet: Couldn't be sure without looking at the code. Doubt if simulation. Probably, programming successive convolutions one step at a time of the simple discrete distributions. (In my browser this link is partially broken). – BruceET Feb 9 '16 at 18:26

It seems to me the spirit of the original context is experimental. In order to get an analytic answer, even in a simple case with only two rolls, you need to specify precisely how the die is loaded.

Here is a simulation in R of an experiment with $n = 3$ rolls of a die loaded so that faces 1 to 6 appear with probabilities $(1/12, 1/12, 1/12, 3/12, 3/12, 3/12)$. In practice, such a bias might be achieved by embedding a heavy weight just below the corner where faces 1, 2, and 3 meet. A million three-roll experiments are simulated. Each pass through the loop simulates one three-roll experiment. The random variable $T$ is the total on the three dice (values between 3 and 18, inclusive).

For the mean and SD of the best-fitting normal distribution, you should be able to find $E(T) = 12.75$ and $SD(T) = ?$, which are approximated in the simulation. (Also, the exact computation is shown for $E(T) = n\sum_{i=1}^6 ip_i.$)

 m = 10^6;  n = 3;  p = c(1, 1, 1, 3, 3, 3)/12
t = numeric(m) # vector to receive totals
for (i in 1:m) {
faces = sample(1:6, n, repl=T, prob = p)
t[i] = sum(faces)  }
mean(t);  sd(t)
## 12.74981  # Approx E(T) = 12.75
## 2.657790  # Approx SD(T)
hist(t, br=(2:18)+.5, prob=T, main="Sum of 3 Rolls", col="wheat", ylim=c(0,.15))
curve(dnorm(x, mean(t), sd(t)), lwd=2, col="blue", add=T)
3*sum((1:6)*p)
## 12.75     # Exact E(T) The normal curve with matching mean and SD is not yet a good fit for only three rolls. Perhaps much better for 20. Using code $mutatis\;mutandis$: Yes 20 looks better. • Might want to browse 'related' pages shown in the right margin for some some interesting combinatorics. – BruceET Feb 8 '16 at 18:05

As an alternative to Matthew Conroy's suggestion to use a computer algebra system, one can also code this with Python using numpy class of Polynomials. The convolution power can be simply obtained by calculating the power of a polynomial. The same example he suggests can be coded as follows:

In : from numpy.polynomial.polynomial import Polynomial

In : p=Polynomial((1/7, 1/7, 1/7, 1/7,1/7, 2/7))

In : p**3  # or alternatively: np.power(p,3)
Out:
Polynomial([ 0.00291545,  0.00874636,  0.01749271,  0.02915452,  0.04373178,
0.06997085,  0.09037901,  0.10495627,  0.11370262,  0.11661808,
0.12244898,  0.09620991,  0.0728863 ,  0.05247813,  0.03498542,
0.02332362], [-1.,  1.], [-1.,  1.])


Of course, one has to work with floats, but otherwise the results agree.

In numpy, the convolution can also be coded with np.convolve , but the application of successive convolutions is more cumbersome.