Re: Normal distribution random numbers

[Mikael Djurfeldt <mikael@djurfeldt.com>]

Hi Zelphir,

random:normal actually uses the Box-Muller-transform. But since it uses 64
bits, we only loose values that would be generated once in 2*10^20. That
is, if we could draw one billion numbers per second, such values would be
drawn once in 7000 years. So, we would start noticing an anomaly after
maybe 100000 years or so.

But maybe we should replace this with some more correct and efficient
algorithm at some point.

Best regards,
Mikael

Den lör 30 maj 2020 22:43Zelphir Kaltstahl <zelphirkaltstahl@posteo.de>
skrev:

> I just realized, that I did not check what Guile implements as
> non-SRFIs. I found:
> https://www.gnu.org/software/guile/manual/html_node/Random.html which
> has `random:normal`! I should have checked that first. Still good to
> know, what a can of worms normal distribution implementation can be.
>
> On 30.05.20 22:21, Zelphir Kaltstahl wrote:
> > Hi Guile Users!
> >
> > I recently wrote a little program involving lots of uniformly
> > distributed random integers. For that I used SRFI-27 and it works fine.
> >
> > Then I thought: How would I get normal distributed random numbers? I
> > don't have a project or program in mind for this, but it struck me, that
> > I do not know, how to get a normal distribution from a uniform
> > distribution. So I dug into the matter …
> >
> > Turns out the math is not really my friend:
> >
> > * https://stackoverflow.com/a/3265174 – OK, if that's true, then don't
> > use Box-Muller-Transform
> > * https://stackoverflow.com/a/86885 – The what? I need to somehow
> > inverse the Gaussian distribution to get a function to calculate normal
> > distributed values from uniformly distributed values? Something like
> > that. Safe to say it is above my current math skills.
> > * The wiki page also does not help me much:
> > https://en.wikipedia.org/wiki/Inverse_transform_sampling Seems too
> > complicated.
> >
> > So I thought: "OK, maybe I can simply copy, how other languages
> > implement it!" The wiki page mentions, that R actually makes use of the
> > inverse thingy. So I set out to look at R source code:
> >
> > * https://github.com/wch/r-source/blob/master/src/nmath/rnorm.c – OK,
> > looks simple enough … Lets see what `norm_rand` is …
> > * https://github.com/wch/r-source/blob/master/src/nmath/snorm.c#L62 –
> > yeah … well … I'm not gonna implement _that_ pile of … Just look at the
> > lines
> > https://github.com/wch/r-source/blob/master/src/nmath/snorm.c#L135-L196
> > what a mess! Not a single comment to help understanding in it. Such a
> > disappointment.
> > * Python also seems to only use an approximation with magic constants:
> > https://github.com/python/cpython/blob/3.8/Lib/random.py#L443
> >
> > So it seems, that there is no easy way to implement it properly with
> > correct tails to the left and right side of the distribution, something
> > clean and not made with mathematical traps built-in. Or is there?
> >
> > I found a post about using 2 normal distributions to do
> > Box-Muller-transform:
> >
> https://www.alanzucconi.com/2015/09/16/how-to-sample-from-a-gaussian-distribution/
> >
> > However, it seems to require a uniform float not integer and it is the
> > Box-Muller-transform, which is said to clamp between -6 and 6 according
> > to the people writing the answers on stackoverflow.
> >
> > So my question is: Is there a good implementation in the Guile universe
> > already? (Or a simple way to implement it?) I don't really need it right
> > now, but I think this thing could be an obstacle for many people without
> > serious math knowledge and it would be good to know, where to find it,
> > should one have need for normal distributed random numbers.
> >
> > Regards,
> > Zelphir
> >
> >
>
>