From: Zelphir Kaltstahl <zelphirkaltstahl@posteo.de>
To: mikael@djurfeldt.com
Cc: guile-user <guile-user@gnu.org>
Subject: Re: Normal distribution random numbers
Date: Thu, 4 Jun 2020 17:08:12 +0200 [thread overview]
Message-ID: <65ab2d2a-d980-d2e1-ce9c-634390eff711@posteo.de> (raw)
In-Reply-To: <CAA2XvwKWk03LMmP-u1b9NBV_w4zE_kSMC-9RMgzEMVQE5b+PEw@mail.gmail.com>
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
> >
> >
>
next prev parent reply other threads:[~2020-06-04 15:08 UTC|newest]
Thread overview: 10+ messages / expand[flat|nested] mbox.gz Atom feed top
2020-05-30 20:21 Normal distribution random numbers Zelphir Kaltstahl
2020-05-30 20:42 ` Zelphir Kaltstahl
2020-06-04 15:03 ` Mikael Djurfeldt
2020-06-04 15:08 ` Zelphir Kaltstahl [this message]
2020-06-04 15:11 ` Mikael Djurfeldt
2020-06-04 15:22 ` Mikael Djurfeldt
2020-05-30 21:30 ` Arne Babenhauserheide
2020-05-30 23:16 ` Zelphir Kaltstahl
-- strict thread matches above, loose matches on Subject: below --
2020-05-31 15:12 tantalum
2020-05-31 19:54 ` Zelphir Kaltstahl
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
List information: https://www.gnu.org/software/guile/
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=65ab2d2a-d980-d2e1-ce9c-634390eff711@posteo.de \
--to=zelphirkaltstahl@posteo.de \
--cc=guile-user@gnu.org \
--cc=mikael@djurfeldt.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).