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