Mathematica giving infinity as a random number!

Question

I am trying to generate several random numbers from normal distribution using Mathematica. The following is the relevant code:

Clear["Global`*"];
SeedRandom[8396]
Sigma = 1;
Dim = 4;

A = RandomReal[NormalDistribution[0, Sigma], Dim, 
   WorkingPrecision -> 6];
B = RandomReal[NormalDistribution[0, Sigma], {Dim, Dim, Dim}, 
   WorkingPrecision -> 6]; 
F = RandomReal[NormalDistribution[0, Sigma], {Dim, Dim, Dim, Dim}, 
   WorkingPrecision -> 6]; 

Now, when I print out B, there are several infinities! B looks as follows:

{{{-1.92981, -1.21226, -0.244472, -0.00362004}, 
  {0.700307, 1.49892, 0.874067, -0.539336}, 
  {-0.643797, 0.356028, \[Infinity], -\[Infinity]}, 
  {-0.130875, -0.379856, 2.20859, 1.65716}}, 
 {{-1.87016, 0.478610, 0.261428, -1.15096}, 
  {-0.0967131, 1.75239, -0.130795, -0.488940}, 
  {-0.771886, 1.04727, -0.499386, 0.180890}, 
  {-0.844033, 0.439520, -0.153382, 0.0686604}}, 
 {{1.75313, -0.917641, -0.222227, 0.746214}, 
  {1.40456, -0.249076, 1.90326, 0.436745}, 
  {1.14792, 0.369685, 0.00756157, 0.407814}, 
  {1.47316, 1.60401, 0.923612, -0.776877}}, 
 {{0.764598, -1.58140, 1.66268, -1.93439}, 
  {-0.520680, 0.227234, 0.908374, -1.11542}, 
  {-2.28119, -1.58329, -0.536070, 0.209272}, 
  {-0.157313, -0.0102335, -0.225685, -0.198312}}}

How is it possible? It shouldn't be because of the WorkingPrecision->6, is it?! This happens to specific RandomSeed, e.g., in this case, 6893. For other RandomSeed, I don't get this problem!

I want to leave WorkingPrecision->6 as it is to avoid larger floating numbers in the computation, unless it is this specific option that's the problem. Any idea?!

Thanks.

Answer 1

So your problem reduces to:

SeedRandom[8396]
RandomReal[NormalDistribution[0, 1], {4, 4, 4}, WorkingPrecision -> 6] 

It does appear to be caused by reducing WorkingPrecision to 6, because it goes away when you stop forcing Mathematica to behave like a bad pocket calculator. I can't see any reason for you to do this... Better to leave WorkingPrecision out altogether, and impose any requirements you have afterwards using N[blah, 6] or any of the other various ways of doing this.

As to why this happens:

I don't know which generator Wolfram now uses for the NormalDistribution, but they are often of the form:

Sqrt[-2 Log[RandomReal[]]] Cos[2Pi RandomReal[]]

RandomReal[] is meant to return a value between 0 and 1. But if you reduce WorkingPrecision to 6, then the help system notes:

RandomReal[spec, WorkingPrecision -> n] yields reals with n-digit precision. Leading or trailing digits in the generated number can turn out to be 0.

and -2 Log[0] yields Infinity.

Not really their fault at all, in my view.

---

Addendum

IN reply to the interesting comment from @george2079 below ... George raises the issue of the likelihood of RandomReal[{0,1}] generating zeroes ... which under normal conditions I believe it should never do.

So, select fabled SeedRandom[42], generate 10 million pseudorandom drawings with WorkingPrecision -> i, and Count how many 0's we get for each i:

Table[
      SeedRandom[42];
      Count[    RandomReal[{0, 1}, {10000000}, WorkingPrecision -> i],  0], 
 {i, 1, 7}]

> {624548, 78050, 9905, 584, 77, 5, 0}

I was going to say that by the time we have WorkingPrecision -> 7, the issue seems to be a non-issue. But, I managed to get a 0 with WorkingPrecision -> 7 by generating 40 million values ... and have now found a 0 with WorkingPrecision -> 8 by generating some more. And after generating about 500 million values, I have now been able to generate a 0 with WorkingPrecision -> 9.