# On random numbers extracted from NormalDistribution by two different methods under the same seed

I've encountered the following situation within various trials.

Why does such a difference arise between bulk sampling and one-by-one sampling?

(*Bulk sampling*)
SeedRandom[1234];
RandomVariate[NormalDistribution[], 10]

(*sampling one-by-one*)
SeedRandom[1234];
Table[RandomVariate[NormalDistribution[]], 10]


Example: For example, consider the following game.

• The player gains points depending on the value of variable x every time.
• x is a random number, and is subject to NormalDistribution [0, 1]
• If x> 0 then the player gets + 1,otherwise, x <0, the player gets -1.

I believed that expectation value of game scores should be equal regardless of how random numbers are generated.

However, the experiment is as follows.

{BlockRandom[SeedRandom[1234];
N@Mean@Table[
Total[If[# > 0, 1, -1] & /@
RandomVariate[NormalDistribution[], 10]], {100}]],
BlockRandom[SeedRandom[1234];
N@Mean@Table[
Total[If[# > 0, 1, -1] & /@
Table[First@RandomVariate[NormalDistribution[], 1],
10]], {100}]]}


{0.18, -0.68}

Histogram[{BlockRandom[SeedRandom[1234];
Table[Total[
If[# > 0, 1, -1] & /@
RandomVariate[NormalDistribution[], 10]], {100}]],
BlockRandom[SeedRandom[1234];
Table[Total[
If[# > 0, 1, -1] & /@
Table[First@RandomVariate[NormalDistribution[], 1],
10]], {100}]]}]


I didn't understand why this happens and found out the phenomenon in this question.

• Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Feb 2 '19 at 15:08

I think this is due to the Box–Muller transform that is employed to generate pairs of normally distributed values. This method is much faster than using the inverse commulative probability function.

The algorithm is essentially as follows:

n = 1000000;
a = Flatten[ReIm[
Times[
Sqrt[-2. Log[RandomReal[{0, 1}, n]]],
Exp[(2. Pi I) RandomReal[{0, 1}, n]]
]
]]; // AbsoluteTiming //First


0.003597

As comparison:

b = RandomVariate[NormalDistribution[], 2 n]; // AbsoluteTiming // First
Histogram[{a, b}]


0.003579

Although I cannot reproduce the precise order used in Mathematica, the timings are very similar. I guess, this supports my hypothesis.

Whenever only a single sampling of the normal distribution is required, the second sampled value is just thrown away. This explains why the results of the bulk sampling equal the results one-by-one sampling, interleaved with other pseudorandom numbers.

• Box-Muller generates two independent Normals. What difference does it make if you drop one of them? None. – wolfies Feb 2 '19 at 12:49
• @wolfies Well, the results of SeedRandom[1234]; RandomVariate[NormalDistribution[], 10] and SeedRandom[1234]; Table[RandomVariate[NormalDistribution[]], 10] are different. But one would not expect that. For example, SeedRandom[1234]; RandomVariate[ExponentialDistribution[1], 10] SeedRandom[1234]; Table[RandomVariate[ExponentialDistribution[1]], 10] return exactly the same result. – Henrik Schumacher Feb 2 '19 at 13:03
• I don't agree that there is any guarantee or requirement that calling 10 pseudorandom numbers from a particular distribution all in one go is required to use the same generator as calling one pseudorandom number at a time. Indeed, if code is well-optimised, one might expect that different algorithms are used to generate multiple values all in one go, rather than single ones. The very same thing happens with the LogNormalDistribution and HalfNormalDistribution, possibly because they use the Normal generator. – wolfies Feb 2 '19 at 13:43
• At the same time, +1 on your very fine (!) sleuthing and explanation as the likely cause and the algorithms being used. – wolfies Feb 2 '19 at 16:48

When I change the distribution from NormalDistribution to UniformDistribution there is no difference in the two results.

SeedRandom[1234];
RandomVariate[UniformDistribution[], 10]


{0.876608, 0.521964, 0.0862234, 0.377913, 0.0116446, 0.927266, 0.543757, 0.479332, 0.245349, 0.759896}

SeedRandom[1234];
Table[First @ RandomVariate[UniformDistribution[], 1], 10]


{0.876608, 0.521964, 0.0862234, 0.377913, 0.0116446, 0.927266, 0.543757, 0.479332, 0.245349, 0.759896}

From this I conclude that the anomalous behavior is caused by the distribution function and not by RandomVariate. This is additional support for Henrik Schumacher's hypothesis.