Generating Uniform Random Numbers

Q: How to generate pure random string?
A: Put a fresh student in front of vi editor and ask him[/her] to quit.
     -- Emiliano Lourbet (@taitooz)

To perform computer-based simulation we need to be able to generate random numbers. Generating random numbers following a uniform distribution are the easiest to generate and are what comes out of the standard programming language "give me a random number" function. Here's a sample Python session:

Uniform random variables are super important because they are the basis from which we generate other random variables, such as binomial, normal, exponential etc.

Discussion

We could generate real random numbers by accessing, for example, noise on the ethernet network device but that noise might not be uniformly distributed. We typically generate pseudorandom numbers that aren't really random but look like they are. From Ross' Simulation book, we see a very easy recursive mechanism (recurrence relation) that generates values in $[0,1)$ using the modulo (remainder) operation:

$x_{i+1} = a x_i$ modulo $m$

That is recursive (or iterative and not closed form) because $x_i$ is a function of a prior value:

$x_1 = ax_0$ modulo $m$
$x_2 = ax_1$ modulo $m$
$x_3 = ax_2$ modulo $m$
$x_4 = ax_3$ modulo $m$
$...$

Ross indicates that the $x_i$ values are in [0,m-1] but setting any $x_i=0$ renders all subsequent $x_i=0$, so we should avoid that. Practically speaking, then, this method returns values in (0,1).

To get random numbers in [0,1) from $x_i$, we use $x_i / m$.

We must pick a value for $a$ and $m$ that make $x_i$ seem random. Ross suggests choosing a large prime number for $m$ that fits in our integer word size, e.g., $m = 2^{31} - 1$, and $a = 7^5 = 16807$.

Initially we set a value for $x_0$, called the random seed (it is not the first random number). Every seed leads to a different sequence of pseudorandom numbers. (In Python, you can set the seed of the standard library by using random.seed([x]).)

Python implementation

Our goal is to take that simple recursive formula and use it to generate uniform random numbers. Function runif01() returns a new random value for every call. Use $m = 2^{31} - 1$, $a = 7^5 = 16807$, and an initial seed of $x_0 = 666$.

Notice that x_i is in the global space not the runif() space.

Let's try it out:

Exercise

Define a new function, runif(a,b) that generates a random number in [a,b) instead of [0,1). Hint: We need to scale and shift a random uniform value in [0,1). Note: You can't use random.random() or any other built-in random number generators for this lab.

def runif(a,b):
    "Return a random value in U(a,b)"
    ...

Exercise

Define a new function, setseed(x), that updates the seed global variable.

def setseed(s):
    "Update the seed global variable but ensure seed > 0"
    ...

This test sequence:

setseed(501)
print runif01()
print runif01()
print runif01()

should generate:

0.00392101099897
0.900431859726
0.558266419712

Random variable density function estimate

Jumping ahead a bit, we can use the histogram plotting example from Manipulating and Visualizing Data as a crude form of density estimation to verify that the distribution of random values is approximately uniform:

In the case of generating pseudorandom numbers, we are interested in the sequence of values generated by the recurrence relation. Now we'll turn our attention to iterative methods that loop until the recurrence relation value converges.