Re: Re: Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51395] Re: [mg51309] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation*From*: DrBob <drbob at bigfoot.com>*Date*: Fri, 15 Oct 2004 02:48:49 -0400 (EDT)*References*: <200410141035.GAA14839@smc.vnet.net>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

The speed difference can be pretty dramatic: Quit a=First@Timing[Table[Random[Real,{0,1000}],{10^6}];] <<RandomReplacement` b=First@Timing[Table[Random[Real,{0,1000}],{10^6}];] b/a 0.375 Second 10.078 Second 26.8747 I can probably live with it, however, in cases where I'm concerned about built-in Random's non-random behavior. Bobby On Thu, 14 Oct 2004 06:35:37 -0400 (EDT), Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > I think I have now fixed the problem reported by Mark Fisher. The new > version of the RandomReplacement now works like this: > > > SeedRandom[5] > > > Timing[First[Table[Random[],{10000}]]] > > > {0.01 Second,0.786599} > > > <<RandomReplacement` > > Timing[First[Table[Random[],{10000}]]] > > > {0.02 Second,0.409222} > > > As you can see there is a loss of speed involved. There may also be > still some as yet undiscovered problems. However, the new package is > available form my sites, though I still should write some documentation > of the latest changes (a few bugs in the previous version have also > been fixed). > > Andrzej Kozlowski > > Andrzej Kozlowski > Chiba, Japan > http://www.akikoz.net/~andrzej/ > http://www.mimuw.edu.pl/~akoz/ > > > > On 12 Oct 2004, at 14:57, Andrzej Kozlowski wrote: > >> It's clearly because of packed arrays: >> >> >> <<Developer` >> >> >> PackedArrayQ[Table[Random[],{249}]] >> >> >> False >> >> >> PackedArrayQ[Table[Random[],{250}]] >> >> >> True >> >> >> >> Obviously Mathematica uses a different method of generating random >> packed arrays. There are various ways of getting round this (e.g. >> constructing lists as joins of lists of length less than 250 etc) but >> one will loose the benefits of packed arrays and with that, presumably, >> a lot of speed. This seems to be something that only WRI can change. >> Note however that not all situations in which one uses Random[] require >> generating this type of lists. >> Moreover, for problems involving the normal distribution it seems to me >> that using the Marsaglia generator gives good results even with the >> built in Random[]. But certainly we should hope that something will >> finally be done do deal with this issue in version 6. >> >> Andrzej >> >> >> >> >> On 11 Oct 2004, at 17:24, DrBob wrote: >> >>> He's right! How the devil does that happen? >>> >>> Here's a test, with Andrzej's package loaded in my Init file. >>> >>> First, with n=250: >>> >>> Quit >>> SeedRandom[5] >>> Table[Random[],{250}]; >>> Last@% >>> >>> 0.107874 >>> >>> Unprotect[Random]; >>> Clear@Random >>> SeedRandom[5] >>> Table[Random[],{250}]; >>> Last@% >>> >>> 0.107874 >>> >>> Now, with n=249: >>> >>> Quit >>> SeedRandom[5] >>> Table[Random[],{249}]; >>> Last@% >>> >>> 0.656099 >>> >>> Unprotect[Random]; >>> Clear@Random >>> SeedRandom[5] >>> Table[Random[],{249}]; >>> Last@% >>> >>> 0.0708373 >>> >>> Bobby >>> >>> On Mon, 11 Oct 2004 01:25:24 -0400 (EDT), Mark Fisher >>> <mark at markfisher.net> wrote: >>> >>>> FYI: I've just a little testing and I find that Mathematica ignors >>>> the user >>>> defined rules for Random in Table[Random[],{n}] when n >= 250. >>>> >>>> Andrzej Kozlowski wrote: >>>>> The problem is not actually with the way Mathematica's >>>>> NormalDistribution but with the uniformly distributed Random[] >>>>> function. >>>>> NormalDistribution itself is generated (in the package >>>>> Statistics`NormalDistribution`) by means of the very classical so >>>>> called Boox-Muller method (although actually the Marsaglia variant >>>>> below works better). You can try downolading my little >>>>> RandomReplacement package form one of my web sites: >>>>> >>>>> http://www.akikoz.net/~andrzej//Mathematica/ in Japan >>>>> >>>>> or >>>>> >>>>> http://www.mimuw.edu.pl/~akoz/Mathematica/ in Poland >>>>> >>>>> The package is based on the post of daniel Lichtblau, which explains >>>>> the nature of the problem: >>>>> >>>>> <http://library.wolfram.com/mathgroup/archive/2000/May/ >>>>> msg00088.html> >>>>> >>>>> >>>>> >>>>> The package is named init.m and can be loaded automatically at the >>>>> start of each Mathematica session. >>>>> However, in fact if you are only concerned with normal distributions >>>>> it >>>>> may be enough to the following. First load the (unmodified) >>>>> NormalDistribution package using >>>>> >>>>> <<Statistics`NormalDistribution` >>>>> >>>>> and then evaluate the code below. It will replace the Box-Muller >>>>> generator by the Marsaglia one. I have found that this is usually >>>>> enough. But if your problem still remains try the RandomReplacement >>>>> package. (Or you can use both of course). Here is the Marsaglia code >>>>> for normal distribution: >>>>> >>>>> >>>>> Statistics`NormalDistribution`Private`normal= >>>>> Compile[{{mu,_Real},{sigma,_Real},{q1,_Real},{q2,_Real}}, >>>>> Module[{va=1.,vb,rad=2.,den=1.}, >>>>> While[rad>=1.,va=2.*Random[]-1.;vb=2.*Random[]-1.; >>>>> rad=va*va+vb*vb]; >>>>> den=Sqrt[-2.*(Log[rad]/rad)]; >>>>> mu+sigma*va*den]]; >>>>> >>>>> >>>>> >>>>> On 5 Oct 2004, at 17:36, Chris Kunkel wrote: >>>>> >>>>> >>>>>> *This message was transferred with a trial version of >>>>>> CommuniGate(tm) >>>>>> Pro* >>>>>> Hello, >>>>>> >>>>>> Recently I have encountered a problem in Mathematica's normal >>>>>> distribution random number generator. The problem arises when I >>>>>> look >>>>>> at the distribution of averages of a list of theses numbers. That >>>>>> is, >>>>>> I generate 1000 numbers and take their average. I do this a number >>>>>> of >>>>>> times and the plot a frequency distribution. Consistently it seems >>>>>> to >>>>>> be skewed positive. Specifically, the number of occurrences less >>>>>> than >>>>>> 3 standard deviations is consistent with the expected number, but >>>>>> the >>>>>> number greater than 3 is always larger (in some cases about twice >>>>>> the >>>>>> expected number). >>>>>> >>>>>> I was wondering if anyone else has noticed this and knows of a fix. >>>>>> Also, does anyone have code for a good quality normal distribution >>>>>> random number generator? >>>>>> >>>>>> Chris >>>>>> >>>>> >>>>> >>>> >>>> >>>> >>>> >>> >>> >>> >>> -- >>> DrBob at bigfoot.com >>> www.eclecticdreams.net >>> >> >> > > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**References**:**Re: Re: Re: normal distribution random number generation***From:*Andrzej Kozlowski <akoz@mimuw.edu.pl>