Re: Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51322] 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:32 -0400 (EDT)*Reply-to*: sean_incali at yahoo.com*Sender*: owner-wri-mathgroup at wolfram.com

I know... It's werid for me... which i wanted to ask you all. when i used the code below I still got the low 24 bin value. In[1]:=<<RandomReplacement` In[2]:= 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]]]]; ListPlot[Log[plotvec],PlotRange->All] ListPlot[plotvec,PlotRange->All] used it exactly as it is shown. but the below works, it seems... SeedRandom[1111111]; vec1 = Table[Random[Integer,1000],{10^6}]; vec2 = Map[If[#<=44,1,0]&, vec1]; ones = Flatten[Position[vec2,1]]; Length[ones]/10.^6 plotvec = Map[Length, Split[Sort[Drop[ones-RotateRight[ones,1],1]]]]; ListPlot[ Log[plotvec], PlotRange->All] ListPlot[plotvec, PlotRange->All] maybe other people can try it and see if it's just me? (I took out the init.m btw) sean --- Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > This is exactly what it is mean to do. The problem > was due to the > usage of the SWB algortithm by the built in Random[] > for generating > random reals. RandomReplace makes Mathematica use > the Wolfran CA random > number generator, which Mathematica uses for > generating random integers > instead of the SWB. This "cures" the problem > reported int hat message. > Have you got any reason to think it does not? > > Andrzej > > > On 13 Oct 2004, at 03:48, sean kim wrote: > > > *This message was transferred with a trial version > of CommuniGate(tm) > > Pro* > > Hi andrej > > > > thank you so much for making your package > available. > > > > I think I may have found another problem... ("may > have > > found" being operative phrase, or non-operative > for > > that matter) > > > > http://groups.google.com/groups? > > > hl=en&lr=&newwindow=1&safe=off&threadm=a5t1fq%245vt%241%40smc.vnet.net& > > > rnum=1&prev=/ > > > groups%3Fhl%3Den%26lr%3D%26newwindow%3D1%26safe%3Doff%26q%3Dchange%2Bth > > > > is%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253 > > > Dcomp.soft-sys.math.mathematica > > > > shows a thread dealing with random reals and > getting a > > much smaller frequency at position 24 than others > > positions. > > > > Daniel posted a work around showing the use of > Random > > Integers to fix the problem, which as far I > > understand, the Random replacement also > implements. > > > > It seems the RandomReplacement doesn't take care > of > > the appearance of the lower than expected > frequency at > > position 24(or interval 24 as the OP had noted). > > > > I must be missing something. > > > > thanks all for any insights. > > > > sean > > > > > > --- Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > > > >> *This message was transferred with a trial > version > >> of CommuniGate(tm) Pro* > >> 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!? > > Yahoo! Mail Address AutoComplete - You start. We > finish. > > http://promotions.yahoo.com/new_mail > > > > __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com