Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51324] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Thu, 14 Oct 2004 06:36:50 -0400 (EDT)
- References: <20041013010600.74710.qmail@web60505.mail.yahoo.com>
- Sender: owner-wri-mathgroup at wolfram.com
Are you using the latest version of the RandomReplacement package? (The
one I posted last night).
The earlier one had a bug and ddi not work with calls to
Random[Real,{a,b}] where a and b are integers. You had to use
Random[Real,N[{a,b}]].
I think (I hope) I fixed that in the latest version.
Andrzej
On 13 Oct 2004, at 10:06, sean kim wrote:
> *This message was transferred with a trial version of CommuniGate(tm)
> Pro*
> I know... It's werid for me... which i wanted to ask
> you all.
>
> when i used the code below I still got the low 24 bin
> value.
>
> In[1]:=<<RandomReplacement`
>
> In[2]:=
> SeedRandom[1111111 ];
> vec1=Table[Random[Real,{0,1000}],{10^6}];
> vec2=Map[If[#<44,1,0]&,vec1];
> ones=Flatten[Position[vec2,1]];
> plotvec=Map[Length,Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];
> ListPlot[Log[plotvec],PlotRange->All]
> ListPlot[plotvec,PlotRange->All]
>
> used it exactly as it is shown.
>
>
> but the below works, it seems...
>
> SeedRandom[1111111];
> vec1 = Table[Random[Integer,1000],{10^6}];
> vec2 = Map[If[#<=44,1,0]&, vec1];
> ones = Flatten[Position[vec2,1]];
> Length[ones]/10.^6
> plotvec =
> Map[Length,
> Split[Sort[Drop[ones-RotateRight[ones,1],1]]]];
> ListPlot[
> Log[plotvec], PlotRange->All]
> ListPlot[plotvec, PlotRange->All]
>
> maybe other people can try it and see if it's just me?
> (I took out the init.m btw)
>
> sean
>
>
>
> --- Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote:
>
>> *This message was transferred with a trial version
>> of CommuniGate(tm) Pro*
>> This is exactly what it is mean to do. The problem
>> was due to the
>> usage of the SWB algortithm by the built in Random[]
>> for generating
>> random reals. RandomReplace makes Mathematica use
>> the Wolfran CA random
>> number generator, which Mathematica uses for
>> generating random integers
>> instead of the SWB. This "cures" the problem
>> reported int hat message.
>> Have you got any reason to think it does not?
>>
>> Andrzej
>>
>>
>> On 13 Oct 2004, at 03:48, sean kim wrote:
>>
>>> *This message was transferred with a trial version
>> of CommuniGate(tm)
>>> Pro*
>>> Hi andrej
>>>
>>> thank you so much for making your package
>> available.
>>>
>>> I think I may have found another problem... ("may
>> have
>>> found" being operative phrase, or non-operative
>> for
>>> that matter)
>>>
>>> http://groups.google.com/groups?
>>>
>>
> hl=en&lr=&newwindow=1&safe=off&threadm=a5t1fq%245vt%241%40smc.vnet.net&
>>
>>> rnum=1&prev=/
>>>
>>
> groups%3Fhl%3Den%26lr%3D%26newwindow%3D1%26safe%3Doff%26q%3Dchange%2Bth
>>
>>>
>>
> is%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253
>>
>>> Dcomp.soft-sys.math.mathematica
>>>
>>> shows a thread dealing with random reals and
>> getting a
>>> much smaller frequency at position 24 than others
>>> positions.
>>>
>>> Daniel posted a work around showing the use of
>> Random
>>> Integers to fix the problem, which as far I
>>> understand, the Random replacement also
>> implements.
>>>
>>> It seems the RandomReplacement doesn't take care
>> of
>>> the appearance of the lower than expected
>> frequency at
>>> position 24(or interval 24 as the OP had noted).
>>>
>>> I must be missing something.
>>>
>>> thanks all for any insights.
>>>
>>> sean
>>>
>>>
>>> --- Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote:
>>>
>>>> *This message was transferred with a trial
>> version
>>>> of CommuniGate(tm) Pro*
>>>> 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:
>>>>
>>>>> *This message was transferred with a trial
>> version
>>>> of CommuniGate(tm)
>>>>> Pro*
>>>>> 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
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>> __________________________________
>>> Do you Yahoo!?
>>> Yahoo! Mail Address AutoComplete - You start. We
>> finish.
>>> http://promotions.yahoo.com/new_mail
>>>
>>
>>
>
>
> __________________________________________________
> Do You Yahoo!?
> Tired of spam? Yahoo! Mail has the best spam protection around
> http://mail.yahoo.com
>