Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: normal distribution random number generation

• To: mathgroup at smc.vnet.net
• Subject: [mg51278] Re: [mg51265] Re: normal distribution random number generation
• From: DrBob <drbob at bigfoot.com>
• Date: Tue, 12 Oct 2004 01:57:41 -0400 (EDT)
• References: <ck0be4$nru$1@smc.vnet.net> <200410110525.BAA05021@smc.vnet.net>
• Reply-to: drbob at bigfoot.com
• Sender: owner-wri-mathgroup at wolfram.com

He's right! How the devil does that happen?

Here's a test, with Andrzej's package loaded in my Init file.

First, with n=250:

Quit
SeedRandom[5]
Table[Random[],{250}];
Last@%

0.107874

Unprotect[Random];
Clear@Random
SeedRandom[5]
Table[Random[],{250}];
Last@%

0.107874

Now, with n=249:

Quit
SeedRandom[5]
Table[Random[],{249}];
Last@%

0.656099

Unprotect[Random];
Clear@Random
SeedRandom[5]
Table[Random[],{249}];
Last@%

0.0708373

Bobby

On Mon, 11 Oct 2004 01:25:24 -0400 (EDT), Mark Fisher <mark at markfisher.net> wrote:

> FYI: I've just a little testing and I find that Mathematica ignors the user
> defined rules for Random in Table[Random[],{n}] when n >= 250.
>
> Andrzej Kozlowski wrote:
>> 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
>>>
>>
>>
>
>
>
>



-- 
DrBob at bigfoot.com
www.eclecticdreams.net

• References:
  • Re: normal distribution random number generation
    • From: Mark Fisher <mark@markfisher.net>