[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
>>
>
>
Prev by Date:
**Re: argMax**
Next by Date:
**Re: Eigenvalues and eigenvectors of a matrix with nonpolynomial elements.**
Previous by thread:
**Re: normal distribution random number generation**
Next by thread:
**Re: Re: Re: Re: normal distribution random number generation**
| |