Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51259] Re: [mg51217] Re: normal distribution random number generation*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sun, 10 Oct 2004 05:52:33 -0400 (EDT)*References*: <ck0ccp$o1u$1@smc.vnet.net> <200410090818.EAA09618@smc.vnet.net> <opsfm2itcmiz9bcq@monster.cox-internet.com>*Sender*: owner-wri-mathgroup at wolfram.com

It's probably best to make a proper package out of this. The function randomSubstitutionFunction should be in the Private` context, invisible to the user while Random[] will be exported, thus replacing the built in Random[]. One can still have the package load at startup, if one choses to. I think I will put this version on my site instead of the earlier fix (with an acknowledgement of all the contributors of course ;-) ) Also, I am going to include a similar short package that will replace the Box-Muller normal random generator in the Statistics`NormalDistribution` package by the Marsaglia version. As I have mentioned before, I have found that when generating normally distributed data this replacement seems to make up for the weakness of the Mathematica Random[] function and is actually a little faster than the Box-Muller normal distribution generator. Andrzej Kozlowski On 10 Oct 2004, at 14:21, DrBob wrote: > > I'm trying to combine that idea with Andrzej Kozlowski's recent fix > for Random, and here's what I came up with: > > Unprotect[Random]; > With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2}, > randomSubstitutionFunction = > Compile[{}, ((Random[Integer, m2] + .5)*m1 + Random[Integer, > m2])*m1]; > Random[] := randomSubstitutionFunction[] > ] > Random[Real, {a_Real, b_Real}] := a + Random[]*(b - a) > Random[Real, b_Real] := Random[Real, {0, b}] > Random[Real] := Random[Real, {0, 1}] > Random[Complex, {a_Complex | a_Real | a_Integer, b_Complex | b_Real | \ > b_Integer}] := a + Random[]*Re[(b - a)] + Random[]*Im[(b - a)]*I > Random[Complex] := Random[Complex, {0, 1 + I}] > Protect[Random]; > > I wanted NOT to use a Global (randomSubstitutionFunction) for the > Compiled function, but I haven't stumbled on a way to manage it. > > Bobby > > On Sat, 9 Oct 2004 04:18:30 -0400 (EDT), Ray Koopman <koopman at sfu.ca> > wrote: > >> Bill Rowe <readnewsciv at earthlink.net> wrote in message >> news:<ck0ccp$o1u$1 at smc.vnet.net>... >>> [...] >>> you will have modified Random to use the Wolfram rule 30 cellular >>> automaton and avoid the subtract with borrow algorithm. The main >>> consequence of this is Random will now be considerably slower. >>> [...] >> >> If time is an issue, you might want to consider generating integers >> on 0...2^n-2 instead of 0...2^n-1. It's always much faster. And if >> you're willing to spend a little of the time you've saved, you can >> add a half and avoid ever having to worry about getting a zero. >> >> In[1]:= ToString[TableForm[Table[With[{m1 = 2^n - 1, m2 = 2^n - 2}, >> {n, First[Timing[Do[Random[Integer,m1],{1*^6}]]]/.Second->1., >> >> First[Timing[Do[Random[Integer,m2],{1*^6}]]]/.Second->1.}], >> {n,2,30}],TableSpacing->{0,2}]] >> >> Out[1]= 2 1.96 1.42 >> 3 2.12 1.5 >> 4 2.38 1.61 >> 5 2.66 1.73 >> 6 2.91 1.86 >> 7 3.16 2. >> 8 3.41 2.1 >> 9 3.68 2.19 >> 10 3.92 2.35 >> 11 4.21 2.56 >> 12 4.5 2.68 >> 13 4.79 2.82 >> 14 5.07 3.02 >> 15 5.34 3.08 >> 16 5.56 3.26 >> 17 5.84 3.38 >> 18 6.09 3.53 >> 19 6.33 3.64 >> 20 6.57 3.77 >> 21 6.84 3.87 >> 22 7.1 4.03 >> 23 7.33 4.2 >> 24 7.63 4.25 >> 25 7.89 4.37 >> 26 8.15 4.56 >> 27 8.4 4.61 >> 28 8.56 4.79 >> 29 8.95 4.95 >> 30 9.16 5.07 >> >> In[2]:= ran1 = With[{m = 2.^-30, m1 = 2^30 - 1}, >> Compile[{},(Random[Integer,m1]*m + Random[Integer,m1])*m]]; >> >> In[3]:= ran2 = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2}, >> Compile[{},(Random[Integer,m2]*m1 + Random[Integer,m2])*m1]]; >> >> In[4]:= ran2h = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2}, >> >> Compile[{},((Random[Integer,m2]+.5)*m1+Random[Integer,m2])*m1]]; >> >> In[5]:= First/@{Timing@Do[ran1[],{1*^5}],Timing@Do[ran2[],{1*^5}], >> Timing@Do[ran2h[],{1*^5}]} >> Out[5]= {2.03 Second, 1.05 Second, 1.08 Second} >> >> >> >> > > > > -- > DrBob at bigfoot.com > www.eclecticdreams.net >

**References**:**Re: normal distribution random number generation***From:*koopman@sfu.ca (Ray Koopman)