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