guile-examples/randomness/custom-distribution.scm

;; (From tantalum on the Guile user list.)

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

;; (Implementation from tantalum on the Guile user list. Changed
;; formatting. Added comments.)


(use-modules
 ;; Only import the procedure ~last~ from srfi-1.
 ((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)
            ;; If ~b~ is already the empty list, meaning there are no further
            ;; arguments given, then a is already the cumulative sum list.
            (list)
            ;; Otherwise continue building the list of cumulative sums. ~a~ is
            ;; already consed in this iteration. The next ~a~ will be the sum of
            ;; the current ~a~ and the next element. The next ~b~ will be the
            ;; tail of the current ~b~. This way we keep summing up in the next
            ;; ~a~ at each iteration.
            (apply cusum
                   (+ a (car b))
                   (cdr b)))))


(define* (random-discrete-f probabilities #:optional (state *random-state*))
  "create a procedure, which returns a random number according to the given
   probabilities of numbers.
   (real ...) [random-state] -> procedure:{-> integer}"
  (let*
      ;; Calculate the cumulative sums of the given probabilities.
      ([cuprob (apply cusum probabilities)]
       ;; Get the last (highest) cumulative sum from the list of cumulative
       ;; sums. This will serve as an upper exclusive bound for the random
       ;; number generation.
       [sum (last cuprob)])
    ;; Return a procedure.
    (lambda ()
      ;; Generate a random positive integer, uniformly distributed, using the
      ;; optionally given random state, maximum excluve upper bound - 1. This
      ;; results in a random integer, which is smaller than the last cumulative
      ;; sum. This means, that there will always be a cumulative sum, which is
      ;; found to be greater than the random integer.
      (let ([deviate (random sum state)])
        ;; Start with index ~a~ being 0 and the current probability
        (let loop ([a 0] [b cuprob])
          ;; Check if the end of the list is reached. If the end of the list is
          ;; reached, return the last index of the list. In theory this case
          ;; should not happen.
          (if (null? b)
              a
              ;; Look for the cumulative sum, for which the generated random
              ;; number ~deviate~ is smaller than the cumulative sum. Then
              ;; return its index. This index is the random number generated by
              ;; the returned random number generating procedure.
              (if (< deviate (car b))
                  a
                  ;; Recur. Increase index by 1 and look in the tail.
                  (loop (+ 1 a) (cdr b)))))))))


;; =====
;; USAGE
;; =====

;; Create a random procedure.
(define random*
  (random-discrete-f (list 1 0 3 1)
                     (random-state-from-platform)))


;; Show a random number generated by the previously created random procedure.
(display (random*))