Likelihood

Random Numbers in MATLAB

Start MATLAB and type the following into the MATLAB prompt

rand

This produces a random number. Repeat several times. Continue typing the material in boxes into the MATLAB prompt.

help rand

Returns help for the rand command. We'll discuss this help file. Try

rand(3)
rand(3,2)
rand(3,2,2)

These commands allow you to create arrays of random numbers without loops. Compare

A = rand(10000,1);

With

for i = 1:10000
B(i,1) = rand;
end

Note however there is really two issues with the second method: first, we could have avoided the loop, and second, the we successively built a larger and larger array without first initializing it. Now that B is defined, repeating the second set of command will go much faster, but still not as fast as avoiding the loop.

Note semicolons suppress a command's insistence on displaying the result on the screen. Try it without the semicolon:

A = rand(10000,1)

Finally, if you are running a model that uses random numbers it is useful to be able to reproduce simulations. This can be done because computer random numbers are really pseudo-random numbers. The number returned the the generator is completely determined by the "seed" used to initialize the generator and the number of times that the generator has been used since the last initialization. Try:

rand('state',0);
rand
rand
rand
rand('state',1);
rand
rand
rand

The command "clock" can be used to set the seed in "random" way. Use "up arrow" to repeat a command.

clock
clock
clock
sum(100*clock)
sum(100*clock)
sum(100*clock)
rand('state',sum(100*clock));
rand
rand
rand('state',sum(100*clock));
rand
rand

If you want a "random" seed and still be able to reproduce your simulations use:

s = sum(100*clock);
rand('state',s);
rand
rand
rand('state',s)
rand
rand

Probability

We want to use the command rand to roll a die. The command rand returns a number between 0 and 1. So 6*rand returns a number between 0 and 6. The command ceil rounds numbers up.

 ceil(6*rand)

Repeat this command several times. But now let's create a new command. Open the Editor and type

function d = die
% DIE Roll a die
% d = die returns a random integer between 1 and 6, inclusive

d = ceil(6*rand);

Save the file as die.m then type

die
die
die

For the next step it is useful to get more than one throw of the die. Save as dice.m then change

function d = dice(n)
% DICE Rolls dice
% d = dice(n) returns n random integers between 1 and 6, inclusive

d = ceil(6*rand(n,1));

Now some statistics:

mean(dice(100))
mean(dice(1000))
mean(dice(10000))
mean(dice(100000))

• Homework Problem 1: Make a guess as to whether the mean converges as the number of dice increases. Assuming you do not know the answer, what numerical experiments would increase your confidence in this regard?

Now try the same thing for

median
mode
std
min
max

The command std is the standard deviation.

• Homework Problem 2: Make guesses as to whether each of these statistics converges as the number of dice increases.