Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: Re: normal distribution random number generation

• To: mathgroup at smc.vnet.net
• Subject: [mg51309] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Thu, 14 Oct 2004 06:35:37 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

I think I have now fixed the problem reported by Mark Fisher. The new  
version of the RandomReplacement now works like this:


SeedRandom[5]


Timing[First[Table[Random[],{10000}]]]


{0.01 Second,0.786599}


<<RandomReplacement`

Timing[First[Table[Random[],{10000}]]]


{0.02 Second,0.409222}


As you can see there is a loss of speed involved. There may also be  
still some as yet undiscovered problems. However, the new package is  
available form my sites, though I still should write some documentation  
of the latest changes (a few bugs in the previous version have also  
been fixed).

Andrzej Kozlowski

Andrzej Kozlowski
Chiba, Japan
http://www.akikoz.net/~andrzej/
http://www.mimuw.edu.pl/~akoz/



On 12 Oct 2004, at 14:57, Andrzej Kozlowski wrote:

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

• Follow-Ups:
  • Re: Re: Re: Re: normal distribution random number generation
    • From: DrBob <drbob@bigfoot.com>