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)
- Re: normal distribution random number generation