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:22:35 +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="119240"; 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:35:47 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 1jgruD-000UcB-PZ for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 17:35:45 +0200 Original-Received: from localhost ([::1]:48364 helo=lists1p.gnu.org) by lists.gnu.org with esmtp (Exim 4.90_1) (envelope-from ) id 1jgruC-0001mk-Q4 for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 11:35:44 -0400 Original-Received: from eggs.gnu.org ([2001:470:142:3::10]:44262) by lists.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.90_1) (envelope-from ) id 1jgrhh-00033X-J1 for guile-user@gnu.org; Thu, 04 Jun 2020 11:22:50 -0400 Original-Received: from mail-vs1-f41.google.com ([209.85.217.41]:39989) by eggs.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_128_GCM_SHA256:128) (Exim 4.90_1) (envelope-from ) id 1jgrhf-00033F-FN for guile-user@gnu.org; Thu, 04 Jun 2020 11:22:49 -0400 Original-Received: by mail-vs1-f41.google.com with SMTP id u17so3732735vsu.7 for ; Thu, 04 Jun 2020 08:22:47 -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=PDZj6CjpzMRxncIYp3rVP1aXUiziTtkzeMNAptz0+8U=; b=NPrmromXnIO16MGqO0RQQYqaniEnv3IsiDC5paKN7IDjopXoQkl+8H83HvvRJ6EzVC 9LvLtg/rrI+cFZYKPuGLJwqLpn1mhAf1TFWlG5+ysBMZzSwlUqBidoyRPKRY4mMBfKyS ihezrfX2fPc4pL0FwrltcmNvcwdKS0GIRXspj27LQx/9DtwyIvY4Z52Hs6vlQFs6Yw5/ kKKpK+EtwuRgr92R20EK/oSilNTtpba4MfACiAHyb067ScHoxXFDh5hCzPqszG8nnnPl bcv49H9nRZMhvtlehcFEBVyIwvsuFGZaMOG+eEPlLQ8SkI4WHWjd5IcMvE5uRoEK36zX Rwuw== X-Gm-Message-State: AOAM531GRAiXkoCuLU0k7lQL7dZhbUENCCrGgZBMMgxv9TwRVGj2wLhp GmusHHYDdHC+8WQxSn/B+NZ/YtwEaGMVu+s2h4g= X-Google-Smtp-Source: ABdhPJzFqWac+4BBj8weo8viFwuwlAd5YCPUvjnTLDKsZRhtTOZcZON6ApozxnZ1LVsL9ozk3G9Sghpx6vW+TYN+Hlc= X-Received: by 2002:a67:ec11:: with SMTP id d17mr3819611vso.101.1591284166333; Thu, 04 Jun 2020 08:22:46 -0700 (PDT) In-Reply-To: Received-SPF: pass client-ip=209.85.217.41; envelope-from=mdjurfeldt@gmail.com; helo=mail-vs1-f41.google.com X-detected-operating-system: by eggs.gnu.org: First seen = 2020/06/04 11:22:46 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:16549 Archived-At: (Maybe I was sloppy now. An anomaly might actually be detected sooner due to increasing distances in the discretization of values in the tail, but I still don't think this is something to worry about.) Den tors 4 juni 2020 17:11Mikael Djurfeldt skrev: > 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, >> that it makes use of that Box-Muller-transform, but wasn't sure, so than= ks >> 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, a= s >> 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. T= hat >> is, if we could draw one billion numbers per second, such values would b= e >> 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 fin= e. >>> > >>> > 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 =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, t= hen don't >>> > use Box-Muller-Transform >>> > * https://stackoverflow.com/a/86885 =E2=80=93 The what? I need to som= ehow >>> > inverse the Gaussian distribution to get a function to calculate norm= al >>> > 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 t= he >>> > 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-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, somethi= ng >>> > 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-di= stribution/ >>> > >>> > However, it seems to require a uniform float not integer and it is th= e >>> > Box-Muller-transform, which is said to clamp between -6 and 6 accordi= ng >>> > to the people writing the answers on stackoverflow. >>> > >>> > So my question is: Is there a good implementation in the Guile univer= se >>> > 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 >>> > >>> > >>> >>>