Chapter 16: Random Numbers and Monte Carlo Simulations

Random numbers are quite useful and interesting to play with. One of the most useful is that of a Monte Carlo simulation which uses random number to solve problems in probability or other fields.

Probability

To fully understand how to use random numbers and how Monte Carlo simulations work, we need some basic probability. When understanding probability there are two types continuous and discrete, which for what we need will come down to using floating point or integers.

Discrete Probability

We’ll first discuss discrete probability. A discrete probability distribution is a finite set (often will be a subset of integers.) For example, let’s use the set

$A={1,2,3,4,5,6}$

and an event, $X$ is subset of these numbers. For example, let $X={1}$, then the probability that this event occurs in the ratio of the number of elements in each set or
$P(X)=\frac{N(X)}{N(A)}$

In this case, $P(X)=1/6$ and that means that the probability that the number 1 comes up is 1/6. Think about this in the case of rolling a die. This says that the probability that 1 comes up is 1/6.

Continuous Probability

Another type of probability distribution is called a continuous distribution and in this case we’ll only consider a type of continuous probability distribution called a uniform distribution.

Let’s consider a set $A={x \; | \; 0 \leq x \leq 1}$ or all real numbers between 0 and 1. Events are still subsets of the set $A$, however the probability that events occur is the fraction of the set.

For example, if $X={x \; | \; 1/3 \leq x \leq 1/2}$, then $P(X)$ is the fraction
$P(X) = \frac{\frac{1}{2}-\frac{1}{3}}{1-0}=\frac{1}{6}$

Pseudo Random Number Generator

The Wikipedia Pseudo-random Number Generator page gives an overview of the subject and a lot of technical details. In short, a truly random number on a computer is very difficult to generate and generally not necessary because a pseudorandom number is sufficient. I also won’t try to explain this in general term since you’d need a significant mathematical background.

A pseudo-random number generator is a function that produces a sequence of numbers that act like random numbers. Let’s examine what this means in terms of the discrete probability with set $A={1,2,3,4,5,6}$.

If a pseudo-random number generator produces a sequence from this set then the sequence should have the following property:

• If the event $X={i}$ for any $i$ between 1 and 6, then $P(X)\approx 1/6$. And by approximately, as the sequence gets larger, the approximation becomes closer to 1/6.

Is this enough? No, the sequence
$${1,2,3,4,5,6,1,2,3,4,5,6,\ldots}$$

satisfies the above property, but I don’t think anyone would consider this random. Another property would be:

• If we know the sequence ${a_1,a_2,a_3,\ldots,a_n}$ then we can’t predict the next number $a_{n+1}$.

This is obviously violated in the sequence above.

A little more technical definition of a sequence of pseudo-random numbers Let $(a_1,a_2,a_3, \ldots,$ be a random sequence. Typically we mean the following properties need to hold:

1. any number in the range 1 to 6 is equally likely to occur.
2. Take N random numbers and let $s_n$ be the number of times the number $n$ occurs. The fraction $s_n/N$ should go to 1/6 in the limit as $n\rightarrow \infty$.
3. Knowing the sequence $s_1, s_2, \ldots, s_k$ does not allow us to predict $s_{k+1}$

If instead we use floating point numbers, there are a few different properties. Assume that the floating point number is in the range $0 \lt a_k \lt 1$. Then the sequence $(a_1,a_2,a_3, \ldots,)$ is random if

Let $s_{[c,d]}=(a_1,a_2,a_3, \ldots)$ be the total of the numbers that satisfy $c \lt a_k \lt d$. Then as the number of random numbers approach infinity $s_{[c,d]}=d-c$

Pseudorandom Numbers in Maple

Maple has some nice features to generate random numbers of many distributions. Look in the RandomTools package for details.

Here’s an example that would be like rolling a die 100 times. First, load the

with(RandomTools):

And the calling

rolls:=Generate(list(integer(range = 1 .. 6), 100))

creates a list of length 100 of random numbers of integers between 1 and 6. (Note: put your cursor back up on the command and hit return to see that you get different sets of numbers each time.)

This list then acts like rolling a single die 100 times. We now want to check how many times the number 1 was rolled. There are a few command in the ListTools package that can help us with these.

with(ListTools)

and now we can type

[SearchAll(1,rolls)]

returns a location where all of the 1s are in the list of 1. We don’t really care where these came up, but typing

nops((#))

where (#) is the line number of the list of 1s above. When I ran this I got 21, but your mileage may vary due to the nature of randomization.

Histograms

A nice visual for rolling the die is that of a Histogram, which is a bar chart of the frequency of rolling each number. The Histogram command of the Statistics package is where we’ll turn:

Histogram(rolls)

gives the following:

and again, your plot will look different or rerunning the lines that generated the rolls will change the plot.

Rolling 2 dice

There are many cases where rolling 2 dice and summing the values is interesting. Although there are many ways to do this, consider

roll2 := [Generate(list(integer(range = 1 .. 6), 100)), Generate(list(integer(range = 1 .. 6), 100))]

returns a list of pairs (as a list) of dice as two lists.

We can then add each pair of dice as

sum2dice := add(i, i in roll2)

which returns a list of the sums of the dice?

Question: Since this returns numbers between 2 and 12, why didn’t we just create a list of length 100 which takes on values between 2 and 12?

What is the probability that we rolled an odd number?

We can do this by the following:

SelectFirst(100, n -> n mod 2 = 1, sum2dice)

where SelectFirst is in the ListTools package. This function returns the first 100 elements of sum2dice that satisfy the function n -> n mod 2 = 1 or odd. If we do nops around this function, we get 57 (again this can vary.)

Let’s Histogram it!!

And if we create a histogram of this with

Histogram(sum2dice)

one example of this is:

Exercise

1. Rerun the commands above however changing the number of rolls of the pairs of dice to 1000. (Note: you probably want to hide the actual results with a colon (:) )

2. What is the probability that you roll a 7. (This is a important roll in craps, a dice game) Compare your result to the known value of the probability if you have background in this (or you can perform an internet search.)

3. Produce a histogram of your results.

Pseudorandom Floating-Point (and other types of ) numbers

Another important use of simulation involves using floating-point numbers as random numbers. The Generate function can do this easily using the following:

Generate(list(float(range = 0 .. 1), 10))

generates a list length 10 of floating points in the range 0 to 1.

Exercise

Try to generate a list of 10 random numbers of type rational in the range from 0 to 5 with a denominator of 24.

Calculating $\pi$ using pseudo random numbers

• Buffon’s Needle Experiment (18th C.)
• Circle in the Square.

Consider the square ${ (x,y)\; | \; 0\leq x\leq 1, 0 \leq y \leq 1}$ and the quarter circle that falls within $x^{2}+y^{2} \leq 1$. The area within the circle is the fraction of the cirle in the square or $\pi/4$.

We can use this to estimate a value of $\pi$.

1. Randomly choose points $x$ and $y$ in the range [0,1].
2. Count all points that fall within the circle.
3. Find the fraction of the points in the circle.
4. Since this fraction should be $\pi/4$, multiply by 4 to estimate $\pi$.

We can do this using Maple as the following:

pts := Matrix(100, 2, Generate(float(range = 0 .. 1), makeproc = true))

will create a 100 by 2 matrix of random numbers between 0 and 1. Then we can calculate the distance-squared of each row of the matrix:

dist2 := [seq(add(x^2, x in pts[n, ..]), n = 1 .. 100)]

and we can count up the number of points that are inside the circle of radius 1 by

inside_circle:=nops(SelectFirst(100, proc (x) options operator, arrow; evalb(x < 1.0) end proc, dist2))

The fraction is this over 100 and then the value of $\pi$ is this times 4.

4*inside_circle/100

And we I ran this I got 2.8

Exercise

Repeat this with 1000, 10_000 or more points.