Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: normal distribution random number generation

• To: mathgroup at smc.vnet.net
• Subject: [mg51283] Re: [mg51265] Re: normal distribution random number generation
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Tue, 12 Oct 2004 01:57:46 -0400 (EDT)
• References: <ck0be4$nru$1@smc.vnet.net> <200410110525.BAA05021@smc.vnet.net> <opsfo5nlc8iz9bcq@monster.cox-internet.com>
• Sender: owner-wri-mathgroup at wolfram.com

It's clearly because of packed arrays:


<<Developer`


PackedArrayQ[Table[Random[],{249}]]


False


PackedArrayQ[Table[Random[],{250}]]


True



Obviously Mathematica uses a different method of generating random 
packed arrays. There are various ways of getting round this (e.g. 
constructing lists as joins of lists of length less than 250 etc) but 
one will loose the benefits of packed arrays and with that, presumably, 
a lot of speed. This seems to be something that only WRI can change.
Note however that not all situations in which one uses Random[] require 
generating this type of lists.
Moreover, for problems involving the normal distribution it seems to me 
that using the Marsaglia generator gives good results even with the 
built in Random[]. But certainly we should hope that something will 
finally be done do deal with this issue in version 6.

Andrzej




On 11 Oct 2004, at 17:24, DrBob wrote:

> 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>