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.

As this is your seventh question, it is about time you learn how to properly format your code. Which is given within the first item of the FAQ. — rcollyer, Apr 29 '13 at 17:01
I fixed up the formatting of the lists to make them more readable. Also, I encourage you using the revisions link to look at how I set up the inline formatting. — rcollyer, Apr 29 '13 at 17:17
@rcollyer, thanks for that. I will try to follow it more carefully from now on. — dbm, Apr 29 '13 at 17:18
I'm just trying to encourage effort. I will willingly (if not happily) fix up bad formatting, but I'd like to see an attempt at formatting, first. So, thanks. — rcollyer, Apr 29 '13 at 17:22

wolfies · Answer 1 · 2013-04-29T22:34:02.193

28

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.

edited Apr 29 '13 at 22:34

answered Apr 29 '13 at 17:18

wolfies

8,722
1
25
54

I still don't understand how Infinity can be an output of a random number generator (is the number bigger than 10^6). I really think an infinity won't be an output of a random number generator. Anyway, if I have to get rid off WorkingPrecision, then after that I have an expression as below: e={2.0799361919940695x[1] + 3.3534325557330327 x[1]^2 - 4.335179297091139` x[1] x[2]}
How would you chop off the coefficients to 6 digits from such an expression?
– dbm Apr 29 '13 at 17:30
If I use N[e,6], then it also changes the indexes on x's, i.e., it shows x[1.00000] etc. instead of x[1]. – dbm Apr 29 '13 at 17:33
That's a fairly different question. Look up PaddedForm, Numberform or even AccountingForm if you want to micromanage how numbers are displayed. – bill s Apr 29 '13 at 17:35
@bill s, Ok. I will have a look at them. However, for RandomReal function, Mathematica's help pages gives the following example: RandomReal[{-1, 1}, 5, WorkingPrecision -> 5]
So, WorkingPrecision should indeed work properly as I had expected it to work!
– dbm Apr 29 '13 at 17:48
Sure, but working precision is a poor way of controlling how your numbers are printed out. – bill s Apr 29 '13 at 17:50
5

@bill s, I really think that this is either a bug or at the least it should have been mentioned in the 'Possible Issues' of Mathematica's relevant help pages. It is a reasonable requirement: I just wanted to generate a few random numbers with a certain precision. – dbm Apr 29 '13 at 18:02
@dbm > RandomReal[{-1, 1}, 5, WorkingPrecision -> 5] ---> well, that's the jackpot ... that's how your 6 digit precision enters ... via RandomReal ... which the Help system notes can generate a 0 ... see updated answer above. – wolfies Apr 29 '13 at 18:24
5

Why not just generate the numbers at the original precision, and then use SetPrecision on the result, rather than making RandomReal do it? – Guillochon Apr 29 '13 at 18:25
@wolfies, Got it! Thanks! Guillochon, I used SetPrecision, but then I have a different issue: http://mathematica.stackexchange.com/questions/24318/setprecision-with-indexed-variables – dbm Apr 29 '13 at 18:33
@Guillochon I suppose the real question is, why doesn't RandomReal do this already? In addition to giving rise to obscure bugs such as the one noted here, it seems like a complete waste of time for it to work at arbitrary precision (below machine precision) rather than just using machine reals and truncating the numbers afterward. – Oleksandr R. Apr 29 '13 at 18:33
2

+1 for diagnostic effort and clarity in writing. – RBerteig Apr 29 '13 at 19:02
4

Note SeedRandom[8396] ; RandomReal[{0, 1}, {20}, WorkingPrecision -> 6] gives a zero.. in fact the chance of an exact zero is about 1 in 10^6 .. indicating they are truncating to 6 places, not using 6 digits precision – george2079 Apr 29 '13 at 20:48
2

Actually, you are using a form of RandomReal that is not currently documented, but using RandomVariate gives the same error. But I think you are right, in RandomVariate's documentation, WorkingPrecision is wrongly documented. WorkingPrecision is a common option to many functions, and as you can see in it's own documentation page, it sets the precision used for internal computations and not the result's precision. In your example clearly 6 digits of precision in its internal computations aren't enough to guarantee 6 digits of output's precision, or no Infinity should be generated – Rojo Apr 30 '13 at 11:14

Mathematica giving infinity as a random number!

1 Answers1