iconEuler Home

Sum of m Dice Throws

It is a very old exercise to compute the probability to get a given sum s in a given number m of dice throws. E.g., for two throws the combinations, which yield the sum of 7 are

Sum of Dice Throws

The total number of combinations is 1/36. This yields chance of 1/6 to get the sum of 7.

We can simulate that in Euler. For this, we take a matrix of nx2 dice throws, sum the rows, make that a row vector, and count the multiplicities of 2 to 12.

>n=1000000; D=sum(randint(n,2,6))'; c=getmultiplicities(2:12,D); c/n
[0.02805,  0.055958,  0.083095,  0.111181,  0.139039,  0.166,
0.139126,  0.111384,  0.083124,  0.055471,  0.027572]

The 6th value belongs to the sum 7, and it should be close to 1/6.

>(c/n)[6], 1/6
0.166
0.166666666667

We get the following experimental fractions for each sum from our sample.

>columnsplot(c/n,lab=2:12):

Sum of Dice Throws

Exact Frequencies

Our first approach to get the exact frequencies of sums in m dice throws involves a function, which creates all possible dice throws and calls another function f$ for each. This is the same approach as in the permutation utilities of Euler. See

Permutations

Our function is recursive.

>function forAllRek (f$,n,x,k) ...
 ## create all combinations of 1 to n in
 ## x[1:k]. For each combination, call f$(x).
 
 if k>1 then
    loop 1 to n
        x[k]=#;
        forAllRek(f$,n,x,k-1);
    end
 else
    loop 1 to n
        x[1]=#;
        f$(x);
    end;
 endif;
 endfunction

We call this starting with k=m.

>function forAll (f$,n,m) ...
 x=zeros(1,m);
 forAllRek(f$,n,x,m);
 endfunction

Let us test by simply printing the combination.

>function pr (x) ...
 x,
 endfunction
>forAll("pr",6,2)
[1,  1]
[2,  1]
[3,  1]
[4,  1]
[5,  1]
[6,  1]
[1,  2]
[2,  2]
[3,  2]
[4,  2]
[5,  2]
[6,  2]
[1,  3]
[2,  3]
[3,  3]
[4,  3]
[5,  3]
[6,  3]
[1,  4]
[2,  4]
[3,  4]
[4,  4]
[5,  4]
[6,  4]
[1,  5]
[2,  5]
[3,  5]
[4,  5]
[5,  5]
[6,  5]
[1,  6]
[2,  6]
[3,  6]
[4,  6]
[5,  6]
[6,  6]

That worked.

Now let us increase a vector v[s] by 1 for each combination x, which has the sum s. We use a global vector v$$ for this.

For simplicity, we take the dice values as 0 to n-1. So we use x-1 instead of x.

>function sumup (x) ...
 s=sum(x-1)+1;
 v$$[s]=v$$[s]+1;
 endfunction

The following function initializes v$$ and does the job.

>function getsums (n,m) ...
 v$$=zeros(1,m*(n-1)+1);
 forAll("sumup",n,m);
 return v$$;
 endfunction

Here are the number of times, each sum from 0 to 10 appeared. We know that the 6th values must be 6.

>getsums(6,2)
[1,  2,  3,  4,  5,  6,  5,  4,  3,  2,  1]

We can plot this, carefully labelling from 2 to 12.

>k=2; columnsplot(getsums(6,k),lab=k:6*k):

Sum of Dice Throws

Here is the result for three dice throws.

>k=3; columnsplot(getsums(6,k),lab=k:6*k):

Sum of Dice Throws

For 6 dice throws, we can get numbers from 6 to 36.

>k=6; columnsplot(getsums(6,k),lab=(0:5*k)+k,width=0.1):

Sum of Dice Throws

The result looks normally distributed. In fact, the central limit theorem tells us, that a normalized sum of equal random variables converges to the 0-1-normal Gauss distribution.

Here, we have the following numbers.

>z=1:6
[1,  2,  3,  4,  5,  6]

The mean value is 3.5.

>sum(z)/6
3.5

The variance of one dice throw is 35/12.

>fraction sum((z-3.5)^2)/6
35/12

Now, we can plot the distribution and add the correct normal distribution with mean k*3.5, and variance k*35/12.

>k=6; v=getsums(6,k); ...
 plot2d((k:6*k)-0.5,v/sum(v),>bar); ...
 plot2d("qnormal(x,k*3.5,sqrt(k*35/12))",>add,color=red):

Sum of Dice Throws

Faster

Of course, the recursive function is not the most speedy way to get the frequencies. However, it requires a bit of thinking to get a faster algorithm.

The idea is the following.

Assume, we know the frequencies for m-1 dice throws. Then a sum of s in m dice throws can occur, if the previous m-1 dice throws had sums s-5 to s, and correspondingly the current dice throw was 0 to 5. So

Sum of Dice Throws

which is a nice iterative algorithm. Euler lets us take the sum of 6 adjacent values in a vector by folding with [1,1,1,1,1,1]. We need to extend the vector with 0 to the left and right before the folding, which is correct, since those numbers are impossible to get in m-1 dice throws.

>function getv (n,m=6) ...
 v=ones(1,m);
 for k=2 to n;
    v=fold(zeros(1,m-1)|v|zeros(1,m-1),ones(1,m));
 end;
 return v;
 endfunction

The result is the same, but many times faster.

>getv(6)
[1,  6,  21,  56,  126,  252,  456,  756,  1161,  1666,  2247,  2856,
3431,  3906,  4221,  4332,  4221,  3906,  3431,  2856,  2247,  1666,
1161,  756,  456,  252,  126,  56,  21,  6,  1]

Euler Home