Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
5  
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
5  
SeedRandom is intended to be used with 42 :) –  belisarius Apr 29 '13 at 17:08
1  
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
3  
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
add comment

1 Answer

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.

share|improve this answer
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
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
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
1  
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
show 8 more comments

Your Answer

 
discard

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.