Re: Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51309] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Thu, 14 Oct 2004 06:35:37 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

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 >> > >

**Follow-Ups**:**Re: Re: Re: Re: normal distribution random number generation***From:*DrBob <drbob@bigfoot.com>