Re: Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51395] Re: [mg51309] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: DrBob <drbob at bigfoot.com>
- Date: Fri, 15 Oct 2004 02:48:49 -0400 (EDT)
- References: <200410141035.GAA14839@smc.vnet.net>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
The speed difference can be pretty dramatic:
Quit
a=First@Timing[Table[Random[Real,{0,1000}],{10^6}];]
<<RandomReplacement`
b=First@Timing[Table[Random[Real,{0,1000}],{10^6}];]
b/a
0.375 Second
10.078 Second
26.8747
I can probably live with it, however, in cases where I'm concerned about built-in Random's non-random behavior.
Bobby
On Thu, 14 Oct 2004 06:35:37 -0400 (EDT), Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote:
> 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
>>>
>>
>>
>
>
>
>
--
DrBob at bigfoot.com
www.eclecticdreams.net
- References:
- Re: Re: Re: normal distribution random number generation
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Re: Re: normal distribution random number generation