Random Number Generation in Mathematica 6

Andrzej Kozlowski

Tokyo Denki University, Chiba, Japan and Warsaw University, Warsaw, Poland

Abstract

We describe the main aspects of the new mechanism of random number generation in Mathematica 6.  After discussing some problems in the design of the random number generation in earlier versions of Mathematica, we show how they were overcome in  versions 6, and explain why the flexibility and extensibility of the new design should prevent such problems arising again in the future. We illustrate the new flexibility with an example from mathematical finance, by computing the value of a "european call option" by Monte-Carlo simulation using different built-in generators with variable parameters, and the new "extensibility" by  defining a  random number generator for a a finite non-uniform distribution.

(Pseudo-) Random Number Generators

Mathematica 6 comes with a completely new set of tools for generating and using pseudo-random sequences of numbers. These may not be the most celebrated or spectacular of the new features but their importance is difficult to  overstate.  Pseudo-random numbers ("pseudo" -to distinguish them from genuine "random numbers", which can, at best, be generated only by certain physical processes), to which from now on we shall refer to as simply as "random numbers"  or RNs, have played  a huge role in many areas of science, mathematics and computer science ever since the invention by Enrico Fermi, Stanislaw Ulam, Nicholas Metropolis and John von Neumann  of what is now known as the "Monte-Carlo" method of integration. In order to apply this technique John Von Neumann invented  the first RNG (Random Number Generator) - the "mid square method". By today's  standards it was of poor quality  but it served its purpose and gave rise to a whole new area of study.

As von Neumann pointed out in an often quoted remark: "there are no such things as random numbers - there are only methods to generate random numbers".  Each such method generates random numbers according to some statistical distribution. Usually one particular distribution, the Uniform Distribution,  is considered as basic; once we have a RNG which approximates the Uniform Distribution we can define RNGs which approximate any other distribution by means o several well known techniques. 

Since a RNG approximates a statistical distribution we can always ask how close is the approximation to the  statistical distribution it is supposed to approximate. A good uniform RNG is one that produces sequences without any "apparent patterns", or more exactly: the patterns that do occur do so not significantly more  or less often than one would expect in the case of genuine uniform distribution. There are a variety of well known statistical tests that are used to assess the quality of RNGs. Many are discussed in Knuth's book [Kn].  One of the most popular  test suites is  George Marsaglia's "Diehard" battery of tests (the pun is deliberate, see [Ma]).  It should always be remembered that even the best RNG  cannot be expected to pass all possible tests.  As Donald Knuth wrote in 1969: "Every random number generator will fail in at least one application". To what extend this is true about some of the modern RNGs may be argued about (in particular see the discussion of RNGs in [Wo2]) but it is clear that rather than aiming for an ideal RNG it is more important to have a number of RNGs that perform well under well understood conditions, and which offer us enough flexibility so that we can choose the best one for the task at hand. Mathematica 6 represents a big advance in this direction over the previous versions.

RNGs in Mathematica  prior to version 6.

Random Number generation has from the beginning played a somewhat special role in Mathematica.  The very first version of Mathematica used two RNGs. One of them  was, at that time, a  new and remarkable non-linear  RNG - which was the first important practical application of Stephen Wolfram's work on Cellular Automata (nowadays usually referred to as NKS [Wo2]). In [Wo1] it was first observed that a simple 1-dimensional CA with neighbourhood rule of size 3, CA30 (so named  because of the decimal value of its eight-row truth table) makes it possible to construct a non-linear RNG (for the uniform distribution) that performed better on tests then  linear feedback shift-register RNGs popular at the time. In fact, the Wolfram CA30 RNG was able to pass all the tests in the  Marsaglia Diehard test suite. The other RNG also used by Mathematica was a well known and, at that time, highly regarded linear congruential RNG  of Marsaglia and Zaman,  which used a "lagged" version of the "subtract with borrow" (SWB) algorithm.  Unfortunately, in versions of Mathematica prior to 6, these generators were not available as two separate, user selectable, methods of generating RNs  but rather Mathematica used one  or the other RNG or a combination of the two depending on what sort of random numbers were required. For example, to generate small integers the CA30 generator was used, to generate machine reals the SWB was used and to generate large integers both were used together. (The precise description of the approach was given by Daniel Lichtblau in [Li] ). 
   As computers became faster more RNs were needed for sophisticated applications and RNGs that were once regarded as satisfactory were found wanting. This is what happened to the Marsaglia Zaman SWB RNG. In 1995 this RNG was found to fail Marsaglia's gap test (one of Marsaglia's Diehard suite), when the gap exceeds the "long lag".  Unfortunately the lags in the Mathematica implementation of the SWB algorithm were rather small and from around the year 2000 posts began to appear on the MathGroup from people who had discovered something "fishy" in the results of simulations they run which depended on the Mathematica Random[] function. In response to this Daniel Lichtblau posted a partial fix on the MathGroup [Li]. The fix used the fact that all known problems were traced to the SWB RNG while the CA30 RNG, which had passed all the Diehard tests with flying colors. Hence by replacing the use of the SWB RNG whenever possible by the CA30 one could bypass the problem. 
  The fix, however, was only a stop-gap measure and a far from satisfactory one. Firstly, it resulted in a very substantial loss in performance. Secondly, it fixed only one (albeit the most common) problem; the one that occurred when machine random reals where generated. Another  problem, which occurred with large random integers, remained. Thirdly, there was no guarantee that, as computer speeds progressed further, new kind of problems affecting perhaps even the CA30 would not be found or at least new RNGs found that perform better than the built-in ones. In fact, there are now known CA based RNGs  that perform better than the  CA30. Thus, only by allowing the user a choice of an RNG to fit his purpose and perhaps the ability to add  new RNGs could one solve the existing problems with enough flexibility to cope with future ones. Such a solution obviously required a fundamentally new approach, which is why it took longer to fix the problems  than many expected it to take. Thanks to that, however, Mathematica 6 now offers practically everything that one could wish for in this respect.

RNGs in Mathematica 6.

Adopting a radically new approach meant, however, that it was impossible to preserve at the same time the old function names and compatibility with earlier versions of Mathematica. Thus it was decided to keep the old function Random[] with its problematic legacy behaviour and add a new collection of functions with names identifying the kind of random number that they generated. Thus,  the old single "form" Random[type, range], where type could be Integer, Real or Complex,  has been superseded by three new functions RandomReal[range], RandomInteger[range], RandomComplex[range] . Furthermore, these new functions can now be used to generate random arrays very efficiently :

RandomInteger[{1,6},{10^3,10^3}];//Timing

which is much faster than

Array[RandomInteger[{1,6}]&,{10^3,10^3}];//Timing

This use of RandomInteger, RandomReal, etc, replaces the function RandomArray from the now obsolete Statistics package. Just as RandomArray did, the new RNGs work with built in distributions other than the Uniform distribution. Here we use this to simulate 20 paths of a random walk approximation to the standard Wiener Process (Brownian motion).

With[{m=500,n=20},ListPlot[Accumulate/@Map[Prepend[#,0]&,RandomReal[NormalDistribution[0,N[Sqrt[1/m]]],{n,m}]],Joined->True,DataRange->{0,1}]]

Now can now demonstrate one of the problems with the original Random function (this one concerns generation of large random integers) and verify that it has been fixed in Version 6. Here is one manifestation of the problem as posted in a message to the MathGroup in April 2004 by Veit Elser (with a slight change in the code)

ls=ConstantArray[0,25];ls[[1]]=1;ls[[21]]=-1;ls[[25]]=1;
r1=Table[Random[Integer,{1,2^62}],{10^5}];
r2=RandomInteger[{1,2^62},{10^5}];

Length[Union[ListCorrelate[ls,r1]]]

Here we see that the problem is fixed in version 6:

Length[Union[ListCorrelate[ls,r2]]]

We can see a striking difference between the two answers, of which the first is clearly far too small. 
Mathematica 6.0 takes also a large step towards avoiding such pitfalls in the future by allowing the user a choice of several built-in RNGs and even the ability to define ones own.  Mathematica now has six built-in RNGs (one is only available on X86 systems). Exactly half of  these  are based on cellular automata; a somewhat unusual situation compared to other programs, but not surprising in view of Mathematica's origins and history. The default generator is no longer  the "Wolfram favorite"  CA30  but a 5-neighbour "ExtendedCA" RNG, which has been shown to produce  higher degree of randomness faster than the CA30 generator. 30CA is still available under the name "Rule30CA", as is another CA based RNG- "Rule50025CA". There are also several popular non-CA based generators: the "MersenneTwister" - a well known "number theoretic" generator due to M. Matsumoto and T. Nishimura,  a user customizable family of Linear Congruential RNGs, several RNG-s from Intel's MKL library and, for the sake of compatibility, the "Legacy" RNGs used in all earlier versions of Mathematica.

An example from Mathematical Finance

We will illustrate the use of these RNGs by using an example from modern stochastic finance: computing the value of the so called "European Call Option" in the "Black-Scholes model". In this case we have an exact analytic formula for the value of the option given by:

BSEuropeanCall[s_,k_,σ_,r_,d_,t_]:=(1/2) (E^(-d t)) s (erf((((((Sqrt[t]) σ)/2)+(((r-d) t+log((s/k)))/((Sqrt[t]) σ)))/(Sqrt[2])))+1)-(1/2) (E^(-r t)) k (erf((((((r-d) t+log((s/k)))/((Sqrt[t]) σ))-(((Sqrt[t]) σ)/2))/(Sqrt[2])))+1)

where s denotes the current value of the stock, k - the strike price, σ the volatility of the stock, r -the risk-free rate of interest, d - the dividend paid by the stock and t the time to expiry of the option. An approximate value of the option can be computed by means of Monte-Carlo simulation. This can be used as a basis of another "test" of a suitability of a random number generator. Better RNGs should produce more accurate answers (actually things are not quite as simple - see below). To see how this works we will compute the value of an European Call option using four RNGs - the default "ExtendedCA" RNG, a version of the "ExtendedCA", with optional parameter values chosen so as to reduce the quality of the randomness and improve speed,  a once popular linear congruential generator known as RANDU and finally a linear congruential generator with parameters chosen so as to make it quite unsuitable for our  purpose.

MCEuropeanCall[s_,k_,σ_,r_,d_,t_,n_]:=(BlockRandom[Module[{B=RandomReal[NormalDistribution[0,1],n]},(E^(-r t)) Mean[Max[0,#-k]&/@(s (E^((-((σ^2)/2)-d+r) t+σ (Sqrt[t]) B)))]]])

MCEuropeanCall1[s_,k_,σ_,r_,d_,t_,n_]:=(SeedRandom[Method->{"ExtendedCA","Skip"->1}];MCEuropeanCall[s,k,σ,r,d,t,n])

MCEuropeanCall2[s_,k_,σ_,r_,d_,t_,n_]:=(SeedRandom[Method->{"Congruential","Multiplier"->2^16+3,"Increment"->0,"Modulus"->2^31}];MCEuropeanCall[s,k,σ,r,d,t,n])

MCEuropeanCall3[s_,k_,σ_,r_,d_,t_,n_]:=(SeedRandom[Method->{"Congruential","Multiplier"->11,"Increment"->0,"Modulus"->63}];MCEuropeanCall[s,k,σ,r,d,t,n])

In the above code we see one modified and one new Mathematica 6 function related to RNGs. The modified version of SeedRandom now admits a Method option, which can be used to set the RNG to be used in subsequent calls to Random functions. The new function BlockRandom localizes all the RNGs, so that subsequent calls to Random functions are not affected by any Random calls enclosed within a block.

We first compute the exact value of an European call option using the Black-Scholes formula:

BSEuropeanCall[23,23,0.3,0,0.1,.3]

Now the approximate value using a Monte-Carlo simulation with the default RNG, using a relatively small number of simulations:

MCEuropeanCall[23,23,0.3,0,0.1,.3,10^4]//Timing

Because of our use of BlockRandom repeating the calculation will return exactly the same value:

MCEuropeanCall[23,23,0.3,0,0.1,.3,10^4]//Timing

This will no longer be the case once we perform a Random calculation outside BlockRandom:

RandomInteger[];

MCEuropeanCall[23,23,0.3,0,0.1,.3,10^4]//Timing

We now use a larger sample and compare our RNGs. First, the default "ExtendedCA". This generator evolves a state vector consisting of 0 and 1, so that an element at a certain position in the new state vector depends on five neighbors of the element in the previous state. For higher randomness a certain number of elements (by default 4) are skipped.

MCEuropeanCall[23,23,0.3,0,0.1,.3,10^7]//Timing

The answer obtained is  accurate enough for almost all practical needs. The second generator is the same but skips only 1 cell. This should result in somewhat higher speed and lower quality of randomness.

MCEuropeanCall1[23,23,0.3,0,0.1,.3,10^7]//Timing

As expected, the second answer is slightly less accurate and the computation time slightly shorter but the difference is not statistically significant. 
The   RANDU generator, which was implemented on early IBM mainframe systems,  nowadays has a bad reputation and fails several randomness tests, but in this particular case  it seems to only a little less accurate  than the high-quality default generator.

MCEuropeanCall2[23,23,0.3,0,0.1,.3,10^7]//Timing

Finally we try our bad RNG. It is again a linear congruential, but with parameters chosen so that RNs have a very short period. This makes the RNG completely unsuitable for this kind of computation.

MCEuropeanCall3[23,23,0.3,0,0.1,.3,10^7]//Timing

At this point is important not to forget to re-seed the RNG to the default, otherwise the bad congruential generator will continue to be used in subsequent calculations!

SeedRandom[Method->"ExtendedCA"]

Monte-Carlo simulations of this kind do not actually test the quality of "randomness" of the RNG. In fact, sequences of so called "low-discrepancy numbers" (or "quasi-random numbers"),  which are completely "regular",  tend to produce more accurate answers in multi-dimensional Monte-Carlo simulations  than pseudo-random numbers , which attempt to imitate true randomness.  There  are a number of known sequences of  low-discrepancy numbers defined and studied by number-theoretic ( rather than statistical) methods. Two of these, "the Sobol numbers" and  "the Niederreiter numbers" are available in Mathematica 6 in the guise of RNGs (with the option Method ->"MKL") as part of the Intel MKL library. However, as this option is only available on Intel-based systems I will not discuss it here.

Adding user defined RNGs and Distributions

Mathematica 6 has a mechanism for adding user defined RNGs (described in the "Random Number Generation" tutorial - one of the many excellent tutorials included as port of the new online Documentation). This mechanism is, however, somewhat complicated and given the large selection of built-in RNG's it seems unlikely that it will be widely used. There is, however, another extensibility in Mathematica 6 Random Number generation that  surely will be welcomed by many. This is the ability to make RNG with user defined distributions from the built-in (or user defined) uniform ones. We give just one example that has been in the past requested by posters to the MathGroup. This defines a Finite Discrete distribution (FD) which corresponds to a random variable taking a list of values  vals with certain, unequal, probabilities probs. (The case of equal probabilities is handled by the default Uniform Distribution). Here is the definition of the distribution FD:

FD/:Random`DistributionVector[FD[probs_,vals_],m_,___]:=RandomChoice[probs->vals,m]

The function FD takes two parameters. The first is a list of probabilities (or weights) , the second a finite number of values, which have to be numbers (in Mathematica's sense, that is NumberQ must return True, so symbols like π or e are not allowed).

Suppose, for example, we toss an unfair coin which comes up Heads with probability (1/3) and Tails with probability (2/3). Assigning the number 1 to Heads and 0 to Tails, we can simulate this as follows:

RandomInteger[FD[{1/3,2/3},{1,0}],{6}]

Note that, in this case, we could just as well have used RandomReal or RandomComplex:

RandomReal[FD[{1/3,2/3},{1,0}],{6}]

We generate a list of 100000 ones or zeros from this distribution:

ls=RandomInteger[FD[{1/3,2/3},{0,1}],10^6];

We check that we got correct frequencies:

N[Tally[ls][[All,2]]/10^6]

Of course, if we use a "bad"  RNG the frequencies will be far from correct:

ls1=BlockRandom[SeedRandom[Method->{"Congruential","Multiplier"->11,"Increment"->0,"Modulus"->63}];RandomInteger[FD[{1/3,2/3},{0,1}],10^6]];

N[Tally[ls1][[All,2]]/10^6]

SeedRandom[Method->"ExtendedCA"]

Conclusion

Mathematica 6 fixes long standing problems as well as bringing many  powerful new features to random number generation. Mathematica users have had to wait for his a long time, but now, in this very important area, with numerous applications in mathematics, science, finance etc., they have a tool that is better than most and second to none.

References

[Kn] Knuth D. E. (1981) The Art of Computer Programming Vol 2:Seminumerical Algorithms. Addison-Wesley
[Li] Lichtblau D. <http://forums.wolfram.com/mathgroup/archive/2000/May/msg00088.html>
[Ma] Marsaglia G. Random Number Generators, Jour. of Modern Statistical Methods, May 2003, vol. 2, No. 1, 2-13. 
[Wo1] Wolfram S. Random Sequence Generation by Cellular Automata, Advances in Applied Mathematics, vol. 7, p 123-169. June 1986
[Wo2] Wolfram S. A New Kind of Science, Wolfram Media, Inc. 2002.