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:03:30 +0200 Message-ID: References: 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="104628"; 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:04:02 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 1jgrPW-000R4r-F7 for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 17:04:02 +0200 Original-Received: from localhost ([::1]:43942 helo=lists1p.gnu.org) by lists.gnu.org with esmtp (Exim 4.90_1) (envelope-from ) id 1jgrPV-0000FV-F9 for guile-user@m.gmane-mx.org; Thu, 04 Jun 2020 11:04:01 -0400 Original-Received: from eggs.gnu.org ([2001:470:142:3::10]:42280) by lists.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.90_1) (envelope-from ) id 1jgrPE-0000FJ-Uk for guile-user@gnu.org; Thu, 04 Jun 2020 11:03:44 -0400 Original-Received: from mail-vk1-f170.google.com ([209.85.221.170]:39849) by eggs.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_128_GCM_SHA256:128) (Exim 4.90_1) (envelope-from ) id 1jgrPD-0006UR-34 for guile-user@gnu.org; Thu, 04 Jun 2020 11:03:44 -0400 Original-Received: by mail-vk1-f170.google.com with SMTP id u15so1422321vkk.6 for ; Thu, 04 Jun 2020 08:03:42 -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=6smmjw54iEQwjZr4KaffKRu/Upyj8HwFUb/Zr3L2B2w=; b=RQiyLm46qGNkxjyFAh3s6UkR4HwobSv4pye5qY8AdZhRLg2tE4Kh5yA41QKSavnjpR RMGlDhjsGTG1eTWLPHXdsrGidaBadrRWzkkgz9Ck17ltxDjPRhrY7UjvUtcMEpszYNc7 /Hiq494NHO9h37U6SMCQqevZ6AyWGGqq2IBfmF4WhnjVoMAASVQAcq/IxVGliKS4EO0V TFydrNQkR3VZwjkFewCRySMR0YiBQnN7+fKmsb7z5Bk/6KMSagei5cRHiN1wyRccTlte yji2ohcM7EQyiSHKnAWQ6uAkeYN7cAu3CTeV3SD1/TN3VeuPB9BQ4MwO4Od2R6KP9ost Ow8g== X-Gm-Message-State: AOAM532e734FkZJlfgkyF+YVlRIsmYk6dod3lbVfmZC2swrdWxVYhwUK NDB/pqT5wf+W4hNJcaWb+t3Ryll/j+TPu8V5k9axhg== X-Google-Smtp-Source: ABdhPJxw32B1U0VG4ltTLh22vwClDmkPvbC3+mTO6AGeVoEj3Q5IiUNTRGUStwcdecoDgSlZA+PWy83vI3mfVR4myGY= X-Received: by 2002:a1f:4185:: with SMTP id o127mr1438416vka.45.1591283021715; Thu, 04 Jun 2020 08:03:41 -0700 (PDT) In-Reply-To: Received-SPF: pass client-ip=209.85.221.170; envelope-from=mdjurfeldt@gmail.com; helo=mail-vk1-f170.google.com X-detected-operating-system: by eggs.gnu.org: First seen = 2020/06/04 11:03:42 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:16546 Archived-At: 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=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, tha= t > > 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, the= n don't > > use Box-Muller-Transform > > * https://stackoverflow.com/a/86885 =E2=80=93 The what? I need to someh= ow > > 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 =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, 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-dist= ribution/ > > > > 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 righ= t > > now, but I think this thing could be an obstacle for many people withou= t > > 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 > > > > > >