Introduction to Random Numbers

Math 245 - Mathematical Models
Goal: Understand how to create computer simulations which incorporate randomness into the model.

6.1 Incorporating Random Numbers

When you are incorporating randomness into your model, you will the command RandomInteger or RandomReal.  The time when RandomInteger is useful is when you are simulating a situation in which there are a fixed number of options, each of which occurs with equal probability.  Otherwise, use RandomReal.  In both cases, the first input gives the range in which to choose the random number; a second input tells how many of such numbers to make, if more than one is desired.  If no input is given at all, the default range is [0,1].

RandomInteger[]
RandomInteger[{-4,4}]
RandomInteger[{1,10},40]
RandomReal[]
RandomReal[0.5]
RandomReal[0.5,10]

A common method of simulating a certain probability of occurrence is to incorporate the RandomReal command.  For example, if we want to simulate an occurrence with probability 0.075, take a random real number, and an If statement to see if the number is less than 0.075.  In the following example, trial is a random number between 0 and 1; and success is 1 in the 7.5% of the time when trial is less than 0.075.

trial=RandomReal[]
success=If[trial≤0.075,1,0]

One might evaluate the above statement multiple times, but why not have Mathematica do it for you.  First, let's create a list where the first entry is the random number and the second entry is whether the occurrence is a success or not.

{trial=RandomReal[],success=If[trial≤0.075,1,0]}

Now simply wrap this statement into a Table to generate as many instances as you desire to be convinced that the If statement works as you expect.

Table[{trial=RandomReal[],If[trial≤0.075,1,0]},{20}]

Notice that a variable is not needed in this Table command :  We simply did Table[expr, {20}] to iterate the expression 20 times.

To see this list in a presentable format, use the TableForm command.  Here the "%" operator references the most recent output.

TableForm[%]

Alternatively, we might not have generated the variable trial; the two steps can be combined into one:

If[RandomReal[]≤0.075,1,0]

and subsequently wrapped in a Table command to simulate 500 occurrences, each with probability 7.5%.

Table[If[RandomReal[]≤0.075,1,0],{500}]

How many 1's were there?  Total the entries in the list.  Then see how close we were to 7.5%:  (Note again the use of the percent operator.)

Total[%]
N[%/500]

The downside to the % operator is that you are never sure if the last command you evaluated is the command directly above.  For instance, re-evaluate this block of code.
(Basically, using % is sloppy coding, but quite useful in many cases.)

Questions from Section 6.1

How might you model a success of 10% using RandomInteger instead of RandomReal?  How about 7.5%?

Practice using the % operator.  
First, define and evaluate some list of even numbers.  
Then divide the list by 2 to divide each entry by 2.  
Then find the sum of these entries in the list.  
Then divide the sum by the number of entries to find the average of the entries.

6.2 Flipping Coins and Rolling Dice

Fair coins

Fair coins and fair dice fit nicely into the "situation in which there are a fixed number of options, each of which occurs with equal probability" mold, so we should use the RandomInteger command.  There are two ways we can approach this---have in our head an interpretation of "1"="Head" and "0"="Tail", or get Mathematica to make this interpretation explicit.  Here is how we would flip one coin and how we would generate 100 flips in each method and determine how many heads appeared.

Method 1, where we use 0’s and 1’s:

RandomInteger[]

Table[RandomInteger[],{100}]

Total[%]

Method 2, where we use “Head” and “Tail”:

If[RandomInteger[]==1,"Head","Tail"]

Table[If[RandomInteger[]==1,"Head","Tail"],{100}]

We can use the Tally command if we want to see how many of each is generated in the second case.

Tally[%]

Fair dice

Rolling fair dice is similar:

RandomInteger[{1,6}]

rolls=Table[RandomInteger[{1,6}],{100}]

Tally[rolls]

Here we might like to have our data sorted; use the Sort command.

Sort[%]

Or we might like a visual representation of this data; use the Histogram command.

Histogram[rolls]

Biased dice

When the dice or coins are biased (faces have unequal probabilities of occurring), we have to fall back on the RandomReal command.  For instance, if we have a four-sided die, where the four sides (A, B, C, and D) come up with probabilities 1/2, 1/4, 1/8, and 1/8, respectively.  We can set up a For loop to generate a random real number between 0 and 1, and increase counters depending on where in the range 0 to 1 the real number falls

aCount=0;bCount=0;cCount=0;dCount=0;
For[i=1,i≤20,i++,
roll=RandomReal[];
If[0≤roll<1/2,aCount++];
If[1/2≤roll<3/4,bCount++];
If[3/4≤roll<7/8,cCount++];
If[7/8≤roll≤1,dCount++];
]
{aCount,bCount,cCount,dCount}

In order to plot this list of data in a Histogram type style, use the command BarChart.

BarChart[%]

A slicker way to do this is to use the Which command, which allows you to nest If statements.  Which should be thought of as a sequence of pairs of tests and outputs.  The command will check to see if the first test is true.  If so, the command will evaluate the first output.  Otherwise, it will check on the second test, and proceed down the line successively.  (Use ? Which to learn more.)

In the following code, I use "1" to represent "A", "2" to represent "B", and so on.  Then I work it into a Table command in order to compile statistics as I did above.

roll=RandomReal[]
Which[roll≤1/2,1,roll≤3/4,2,roll≤7/8,3,roll≤1,4]

rollList=Table[roll=RandomReal[];Which[roll≤1/2,1,roll≤3/4,2,roll≤7/8,3,roll≤1,4],{20}]

Tally[rollList]

Histogram[rollList]

Questions from Section 6.2

Why is the number of Heads generated using method one different from the number of Heads generated using method two?

How would you simulate rolling a six-sided die where the each of the first three sides come up three times as often as each of the last three sides?

6.3 Using Monte Carlo Methods to Approximate Area

Idea behind the method:  Given a region of unknown area A that we surround by a box [(Subscript[x, 1]),(Subscript[x, 2])]×[(Subscript[y, 1]),(Subscript[y, 2])] with area B=((Subscript[x, 2])-(Subscript[x, 1]))((Subscript[y, 2])-(Subscript[y, 1])), if we randomly choose a point from the box, we expect the point will fall in the region with probability A/B.  So if we take N random points in the box and M fall inside the region, then we would expect that A/B≈M/N.  Therefore, we can approximate A≈B M/N.

The RegionPlot Command

To plot a region defined by a set of inequalities, use the RegionPlot command.  The first input tells Mathematica which points are valid points.  For example, in the following example, there are three inequalities, x≥0, y≥0, and x+y<1.5, separated by && because points must satisfy inequality 1 AND inequality 2 AND inequality 3.  The second and third inputs are the ranges for x and y.

RegionPlot[x≥0&&y≥0&&x+y<1.5,{x,-2,2},{y,-2,2}]

Approximating area, a first example

Given a region defined by inequalities, such as 0≤y≤Sin[x], then we can plot the region using the command RegionPlot.

region=RegionPlot[0≤y≤Sin[x],{x,0,Pi},{y,0,1}]

If we want to generate random points in the rectangle [0,π]×[0,1], then we need to generate an x-coordinate and a y-coordinate.
An x-coordinate would be a random real number between 0 and π.
A y-coordinate would be a random real number between 0 and 1.
We need to determine if the y-coordinate is less than or equal to the sine of the x-coordinate.
The following command will test this and output 1 if the point is in the region, and output 0 if not.

xval=RandomReal[{0,1}];
yval=RandomReal[{0,Pi}];
If[xval≤Sin[yval],1,0]

As before, wrap it in a Table command to simulate many points, and sum the data to calculate how many points are in the region.

randPoints=Table[
xval=RandomReal[{0,1}];
yval=RandomReal[{0,Pi}];
If[xval≤Sin[yval],1,0],{500}];
inRegion=Total[randPoints]

The area of the region is approximated by:

N[(Pi)*inRegion/500]

The actual area is:

Integrate[Sin[x],{x,0,Pi}]

Approximating area, a second example

Here is another region, defined by the inequalities |x|+|y|≤1.

region=RegionPlot[Abs[x]+Abs[y]≤1,{x,-1,1},{y,-1,1}]

xval=RandomReal[{-1,1}];
yval=RandomReal[{-1,1}];
If[Abs[xval]+Abs[yval]≤1,1,0]

Total[Table[
xval=RandomReal[{-1,1}];
yval=RandomReal[{-1,1}];
If[Abs[xval]+Abs[yval]≤1,1,0],{500}]]

Let's show a histogram of 1000 of these 500-point trials.

trials=Table[Total[Table[If[Abs[RandomReal[{-1,1}]]+Abs[RandomReal[{-1,1}]]≤1,1,0],{500}]],{1000}];
Histogram[trials]
N[Mean[trials]]

Questions for Section 6.3

What is the true area of the second region?  What is the approximated area of the second region?

What from the second example changes if you instead choose random points from the rectangle [-1,2]x[-2,1]?

How would you plot some random points in the box?  [Hint: Use a Table command to create a list of (x,y) coordinate pairs.]