Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51267] Re: [mg51217] Re: normal distribution random number generation
- From: DrBob <drbob at bigfoot.com>
- Date: Mon, 11 Oct 2004 01:25:26 -0400 (EDT)
- References: <ck0ccp$o1u$1@smc.vnet.net> <200410090818.EAA09618@smc.vnet.net> <opsfm2itcmiz9bcq@monster.cox-internet.com> <95492AA4-1AD1-11D9-BEB1-000A95B4967A@mimuw.edu.pl>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
I really didn't contribute enough to earn a mention, Andrzej. I'm still scratching my head trying to eliminate the named helper function. Bobby On Mon, 11 Oct 2004 00:32:23 +0900, Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > 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 >> > > > > -- DrBob at bigfoot.com www.eclecticdreams.net
- References:
- Re: normal distribution random number generation
- From: koopman@sfu.ca (Ray Koopman)
- Re: normal distribution random number generation