Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: normal distribution random number generation

• To: mathgroup at smc.vnet.net
• Subject: [mg51143] Re: [mg51113] normal distribution random number generation
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Wed, 6 Oct 2004 04:34:18 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

The problem is not actually with the way Mathematica's 
NormalDistribution but with the uniformly distributed Random[] 
function.
NormalDistribution itself is generated (in the package 
Statistics`NormalDistribution`) by means of the very classical so 
called Boox-Muller method (although actually the Marsaglia variant 
below works better). You can try downolading my little 
RandomReplacement package form one of my web sites:

http://www.akikoz.net/~andrzej//Mathematica/  in Japan

or

http://www.mimuw.edu.pl/~akoz/Mathematica/   in Poland

The package is based on the post of daniel Lichtblau, which explains 
the nature of the problem:

<http://library.wolfram.com/mathgroup/archive/2000/May/msg00088.html>



The package is named init.m and can be loaded automatically at the 
start of each Mathematica session.
However, in fact if you are only concerned with normal distributions it 
may be enough to the following. First load the (unmodified)  
NormalDistribution package using

<<Statistics`NormalDistribution`

and then evaluate the code below. It will replace the Box-Muller 
generator by the Marsaglia one. I have found that this is usually 
enough. But if your problem still remains try the RandomReplacement 
package. (Or you can use both of course). Here is the Marsaglia code 
for normal distribution:


Statistics`NormalDistribution`Private`normal=
     Compile[{{mu,_Real},{sigma,_Real},{q1,_Real},{q2,_Real}},
       Module[{va=1.,vb,rad=2.,den=1.},
         While[rad>=1.,va=2.*Random[]-1.;vb=2.*Random[]-1.;
           rad=va*va+vb*vb];
         den=Sqrt[-2.*(Log[rad]/rad)];
         mu+sigma*va*den]];



On 5 Oct 2004, at 17:36, Chris Kunkel wrote:

> *This message was transferred with a trial version of CommuniGate(tm) 
> Pro*
> Hello,
>
> Recently I have encountered a problem in Mathematica's normal
> distribution random number generator.  The problem arises when I look
> at the distribution of averages of a list of theses numbers.  That is,
> I generate 1000 numbers and take their average.  I do this a number of
> times and the plot a frequency distribution.  Consistently it seems to
> be skewed positive.  Specifically, the number of occurrences less than
> 3 standard deviations is consistent with the expected number, but the
> number greater than 3 is always larger (in some cases about twice the
> expected number).
>
> I was wondering if anyone else has noticed this and knows of a fix.
> Also, does anyone have code for a good quality normal distribution
> random number generator?
>
> Chris
>