Discrete probablility simulations, Grinstead and Snell section 1.1
Nancy K. Stanton
code based in part on Peter Doyle's code
Contents
Random numbers in MATLAB
MATLAB has the command rand to generate random numbers between 0 and 1.
rand
ans =
0.9575
rand(1,n) will produce an 
component vector of random numbers, so to get 5 random numbers, I could give the command
rand(1,5)
ans =
0.9649 0.1576 0.9706 0.9572 0.4854
To get a table of 20 random numbers like the one on p. 3 which has 5 rows of 4 numbers, I can give the command
rand(5,4)
ans =
0.8003 0.9595 0.6787 0.1712
0.1419 0.6557 0.7577 0.7060
0.4218 0.0357 0.7431 0.0318
0.9157 0.8491 0.3922 0.2769
0.7922 0.9340 0.6555 0.0462
Simulating cointossing
I can simulate tossing a fair coin by generating a random number and calling it heads if it is less than 1/2 and tails if it is greater. The MATLAB function cointosses.m simulates this. To see how it works, I can give the command doc cointosses which will pop up the help for the command in a separate window.
doc cointosses
So, if I want to toss a fair coin 20 times, I give the command
cointosses(20,1/2)
The result of tossing a coin 20 times is HTTHTHHHTTHHHTTTHTTH The proportion of heads is 0.5
To see the code, I can give the command
type cointosses
function cointosses(n,p)
% cointosses(n,p) simulates tossing a coin n times, with probability p of heads
A = 0;
heads = 0;
for i=1:n
if rand<=p
A = [A,'H'];
heads = heads+1;
else A = [A,'T'];
end
end
disp(['The result of tossing a coin ', num2str(n) ' times is ', num2str(A(1,2:n+1))])
disp(['The proportion of heads is ', num2str(heads/n)])
Let's see what happens if we toss the coin 1000 times and repeat this a 5 times. We probably don't want to list all 1000 outcomes, just the proportion, which we can do with cointosses2 which is the same as cointosses except that it only displays the proportion of heads.
cointosses2(1000,1/2)
The proportion of heads is 0.516
cointosses2(1000,1/2)
The proportion of heads is 0.482
cointosses2(1000,1/2)
The proportion of heads is 0.479
cointosses2(1000,1/2)
The proportion of heads is 0.481
cointosses2(1000,1/2)
The proportion of heads is 0.505
de Mere's first experiment
Is it reasonable to think that in 4 rolls of a die a 6 will turn up? The command demere1(n) simulates this experiment n times.
demere1(100)
The proportion of trials with a 6 in the first 4 rolls is 0.51
demere1(1000)
The proportion of trials with a 6 in the first 4 rolls is 0.514
demere1(1000)
The proportion of trials with a 6 in the first 4 rolls is 0.497
demere1(1000)
The proportion of trials with a 6 in the first 4 rolls is 0.525
demere1(1000)
The proportion of trials with a 6 in the first 4 rolls is 0.515
de Mere's second experiment
We'll try this 5 times with 10000 trials, 24 rolls per trial.
demere2(10000,24)
The proportion of trials with a double 6 in the first 24 is 0.4974
demere2(10000,24)
The proportion of trials with a double 6 in the first 24 is 0.4781
demere2(10000,24)
The proportion of trials with a double 6 in the first 24 is 0.4847
demere2(10000,24)
The proportion of trials with a double 6 in the first 24 is 0.4999
demere2(10000,24)
The proportion of trials with a double 6 in the first 24 is 0.4869
Now we'll do it with 25 rolls, which is the default for demere2.
demere2(10000)
The proportion of trials with a double 6 in the first 25 is 0.5059
demere2(10000)
The proportion of trials with a double 6 in the first 25 is 0.502
demere2(10000)
The proportion of trials with a double 6 in the first 25 is 0.4994
demere2(10000)
The proportion of trials with a double 6 in the first 25 is 0.5037
demere2(10000)
The proportion of trials with a double 6 in the first 25 is 0.506
The commands cointosses, cointosses2, demere1 and demere2 are not built-in MATLAB commands, so if you want to use them, you'll have to download the m-files and put them in your MATLAB path.