Normal distribution random numbers

[Zelphir Kaltstahl <zelphirkaltstahl@posteo.de>]

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