﻿ Euler Math Toolbox - Tutorials

# 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 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), 1/6
```
```0.166
0.166666666667
```

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

```>columnsplot(c/n,lab=2:12):
``` # 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

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=#;
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):
``` Here is the result for three dice throws.

```>k=3; columnsplot(getsums(6,k),lab=k:6*k):
``` 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):
``` 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); ...
``` # 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 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