Re: Normal distribution random numbers

[Zelphir Kaltstahl]

Hi Mikael!

Thanks for putting that into perspective and giving some numbers!

When I looked at the code of Guile for random:normal, I also guessed,
that it makes use of that Box-Muller-transform, but wasn't sure, so
thanks for confirming that as well.

So basically the tails are wrong, but to draw a number in the area where
the tails are wrong is so unlikely, that it would take that much time,
as stated in your number example, if I understand this correctly(?/.)

Regards,

Zelphir

On 04.06.20 17:03, Mikael Djurfeldt wrote:
> 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
> <mailto: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
>     >
>     >
>