Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}.

Plotting the successive pairs, I get this:

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?

share|cite|improve this question
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. 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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.