Re: Re: Re: normal distribution random number generation
- To: mathgroup at smc.vnet.net
- Subject: [mg51316] 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:04 -0400 (EDT)
- Reply-to: sean_incali at yahoo.com
- Sender: owner-wri-mathgroup at wolfram.com
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%2Bthis%2Bnumber%2Bto%2Bget%2Banother%2Brandom%2Bsequence%26meta%3Dgroup%253Dcomp.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:
> 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