Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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
>


  • Prev by Date: Re: cross-product in cylindrical problem
  • Next by Date: Re: Re: symbol replace
  • Previous by thread: Re: Re: normal distribution random number generation
  • Next by thread: Re: Re: normal distribution random number generation