Troubling Bug involving RandomReal, NestList, and Table

[jwme...@gmail.com]

When making tables of 250 or more directions, Mathematica seems to
forget how to do arithmetic. I realize these "random" directions are
not properly uniform, but that's beside the point.

$Version
> "6.0 for Mac OS X x86 (32-bit) (April 20, 2007)"

randomDirection3[] :=
NestList[Sqrt[1 - #^2] &, RandomReal[], 1]

ListPlot[Table[randomDirection3[], {249}], AspectRatio -> Automatic]
ListPlot[Table[randomDirection3[], {250}], AspectRatio -> Automatic]

> [Graphs you'll need to generate for yourself]

Variance[Table[Norm[randomDirection3[]], {249}]]
Variance[Table[Norm[randomDirection3[]], {250}]]

> 1.49104*10^-33
> 0.0513843

What!?!

Jason Merrill

[Mark Fisher]

This is bad. It looks to me like the bug is located in compiling
randomDirection3. (I'm guessing Table compiles its first argument when
the length of the list is 250 or greater.) Consider this:

rd = Compile[{}, NestList[Sqrt[1 - #^2] &, RandomReal[], 1]]

rd[] returns a pair of numbers that don't lie on the unit circle.

I'm using 6.0.1. Does anyone know if this is fixed in 6.0.2?

Hey Dan and Rob! What gives?

--Mark

[Jens-Peer Kuska]

Hi,

this is a bug when Mathematica "Compile" the expression
in Table[] try

opts = "CompileOptions" /. SystemOptions[];
opts = opts /. ( "TableCompileLength" -> _) :> (
"TableCompileLength" -> Infinity);
SetSystemOptions["CompileOptions" -> opts];

and the bug disappear.

Regards
Jens

[Darren Glosemeyer]

jwme...@gmail.com wrote:
> When making tables of 250 or more directions, Mathematica seems to
> forget how to do arithmetic. I realize these "random" directions are
> not properly uniform, but that's beside the point.
>
> $Version
>
>> "6.0 for Mac OS X x86 (32-bit) (April 20, 2007)"
>>
>
> randomDirection3[] :=
> NestList[Sqrt[1 - #^2] &, RandomReal[], 1]
>
> ListPlot[Table[randomDirection3[], {249}], AspectRatio -> Automatic]
> ListPlot[Table[randomDirection3[], {250}], AspectRatio -> Automatic]
>
>
>> [Graphs you'll need to generate for yourself]
>>
>
> Variance[Table[Norm[randomDirection3[]], {249}]]
> Variance[Table[Norm[randomDirection3[]], {250}]]
>
>
>> 1.49104*10^-33
>> 0.0513843
>>
>
> What!?!
>
> Jason Merrill
>

This appears to be Table compilation problem. Note that problems occur
when the table size hits the default length for compiling tables.

In[1]:= "TableCompileLength" /. ("CompileOptions" /.
Developer`SystemOptions[])

Out[1]= 250

I've passed this along to the appropriate people in the company so it
can be investigated further.

One possible solution would be to increase the value of
"TableCompileLength" (I'll show this below), but I would suggest instead
writing a function to get n values which gets all RandomReal values in
one call. This avoids the problem and will be faster because only one
call to RandomReal is needed to get n random values.

In[2]:= randomDirection[n_] :=
Map[NestList[Sqrt[1 - #^2] &, #, 1] &, RandomReal[1, n]]

In[3]:= Timing[randomDirection[10^6];]

Out[3]= {1.984, Null}

The following shows the approach of increasing the value of
"TableCompileLength" and obtaining values via a Table of
randomDirection3[] values.

In[4]:= Developer`SetSystemOptions[
"CompileOptions" -> {"TableCompileLength" -> Infinity}];

In[5]:= randomDirection3[] := NestList[Sqrt[1 - #^2] &, RandomReal[], 1]

In[6]:= Timing[Table[randomDirection3[], {10^6}];]

Out[6]= {7.047, Null}

Plots for points from both of the approaches above appear as expected.

Darren Glosemeyer
Wolfram Research

[Szabolcs Horvát]

Here's my guess about what's wrong:

Functions like Table compile their arguments when used with lists larger
than a certain size. This size is 250 by default for Table, but it can
be changed like this:

SetSystemOptions["CompileOptions" -> "TableCompileLength" -> 1000]

(Use SystemOptions[] and SystemOptions["CompileOptions"] to list
available options. Unfortunately these things are not very well
documented.)

When randomDirection3[] is compiled, something goes wrong. We get the
same incorrect result if we explicitly compile the function:

cp = Compile[{}, NestList[Sqrt[1 - #^2] &, RandomReal[], 1]]

and use cp[] instead of randomDirection3[].

It seems that Compile inlines RandomReal[], disregarding the fact that
RandomReal[] may return different values on subsequent calls, and the
compiled function is equivalent to

{RandomReal[], Sqrt[1 - RandomReal[]^2]}

One workaround that you could use (without giving up compilation) is this:

cp = Compile[{}, With[{r = Random[]}, {r, Sqrt[1 - r^2]}]]

(But I am sure you already found many other ways around the problem.)

Szabolcs

[jwme...@gmail.com]

On Mar 7, 2:48 am, Szabolcs Horv=E1t <szhor...@gmail.com> wrote:

> jwmerr...@gmail.com wrote:
> > When making tables of 250 or more directions, Mathematica seems to
> > forget how to do arithmetic. I realize these "random" directions are
> > not properly uniform, but that's beside the point.

> > randomDirection3[] :=

> > NestList[Sqrt[1 - #^2] &, RandomReal[], 1]

> > Variance[Table[Norm[randomDirection3[]], {249}]]

> > Variance[Table[Norm[randomDirection3[]], {250}]]
>
> >> 1.49104*10^-33
> >> 0.0513843

> One workaround that you could use (without giving up compilation) is this:=

>
> cp = Compile[{}, With[{r = Random[]}, {r, Sqrt[1 - r^2]}]]
>
> (But I am sure you already found many other ways around the problem.)
>
> Szabolcs

Thanks to all for your explanations and advice. My eventual solution
(pilfered from elsewhere on the internet), was

randomUnitVector[dimensions_] :=
Normalize[RandomReal[NormalDistribution[], dimensions]]

Regards,

JM