# Bays-Durham Shuffling

Reading "Random number generation and Monte Carlo methods" by James E. Gentle, I encountered an example concerning the Bays-Durham Shuffling algorithm whose result I can't reproduce.

Basically, such an example is written in Random Number Generation - Durham Shuffling of Uniform Deviates. For the full text you can go to Random number generation and Monte Carlo methods.

Well, I'm trying to reproduce Figure 1.7 for the linear congruential generator of $$x_{i} \equiv 3x_{i-1} \text{mod} 31$$. In order to do so, I wrote the following piece of code in Mathematica:

Defines initial variables for the generator:

a = 3; mod = 31; x = 9; n = 0; m = {};


Produce a list with the initial set of numbers (In this case: 27, 19, 26, ... <<27>> ... 1, 3, 9):

While[n < 30, AppendTo[m, Mod[a x, mod]]; x = Mod[a x, mod]; n++];


The algorithm generates the first sequence of numbers in a slightly different way using the input value k, so in order to account for this, I had to use this While for the first sequence:

k = 8; i = 1; y = m[[k + 1]];

While[i == 1, choose = Mod[y, 8] + 1; y = l[[choose]];


Then, for the following sequences:

i = i + 1; m[[choose]] = k; k = k + 2];

While[i < 10000, choose = Mod[y, 8] + 1; y = m[[choose]];

i = i + 1; k++;


Given that the congruence relation is congruence modulo 31, the sequence will repeat itself after 30 digits. To avoid looking for the non-existent element m[[31]], I added the conditional:

  If[k < 30, m[[choose]] = m[[k]], k = 1; m[[choose]] = m[[k]]]];


The book says:

'By continuing in this manner to yield 10,000 deviates and plotting the successive pairs, we get Figure 1.7'.

Well, as you can see, the procedure takes places 10,000 times. The output is:

{4, 4, 21, 18, 2, 4, 21, 3, 13, 8, 24, 10, 30, 28, 22, 4, 12, 5, 15, 14, 11, 2, 6, 18, 23, 7, 21, 1, 3, 9}.


Unfortunately, even if some of the points look like the ones in Figure 1.7, obviously those images are not the same. Actually, Figure 1.7 contains 122 points (assuming I didn't miscount), however, given the generator, I don't know how you can produce more than 30 points. Actually, the author writes:

... Remember that there are only 30 different values, however.

Am I missing something? Any ideas?

• Your first two code lines are more compactly written as m=Rest[NestList[Mod[3 #, 31] &, 9, 30]]. I don't have enough time to thoroughly look at your problem now but will get back to you later. In any event, LCGs are supposed to show those "planes" noted by Marsaglia, and should certainly show up for so tiny a modulus. – J. M. is a poor mathematician Oct 21 '10 at 1:14
• You are right. It is a nice change, although not so intuitive. Sure, those parallel hyperplanes (in this case are simply lines) are shown, however, the result of performing the shuffling is to reduce the lattice structure. I haven't read Marsaglia's paper ('Random numbers fall mainly in the planes') but I think that those planes are a necessity for '... overlapping sequences of the output of a congruential generator with modulus m...'. By the way, thank you for the help. – Robert Smith Oct 21 '10 at 3:03

I know this question just celebrated its 8th birthday, sorry. I always think that better late than never.

The range of the output has nothing to do with the period. You can have a truly random generator that produces bits one by one, and it won't be periodic, even if its output only has 2 different values. In your case, the range of the output is from 1 to 30, rather than from 0 to 1, but the idea is the same.

Now, the output of the LCG does have a period of 30, but the output of the shuffle does not. Numbers that enter the Bays-Durham array aren't extracted in the same order they entered, but in a pseudo-random order; a particular number in the array may be delayed for quite some draws compared to others. This changes the periodicity: while the LCG's state can only hold 30 possible values, and hence its period can't be more than 30, the Bays-Durham array and internal index provide additional state that alters the period, without changing the output range nor significantly altering the probability of any particular number. An upper bound for the period of this combined generator is given by the number of its possible states, which if my numbers aren't wrong, is $$30^9·8 = 1.57464·10^{14}$$ (eight numbers from 1 to 30 in the array, one number from 1 to 30 in the LCG and 8 possible values for the internal index); that's approximately 47.16 bits of information. I don't know if it will achieve that period, though, but that upper limit suggests that it's not expected to be short.

As another example, if the input to the shuffle is the sequence $$1, 2, 3, ..., 29, 30, 1, 2, 3, ...$$ which obviously has period 30, the first 100 numbers of the output are:

3, 1, 11, 10, 13, 4, 12, 14, 15, 18, 5, 2, 16, 20, 6, 21, 24, 7, 25, 27, 29, 8, 28, 1, 22, 26, 30, 3, 4, 8, 2, 9, 17, 23, 5, 10, 12, 16, 13, 19, 18, 20, 14, 19, 21, 22, 25, 6, 15, 23, 26, 27, 1, 11, 17, 24, 2, 3, 7, 28, 7, 9, 4, 8, 11, 12, 15, 29, 10, 16, 5, 14, 17, 20, 30, 18, 23, 24, 6, 21, 27, 28, 25, 1, 13, 22, 30, 2, 4, 8, 29, 7, 10, 19, 26, 3, 9, 13, 5, 12

No sign of periodicity there, as you can see, even though the input is anything but random. The number 3, for example, appears 4 times, in positions 1, 28, 58 and 96. It's followed by 1 the first time, by 4 the second, by 7 the third, and by 9 the fourth.

If you draw the output in pairs, you can of course plot them in a 30×30 square, because you will have lots of (pseudo)independent pairs. You're not limited to 15 different pairs; you would if you took the period-30 output of the LCG.

This also explains why Bays-Durham breaks the "hyperplanarity" of the LCG: the numbers that were consecutive are no longer consecutive, but instead separated by a (pseudo)random distance, therefore if we try to plot pairs of numbers, we won't be plotting numbers that are consecutive in the LCG.

The above series uses Knuth's method to derive the index, which applied to this case, is $$j = \lfloor\frac{8·(X_i-1)}{30}\rfloor$$ (as a zero-based index). But the book uses a different method: $$j = X_i\bmod 8$$, adding 1 to make it 1-based. With that method, I get the following plot out of 5,000 pairs:

I don't know how it compares to either of the images because I can't see them.

• Thank you for your response. It has been a while, so I need to read again this chapter to remember some details. I was able to obtain Figure 1.7: imgur.com/a/7pOmdzq for the comparison. It looks different, although to some extent, this might be normal result with a different method. – Robert Smith Oct 28 '18 at 7:05
• Thanks for the image. Wow, my plot looks nothing like that, and what's more, it's not supposed to look like that. I was surprised by the near absence of consecutive pairs with the same value (the diagonal), but I wasn't surprised in the least that the points almost filled the 30×30 square otherwise. That's why 5,000 points, according to the wikibook, sounded excessive to me. I expected the image in the book to be similar to mine. I wonder if they had a bug in the implementation they used to make the plot. – Pedro Gimeno Oct 28 '18 at 15:26
• It is possible. I'm also confused about the differences in both plots. – Robert Smith Oct 30 '18 at 18:06