unofficial mirror of guile-user@gnu.org 
 help / color / mirror / Atom feed
* Normal distribution random numbers
@ 2020-05-30 20:21 Zelphir Kaltstahl
  2020-05-30 20:42 ` Zelphir Kaltstahl
  2020-05-30 21:30 ` Arne Babenhauserheide
  0 siblings, 2 replies; 10+ messages in thread
From: Zelphir Kaltstahl @ 2020-05-30 20:21 UTC (permalink / raw)
  To: Guile User

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




^ permalink raw reply	[flat|nested] 10+ messages in thread
* Re: Normal distribution random numbers
@ 2020-05-31 15:12 tantalum
  2020-05-31 19:54 ` Zelphir Kaltstahl
  0 siblings, 1 reply; 10+ messages in thread
From: tantalum @ 2020-05-31 15:12 UTC (permalink / raw)
  To: guile-user

surely not the ideal way to generate numbers with a normal distribution, 
but there is a way to use custom probabilities from a list, which i 
think is nice to know.
it works like this:

have a list of probabilities. can be as long as you want. i think that 
is called probability density.
->
(1 0 3 1)

create the cumulative sums for this list. that is, sums like this
   (a b c ...) -> (a (+ a b) (+ a b c) ...)
i think that is called cumulative distribution.
->
(1 1 4 5)

create a random number up to the largest sum
->
(random 5)

return the first index of the list of cumulative sums that is greater 
than the random number.
given the distribution above, you would see index 2 a lot, never index 
1, and index 0 and 3 rarely.

~~~
(use-modules ((srfi srfi-1) #:select (last)))

(define (cusum a . b)
   "calculate cumulative sums from the given numbers.
    (a b c ...) -> (a (+ a b) (+ a b c) ...)"
   (cons a (if (null? b) (list) (apply cusum (+ a (car b)) (cdr b)))))

(define* (random-discrete-f probabilities #:optional (state 
*random-state*))
   "(real ...) [random-state] -> procedure:{-> integer}"
   (let* ((cuprob (apply cusum probabilities)) (sum (last cuprob)))
     (lambda ()
       (let ((deviate (random sum state)))
         (let loop ((a 0) (b cuprob))
           (if (null? b) a (if (< deviate (car b)) a (loop (+ 1 a) (cdr 
b)))))))))

(define random* (random-discrete-f (list 1 0 3 1) 
(random-state-from-platform)))
(display (random*))
~~~




^ permalink raw reply	[flat|nested] 10+ messages in thread
end of thread, other threads:[~2020-06-04 15:22 UTC | newest]

Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
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
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

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).