Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51322] Re: [mg51283] Re: [mg51265] Re: normal distribution random number generation
- From: sean kim <sean_incali at yahoo.com>
- Date: Thu, 14 Oct 2004 06:36:32 -0400 (EDT)
- Reply-to: sean_incali at yahoo.com
- Sender: owner-wri-mathgroup at wolfram.com
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 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