B Random Number Generators
The purpose of this appendix is to give a brief overview of random number generation. Random number generators are a vital tool for cryptography, statistical sampling, video gaming, gambling and lotteries. Our application of interest is the production of artificial data for simulation experiments. Both Gentle (2003) and Eubank and Kupresanin (2012) can be consulted for more thorough treatments of that topic, and Ismay (2013) for methods that test the independence of parallel streams.
The discussion in this and the subsequent appendix requires a little background. First, there is the concept of a bit. Specifically, every nonnegative integer, \(x\), has a binary representation as a sum of powers of 2. This looks like \[\begin{equation} x = \sum_{k = 0}^mb_k \cdot 2^k \label{eq:binrep} \end{equation}\] for some non-negative integer \(m\) and integers \(b_0, \ldots, b_m\) having values of either 0 or 1. The \(b_k\) are called bits and knowing or manipulating the bits is equivalent to knowing or working with the number itself. Computers store and work with these types of bit-wise representations. Typical values for \(m\) in (??) are 31 and 63.
We can illustrate these ideas using Python. In that regard, a Python function that returns the bitwise representation of a nonnegative integer is bin.
One may check, for example, that binary forms for the integers 3, 15 and 203 are
## '0b11'
## '0b1111'
## '0b11001011'Note that bin returns a text string prefaced by 0b to indicate that it is in binary form. To be clear on what these numbers mean observe that
\[\begin{align*}
3 &= 2^1 + 2^0 \\
15 &= 2^3 + 2^2 + 2^1 + 2^0 \\
203 &= 2^7 + 2^6 + 2^3 + 2^1 + 2^0
\end{align*}\]
To manipulate integers at the bit-wise level we can use various operators that include AND (&), OR (|), the exclusive OR, XOR, denoted by ^ in Python, and the right and left shift functions >> and <<. The bit-wise AND compares two numbers bit by bit, returning a new number wherein bits are set to 1 only if both corresponding bits in the input numbers are 1: e.g.,
## '0b111'
## '0b11001001'
## '0b1'The & operator appears in masking where one wants to enforce a specific bit size for a binary value. For example, X & 255 returns an 8-bit output and is equivalent to X % 256
.
The bit-wise OR operator returns a 1 if either bit is 1 and 0 otherwise while the exclusive OR compares two numbers to produce a new number wherein if the bits being compared are different the result is 1 and is 0 otherwise. These produce results like
## '0b1111'
## '0b100001'
## '0b101111'
## '0b101110'The right shift operator moves bits to the right by a specified value. Bits on the far right fall off (are discarded). For example,
## '0b100101101'
## '0b100101'The left shift moves bits to the left a specified number of slots as in
## '0b100100111'
## '0b100100111000'Another bit-wise operation that appears in random number generation is the XOR shift which takes the form
for integers a and b. The steps involved in the computation are illustrated by
## '0b11110000'
## '0b111100'
## '0b11001100'
## '0b11001'
## '0b11001'As demonstrated in this simple example the xorshift operation mixes bits in a way that allows information from high-order bits to flow down into lower-order slots.
The xorshift or shift-register random number generator suggested by Marsaglia (2003) takes the form
for \(a, b, c\) positive integers. Here ^= performs a bitwise XOR operation on the right hand side and then assigns the result back to that same variable. Marsaglia’s original xorshift generator is very fast. But, the plain three-shift versions do not pass the stringent BigCrush tests in the TestU01 suite of performance tests for random number generators in L’Ecuyer and Simard (2007) without additional scrambling. Xorshift operations are used as building blocks in more sophisticated generators, including the PCGs discussed in the next Appendix.
B.1 Congruential Generators
The once pervasive linear congruential generator (LCG) has the form \[\begin{equation} state_{t+1} = (a \ast state_t + c) \bmod m \tag{B.1} \end{equation}\] for integers \(m > 0, 0 < a < m, 0 \leq c\) and \(0 < state_0 < m\) the seed that starts the recursion. The period of a generator is the number of values in the sequence (B.1) until the recursion starts to repeat. In this instance it is at most \(m\).
The Hull-Dobell theorem, e.g., Knuth (1997), has the consequence that when \(c \neq 0\) the LCG will have full period if and only if (i) \(m\) and \(c\) are coprime, (ii) \(a - 1\) is divisible by all prime factors of \(m\), and (iii) \(a - 1\) is divisible by 4 when \(m\) is divisible by 4.
A case of interest for the next appendix uses \(m = 2^{64}\) and the multiplier suggested by O’Neill (2014): \(a = 47026247687942121848144207491837523525\). Since the only prime factor of \(2^{64}\) is 2 and \(a - 1\) is divisible by 4 all that is needed for a full period is that \(c\) be odd.
Congruential generators are fast and easy to implement. Everything looks fine at first glance, but the problem is in the bits. The evolution of each bit is not equally random and low order bits are weak. In particular, if the generator has full period the lowest bit will simply alternate between 0 and 1. This is certainly not random. The next-lowest bits are also weak, although less obviously so. In general, lower bits have shorter periods and simpler structure, while higher-order bits tend to behave more randomly.
The entire LCG is linear modulo \(m\). This creates a strong algebraic structure. For example, in higher dimensions, successive triples tend to fall on a relatively small number of planes. This is the well-known lattice structure of congruential generators that, in particular, causes them to fail the BigCrush tests in TestU01.
The inherent flaw of starting the same generator with different seeds to produce parallel streams can be easily seen here. Let \(x_0\) and \(y_0\) be two different seeds that are used to initialize a congruential generator \(x_t = ax_{t - 1} + c \mod m\). Then, \(y_t - x_t = a (y_{t-1} - x_{t-1}) = a^t(y_0 - x_0)\) which means that \(y_t = x_t + a^t (y_0 - x_0)\) and the two streams are linear transformations of each other. Perhaps more telling is that if the generator has full period there will be a \(k\) such that \(y_t = x_{t + k}\) for all \(t\). To see this note that the statement is clearly true for each fixed \(t\). So, in particular, for \(t = 0\) there is a \(k_0\) with \(y_0 = x_{k_0}\). Thus, \(y_1 = ay_0 + c = ax_{k_0} + c = x_{k_0 + 1}\). Proceeding by induction we find that \(y_t = x_{t + k_0}\) for all \(t\). The \(x\) and \(y\) streams are therefore just shifted versions of each other.
A related question concerns two congruential state sequences generated with the same multiplier but different increments. Although the sequences may visit the same individual states, they cannot be shifted copies of one another. Indeed, if one sequence agreed with a shifted version of the other at every index, comparison of their recursions would imply that their increments were equal modulo the modulus. This property allows generators such as PCG in the next appendix to define distinct parallel streams by assigning different admissible increments to different processes.
B.2 Other Generators
Another generator that is still widely used in statistical computing is the Mersenne Twister developed by Matsumoto and Nishimura (1998). It is a generalized feedback shift register generator with the addition of a twist transformation. It has a very long period of \(2^{19937}\) but is known to fail some of the TestU01 evaluations.
Combined multiple recursive generators such as the L’Ecuyer et al. (2007) MRG32k3a are especially important for simulation. This generator passes the TestU01 tests, has a period of roughly \(2^{191}\) and is designed to support independent streams and substreams for use in parallel number generation. It is readily available in R through the parallel package but is not part of Python’s standard random-number infrastructure.
A different approach is represented by counter-based generators such as Threefry and Philox from Salmon et al. (2011). These generators do not advance by a recurrence in the usual sense. Instead, they apply a deterministic transformation to a counter. A commonly used version of Threefry has a period of \(2^{128}\) and produces \(2^{64}\) independent streams. Philox has period \(2^{256}\) and produces at least \(2^{64}\) non-overlapping number streams. These generators are available in R through the dqrng package. In Python Philox can be found in NumPy.random with Threefry in the randomgen module.
B.3 The tale of RANDU
The infamous RANDU is an LCG used on IBM mainframes in the 1960s and 1970s that became the default random number generator in some widely used software. It provides a notorious example of a poorly chosen linear congruential generator whose outputs, when viewed in three dimensions, lie on a small number of planes rather than filling the space uniformly. Its widespread use led to unreliable simulation results and motivated the development of more rigorous tests for random number generators such as TestU01.
RANDU has the form \[ state_{t+1} = (65539\, state_t) \bmod 2^{31}. \] The generator has period \(2^{29}\) but has fundamental flaws. The generated sequence satisfies a linear relation among triples:
\[ state_{t+2} = 6\cdot state_{t+1} - 9\cdot state_t \pmod{2^{31}} \]
This has the consequence that when you plot the triples \((state_t, state_{t+1}, state_{t+2})\) instead of filling the cube, they fall on only 15 parallel planes. This was discovered by Marsaglia (1968).
Figure B.1 makes RANDU’s failure visible. The trick is that a single viewpoint is not convincing: from most angles the points look like a fuzzy cube. The structure only jumps out when the cube is rotated. The figure shows the two rotations that make it clearest: a steep one, and a nearly edge-on one, where the points separate into about fifteen thin parallel sheets with empty bands between. A good generator would fill the cube uniformly from both angles.
Figure B.1: RANDU triples viewed from two rotations. The edge-on view reveals approximately fifteen parallel planes.