From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.io!.POSTED.ciao.gmane.io!not-for-mail From: Mikael Djurfeldt Newsgroups: gmane.lisp.guile.user Subject: Re: Normal distribution random numbers Date: Thu, 4 Jun 2020 17:11:07 +0200 Message-ID: References: <65ab2d2a-d980-d2e1-ce9c-634390eff711@posteo.de> Reply-To: mikael@djurfeldt.com Mime-Version: 1.0 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable Injection-Info: ciao.gmane.io; posting-host="ciao.gmane.io:159.69.161.202"; logging-data="10856"; mail-complaints-to="usenet@ciao.gmane.io" Cc: guile-user To: Zelphir Kaltstahl Original-X-From: guile-user-bounces+guile-user=m.gmane-mx.org@gnu.org Thu Jun 04 17:12:11 2020 Return-path: Envelope-to: guile-user@m.gmane-mx.org Original-Received: from lists.gnu.org ([209.51.188.17]) by ciao.gmane.io with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.92) (envelope-from ) id 1jgrXO-0002j4-Te for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 17:12:11 +0200 Original-Received: from localhost ([::1]:48566 helo=lists1p.gnu.org) by lists.gnu.org with esmtp (Exim 4.90_1) (envelope-from ) id 1jgrXN-0003gY-VM for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 11:12:09 -0400 Original-Received: from eggs.gnu.org ([2001:470:142:3::10]:43278) by lists.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.90_1) (envelope-from ) id 1jgrWc-00030K-Lw for guile-user@gnu.org; Thu, 04 Jun 2020 11:11:22 -0400 Original-Received: from mail-vs1-f51.google.com ([209.85.217.51]:36275) by eggs.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_128_GCM_SHA256:128) (Exim 4.90_1) (envelope-from ) id 1jgrWb-0000yz-2z for guile-user@gnu.org; Thu, 04 Jun 2020 11:11:22 -0400 Original-Received: by mail-vs1-f51.google.com with SMTP id j13so3727187vsn.3 for ; Thu, 04 Jun 2020 08:11:20 -0700 (PDT) X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:mime-version:references:in-reply-to:reply-to :from:date:message-id:subject:to:cc; bh=2qKZPVXg2Ey0e+WyeKDWLmtLp05Bji3ufAiDiyd9SSE=; b=il+BIO4h2EL+HbOPnil7EUteN+gFGPAn6K7u9sOu5YSJO9jl9LUsa1pT6KHxSyybNq 0B3fspNoCDdDvFmkB57Xb4dYxAqTRxgV6tg8/etrTNpPYOqjALIhWwZKnq3ZEZNypfGm Du8D1ZC1vd9TOPX4jJrDIjXEieWYVmQQCOzmQcGVaAwEJxkutTo/XOTuQrvy/MYMrMOB ASS9Qcttto/uwmYUxGqRjjKBxsQPzNoa3cvT3gFqwXxNVV6eOaMW/GLG1nQ0YIMZvA1l Vn7lkKXEL+HpTXlBQjwDoTtEsj0LbSh/uFLN2jUoWHNxtuaJp4x5lHgt/g80JtO01Sro sT8g== X-Gm-Message-State: AOAM532pnCqvjJncx8vJk8GCxuPqIZ7oADWKjTVUaIwYMQAxa1mgLQd0 6Dvz4N1BD25ODe5KSiW/B7olHnIizcsAxcFhyTY= X-Google-Smtp-Source: ABdhPJx186q1gT1r/ixWkmuN8nMakiQKBAkoDWc/RhSqKFHUdYrYm0R+lUmWeUgVfKnPJzw7KB9MHXgIbGr5/1ZRutY= X-Received: by 2002:a67:ec11:: with SMTP id d17mr3772345vso.101.1591283479777; Thu, 04 Jun 2020 08:11:19 -0700 (PDT) In-Reply-To: <65ab2d2a-d980-d2e1-ce9c-634390eff711@posteo.de> Received-SPF: pass client-ip=209.85.217.51; envelope-from=mdjurfeldt@gmail.com; helo=mail-vs1-f51.google.com X-detected-operating-system: by eggs.gnu.org: First seen = 2020/06/04 11:11:20 X-ACL-Warn: Detected OS = Linux 2.2.x-3.x [generic] [fuzzy] X-Spam_score_int: -15 X-Spam_score: -1.6 X-Spam_bar: - X-Spam_report: (-1.6 / 5.0 requ) BAYES_00=-1.9, FREEMAIL_FORGED_FROMDOMAIN=0.001, FREEMAIL_FROM=0.001, HEADER_FROM_DIFFERENT_DOMAINS=0.249, HTML_MESSAGE=0.001, RCVD_IN_DNSWL_NONE=-0.0001, RCVD_IN_MSPIKE_H2=-0.001, SPF_PASS=-0.001, URIBL_BLOCKED=0.001 autolearn=_AUTOLEARN X-Spam_action: no action X-Content-Filtered-By: Mailman/MimeDel 2.1.23 X-BeenThere: guile-user@gnu.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: General Guile related discussions List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-user-bounces+guile-user=m.gmane-mx.org@gnu.org Original-Sender: "guile-user" Xref: news.gmane.io gmane.lisp.guile.user:16548 Archived-At: Yes Den tors 4 juni 2020 17:08Zelphir Kaltstahl skrev: > 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, tha= t > 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 6= 4 > 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=C3=B6r 30 maj 2020 22:43Zelphir Kaltstahl > 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, th= at >> > I do not know, how to get a normal distribution from a uniform >> > distribution. So I dug into the matter =E2=80=A6 >> > >> > Turns out the math is not really my friend: >> > >> > * https://stackoverflow.com/a/3265174 =E2=80=93 OK, if that's true, th= en don't >> > use Box-Muller-Transform >> > * https://stackoverflow.com/a/86885 =E2=80=93 The what? I need to some= how >> > inverse the Gaussian distribution to get a function to calculate norma= l >> > 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 th= e >> > inverse thingy. So I set out to look at R source code: >> > >> > * https://github.com/wch/r-source/blob/master/src/nmath/rnorm.c =E2=80= =93 OK, >> > looks simple enough =E2=80=A6 Lets see what `norm_rand` is =E2=80=A6 >> > * https://github.com/wch/r-source/blob/master/src/nmath/snorm.c#L62 = =E2=80=93 >> > yeah =E2=80=A6 well =E2=80=A6 I'm not gonna implement _that_ pile of = =E2=80=A6 Just look at the >> > lines >> > https://github.com/wch/r-source/blob/master/src/nmath/snorm.c#L135-L19= 6 >> > 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, somethin= g >> > 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-dis= tribution/ >> > >> > 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 accordin= g >> > to the people writing the answers on stackoverflow. >> > >> > So my question is: Is there a good implementation in the Guile univers= e >> > already? (Or a simple way to implement it?) I don't really need it rig= ht >> > now, but I think this thing could be an obstacle for many people witho= ut >> > 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 >> > >> > >> >>