Re: Re: normal distribution random number generation

*To*: mathgroup at smc.vnet.net*Subject*: [mg51266] Re: [mg51217] Re: normal distribution random number generation*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Mon, 11 Oct 2004 01:25:25 -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

I have put on both my web sites a new version of the Random[] fix incorporating the ideas of Ray and Bob. The package is now called RadnomReplacement.m and is a proper Mathematica package, which defines its own context: RandomReplacement`. When loaded, it will replace the built in Random[] function with one based on the Wolfram CA algorithm. I have also put on the same page another "package' called MarsagliaNormal.m which changes the definition of the normal distribution generator in the Statistics`NormalDistribution` package to that of Marsaglia's. This is not a "true package", in that it does not define any new contexts, and has to be read in using Get and its full name including the extension *.m. It can be read in after or before the Statistics`NormalDistribution` package. In the latter case it will automatically load the package and seamlessly redifine the normal generator. Doing this seems to have two advantages. It seems to be marginally faster than the Box-Muller generator, though this may be hardware dependent. It also seems to "eliminate' the problem caused by "insufficient randomness' of Random[] so if one is dealing only with normally distributed quantities it may not be necessary to read in the RandomReplacement package. Andrzej 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)