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