Lazy generator for Gaussian random samples using efficient box-muller algorithm

randn.clj

(ns statistics.lazy.rand 
  (:require [clojure.math :as math]))

(defn rand-gaussian
  "Get a random sample from the standard Gaussian distribution. 
  Optionall specify the mean m and the standard deviation sd. E.g.:
  (rand-gaussian) # => 0.1324...
  (rand-gaussian 5 0.1) # => 5.3397...
  "
  ([] (rand-gaussian 0 1))
  ([m sd]
    (defn scale [x] (+ m (* sd x)))
    (def p (rand))
    (def q (rand))

  ; We use the Box-Muller transform
  (let [rho (clojure.math/sqrt (* -2 (clojure.math/log q)))
        theta (* 2 clojure.math/PI p)
        _box    (* rho (clojure.math/cos theta))
        _muller (* rho (clojure.math/sin theta))
        box    (scale _box)
        muller (scale _muller)]

    ; (clojure.pprint/pprint {:mean m :sd sd 
    ;                          :_box _box :_muller _muller 
    ;                          :box box :muller muller})

    ; tail call
    (lazy-seq (cons box 
              (cons muller 
              (rand-gaussian m sd)))))))

(defn sample-n
  "Generate n samples based on the random sampler f. E.g.
  (sample-n |(rand-int 0 3) 4) # => @[0 1 2 0]
  (sample-n |(rand-uniform) 4)
  "
  [f n]
  (take n ; (generate [_ :iterate true])
            (f)))


(take 4  (rand-gaussian)) ; =>
; (-1.100044284767874
;  -0.8147643792371908
;  2.09931093798589
;  -0.8900201976621036)

(rand-gaussian) ; generates infinitely