Re: Troubling Bug involving RandomReal, NestList, and Table
- To: mathgroup at smc.vnet.net
- Subject: [mg86286] Re: Troubling Bug involving RandomReal, NestList, and Table
- From: Szabolcs Horvát <szhorvat at gmail.com>
- Date: Fri, 7 Mar 2008 02:32:03 -0500 (EST)
- Organization: University of Bergen
- References: <fqo8tn$t2o$1@smc.vnet.net>
jwmerrill at 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!?!
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