Problem Session,09/25/98: De Méré's Problem

Example:
De Méré's problem is whether or not it is more likely to get at least one double six in 24 throws of a pair of dice or to get at least one six in 4 throws of a single die?

We are going to ask the computer to help `solve' de Mere's problem:
(later we will do the mathematics necessary)
We define the experiment 1: toss of a die 4 times, success occurs if the die comes up with a six at least once. Experiment 2: toss 2 dice 24 times, success occurs if the dice come up with a double six at least once.

We want to compute the probability of success in each case as a long term frequency.

First we show how each experiment can be run by one command in matlab:

Experiment 1:

>>exp1=unidrnd(6,1,4)
exp1 =3     3     3     4

We code it as a success if there is a 6, success is a binary variable which is 1 if a success occurs and 0 if not.

>>sucess=(max(exp1)==6)
sucess =
     0

Now in order to run this a large number of times, we are going to loop and repaet these instructions many time and stock them in a vector, result1 that we initialize to be all zeros at first:

>>result1=zeros(1,10)
result1 =
     0     0     0     0     0     0     0     0     0     0
>>result1=zeros(1,1000)

>>for i =(1:1000)
exp1=unidrnd(6,1,4);
result1(i)=(max(exp1)==6);
end

>>sum(result1)/1000
ans =
    0.5220

After having created the simulation line by line we would like to stock it for future use, this is the procedure:

1.

You can open your favorite editor (say emacs) from within matlab by typing at the prompt:

>>!emacs demere1.m

Or if you in a windowing system, you can open a window at the same time by:

>>!emacs demere1.m &

(You may also open up another window, and edit from there, careful that you are editing under your current directory or matlab won't see the function).

2.

function out=demere1(nsimul)
%DeMere's first experiment: one dice
%nsimul is the input=number of simulations
%This function outputs the proportion of success

result=zeros(1,nsimul);
     for i=1:nsimul
     exp1=unidrnd(6,1,4);
     
result(i)=(max(exp1)==6);
end
out=sum(result)/nsimul;

3.

Once you have saved your function, you can just call it with matlab by typing:

>>demere1(10)
ans =
    0.6000
>>demere1(1000)
>> demere1(100)
ans =
    0.5500

Experiment 2:
Here is the function to try for the second experiment, same procedure but this time save the file as demere2.m.

function out=demere2(nsimul)
result2=zeros(1,nsimul);
     for i=(1:nsimul)
      twod=sum(unidrnd(6,2,24));
      result2(i)=(sum(twod==12)>0);
     end
out=sum(result2)/nsimul;

Now try the function out with:

>> demere2(100)
ans =
    0.5200
>> demere2(10000)  % Take a break.....
ans =
    0.4859