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.

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.

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.

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}\]

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:

- any number in the range 1 to 6 is equally likely to occur.
- 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 \).
- 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\)

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.

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.

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?

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.)

And if we create a histogram of this with

`Histogram(sum2dice)`

one example of this is:

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 (:) )

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.)

Produce a histogram of your results.

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.

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

- 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\).

- Randomly choose points \(x\) and \(y\) in the range [0,1].
- Count all points that fall within the circle.
- Find the fraction of the points in the circle.
- 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

Repeat this with 1000, 10_000 or more points.