Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51323] 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:38 -0400 (EDT)
- Reply-to: sean_incali at yahoo.com
- Sender: owner-wri-mathgroup at wolfram.com
i just took a look at Part[plotvec, 24]
the respective values, for me, are,
375 and 686.
so i did a test with the following.
in the fresh kernel, I get exactly the same thing as
you have.
but I decided to evaluate the code
In[5]:=
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]]]];
Part[plotvec, 24]
ListPlot[Log[plotvec],PlotRange->All]
ListPlot[plotvec,PlotRange->All]
Out[10]=
375
and then
reevaluate the codes
In[13]:= SeedRandom[5]
In[14]:= Timing[First[Table[Random[],{10000}]]]
Out[14]= {0.01 Second,0.786599}
In[15]:= <<RandomReplacement`
Timing[First[Table[Random[],{10000}]]]
Out[16]= {0.01 Second,0.611425}
and I get the different value than before(one with the
fresh kernel)
i wonder what's going on here?
sean
--- 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:
>
> > *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!?
Read only the mail you want - Yahoo! Mail SpamGuard.
http://promotions.yahoo.com/new_mail