iconEuler Examples

Dice Throws

by R. Grothmann

This notebook shows some Monte Carlo methods in Euler. First we want
to simulate dice throws until we get a 6. To do this, we simply throw
20 times, and determine the number of throws necessary to get a 6.

The trick is to use nonzeros(k==6) to get the indices of all elements
of k, which are equal to 6. Then we take the minimum. Thus we can
avoid a loop.
>k=intrandom(1,20,6); min(nonzeros(k==6))
3
Let us see, if a loop would be faster. So we write a loop, which
does the same.
>function test ...
k=0;
repeat
  k=k+1; 
  if intrandom(6)==6 then return k; endif;
end;
endfunction
Then we time it.
>t=time; loop 1 to 100000; test; end; time-t,
5.413
Try the other method above.
>t=time; loop 1 to 100000; k=intrandom(1,20,6); min(nonzeros(k==6)); end; time-t,
2.917
Indeed the matrix language is a bit faster.

So what is the average waiting time until we get a 6?
>function waitforsix (m=10000) ...
r=zeros(1,m);
loop 1 to m;
  k=intrandom(1,50,6); r[#]=min(nonzeros(k==6));
end;
return r
endfunction
>r=waitforsix(); mean(r),
6.0581
Let us show the distribution of the data.
>plot2d(r,distribution=1):

Average Number of necessary Dice Throws

Of course, the exact answer is 6.

One direct way to show this is to compute the expected value of n,
where the event n occurs, if there is a 6 at the n-th throw, and no 6
before that.

With the default methods, Maxima cannot evaluate the following sum,
which would be our result.
>su &= sum(n*(5/6)^(n-1)*1/6,n,1,inf)|simpsum
                            inf
                            ====     n - 1
                            \     n 5
                             >    --------
                            /         n
                            ====     6
                            n = 1

There is the package simplify_sum, which can handle many sums.
>&load(simplify_sum); &simplify_sum(su)
                                  6

Alternatively, we can use the following trick.
>&assume(abs(x)<1); &sum(x^n,n,1,inf)|simpsum
                                  x
                                -----
                                1 - x

Then we have
>&diff(x^n,x)
                                  n - 1
                               n x

Thus, we get our result in the following way.
>&diffat(x/(1-x),x=5/6)*1/6
                                  6

Let us add the plot of the distribution of the waiting time to our
experimental waiting time.
>n=1:30; plot2d((5/6)^(n-1)*1/6,add=1,thickness=2):

Average Number of necessary Dice Throws

Our next problem is to determine the average waiting time, until
we get all dices at least once.

We setup a Monte Carlo simulation of one try. To determine, when
we have all dices for the first time, we take the maximum of the
indices of the first time appearance of each number.
>function tryone (n:index=100) ...
k=intrandom(1,n,6);
r=0;
loop 1 to 6
  r=max(r,min(nonzeros(k==#)));
end;
return r
endfunction
>tryone
25
Now we do that 10000 times, and denote the results in a vector.
>m=10000; s=zeros(1,m); ...
 loop 1 to m; s[#]=tryone(); end; ...
 plot2d(s,distribution=1):

Average Number of necessary Dice Throws

The following is the mean.
>mean(s)
14.7676
What is the exact answer?

To answer this, assume you have already 4 different dice throws.
How long will you have to wait for the fifth? Using the same reasoning
as above, you get 1/p, where p=2/6 is the chance to throw any
a new dice value.

So all we have to do is add these expectation times starting with
p=5/6. Note we have to add 1 for our first dice throw.
>1+sum(6/(1:5))
14.7

Examples Homepage