Technical Note on MATLAB^®’s Random Number Generator

Any Monte Carlo simulation code by definition uses numbers that come from some probability distribution. Some codes use a lot of random numbers; a few of the codes in this book use millions of them. MATLAB^®’s random number generator provides a convenient software source for these numbers. MathWorks, the creator of MATLAB^®, has gone through a series of designs for its generators, and the latest version uses a highly sophisticated combination of shift register/bit manipulation processes that, unlike earlier generators, require no multiplication or division operations. The latest MATLAB^® generator, then, is very fast. At least as important, however, is that such a design results in a stupendously long period, where the period is the number of numbers a generator can produce before it starts repeating.

The fact that a software generator is certain to eventually repeat is due to the fact that the generator is what computer scientists call a finite-state machine. If a digital entity is built from n elementary two-state elements (that is, from n binary elements that are each either 0 or 1), then there can be at most 2ⁿ distinct combinations of 0s and 1s (each such combination is called a state of the entity). To use repeating random numbers in a simulation means the numbers aren’t random anymore, a situation that invalidates the underlying, fundamental assertion that the simulation somehow emulates the random process being simulated.

So, having a large period is an essential feature of a good random number generator. MATLAB^®’s random number generator period is so huge that, even if producing numbers at the impressive rate of a billion per second beginning with the birth of the universe in the moment of the Big Bang, we would today have observed only an infinitesimally tiny fraction of the period. No simulation code you or anybody else will ever write will even begin to exhaust the ability of MATLAB^® to provide the code with new random numbers.

The MATLAB^® command rand produces a number from a distribution uniform from 0 to 1. If you want a number from a distribution uniform from a to b, where b > a, then a + (b − a)*rand does the job. MATLAB^® also provides the command randn that produces a number from a normal (bell-shaped) distribution with zero mean (m = 0) and unit standard deviation (σ = 1), that is, from a distribution with density function

If you want a number from a normal distribution with mean m and standard deviation σ, then m + σ*randn does the job. If you want a number from a distribution other than a uniform one or a normal one, then you’ll have to write some additional code of your own: see my book Digital Dice (Princeton 2008), pp. 252–254, for an example on one way to do that for the case of an exponential distribution.

When you start a new MATLAB^® session the generator is automatically put into a predetermined initial state; then the generator transitions from its present state to the next state as each new number is produced. The details of the state transitions are determined by the specific design of the generator, the details of which you don’t need to know. This means that if you run one of the codes in this book multiple times in the same session, the first run uses one set of numbers from the generator, the second run uses a different set, and so on. That means that while the results produced by the code from one run to the next will be close (hopefully!), they won’t be exactly the same.

When you are first writing a simulation code, getting different answers every time you run the code can complicate the process of debugging; is the code giving different results because its logic is in error (this is bad) or because it’s getting different numbers from the generator (this is okay)? So, at the start of writing a code, it’s often desirable to be able to reinitialize the generator to the same state, every time the code is run during the same MATLAB^® session. This can be done by writing, somewhere in the code before the first use of the generator, the command rand(‘state’,0).

Once the debugging process is completed, however, you’ll want every new run of the code from session to session to use random numbers it has never seen before. After all, why bother running a simulation code with the same generator numbers as before—you’ll just get the same answer as before! So, how do you achieve that when I’ve just told you that MATLAB^® automatically forces the generator into the same initial state at the start of every new session? There is a simple way to do this using the clock command.

MATLAB^®, residing in your computer, has access to your computer’s knowledge of the time and uses that access to create and continually update the six-element vector clock with the format [year, month, day, hour, minute, second]. For example, I am typing this at one minute and a few seconds past 1 p.m. on the afternoon of February 3, 2012, and “right now” the clock vector is [2012, 2, 3, 13, 1, 9.5]. Since the march of time is inexorably unidirectional, we can use this vector to generate a new, unique (almost always) initial state every time we run a simulation code. This is done by writing, somewhere in your code before the first use of the generator, the command rand(‘state’,100*sum(clock)). The command sum adds the six elements of the clock vector, and the multiplication by 100 results in an integer (for the above vector the integer is 204050) to set the initial state of the generator. A different time, a different clock vector, a different (usually) integer, and so a different initial state for the generator and a different stream of numbers each time you run your code.

You’ll see in the codes given in this book that I haven’t bothered to do anything about the initial state of the random number generator. The results I give from running the codes in this book are typical—but were produced by whatever numbers came out of the generator when I happened to run the simulations. If you run these codes on your machine (if you have MATLAB^®) you’ll almost certainly get results that are near mine, but not exactly the same. It is, in fact, a singular feature of Monte Carlo computer codes to get something new almost every time!

When a code runs hundreds of thousands (or even millions) of simulations and then computes the average of all those individual simulations, the results are generally ones having a lot of digits. Do all of those digits really have meaning? Almost certainly not, but nonetheless I have reported all the digits the codes in this book have generated, even though probably only the first three (or maybe four, if we’re lucky) are correct. You should be skeptical of digits beyond the fourth decimal place produced by the codes in this book. If you need more digits, however, theoretical analyses of the underlying mathematics of the Monte Carlo method show that the statistical error decreases as , where N is the number of simulations.

For example, increasing N from 10,000 to 1,000,000 (a factor of 100) should reduce the error in the code’s estimates of whatever parameters are being studied by a factor of about . A specific illustration of this behavior is given in my book Digital Dice, pp. 11–15, along with an example of a quite different approach (called variance reduction) on pp. 223–228.