Re: Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51396] Re: [mg51309] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Fri, 15 Oct 2004 02:48:52 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
I know. I now have some ideas about how to make it a bit faster, but
need to test them first. But there is no way it can ever be as fast as
the built-in Random. I do hope the decision makers at WRI take this
issue more seriously and do something about it in Mathematica 6. As
computers grow faster the problem is becoming more visible and the
number of users affected by it seems to be growing.
Andrzej
On 15 Oct 2004, at 14:57, DrBob wrote:
> 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
>