       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:= 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= 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:= ran1 = With[{m = 2.^-30, m1 = 2^30 - 1},
>>         Compile[{},(Random[Integer,m1]*m + Random[Integer,m1])*m]];
>>
>> In:= ran2 = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2},
>>         Compile[{},(Random[Integer,m2]*m1 + Random[Integer,m2])*m1]];
>>
>> In:= ran2h = With[{m1 = 1/(2.^30 - 1.), m2 = 2^30 - 2},
>>
>> Compile[{},((Random[Integer,m2]+.5)*m1+Random[Integer,m2])*m1]];
>>
>> In:= First/@{Timing@Do[ran1[],{1*^5}],Timing@Do[ran2[],{1*^5}],
>>                                          Timing@Do[ran2h[],{1*^5}]}
>> Out= {2.03 Second, 1.05 Second, 1.08 Second}
>>
>>
>>
>>
>
>
>
> --
> DrBob at bigfoot.com
> www.eclecticdreams.net
>

```

• Prev by Date: Re: Generating r-combination or r-permutations with replacement.
• Next by Date: Subscript problem
• Previous by thread: Re: Re: normal distribution random number generation
• Next by thread: Re: Re: normal distribution random number generation