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 >