Re: Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51324] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Thu, 14 Oct 2004 06:36:50 -0400 (EDT)*References*: <20041013010600.74710.qmail@web60505.mail.yahoo.com>*Sender*: owner-wri-mathgroup at wolfram.com

Are you using the latest version of the RandomReplacement package? (The one I posted last night). The earlier one had a bug and ddi not work with calls to Random[Real,{a,b}] where a and b are integers. You had to use Random[Real,N[{a,b}]]. I think (I hope) I fixed that in the latest version. Andrzej On 13 Oct 2004, at 10:06, sean kim wrote: > *This message was transferred with a trial version of CommuniGate(tm) > Pro* > 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 message was transferred with a trial version >> of CommuniGate(tm) Pro* >> 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 >

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

**Mathematica Symbolic Toolbox for MATLAB--Version 2.0**

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

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