Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51323] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: sean kim <sean_incali at yahoo.com>
- Date: Thu, 14 Oct 2004 06:36:38 -0400 (EDT)
- Reply-to: sean_incali at yahoo.com
- Sender: owner-wri-mathgroup at wolfram.com
i just took a look at Part[plotvec, 24] the respective values, for me, are, 375 and 686. so i did a test with the following. in the fresh kernel, I get exactly the same thing as you have. but I decided to evaluate the code In[5]:= SeedRandom[1111111 ]; vec1=Table[Random[Real,{0,1000}],{10^6}]; vec2=Map[If[#<44,1,0]&,vec1]; ones=Flatten[Position[vec2,1]]; plotvec=Map[Length, Split[Sort[Drop[ones-RotateRight[ones,1],1]]]]; Part[plotvec, 24] ListPlot[Log[plotvec],PlotRange->All] ListPlot[plotvec,PlotRange->All] Out[10]= 375 and then reevaluate the codes In[13]:= SeedRandom[5] In[14]:= Timing[First[Table[Random[],{10000}]]] Out[14]= {0.01 Second,0.786599} In[15]:= <<RandomReplacement` Timing[First[Table[Random[],{10000}]]] Out[16]= {0.01 Second,0.611425} and I get the different value than before(one with the fresh kernel) i wonder what's going on here? sean --- 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: > > > *This message was transferred with a trial version > of CommuniGate(tm) > > Pro* > > 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 > >> > > > > > > __________________________________ Do you Yahoo!? Read only the mail you want - Yahoo! Mail SpamGuard. http://promotions.yahoo.com/new_mail