4
$\begingroup$

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$)

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – V-X Feb 8 '16 at 15:17
  • 1
    $\begingroup$ 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. $\endgroup$ – Ramon Crehuet Feb 8 '16 at 16:14
6
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – BruceET Feb 8 '16 at 18:21
  • 1
    $\begingroup$ 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. $\endgroup$ – BruceET Feb 8 '16 at 18:54
  • $\begingroup$ 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... $\endgroup$ – Ramon Crehuet Feb 9 '16 at 9:09
  • 1
    $\begingroup$ @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). $\endgroup$ – BruceET Feb 9 '16 at 18:26
1
$\begingroup$

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)

enter image description here

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. enter image description here

$\endgroup$
  • $\begingroup$ Might want to browse 'related' pages shown in the right margin for some some interesting combinatorics. $\endgroup$ – BruceET Feb 8 '16 at 18:05
1
$\begingroup$

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 [14]: from numpy.polynomial.polynomial import Polynomial

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

In [16]: p**3  # or alternatively: np.power(p,3)
Out[16]: 
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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.